UPP  11.0.0
 All Data Structures Files Functions Variables Pages
CLDRAD.f
Go to the documentation of this file.
1 
79  SUBROUTINE cldrad
80 
81 !
82  use vrbls4d, only: dust,suso, salt, soot, waso,no3,nh4
83  use vrbls3d, only: qqw, qqr, t, zint, cfr, qqi, qqs, q, ext, zmid,pmid,&
84  pint, duem, dusd, dudp, duwt, dusv, ssem, sssd,ssdp,&
85  sswt, sssv, bcem, bcsd, bcdp, bcwt, bcsv, ocem,ocsd,&
86  ocdp, ocwt, ocsv, sca, asy,cfr_raw
87  use vrbls2d, only: cldefi, cfracl, avgcfracl, cfracm, avgcfracm, cfrach,&
88  avgcfrach, avgtcdc, ncfrst, acfrst, ncfrcv, acfrcv, &
89  hbot, hbotd, hbots, htop, htopd, htops, fis, pblh, &
90  pbot, pbotl, pbotm, pboth, cnvcfr, ptop, ptopl, &
91  ptopm, ptoph, ttopl, ttopm, ttoph, pblcfr, cldwork, &
92  aswin, auvbin, auvbinc, aswout,alwout, aswtoa, &
93  rlwtoa, czmean, czen, rswin, alwin, alwtoa, rlwin, &
94  sigt4, rswout, radot, rswinc, aswinc, aswoutc, &
95  aswtoac, alwoutc, aswtoac, avisbeamswin, &
96  avisdiffswin, aswintoa, aswtoac, airbeamswin, &
97  airdiffswin, dusmass, dusmass25, ducmass, ducmass25, &
98  alwinc, alwtoac, swddni, swddif, swdnbc, swddnic, &
99  swddifc, swupbc, lwdnbc, lwupbc, swupt, &
100  taod5502d, aerssa2d, aerasy2d, mean_frp, ebb, hwp, &
101  lwp, iwp, avgcprate, &
102  dustcb,sscb,bccb,occb,sulfcb,dustpm,sspm,aod550, &
103  du_aod550,ss_aod550,su_aod550,oc_aod550,bc_aod550, &
104  pwat,dustpm10,maod,no3cb,nh4cb,aqm_aod550
105  use masks, only: lmh, htm
106  use params_mod, only: tfrz, d00, h99999, qcldmin, small, d608, h1, rog, &
107  gi, rd, qconv, abscoefi, abscoef, stbol, pq0, a2, &
108  a3, a4
109  use ctlblk_mod, only: jsta, jend, spval, modelname, grib, cfld,datapd, &
110  fld_info, avrain, theat, ifhr, ifmin, avcnvc, &
111  tclod, ardsw, trdsw, ardlw, nbin_du, trdlw, im, &
112  nbin_ss, nbin_oc,nbin_bc,nbin_su,nbin_no3,dtq2, &
113  jm, lm, gocart_on, gccpp_on, nasa_on, me, rdaod, &
114  ista, iend,aqf_on
115  use rqstfld_mod, only: iget, id, lvls, iavblfld
116  use gridspec_mod, only: dyval, gridtype
117  use cmassi_mod, only: trad_ice
118  use machine_post, only: kind_phys
119  use upp_physics, only: calrh, calcape
120 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121  implicit none
122 !
123 ! SET CELSIUS TO KELVIN CONVERSION.
124  REAL,PARAMETER :: c2k=273.15, ptop_low=64200., ptop_mid=35000., &
125  ptop_high=15000.
126 !
127 ! DECLARE VARIABLES.
128 !
129 ! LOGICAL,dimension(im,jm) :: NEED
130  INTEGER :: lcbot,lctop,jc,ic !bsf
131  INTEGER,dimension(ista:iend,jsta:jend) :: ibott, ibotcu, ibotdcu, ibotscu, ibotgr, &
132  itopt, itopcu, itopdcu, itopscu, itopgr
133  REAL,dimension(im,jm) :: grid1
134  REAL,dimension(ista:iend,jsta:jend) :: grid2, egrid1, egrid2, egrid3, &
135  cldp, cldz, cldt, cldzcu
136  REAL,dimension(lm) :: rhb, watericetotal, pabovesfc
137  REAL :: watericemax, wimin, zcldbase, zcldtop, zpbltop, &
138  rhoice, coeffp, exponfp, const1, cloud_def_p, &
139  pcldbase, rhoair, vovermd, concfp, betav, &
140  vertvis, tx, tv, pol, esx, es, e, zsf, zcld, frac
141  integer nfog, nfogn(7),npblcld,nlifr, k1, k2, ll, ii, ib, n, jj, &
142  numr, numpts
143  real,dimension(lm) :: cldfra, cfr_layer_sum
144  real :: ceiling_thresh_cldfra, cldfra_max, &
145  zceil, zceil1, zceil2, previous_sum, &
146  ceil_min, ceil_neighbor
147 
148  real,dimension(im,jm) :: ceil
149 
150 ! B ZHOU: For aviation:
151  REAL, dimension(ista:iend,jsta:jend) :: tcld, ceiling
152  real cu_ir(lm), q_conv !bsf
153 !jw
154  integer i,j,l,k,ibot,itclod,lbot,ltop,itrdsw,itrdlw, &
155  llmh,itheat,ifincr,itype,itop,num_thick
156  real dpbnd,rrnum,qcld,rsum,tlmh,factrs,factrl,dp, &
157  opdepth, tmp,qsat,rhum,tcext,delz,dely,dy_m
158 !
159  real full_cld(im,jm) !-- Must be dimensioned for the full domain
160  real, allocatable :: full_ceil(:,:), full_fis(:,:)
161 !
162  real dummy(ista:iend,jsta:jend)
163  integer idummy(ista:iend,jsta:jend)
164  real full_dummy(im,jm)
165 !
166 ! --- Revision added for GOCART ---
167 
168 !! GOCART aerosol optical data from GSFC Mie code calculations,
169 !! mapped to 7 channels: 0.34, 0.44, 0.55, 0.66, 0.86, 1.63, 11.1 micron
170 !! data wvnum1/0.338, 0.430, 0.545, 0.62, 0.841, 1.628, 11.0/
171 !! data wvnum2/0.342, 0.450, 0.565, 0.67, 0.876, 1.652, 11.2/
172 
173  integer, parameter :: krhlev = 36 ! num of rh levels for rh-dep components
174  integer, parameter :: kcm1 = 5 ! num of rh independent aer species
175  integer, parameter :: kcm2 = 6 ! num of rh dependent aer species
176  integer, parameter :: nbdsw = 7 ! total num of sw bands
177  integer, parameter :: noaer = 20 ! unit for LUTs file
178  CHARACTER :: aerosolname(kcm2)*4, aerosolname_rd*4, aerosol_file*30
179  CHARACTER :: aername_rd*4, aeropt*3
180 
181 ! - aerosol optical properties: mass extinction efficiency
182  REAL, ALLOCATABLE :: extrhd_du(:,:,:), extrhd_ss(:,:,:), &
183  & extrhd_SU(:,:,:), extrhd_BC(:,:,:), &
184  & extrhd_OC(:,:,:), extrhd_NI(:,:,:)
185 
186 ! - aerosol optical properties: mass scattering efficienc
187  REAL, ALLOCATABLE :: scarhd_du(:,:,:), scarhd_ss(:,:,:), &
188  & scarhd_SU(:,:,:), scarhd_BC(:,:,:), &
189  & scarhd_OC(:,:,:), scarhd_NI(:,:,:)
190 
191 ! - aerosol optical properties: asymmetry factor
192  REAL, ALLOCATABLE :: asyrhd_du(:,:,:), asyrhd_ss(:,:,:), &
193  & asyrhd_SU(:,:,:), asyrhd_BC(:,:,:), &
194  & asyrhd_OC(:,:,:), asyrhd_NI(:,:,:)
195 
196 ! - aerosol optical properties: single scatter albedo
197  REAL, ALLOCATABLE :: ssarhd_du(:,:,:), ssarhd_ss(:,:,:), &
198  & ssarhd_SU(:,:,:), ssarhd_BC(:,:,:), &
199  & ssarhd_OC(:,:,:), ssarhd_NI(:,:,:)
200 
201 ! --- aerosol optical properties mapped onto specified spectral bands
202 ! - relative humidity independent aerosol optical properties: du
203  real (kind=kind_phys) :: extrhi(kcm2,nbdsw) ! extinction coefficient
204 
205 ! - relative humidity dependent aerosol optical properties: oc, bc, su, ss001-005
206  real (kind=kind_phys) :: extrhd(krhlev,kcm2,nbdsw) ! extinction coefficient
207 !
208  REAL,dimension(ista:iend,jsta:jend) :: p1d,t1d,q1d,egrid4
209 ! REAL, allocatable :: RH3D(:,:,:) ! RELATIVE HUMIDITY
210  real, allocatable:: rdrh(:,:,:)
211  integer, allocatable :: ihh(:,:,:)
212  REAL :: rh3d, drh0, drh1, ext01, ext02,sca01,asy01
213  INTEGER :: ih1, ih2,naero
214  INTEGER :: ios, indx, issam, isscm, isuso, iwaso, isoot, nbin
215  REAL :: ccdry, ccwet, ssam, sscm
216  REAL,dimension(ista:iend,jsta:jend) :: aod_du, aod_ss, aod_su, aod_oc, aod_bc, aod_ni, aod
217  REAL,dimension(ista:iend,jsta:jend) :: sca_du, sca_ss, sca_su, sca_oc,sca_bc, sca_ni,sca2d
218  REAL,dimension(ista:iend,jsta:jend) :: asy_du, asy_ss, asy_su, asy_oc, asy_bc,asy_ni,asy2d
219  REAL,dimension(ista:iend,jsta:jend) :: angst, aod_440, aod_860 ! FORANGSTROM EXPONENT
220  REAL :: ang1, ang2
221  INTEGER :: indx_ext(kcm2), indx_sca(kcm2)
222  LOGICAL :: laeropt, lext, lsca, lasy
223  LOGICAL :: laersmass
224  REAL, allocatable :: fpm25_du(:),fpm25_ss(:)
225  REAL, allocatable, dimension(:,:) :: rhosfc, smass_du_cr,smass_du_fn, &
226  & smass_ss_cr, smass_ss_fn, smass_oc,smass_bc, &
227  & smass_su, smass_cr, smass_fn
228  real (kind=kind_phys), dimension(KRHLEV) :: rhlev
229  data rhlev(:)/ .0, .05, .10, .15, .20, .25, .30, .35, &
230  & .40, .45, .50, .55, .60, .65, .70, .75, &
231  & .80, .81, .82, .83, .84, .85, .86, .87, &
232  & .88, .89, .90, .91, .92, .93, .94, .95, &
233  & .96, .97, .98, .99/
234 !
235  data aerosolname /'DUST', 'SALT', 'SUSO', 'SOOT', 'WASO', 'NITR'/
236 ! INDEX FOR TOTAL AND SPECIATED AEROSOLS (DU, SS, SU, OC, BC)
237  data indx_ext / 610, 611, 612, 613, 614, 615 /
238  data indx_sca / 651, 652, 653, 654, 655, 687 /
239  logical, parameter :: debugprint = .false.
240  logical :: model_pwat
241 !
242 !
243 !*************************************************************************
244 ! START CLDRAD HERE.
245 !
246 !*** BLOCK 1. SOUNDING DERIVED FIELDS.
247 !
248 ! ETA SURFACE TO 500MB LIFTED INDEX. TO BE CONSISTENT WITH THE
249 ! LFM AND NGM POSTING WE ADD 273.15 TO THE LIFTED INDEX
250 ! GSM WILL NOT ADD 273 TO VALUE FOR RAPID REFRESH TO BE
251 ! CONSISTENT WITH RUC
252 !
253 ! THE BEST (SIX LAYER) AND BOUNDARY LAYER LIFTED INDICES ARE
254 ! COMPUTED AND POSTED IN SUBROUTINE MISCLN.
255 !
256  IF (iget(030)>0.OR.iget(572)>0) THEN
257 !$omp parallel do private(i,j)
258  DO j=jsta,jend
259  DO i=ista,iend
260  egrid1(i,j) = spval
261  ENDDO
262  ENDDO
263 !
264  CALL otlift(egrid1)
265 !
266  IF(modelname == 'RAPR') THEN
267 !$omp parallel do private(i,j)
268  DO j=jsta,jend
269  DO i=ista,iend
270  IF(egrid1(i,j) < spval) grid1(i,j) = egrid1(i,j)
271  ENDDO
272  ENDDO
273  ELSE
274 !$omp parallel do private(i,j)
275  DO j=jsta,jend
276  DO i=ista,iend
277  IF(egrid1(i,j) < spval) grid1(i,j) = egrid1(i,j) + tfrz
278  ENDDO
279  ENDDO
280  ENDIF
281 !
282  if(iget(030) > 0) then
283  if(grib == "grib2" )then
284  cfld = cfld+1
285  fld_info(cfld)%ifld = iavblfld(iget(030))
286 !$omp parallel do private(i,j,ii,jj)
287  do j=1,jend-jsta+1
288  jj = jsta+j-1
289  do i=1,iend-ista+1
290  ii = ista+i-1
291  datapd(i,j,cfld) = grid1(ii,jj)
292  enddo
293  enddo
294  endif
295  endif
296 !for GFS
297  if(iget(572) > 0) then
298  if(grib == "grib2" )then
299  cfld = cfld+1
300  fld_info(cfld)%ifld = iavblfld(iget(572))
301 ! where(GRID1 /= SPVAL) GRID1 = GRID1-TFRZ
302 !$omp parallel do private(i,j,ii,jj)
303  do j=1,jend-jsta+1
304  jj = jsta+j-1
305  do i=1,iend-ista+1
306  ii = ista+i-1
307  if (grid1(ii,jj) /= spval) grid1(ii,jj) = grid1(ii,jj) - tfrz
308  datapd(i,j,cfld) = grid1(ii,jj)
309  enddo
310  enddo
311  endif
312  endif
313 
314  ENDIF
315 !
316 ! SOUNDING DERIVED AREA INTEGRATED ENERGIES - CAPE AND CIN.
317 ! THIS IS THE SFC-BASED CAPE/CIN (lowest 70 mb searched)
318 !
319 ! CONVECTIVE AVAILABLE POTENTIAL ENERGY.
320  IF ((iget(032) > 0))THEN
321 ! dong add missing value for cape
322  grid1 = spval
323  IF ( (lvls(1,iget(032))>0) )THEN
324  itype = 1
325  dpbnd = 10.e2
326  dummy = 0.
327  idummy = 0
328  CALL calcape(itype,dpbnd,dummy,dummy,dummy,idummy,egrid1,egrid2, &
329  egrid3,dummy,dummy)
330 !$omp parallel do private(i,j)
331  DO j=jsta,jend
332  DO i=ista,iend
333  IF(fis(i,j) < spval) grid1(i,j) = egrid1(i,j)
334  ENDDO
335  ENDDO
336  CALL bound(grid1,d00,h99999)
337  if(grib == "grib2" )then
338  cfld = cfld+1
339  fld_info(cfld)%ifld = iavblfld(iget(032))
340 !$omp parallel do private(i,j,ii,jj)
341  do j=1,jend-jsta+1
342  jj = jsta+j-1
343  do i=1,iend-ista+1
344  ii=ista+i-1
345  datapd(i,j,cfld) = grid1(ii,jj)
346  enddo
347  enddo
348  endif
349  END IF
350  END IF
351 !
352 ! CONVECTIVE INHIBITION.
353  IF ((iget(107) > 0))THEN
354 ! dong add missing value for cin
355  grid1 = spval
356  IF ( (lvls(1,iget(107)) > 0) )THEN
357  IF ((iget(032) > 0))THEN
358  IF ( (lvls(1,iget(032)) > 0) )THEN
359 !$omp parallel do private(i,j)
360  DO j=jsta,jend
361  DO i=ista,iend
362  IF(fis(i,j) < spval) grid1(i,j) = - egrid2(i,j)
363  ENDDO
364  ENDDO
365  END IF
366  ELSE
367  itype = 1
368  dpbnd = 10.e2
369  dummy = 0.
370  idummy = 0
371  CALL calcape(itype,dpbnd,dummy,dummy,dummy,idummy,egrid1,egrid2, &
372  egrid3,dummy,dummy)
373 !$omp parallel do private(i,j)
374  DO j=jsta,jend
375  DO i=ista,iend
376  IF(fis(i,j) < spval) grid1(i,j) = - egrid2(i,j)
377  ENDDO
378  ENDDO
379  END IF
380  CALL bound(grid1,d00,h99999)
381 !$omp parallel do private(i,j)
382  DO j=jsta,jend
383  DO i=ista,iend
384  IF(fis(i,j) < spval) grid1(i,j) = - grid1(i,j)
385  ENDDO
386  ENDDO
387  if(grib == "grib2" )then
388  cfld = cfld+1
389  fld_info(cfld)%ifld = iavblfld(iget(107))
390 !$omp parallel do private(i,j,ii,jj)
391  do j=1,jend-jsta+1
392  jj = jsta+j-1
393  do i=1,iend-ista+1
394  ii=ista+i-1
395  datapd(i,j,cfld) = grid1(ii,jj)
396  enddo
397  enddo
398  endif
399  END IF ! end for lvls(107)
400  END IF ! end of iget(107)
401 
402 !!!=======================================================================
403 !
404 ! TOTAL COLUMN PRECIPITABLE WATER (SPECIFIC HUMIDITY).
405  IF (iget(080) > 0) THEN
406 ! dong
407  grid1 = spval
408  model_pwat = .false.
409  DO j=jsta,jend
410  DO i=ista,iend
411  IF(abs(pwat(i,j)-spval)>small) THEN
412  model_pwat = .true.
413  exit
414  ENDIF
415  END DO
416  END DO
417  IF (model_pwat) THEN
418  DO j=jsta,jend
419  DO i=ista,iend
420  grid1(i,j) = pwat(i,j)
421  END DO
422  END DO
423  ELSE
424  CALL calpw(grid1(ista:iend,jsta:jend),1)
425  DO j=jsta,jend
426  DO i=ista,iend
427  IF(fis(i,j) >= spval) grid1(i,j)=spval
428  END DO
429  END DO
430  ENDIF
431  CALL bound(grid1,d00,h99999)
432  if(grib == "grib2" )then
433  cfld = cfld + 1
434  fld_info(cfld)%ifld = iavblfld(iget(080))
435 !$omp parallel do private(i,j,ii,jj)
436  do j=1,jend-jsta+1
437  jj = jsta+j-1
438  do i=1,iend-ista+1
439  ii=ista+i-1
440  datapd(i,j,cfld) = grid1(ii,jj)
441  enddo
442  enddo
443  endif
444  ENDIF
445 !
446 ! E. James - 8 Dec 2017
447 ! TOTAL COLUMN AOD (TAOD553D FROM HRRR-SMOKE)
448 !
449  IF (iget(735) > 0) THEN
450  IF (modelname == 'RAPR' .OR. modelname == 'FV3R') THEN
451  CALL calpw(grid1(ista:iend,jsta:jend),19)
452  CALL bound(grid1,d00,h99999)
453  ENDIF
454  if(grib == "grib2" )then
455  cfld = cfld + 1
456  fld_info(cfld)%ifld = iavblfld(iget(735))
457 !$omp parallel do private(i,j,ii,jj)
458  do j=1,jend-jsta+1
459  jj = jsta+j-1
460  do i=1,iend-ista+1
461  ii=ista+i-1
462  datapd(i,j,cfld) = grid1(ii,jj)
463  enddo
464  enddo
465  endif
466  ENDIF
467 !
468 ! E. James - 8 Dec 2017
469 ! TOTAL COLUMN FIRE SMOKE (tracer_1a FROM HRRR-SMOKE)
470 !
471  IF (iget(736) > 0) THEN
472  CALL calpw(grid1(ista:iend,jsta:iend),18)
473  CALL bound(grid1,d00,h99999)
474  if(grib == "grib2" )then
475  cfld = cfld + 1
476  fld_info(cfld)%ifld = iavblfld(iget(736))
477 !$omp parallel do private(i,j,ii,jj)
478  do j=1,jend-jsta+1
479  jj = jsta+j-1
480  do i=1,iend-ista+1
481  ii=ista+i-1
482  datapd(i,j,cfld) = grid1(ii,jj)
483  enddo
484  enddo
485  endif
486  ENDIF
487 !
488 ! TOTAL COLUMN DUST
489 !
490  IF (iget(741) > 0) THEN
491  CALL calpw(grid1(ista:iend,jsta:iend),22)
492  CALL bound(grid1,d00,h99999)
493  if(grib == "grib2" )then
494  cfld = cfld + 1
495  fld_info(cfld)%ifld = iavblfld(iget(741))
496 !$omp parallel do private(i,j,ii,jj)
497  do j=1,jend-jsta+1
498  jj = jsta+j-1
499  do i=1,iend-ista+1
500  ii=ista+i-1
501  datapd(i,j,cfld) = grid1(ii,jj)
502  enddo
503  enddo
504  endif
505  ENDIF
506 !
507 ! TOTAL COLUMN COARSEPM
508 !
509  IF (iget(1011) > 0) THEN
510  CALL calpw(grid1(ista:iend,jsta:iend),23)
511  CALL bound(grid1,d00,h99999)
512  if(grib == "grib2" )then
513  cfld = cfld + 1
514  fld_info(cfld)%ifld = iavblfld(iget(1011))
515 !$omp parallel do private(i,j,ii,jj)
516  do j=1,jend-jsta+1
517  jj = jsta+j-1
518  do i=1,iend-ista+1
519  ii=ista+i-1
520  datapd(i,j,cfld) = grid1(ii,jj)
521  enddo
522  enddo
523  endif
524  ENDIF
525 !
526 ! TOTAL COLUMN CLOUD WATER
527  IF (iget(200) > 0 .or. iget(575) > 0) THEN
528  grid1 = spval
529  grid2 = spval
530  IF (modelname == 'RAPR') THEN
531  DO j=jsta,jend
532  DO i=ista,iend
533  IF(lwp(i,j) < spval) grid1(i,j) = lwp(i,j)/1000.0 ! use WRF-diagnosed value
534  ENDDO
535  ENDDO
536  ELSE
537  CALL calpw(grid1(ista:iend,jsta:jend),2)
538  IF(modelname == 'GFS')then
539 ! GFS combines cloud water and cloud ice, hoping to seperate them next implementation
540  CALL calpw(grid2(ista:iend,jsta:jend),3)
541 !$omp parallel do private(i,j)
542  DO j=jsta,jend
543  DO i=ista,iend
544  IF(grid1(i,j)<spval.and.grid2(i,j)<spval)THEN
545  grid1(i,j) = grid1(i,j) + grid2(i,j)
546  ELSE
547  grid1(i,j) = spval
548  ENDIF
549  ENDDO
550  ENDDO
551  END IF ! GFS
552  END IF ! RAPR
553 
554  CALL bound(grid1,d00,h99999)
555  if(iget(200) > 0) then
556  if(grib == "grib2" )then
557  cfld = cfld + 1
558  fld_info(cfld)%ifld = iavblfld(iget(200))
559 !$omp parallel do private(i,j,ii,jj)
560  do j=1,jend-jsta+1
561  jj = jsta+j-1
562  do i=1,iend-ista+1
563  ii=ista+i-1
564  datapd(i,j,cfld) = grid1(ii,jj)
565  enddo
566  enddo
567  endif
568  endif
569  if(iget(575) > 0) then
570  if(grib == "grib2" )then
571  cfld = cfld + 1
572  fld_info(cfld)%ifld = iavblfld(iget(575))
573 !$omp parallel do private(i,j,ii,jj)
574  do j=1,jend-jsta+1
575  jj = jsta+j-1
576  do i=1,iend-ista+1
577  ii=ista+i-1
578  datapd(i,j,cfld) = grid1(ii,jj)
579  enddo
580  enddo
581  endif
582 
583  endif
584  ENDIF
585 !
586 ! TOTAL COLUMN CLOUD ICE
587  IF (iget(201) > 0) THEN
588  grid1 = spval
589  IF (modelname == 'RAPR') THEN
590  DO j=jsta,jend
591  DO i=ista,iend
592  IF(iwp(i,j) < spval) grid1(i,j) = iwp(i,j)/1000.0 ! use WRF-diagnosed value
593  ENDDO
594  ENDDO
595  ELSE
596  CALL calpw(grid1(ista:iend,jsta:jend),3)
597  END IF
598  CALL bound(grid1,d00,h99999)
599  if(grib == "grib2" )then
600  cfld = cfld + 1
601  fld_info(cfld)%ifld = iavblfld(iget(201))
602 !$omp parallel do private(i,j,ii,jj)
603  do j=1,jend-jsta+1
604  jj = jsta+j-1
605  do i=1,iend-ista+1
606  ii=ista+i-1
607  datapd(i,j,cfld) = grid1(ii,jj)
608  enddo
609  enddo
610  endif
611  ENDIF
612 !
613 ! TOTAL COLUMN RAIN
614  IF (iget(202) > 0) THEN
615  CALL calpw(grid1(ista:iend,jsta:jend),4)
616  CALL bound(grid1,d00,h99999)
617  if(grib=="grib2" )then
618  cfld=cfld+1
619  fld_info(cfld)%ifld=iavblfld(iget(202))
620 !$omp parallel do private(i,j,ii,jj)
621  do j=1,jend-jsta+1
622  jj = jsta+j-1
623  do i=1,iend-ista+1
624  ii=ista+i-1
625  datapd(i,j,cfld) = grid1(ii,jj)
626  enddo
627  enddo
628  endif
629  ENDIF
630 !
631 ! TOTAL COLUMN SNOW
632  IF (iget(203) > 0) THEN
633  CALL calpw(grid1(ista:iend,jsta:jend),5)
634  CALL bound(grid1,d00,h99999)
635  if(grib=="grib2" )then
636  cfld=cfld+1
637  fld_info(cfld)%ifld=iavblfld(iget(203))
638 !$omp parallel do private(i,j,ii,jj)
639  do j=1,jend-jsta+1
640  jj = jsta+j-1
641  do i=1,iend-ista+1
642  ii=ista+i-1
643  datapd(i,j,cfld) = grid1(ii,jj)
644  enddo
645  enddo
646  endif
647  ENDIF
648 !
649 ! SRD
650 ! TOTAL COLUMN GRAUPEL
651  IF (iget(428) > 0) THEN
652  CALL calpw(grid1(ista:iend,jsta:jend),16)
653  CALL bound(grid1,d00,h99999)
654  if(grib=="grib2" )then
655  cfld=cfld+1
656  fld_info(cfld)%ifld=iavblfld(iget(428))
657 !$omp parallel do private(i,j,ii,jj)
658  do j=1,jend-jsta+1
659  jj = jsta+j-1
660  do i=1,iend-ista+1
661  ii=ista+i-1
662  datapd(i,j,cfld) = grid1(ii,jj)
663  enddo
664  enddo
665  endif
666  ENDIF
667 ! SRD
668 
669 ! TOTAL COLUMN CONDENSATE
670  IF (iget(204) > 0) THEN
671  CALL calpw(grid1(ista:iend,jsta:jend),6)
672  CALL bound(grid1,d00,h99999)
673  if(grib=="grib2" )then
674  cfld=cfld+1
675  fld_info(cfld)%ifld=iavblfld(iget(204))
676 !$omp parallel do private(i,j,ii,jj)
677  do j=1,jend-jsta+1
678  jj = jsta+j-1
679  do i=1,iend-ista+1
680  ii=ista+i-1
681  datapd(i,j,cfld) = grid1(ii,jj)
682  enddo
683  enddo
684  endif
685  ENDIF
686 !
687 ! TOTAL COLUMN SUPERCOOLED (<0C) LIQUID WATER
688  IF (iget(285) > 0) THEN
689  CALL calpw(grid1(ista:iend,jsta:jend),7)
690  CALL bound(grid1,d00,h99999)
691  if(grib=="grib2" )then
692  cfld=cfld+1
693  fld_info(cfld)%ifld=iavblfld(iget(285))
694 !$omp parallel do private(i,j,ii,jj)
695  do j=1,jend-jsta+1
696  jj = jsta+j-1
697  do i=1,iend-ista+1
698  ii=ista+i-1
699  datapd(i,j,cfld) = grid1(ii,jj)
700  enddo
701  enddo
702  endif
703  ENDIF
704 !
705 ! TOTAL COLUMN MELTING (>0C) ICE
706  IF (iget(286) > 0) THEN
707  CALL calpw(grid1(ista:iend,jsta:jend),8)
708  CALL bound(grid1,d00,h99999)
709  if(grib=="grib2" )then
710  cfld=cfld+1
711  fld_info(cfld)%ifld=iavblfld(iget(286))
712 !$omp parallel do private(i,j,ii,jj)
713  do j=1,jend-jsta+1
714  jj = jsta+j-1
715  do i=1,iend-ista+1
716  ii=ista+i-1
717  datapd(i,j,cfld) = grid1(ii,jj)
718  enddo
719  enddo
720  endif
721  ENDIF
722 !
723 ! TOTAL COLUMN SHORT WAVE T TENDENCY
724  IF (iget(290) > 0) THEN
725  CALL calpw(grid1(ista:iend,jsta:jend),9)
726  if(grib=="grib2" )then
727  cfld=cfld+1
728  fld_info(cfld)%ifld=iavblfld(iget(290))
729 !$omp parallel do private(i,j,ii,jj)
730  do j=1,jend-jsta+1
731  jj = jsta+j-1
732  do i=1,iend-ista+1
733  ii=ista+i-1
734  datapd(i,j,cfld) = grid1(ii,jj)
735  enddo
736  enddo
737  endif
738  ENDIF
739 !
740 ! TOTAL COLUMN LONG WAVE T TENDENCY
741  IF (iget(291) > 0) THEN
742  CALL calpw(grid1(ista:iend,jsta:jend),10)
743  if(grib=="grib2" )then
744  cfld=cfld+1
745  fld_info(cfld)%ifld=iavblfld(iget(291))
746 !$omp parallel do private(i,j,ii,jj)
747  do j=1,jend-jsta+1
748  jj = jsta+j-1
749  do i=1,iend-ista+1
750  ii=ista+i-1
751  datapd(i,j,cfld) = grid1(ii,jj)
752  enddo
753  enddo
754  endif
755  ENDIF
756 !
757 ! TOTAL COLUMN GRID SCALE LATENT HEATING (TIME AVE)
758  IF (iget(292) > 0) THEN
759  CALL calpw(grid1(ista:iend,jsta:jend),11)
760  IF(avrain > 0.)THEN
761  rrnum = 1./avrain
762  ELSE
763  rrnum = 0.
764  ENDIF
765 !$omp parallel do private(i,j)
766  DO j=jsta,jend
767  DO i=ista,iend
768  IF(grid1(i,j) < spval) grid1(i,j) = grid1(i,j)*rrnum
769  ENDDO
770  ENDDO
771  id(1:25)=0
772  itheat = nint(theat)
773  IF (itheat /= 0) THEN
774  ifincr = mod(ifhr,itheat)
775  ELSE
776  ifincr=0
777  END IF
778  id(19) = ifhr
779  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
780  id(20) = 3
781  IF (ifincr==0) THEN
782  id(18) = ifhr-itheat
783  ELSE
784  id(18) = ifhr-ifincr
785  ENDIF
786  IF(ifmin >= 1)id(18)=id(18)*60
787  IF (id(18)<0) id(18) = 0
788  if(grib=="grib2" )then
789  cfld=cfld+1
790  fld_info(cfld)%ifld=iavblfld(iget(292))
791  if(itheat>0) then
792  fld_info(cfld)%ntrange=1
793  else
794  fld_info(cfld)%ntrange=0
795  endif
796  fld_info(cfld)%tinvstat=ifhr-id(18)
797 !$omp parallel do private(i,j,ii,jj)
798  do j=1,jend-jsta+1
799  jj = jsta+j-1
800  do i=1,iend-ista+1
801  ii=ista+i-1
802  datapd(i,j,cfld) = grid1(ii,jj)
803  enddo
804  enddo
805  endif
806  ENDIF
807 !
808 ! TOTAL COLUMN CONVECTIVE LATENT HEATING (TIME AVE)
809  IF (iget(293) > 0) THEN
810  CALL calpw(grid1(ista:iend,jsta:jend),12)
811  IF(avrain > 0.)THEN
812  rrnum = 1./avcnvc
813  ELSE
814  rrnum = 0.
815  ENDIF
816 !$omp parallel do private(i,j)
817  DO j=jsta,jend
818  DO i=ista,iend
819  IF(grid1(i,j) < spval) grid1(i,j) = grid1(i,j)*rrnum
820  ENDDO
821  ENDDO
822  id(1:25)=0
823  itheat = nint(theat)
824  IF (itheat /= 0) THEN
825  ifincr = mod(ifhr,itheat)
826  ELSE
827  ifincr=0
828  END IF
829  id(19) = ifhr
830  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
831  id(20) = 3
832  IF (ifincr==0) THEN
833  id(18) = ifhr-itheat
834  ELSE
835  id(18) = ifhr-ifincr
836  ENDIF
837  IF(ifmin >= 1)id(18)=id(18)*60
838  IF (id(18)<0) id(18) = 0
839  if(grib=="grib2" )then
840  cfld=cfld+1
841  fld_info(cfld)%ifld=iavblfld(iget(293))
842  if(itheat>0) then
843  fld_info(cfld)%ntrange=1
844  else
845  fld_info(cfld)%ntrange=0
846  endif
847  fld_info(cfld)%tinvstat=ifhr-id(18)
848 !$omp parallel do private(i,j,ii,jj)
849  do j=1,jend-jsta+1
850  jj = jsta+j-1
851  do i=1,iend-ista+1
852  ii=ista+i-1
853  datapd(i,j,cfld) = grid1(ii,jj)
854  enddo
855  enddo
856  endif
857  ENDIF
858 !
859 ! TOTAL COLUMN moisture convergence
860  IF (iget(295)>0) THEN
861  CALL calpw(grid1(ista:iend,jsta:jend),13)
862  if(grib=="grib2" )then
863  cfld=cfld+1
864  fld_info(cfld)%ifld=iavblfld(iget(295))
865  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
866  endif
867  ENDIF
868 !
869 ! TOTAL COLUMN RH
870  IF (iget(312)>0) THEN
871  CALL calpw(grid1(ista:iend,jsta:jend),14)
872  if(grib=="grib2" )then
873  cfld=cfld+1
874  fld_info(cfld)%ifld=iavblfld(iget(312))
875  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
876  endif
877  ENDIF
878 !
879 ! TOTAL COLUMN OZONE
880  IF (iget(299) > 0) THEN
881  CALL calpw(grid1(ista:iend,jsta:jend),15)
882  if(grib=="grib2" )then
883  cfld=cfld+1
884  fld_info(cfld)%ifld=iavblfld(iget(299))
885 !$omp parallel do private(i,j,ii,jj)
886  do j=1,jend-jsta+1
887  jj = jsta+j-1
888  do i=1,iend-ista+1
889  ii=ista+i-1
890  datapd(i,j,cfld) = grid1(ii,jj)
891  enddo
892  enddo
893  endif
894  ENDIF
895 !
896 ! BOTTOM AND/OR TOP OF SUPERCOOLED (<0C) LIQUID WATER LAYER
897  IF (iget(287)>0 .OR. iget(288)>0) THEN
898  DO j=jsta,jend
899  DO i=ista,iend
900  grid1(i,j)=-5000.
901  grid2(i,j)=-5000.
902 !-- Search for the base first, then look for the top if supercooled liquid exists
903  lbot=0
904  lm=nint(lmh(i,j))
905  DO l=lm,1,-1
906  qcld=qqw(i,j,l)+qqr(i,j,l)
907  IF (qcld>=qcldmin .AND. t(i,j,l)<tfrz) THEN
908  lbot=l
909  EXIT
910  ENDIF
911  ENDDO !--- End L loop
912  IF (lbot > 0) THEN
913 !-- Supercooled liquid exists, so get top & bottom heights. In this case,
914 ! be conservative and select the lower interface height at the bottom of the
915 ! layer and the top interface height at the top of the layer.
916  grid1(i,j)=zint(i,j,lbot+1)
917  DO l=1,lm
918  qcld=qqw(i,j,l)+qqr(i,j,l)
919  IF (qcld>=qcldmin .AND. t(i,j,l)<tfrz) THEN
920  ltop=l
921  EXIT
922  ENDIF
923  ENDDO !--- End L loop
924  ltop=min(lbot,ltop)
925  grid2(i,j)=zint(i,j,ltop)
926  ENDIF !--- End IF (LBOT > 0)
927  ENDDO !--- End I loop
928  ENDDO !--- End J loop
929  IF (iget(287)>0) THEN
930  if(grib=="grib2" )then
931  cfld=cfld+1
932  fld_info(cfld)%ifld=iavblfld(iget(287))
933  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
934  endif
935  ENDIF
936  IF (iget(288)>0) THEN
937 !$omp parallel do private(i,j)
938  DO j=jsta,jend
939  DO i=ista,iend
940  grid1(i,j)=grid2(i,j)
941  ENDDO
942  ENDDO
943  if(grib=="grib2" )then
944  cfld=cfld+1
945  fld_info(cfld)%ifld=iavblfld(iget(288))
946 !$omp parallel do private(i,j,ii,jj)
947  do j=1,jend-jsta+1
948  jj = jsta+j-1
949  do i=1,iend-ista+1
950  ii=ista+i-1
951  datapd(i,j,cfld) = grid1(ii,jj)
952  enddo
953  enddo
954  endif
955  ENDIF
956  ENDIF
957 !
958 !
959 ! Convective cloud efficiency parameter used in convection ranges
960 ! from 0.2 (EFIMN in cuparm in model) to 1.0 (Ferrier, Feb '02)
961  IF (iget(197)>0) THEN
962  DO j=jsta,jend
963  DO i=ista,iend
964  grid1(i,j) = cldefi(i,j)
965  ENDDO
966  ENDDO
967  if(grib=="grib2" )then
968  cfld=cfld+1
969  fld_info(cfld)%ifld=iavblfld(iget(197))
970  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
971  endif
972  ENDIF
973 !
974  IF ((modelname=='NMM' .AND. gridtype=='B') .OR. &
975  modelname=='FV3R') THEN
976 !nmmb_clds1
977 !
978 !-- Initialize low, middle, high, and total cloud cover;
979 ! also a method for cloud ceiling height
980 !
981  DO j=jsta,jend
982  DO i=ista,iend
983  cfracl(i,j)=0.
984  cfracm(i,j)=0.
985  cfrach(i,j)=0.
986  tcld(i,j)=0.
987  ENDDO
988  ENDDO
989 !
990 !-- Average cloud fractions over a 10 mi (16.09 km) radius (R),
991 ! approximated by a box of the same area = pi*R**2. Final
992 ! distance (d) is 1/2 of box size, d=0.5*sqrt(pi)*R=14259 m.
993 !
994  if(grib == "grib2" )then
995  dy_m=dyval*0.1112 !- DY_m in m
996  endif
997  dely=14259./dy_m
998  numr=nint(dely)
999  write (*,*) 'numr,dyval,DY_m=',numr,dyval,dy_m
1000  DO l=lm,1,-1
1001  DO j=jsta,jend
1002  DO i=ista,iend
1003  if(cfr(i,j,l)<spval) then
1004  full_cld(i,j)=cfr(i,j,l) !- 3D cloud fraction (from radiation)
1005  else
1006  full_cld(i,j)=spval
1007  endif
1008  ENDDO
1009  ENDDO
1010 ! CALL AllGETHERV(FULL_CLD)
1011  full_dummy=spval
1012  CALL collect_all(full_cld(ista:iend,jsta:jend),full_dummy)
1013  full_cld=full_dummy
1014  DO j=jsta,jend
1015  DO i=ista,iend
1016  numpts=0
1017  frac=0.
1018  DO jc=max(1,j-numr),min(jm,j+numr)
1019  DO ic=max(1,i-numr),min(im,i+numr)
1020 ! if(IC>=1.and.IC<=IM.and.JM>=JSTA.and.JM<=JEND) then
1021  IF(full_cld(ic,jc) /= spval) THEN
1022  numpts=numpts+1
1023  frac=frac+full_cld(ic,jc)
1024  ENDIF
1025 ! else
1026 ! FRAC=spval
1027 ! endif
1028  ENDDO
1029  ENDDO
1030  IF (numpts>0) frac=frac/REAL(numpts)
1031  if(pmid(i,j,l)<spval) then
1032  pcldbase=pmid(i,j,l) !-- Using PCLDBASE variable for convenience
1033  IF (pcldbase>=ptop_low) THEN
1034  cfracl(i,j)=max(cfracl(i,j),frac)
1035  ELSE IF (pcldbase>=ptop_mid) THEN
1036  cfracm(i,j)=max(cfracm(i,j),frac)
1037  ELSE
1038  cfrach(i,j)=max(cfrach(i,j),frac)
1039  ENDIF
1040  tcld(i,j)=max(tcld(i,j),frac)
1041  else
1042  cfracl(i,j)=spval
1043  cfracm(i,j)=spval
1044  cfrach(i,j)=spval
1045  tcld(i,j)=spval
1046  endif
1047  ENDDO ! I
1048  ENDDO ! J
1049  ENDDO ! L
1050 !end nmmb_clds1
1051  ELSEIF (modelname=='GFS') THEN
1052 !Initialize for GLOBAL FV3 which has cluod fraction in range from
1053 !0.0 to 1.0
1054 !
1055 !-- Initialize low, middle, high, and total cloud cover;
1056 ! also a method for cloud ceiling height
1057 !
1058  DO j=jsta,jend
1059  DO i=ista,iend
1060  cfracl(i,j)=0.
1061  cfracm(i,j)=0.
1062  cfrach(i,j)=0.
1063  tcld(i,j)=0.
1064  ENDDO
1065  ENDDO
1066  DO l=lm,1,-1
1067  DO j=jsta,jend
1068  DO i=ista,iend
1069  frac=cfr(i,j,l) !- 3D cloud fraction at model layers
1070  pcldbase=pmid(i,j,l) !-- Using PCLDBASE variable for convenience
1071  IF (pcldbase>=ptop_low) THEN
1072  cfracl(i,j)=max(cfracl(i,j),frac)
1073  ELSE IF (pcldbase>=ptop_mid) THEN
1074  cfracm(i,j)=max(cfracm(i,j),frac)
1075  ELSE
1076  cfrach(i,j)=max(cfrach(i,j),frac)
1077  ENDIF
1078  tcld(i,j)=max(tcld(i,j),frac)
1079  ENDDO ! I
1080  ENDDO ! J
1081  ENDDO ! L
1082  ENDIF
1083 !
1084 !*** BLOCK 2. 2-D CLOUD FIELDS.
1085 
1086 ! GSD maximum cloud fraction in (PBL + 1 km) (J. Kenyon, 8 Aug 2019)
1087  IF (iget(799)>0) THEN
1088 !$omp parallel do private(i,j,k)
1089  DO j=jsta,jend
1090  DO i=ista,iend
1091  grid1(i,j)=0.0
1092  DO k = 1,lm
1093  IF (zmid(i,j,lm-k+1) <= pblh(i,j)+1000.0) THEN
1094  grid1(i,j)=max(grid1(i,j),cfr(i,j,lm-k+1)*100.0)
1095  ENDIF
1096  ENDDO
1097  ENDDO
1098  ENDDO
1099  if(grib=="grib2" )then
1100  cfld=cfld+1
1101  fld_info(cfld)%ifld=iavblfld(iget(799))
1102  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1103  endif
1104  ENDIF
1105 !
1106 ! LOW CLOUD FRACTION.
1107  IF (iget(037) > 0) THEN
1108 !$omp parallel do private(i,j)
1109  DO j=jsta,jend
1110  DO i=ista,iend
1111  IF(cfracl(i,j) < spval) then
1112  grid1(i,j) = cfracl(i,j)*100.
1113  else
1114  grid1(i,j) = spval
1115  endif
1116  ENDDO
1117  ENDDO
1118  if(grib=="grib2" )then
1119  cfld=cfld+1
1120  fld_info(cfld)%ifld=iavblfld(iget(037))
1121 !$omp parallel do private(i,j,ii,jj)
1122  do j=1,jend-jsta+1
1123  jj = jsta+j-1
1124  do i=1,iend-ista+1
1125  ii=ista+i-1
1126  datapd(i,j,cfld) = grid1(ii,jj)
1127  enddo
1128  enddo
1129  endif
1130  ENDIF
1131 !
1132 ! TIME AVERAGED LOW CLOUD FRACTION.
1133  IF (iget(300) > 0) THEN
1134 !$omp parallel do private(i,j)
1135  DO j=jsta,jend
1136  DO i=ista,iend
1137  IF(avgcfracl(i,j) < spval) then
1138  grid1(i,j) = avgcfracl(i,j)*100.
1139  else
1140  grid1(i,j) = spval
1141  endif
1142  ENDDO
1143  ENDDO
1144  id(1:25)=0
1145  itclod = nint(tclod)
1146  IF(itclod /= 0) then
1147  ifincr = mod(ifhr,itclod)
1148  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
1149  ELSE
1150  ifincr = 0
1151  endif
1152 
1153  id(19) = ifhr
1154  IF(ifmin >= 1)id(19)=ifhr*60+ifmin !USE MIN FOR OFF-HR FORECAST
1155  id(20) = 3
1156  IF (ifincr==0) THEN
1157  id(18) = ifhr-itclod
1158  ELSE
1159  id(18) = ifhr-ifincr
1160  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1161  ENDIF
1162  IF (id(18)<0) id(18) = 0
1163  if(grib=="grib2" )then
1164  cfld=cfld+1
1165  fld_info(cfld)%ifld=iavblfld(iget(300))
1166  if(itclod>0) then
1167  fld_info(cfld)%ntrange=1
1168  else
1169  fld_info(cfld)%ntrange=0
1170  endif
1171  fld_info(cfld)%tinvstat=ifhr-id(18)
1172 !$omp parallel do private(i,j,ii,jj)
1173  do j=1,jend-jsta+1
1174  jj = jsta+j-1
1175  do i=1,iend-ista+1
1176  ii=ista+i-1
1177  datapd(i,j,cfld) = grid1(ii,jj)
1178  enddo
1179  enddo
1180  endif
1181  ENDIF
1182 !
1183 ! MIDDLE CLOUD FRACTION.
1184  IF (iget(038) > 0) THEN
1185 ! GRID1=SPVAL
1186 !$omp parallel do private(i,j)
1187  DO j=jsta,jend
1188  DO i=ista,iend
1189  IF(cfracm(i,j) < spval) then
1190  grid1(i,j) = cfracm(i,j)*100.
1191  else
1192  grid1(i,j) = spval
1193  endif
1194  ENDDO
1195  ENDDO
1196  if(grib=="grib2" )then
1197  cfld=cfld+1
1198  fld_info(cfld)%ifld=iavblfld(iget(038))
1199 !$omp parallel do private(i,j,ii,jj)
1200  do j=1,jend-jsta+1
1201  jj = jsta+j-1
1202  do i=1,iend-ista+1
1203  ii=ista+i-1
1204  datapd(i,j,cfld) = grid1(ii,jj)
1205  enddo
1206  enddo
1207  endif
1208  ENDIF
1209 !
1210 ! TIME AVERAGED MIDDLE CLOUD FRACTION.
1211  IF (iget(301) > 0) THEN
1212 !$omp parallel do private(i,j)
1213  DO j=jsta,jend
1214  DO i=ista,iend
1215  IF(abs(avgcfracm(i,j)-spval)>small)THEN
1216  grid1(i,j) = avgcfracm(i,j)*100.
1217  ELSE
1218  grid1(i,j) = spval
1219  END IF
1220  ENDDO
1221  ENDDO
1222  id(1:25)=0
1223  itclod = nint(tclod)
1224  IF(itclod /= 0) then
1225  ifincr = mod(ifhr,itclod)
1226  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
1227  ELSE
1228  ifincr = 0
1229  endif
1230 
1231  id(19) = ifhr
1232  IF(ifmin >= 1)id(19)=ifhr*60+ifmin !USE MIN FOR OFF-HR FORECAST
1233  id(20) = 3
1234  IF (ifincr==0) THEN
1235  id(18) = ifhr-itclod
1236  ELSE
1237  id(18) = ifhr-ifincr
1238  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1239  ENDIF
1240  IF (id(18)<0) id(18) = 0
1241  if(grib=="grib2" )then
1242  cfld=cfld+1
1243  fld_info(cfld)%ifld=iavblfld(iget(301))
1244  if(itclod>0) then
1245  fld_info(cfld)%ntrange=1
1246  else
1247  fld_info(cfld)%ntrange=0
1248  endif
1249  fld_info(cfld)%tinvstat=ifhr-id(18)
1250 !$omp parallel do private(i,j,ii,jj)
1251  do j=1,jend-jsta+1
1252  jj = jsta+j-1
1253  do i=1,iend-ista+1
1254  ii=ista+i-1
1255  datapd(i,j,cfld) = grid1(ii,jj)
1256  enddo
1257  enddo
1258  endif
1259  ENDIF
1260 !
1261 ! HIGH CLOUD FRACTION.
1262  IF (iget(039)>0) THEN
1263 ! GRID1=SPVAL
1264 !$omp parallel do private(i,j)
1265  DO j=jsta,jend
1266  DO i=ista,iend
1267  IF(cfrach(i,j) < spval) then
1268  grid1(i,j) = cfrach(i,j)*100.
1269  else
1270  grid1(i,j) = spval
1271  endif
1272  ENDDO
1273  ENDDO
1274  if(grib=="grib2" )then
1275  cfld=cfld+1
1276  fld_info(cfld)%ifld=iavblfld(iget(039))
1277 !$omp parallel do private(i,j,ii,jj)
1278  do j=1,jend-jsta+1
1279  jj = jsta+j-1
1280  do i=1,iend-ista+1
1281  ii = ista+i-1
1282  datapd(i,j,cfld) = grid1(ii,jj)
1283  enddo
1284  enddo
1285  endif
1286  ENDIF
1287 !
1288 ! TIME AVERAGED HIGH CLOUD FRACTION.
1289  IF (iget(302) > 0) THEN
1290 ! GRID1=SPVAL
1291 !$omp parallel do private(i,j)
1292  DO j=jsta,jend
1293  DO i=ista,iend
1294  IF(avgcfrach(i,j) < spval) then
1295  grid1(i,j) = avgcfrach(i,j)*100.
1296  else
1297  grid1(i,j) = spval
1298  endif
1299  ENDDO
1300  ENDDO
1301  id(1:25)=0
1302  itclod = nint(tclod)
1303  IF(itclod /= 0) then
1304  ifincr = mod(ifhr,itclod)
1305  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
1306  ELSE
1307  ifincr = 0
1308  endif
1309 
1310  id(19) = ifhr
1311  IF(ifmin >= 1)id(19)=ifhr*60+ifmin !USE MIN FOR OFF-HR FORECAST
1312  id(20) = 3
1313  IF (ifincr==0) THEN
1314  id(18) = ifhr-itclod
1315  ELSE
1316  id(18) = ifhr-ifincr
1317  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1318  ENDIF
1319  IF (id(18)<0) id(18) = 0
1320  if(grib=="grib2" )then
1321  cfld=cfld+1
1322  fld_info(cfld)%ifld=iavblfld(iget(302))
1323  if(itclod>0) then
1324  fld_info(cfld)%ntrange=1
1325  else
1326  fld_info(cfld)%ntrange=0
1327  endif
1328  fld_info(cfld)%tinvstat=ifhr-id(18)
1329 !$omp parallel do private(i,j,ii,jj)
1330  do j=1,jend-jsta+1
1331  jj = jsta+j-1
1332  do i=1,iend-ista+1
1333  ii = ista+i-1
1334  datapd(i,j,cfld) = grid1(ii,jj)
1335  enddo
1336  enddo
1337  endif
1338  ENDIF
1339 !
1340 ! TOTAL CLOUD FRACTION (INSTANTANEOUS).
1341  IF ((iget(161) > 0) .OR. (iget(260) > 0)) THEN
1342 ! GRID1=SPVAL
1343  IF(modelname=='NCAR' .OR. modelname=='RAPR')THEN
1344 !$omp parallel do private(i,j)
1345  DO j=jsta,jend
1346  DO i=ista,iend
1347  grid1(i,j) = spval
1348  egrid1(i,j)=0.
1349  do l = 1,lm
1350  egrid1(i,j)=max(egrid1(i,j),cfr(i,j,l))
1351  end do
1352  ENDDO
1353  ENDDO
1354 
1355  ELSE IF (modelname=='NMM'.OR.modelname=='FV3R' &
1356  .OR. modelname=='GFS')THEN
1357  DO j=jsta,jend
1358  DO i=ista,iend
1359 ! EGRID1(I,J)=AMAX1(CFRACL(I,J),
1360 ! 1 AMAX1(CFRACM(I,J),CFRACH(I,J)))
1361 ! EGRID1(I,J)=1.-(1.-CFRACL(I,J))*(1.-CFRACM(I,J))* &
1362 ! & (1.-CFRACH(I,J))
1363  grid1(i,j)=spval
1364  egrid1(i,j)=tcld(i,j)
1365  ENDDO
1366  ENDDO
1367  END IF
1368 !$omp parallel do private(i,j)
1369  DO j=jsta,jend
1370  DO i=ista,iend
1371  IF(abs(egrid1(i,j)-spval) > small) THEN
1372  grid1(i,j) = egrid1(i,j)*100.
1373  tcld(i,j) = egrid1(i,j)*100. !B ZHOU, PASSED to CALCEILING
1374  END IF
1375  ENDDO
1376  ENDDO
1377  IF (iget(161)>0) THEN
1378  if(grib=="grib2" )then
1379  cfld=cfld+1
1380  fld_info(cfld)%ifld=iavblfld(iget(161))
1381 !$omp parallel do private(i,j,ii,jj)
1382  do j=1,jend-jsta+1
1383  jj = jsta+j-1
1384  do i=1,iend-ista+1
1385  ii = ista+i-1
1386  datapd(i,j,cfld) = grid1(ii,jj)
1387  enddo
1388  enddo
1389  endif
1390  ENDIF
1391  ENDIF
1392 !
1393 ! TIME AVERAGED TOTAL CLOUD FRACTION.
1394  IF (iget(144) > 0) THEN
1395 ! GRID1=SPVAL
1396  IF(modelname == 'GFS' .OR. modelname == 'FV3R')THEN
1397 !$omp parallel do private(i,j)
1398  DO j=jsta,jend
1399  DO i=ista,iend
1400  IF(abs(avgtcdc(i,j)-spval) > small) then
1401  grid1(i,j) = avgtcdc(i,j)*100.
1402  else
1403  grid1(i,j) = spval
1404  endif
1405  END DO
1406  END DO
1407 
1408  ELSE IF(modelname == 'NMM')THEN
1409  DO j=jsta,jend
1410  DO i=ista,iend
1411 ! RSUM = NCFRST(I,J)+NCFRCV(I,J)
1412 ! IF (RSUM>0.0) THEN
1413 ! EGRID1(I,J)=(ACFRST(I,J)+ACFRCV(I,J))/RSUM
1414 ! ELSE
1415 ! EGRID1(I,J) = D00
1416 ! ENDIF
1417 !ADDED BRAD'S MODIFICATION
1418  rsum = d00
1419  IF (ncfrst(i,j)<spval.and.acfrst(i,j)<spval)THEN
1420  IF (ncfrst(i,j) > 0) rsum=acfrst(i,j)/ncfrst(i,j)
1421  IF (ncfrcv(i,j) > 0) &
1422  rsum=max(rsum, acfrcv(i,j)/ncfrcv(i,j))
1423  grid1(i,j) = rsum*100.
1424  ELSE
1425  grid1(i,j) = spval
1426  ENDIF
1427  ENDDO
1428  ENDDO
1429  END IF
1430  IF(modelname == 'NMM' .OR. modelname == 'GFS' .OR. &
1431  modelname == 'FV3R')THEN
1432  id(1:25)= 0
1433  itclod = nint(tclod)
1434  IF(itclod /= 0) then
1435  ifincr = mod(ifhr,itclod)
1436  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
1437  ELSE
1438  ifincr = 0
1439  endif
1440 
1441  id(19) = ifhr
1442  IF(ifmin >= 1)id(19)=ifhr*60+ifmin !USE MIN FOR OFF-HR FORECAST
1443  id(20) = 3
1444  IF (ifincr==0) THEN
1445  id(18) = ifhr-itclod
1446  ELSE
1447  id(18) = ifhr-ifincr
1448  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1449  ENDIF
1450  IF (id(18)<0) id(18) = 0
1451  ENDIF
1452  if(grib=="grib2" )then
1453  cfld=cfld+1
1454  fld_info(cfld)%ifld=iavblfld(iget(144))
1455  if(itclod>0) then
1456  fld_info(cfld)%ntrange=1
1457  else
1458  fld_info(cfld)%ntrange=0
1459  endif
1460  fld_info(cfld)%tinvstat=ifhr-id(18)
1461 !$omp parallel do private(i,j,ii,jj)
1462  do j=1,jend-jsta+1
1463  jj = jsta+j-1
1464  do i=1,iend-ista+1
1465  ii = ista+i-1
1466  datapd(i,j,cfld) = grid1(ii,jj)
1467  enddo
1468  enddo
1469  endif
1470  ENDIF
1471 !
1472 ! TIME AVERAGED STRATIFORM CLOUD FRACTION.
1473  IF (iget(139)>0) THEN
1474  IF(modelname /= 'NMM')THEN
1475  grid1=spval
1476  ELSE
1477  DO j=jsta,jend
1478  DO i=ista,iend
1479  IF (ncfrst(i,j)<spval.and.acfrst(i,j)<spval)THEN
1480  IF (ncfrst(i,j)>0.0) THEN
1481  grid1(i,j) = acfrst(i,j)/ncfrst(i,j)*100.
1482  ELSE
1483  grid1(i,j) = d00
1484  ENDIF
1485  ELSE
1486  grid1(i,j) = spval
1487  ENDIF
1488  ENDDO
1489  ENDDO
1490  END IF
1491  IF(modelname=='NMM' .or. modelname=='FV3R')THEN
1492  id(1:25)=0
1493  itclod = nint(tclod)
1494  IF(itclod /= 0) then
1495  ifincr = mod(ifhr,itclod)
1496  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
1497  ELSE
1498  ifincr = 0
1499  endif
1500  id(19) = ifhr
1501  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1502  id(20) = 3
1503  IF (ifincr==0) THEN
1504  id(18) = ifhr-itclod
1505  ELSE
1506  id(18) = ifhr-ifincr
1507  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1508  ENDIF
1509  IF (id(18)<0) id(18) = 0
1510  ENDIF
1511  if(grib=="grib2" )then
1512  cfld=cfld+1
1513  fld_info(cfld)%ifld=iavblfld(iget(139))
1514  if(itclod>0) then
1515  fld_info(cfld)%ntrange=1
1516  else
1517  fld_info(cfld)%ntrange=0
1518  endif
1519  fld_info(cfld)%tinvstat=ifhr-id(18)
1520  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1521  endif
1522  ENDIF
1523 !
1524 ! TIME AVERAGED CONVECTIVE CLOUD FRACTION.
1525  IF (iget(143)>0) THEN
1526  IF(modelname /= 'NMM')THEN
1527  grid1=spval
1528  ELSE
1529  DO j=jsta,jend
1530  DO i=ista,iend
1531  IF (ncfrcv(i,j)<spval.and.acfrcv(i,j)<spval)THEN
1532  IF (ncfrcv(i,j)>0.0) THEN
1533  grid1(i,j) = acfrcv(i,j)/ncfrcv(i,j)*100.
1534  ELSE
1535  grid1(i,j) = d00
1536  ENDIF
1537  ELSE
1538  grid1(i,j) = spval
1539  ENDIF
1540  ENDDO
1541  ENDDO
1542  END IF
1543  IF(modelname=='NMM')THEN
1544  id(1:25)=0
1545  itclod = nint(tclod)
1546  IF(itclod /= 0) then
1547  ifincr = mod(ifhr,itclod)
1548  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
1549  ELSE
1550  ifincr = 0
1551  endif
1552  id(19) = ifhr
1553  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
1554  id(20) = 3
1555  IF (ifincr==0) THEN
1556  id(18) = ifhr-itclod
1557  ELSE
1558  id(18) = ifhr-ifincr
1559  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
1560  ENDIF
1561  IF (id(18)<0) id(18) = 0
1562  ENDIF
1563  if(grib=="grib2" )then
1564  cfld=cfld+1
1565  fld_info(cfld)%ifld=iavblfld(iget(143))
1566  if(itclod>0) then
1567  fld_info(cfld)%ntrange=1
1568  else
1569  fld_info(cfld)%ntrange=0
1570  endif
1571  fld_info(cfld)%tinvstat=ifhr-id(18)
1572  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1573  endif
1574  ENDIF
1575 !
1576 ! CLOUD BASE AND TOP FIELDS
1577  IF((iget(148)>0) .OR. (iget(149)>0) .OR. &
1578  (iget(168)>0) .OR. (iget(178)>0) .OR. &
1579  (iget(179)>0) .OR. (iget(194)>0) .OR. &
1580  (iget(408)>0) .OR. &
1581  (iget(409)>0) .OR. (iget(406)>0) .OR. &
1582  (iget(195)>0) .OR. (iget(260)>0) .OR. &
1583  (iget(275)>0)) THEN
1584 !
1585 !--- Calculate grid-scale cloud base & top arrays (Ferrier, Feb '02)
1586 !
1587 !--- Rain is not part of cloud, only cloud water + cloud ice + snow
1588 !
1589  DO j=jsta,jend
1590  DO i=ista,iend
1591 !
1592 !--- Various convective cloud base & cloud top levels
1593 !
1594 ! write(*,*)' hbot=',hbot(i,j),' hbotd=',hbotd(i,j),'
1595 ! hbots=',hbots(i,j)&
1596 ! ,' htop=',htop(i,j),' htopd=',htopd(i,j),' htops=',htops(i,j),i,j
1597 ! Initilize
1598  ibotcu(i,j) = 0
1599  itopcu(i,j) = 100
1600  ibotdcu(i,j) = 0
1601  itopdcu(i,j) = 100
1602  ibotscu(i,j) = 0
1603  itopscu(i,j) = 100
1604  if (hbot(i,j) /= spval) then
1605  ibotcu(i,j) = nint(hbot(i,j))
1606  endif
1607  if (hbotd(i,j) /= spval) then
1608  ibotdcu(i,j) = nint(hbotd(i,j))
1609  endif
1610  if (hbots(i,j) /= spval) then
1611  ibotscu(i,j) = nint(hbots(i,j))
1612  endif
1613  if (htop(i,j) /= spval) then
1614  itopcu(i,j) = nint(htop(i,j))
1615  endif
1616  if (htopd(i,j) /= spval) then
1617  itopdcu(i,j) = nint(htopd(i,j))
1618  endif
1619  if (htops(i,j) /= spval) then
1620  itopscu(i,j) = nint(htops(i,j))
1621  endif
1622  IF (ibotcu(i,j)-itopcu(i,j) <= 1) THEN
1623  ibotcu(i,j) = 0
1624  itopcu(i,j) = 100
1625  ENDIF
1626  IF (ibotdcu(i,j)-itopdcu(i,j) <= 1) THEN
1627  ibotdcu(i,j) = 0
1628  itopdcu(i,j) = 100
1629  ENDIF
1630  IF (ibotscu(i,j)-itopscu(i,j) <= 1) THEN
1631  ibotscu(i,j) = 0
1632  itopscu(i,j) = 100
1633  ENDIF
1634 ! Convective cloud top height
1635  itop = itopcu(i,j)
1636  IF (itop > 0 .AND. itop < 100) THEN
1637 ! print *, 'aha ', ITOP
1638  ENDIF
1639  IF (itop > 0 .AND. itop <= nint(lmh(i,j))) THEN
1640  cldzcu(i,j) = zmid(i,j,itop)
1641  else
1642  cldzcu(i,j) = -5000.
1643  endif
1644 
1645 ! !
1646  !--- Grid-scale cloud base & cloud top levels
1647  !
1648  !--- Grid-scale cloud occurs when the mixing ratio exceeds QCLDmin
1649  ! or in the presence of snow when RH>=95% or at/above the PBL top.
1650  !
1651  if(modelname == 'RAPR') then
1652  ibotgr(i,j)=0
1653  DO l=nint(lmh(i,j)),1,-1
1654  qcld=qqw(i,j,l)+qqi(i,j,l)+qqs(i,j,l)
1655  IF (qcld >= qcldmin) THEN
1656  ibotgr(i,j)=l
1657  EXIT
1658  ENDIF
1659  ENDDO !--- End L loop
1660  itopgr(i,j)=100
1661  DO l=1,nint(lmh(i,j))
1662  qcld=qqw(i,j,l)+qqi(i,j,l)+qqs(i,j,l)
1663  IF (qcld >= qcldmin) THEN
1664  itopgr(i,j)=l
1665  EXIT
1666  ENDIF
1667  ENDDO !--- End L loop
1668  else
1669  ibotgr(i,j) = 0
1670  zpbltop = pblh(i,j)+zint(i,j,nint(lmh(i,j))+1)
1671  DO l=nint(lmh(i,j)),1,-1
1672  qcld = qqw(i,j,l)+qqi(i,j,l) !- no snow +QQS(I,J,L)
1673  IF (qcld >= qcldmin) THEN
1674  ibotgr(i,j) = l
1675  EXIT
1676  ENDIF
1677 snow_check: IF (qqs(i,j,l)>=qcldmin) THEN
1678  tmp=t(i,j,l)
1679  IF (tmp>=c2k) THEN
1680  qsat=pq0/pmid(i,j,l)*exp(a2*(tmp-a3)/(tmp-a4))
1681  ELSE
1682 !-- Use Teten's formula for ice from Murray (1967). More info at
1683 ! http://faculty.eas.ualberta.ca/jdwilson/EAS372_13/Vomel_CIRES_satvpformulae.html
1684  qsat=pq0/pmid(i,j,l)*exp(21.8745584*(tmp-a3)/(tmp-7.66))
1685  ENDIF
1686  rhum=q(i,j,l)/qsat
1687  IF (rhum>=0.98 .AND. zmid(i,j,l)>=zpbltop) THEN
1688  ibotgr(i,j)=l
1689  EXIT
1690  ENDIF
1691  ENDIF snow_check
1692  ENDDO !--- End L loop
1693  itopgr(i,j) = 100
1694  DO l=1,nint(lmh(i,j))
1695  qcld=qqw(i,j,l)+qqi(i,j,l)+qqs(i,j,l)
1696  IF (qcld >= qcldmin) THEN
1697  itopgr(i,j)=l
1698  EXIT
1699  ENDIF
1700  ENDDO !--- End L loop
1701  endif
1702  !
1703  !--- Combined (convective & grid-scale) cloud base & cloud top levels
1704  IF(modelname == 'NCAR' .OR. modelname == 'RAPR')THEN
1705  ibott(i,j) = ibotgr(i,j)
1706  itopt(i,j) = itopgr(i,j)
1707  ELSE
1708  ibott(i,j) = max(ibotgr(i,j), ibotcu(i,j))
1709 ! if(i==200 .and. j==139)print*,'Debug cloud base 1: ',&
1710 ! IBOTGr(I,J),IBOTCu(I,J),ibott(i,j)
1711  itopt(i,j) = min(itopgr(i,j), itopcu(i,j))
1712  END IF
1713  ENDDO !--- End I loop
1714  ENDDO !--- End J loop
1715  ENDIF !--- End IF tests
1716 !
1717 ! CONVECTIVE CLOUD TOP HEIGHT
1718  IF (iget(758)>0) THEN
1719 
1720  DO j=jsta,jend
1721  DO i=ista,iend
1722  grid1(i,j) = cldzcu(i,j)
1723  ENDDO
1724  ENDDO
1725  if(grib=="grib2" )then
1726  cfld=cfld+1
1727  fld_info(cfld)%ifld=iavblfld(iget(758))
1728  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1729  endif
1730  ENDIF
1731 !
1732 !-------------------------------------------------
1733 !----------- VARIOUS CLOUD BASE FIELDS ----------
1734 !-------------------------------------------------
1735 !
1736 !--- "TOTAL" CLOUD BASE FIELDS (convective + grid-scale; Ferrier, Feb '02)
1737 !
1738  IF ((iget(148)>0) .OR. (iget(178)>0) .OR.(iget(260)>0) ) THEN
1739  DO j=jsta,jend
1740  DO i=ista,iend
1741  ibot=ibott(i,j) !-- Cloud base ("bottoms")
1742  IF(modelname == 'RAPR') then
1743  IF (ibot <= 0) THEN
1744  cldp(i,j) = spval
1745  cldz(i,j) = spval
1746  ELSE IF (ibot <= nint(lmh(i,j))) THEN
1747  cldp(i,j) = pmid(i,j,ibot)
1748  IF (ibot == lm) THEN
1749  cldz(i,j) = zint(i,j,lm)
1750  ELSE
1751  cldz(i,j) = htm(i,j,ibot+1)*t(i,j,ibot+1) &
1752  *(q(i,j,ibot+1)*d608+h1)*rog* &
1753  (log(pint(i,j,ibot+1))-log(cldp(i,j)))&
1754  +zint(i,j,ibot+1)
1755  ENDIF !--- End IF (IBOT == LM) ...
1756  ENDIF !--- End IF (IBOT <= 0) ...
1757  ELSE
1758  IF (ibot>0 .AND. ibot<=nint(lmh(i,j))) THEN
1759  cldp(i,j) = pmid(i,j,ibot)
1760  cldz(i,j) = zmid(i,j,ibot)
1761  ELSE
1762  cldp(i,j) = -50000.
1763  cldz(i,j) = -5000.
1764  ENDIF !--- End IF (IBOT <= 0) ...
1765  ENDIF
1766  ENDDO !--- End DO I loop
1767  ENDDO !--- End DO J loop
1768 ! CLOUD BOTTOM PRESSURE
1769  IF (iget(148)>0) THEN
1770  DO j=jsta,jend
1771  DO i=ista,iend
1772  grid1(i,j) = cldp(i,j)
1773  ENDDO
1774  ENDDO
1775  if(grib=="grib2" )then
1776  cfld=cfld+1
1777  fld_info(cfld)%ifld=iavblfld(iget(148))
1778  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1779  endif
1780  ENDIF
1781 ! CLOUD BOTTOM HEIGHT
1782  IF (iget(178)>0) THEN
1783 !--- Parameter was set to 148 in operational code (Ferrier, Feb '02)
1784  DO j=jsta,jend
1785  DO i=ista,iend
1786  grid1(i,j) = cldz(i,j)
1787  ENDDO
1788  ENDDO
1789  if(grib=="grib2" )then
1790  cfld=cfld+1
1791  fld_info(cfld)%ifld=iavblfld(iget(178))
1792  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1793  endif
1794  ENDIF
1795  ENDIF
1796 
1797 ! GSD CLOUD CEILING ALGORITHMS...
1798 
1799 ! Parameter 408: legacy ceiling diagnostic
1800  IF (iget(408)>0) THEN
1801 !- imported from RUC post
1802 ! -- constants for effect of snow on ceiling
1803 ! Also found in calvis.f
1804  rhoice = 970.
1805  coeffp = 10.36
1806 ! - new value from Roy Rasmussen - Dec 2003
1807 ! exponfp = 0.7776
1808 ! change consistent with CALVIS_GSD.f
1809  exponfp = 1.
1810  const1 = 3.912
1811 
1812  nfog = 0
1813  do k=1,7
1814  nfogn(k) = 0
1815  end do
1816  npblcld = 0
1817 
1818  cloud_def_p = 0.0000001
1819 
1820  DO j=jsta,jend
1821  DO i=ista,iend
1822 !
1823 !- imported from RUC post
1824  cldz(i,j) = spval
1825  pcldbase = spval
1826  zcldbase = spval
1827  watericemax = -99999.
1828  do k=1,lm
1829  ll=lm-k+1
1830  watericetotal(k) = qqw(i,j,ll) + qqi(i,j,ll)
1831  watericemax = max(watericemax,watericetotal(k))
1832  end do
1833 
1834  if (watericemax>=cloud_def_p) then
1835 
1836 ! Cloud base
1837 !====================
1838 
1839 ! --- Check out no. of points with thin cloud layers near surface
1840  do k=2,3
1841  pabovesfc(k) = pint(i,j,lm) - pint(i,j,lm-k+1)
1842  if (watericetotal(k)<cloud_def_p) then
1843 ! --- wimin is watericemin in lowest few levels
1844  wimin = 100.
1845  do k1=k-1,1,-1
1846  wimin = min(wimin,watericetotal(k1))
1847  end do
1848  if (wimin>cloud_def_p) then
1849  nfogn(k)= nfogn(k)+1
1850  end if
1851  end if
1852  end do
1853 
1854 ! Eliminate fog layers near surface in watericetotal array
1855  loop1778 : do k=2,3
1856 ! --- Do this only when at least 10 mb (1000 Pa) above surface
1857 ! if (pabovesfc(k)>1000.) then
1858  if (watericetotal(k)<cloud_def_p) then
1859  if (watericetotal(1)>cloud_def_p) then
1860  nfog = nfog+1
1861  do k1=1,k-1
1862  if (watericetotal(k1)>=cloud_def_p) then
1863  watericetotal(k1)=0.
1864  end if
1865  end do
1866  end if
1867  end if
1868  exit loop1778
1869 ! end if
1870  end do loop1778
1871 
1872 !! At surface?
1873 !commented out 16aug11
1874 ! if (watericetotal(1)>cloud_def_p) then
1875 ! zcldbase = zmid(i,j,lm)
1876 ! go to 3788
1877 ! end if
1878 !! Aloft?
1879  loop371: do k=2,lm
1880  k1 = k
1881  if (watericetotal(k)>cloud_def_p) then
1882  if (k1<=4) then
1883 ! -- If within 4 levels of surface, just use lowest cloud level
1884 ! as ceiling WITHOUT vertical interpolation.
1885  zcldbase = zmid(i,j,lm-k1+1)
1886  pcldbase = pmid(i,j,lm-k1+1)
1887  else
1888 ! -- Use vertical interpolation to obtain cloud level
1889  zcldbase = zmid(i,j,lm-k1+1) + (cloud_def_p-watericetotal(k1)) &
1890  * (zmid(i,j,lm-k1+2)-zmid(i,j,lm-k1+1)) &
1891  / (watericetotal(k1-1) - watericetotal(k1))
1892  pcldbase = pmid(i,j,lm-k1+1) + (cloud_def_p-watericetotal(k1)) &
1893  * (pmid(i,j,lm-k1+2)-pmid(i,j,lm-k1+1)) &
1894  / (watericetotal(k1-1) - watericetotal(k1))
1895  end if
1896  zcldbase = max(zcldbase,fis(i,j)*gi+5.)
1897 
1898 ! 3788 continue
1899 
1900 ! -- consider lowering of ceiling due to falling snow
1901 ! -- extracted from calvis.f (visibility diagnostic)
1902  if (qqs(i,j,lm)>0.) then
1903  tv=t(i,j,lm)*(h1+d608*q(i,j,lm))
1904  rhoair=pmid(i,j,lm)/(rd*tv)
1905  vovermd = (1.+q(i,j,lm))/rhoair + qqs(i,j,lm)/rhoice
1906  concfp = qqs(i,j,lm)/vovermd*1000.
1907  betav = coeffp*concfp**exponfp + 1.e-10
1908  vertvis = 1000.*min(90., const1/betav)
1909  if (vertvis < zcldbase-fis(i,j)*gi ) then
1910  zcldbase = fis(i,j)*gi + vertvis
1911  loop3741: do k2=2,lm
1912  k1 = k2
1913  if (zmid(i,j,lm-k2+1) > zcldbase) then
1914  pcldbase = pmid(i,j,lm-k1+2) + (zcldbase-zmid(i,j,lm-k1+2)) &
1915  *(pmid(i,j,lm-k1+1)-pmid(i,j,lm-k1+2) ) &
1916  /(zmid(i,j,lm-k1+1)-zmid(i,j,lm-k1+2) )
1917  exit loop3741
1918  endif
1919  end do loop3741
1920  end if
1921  end if
1922  exit loop371
1923  endif
1924  end do loop371
1925  endif
1926 
1927 !new 15 aug 2011
1928  cldz(i,j) = zcldbase
1929  cldp(i,j) = pcldbase
1930 
1931 ! --- Now, do a PBL cloud check.
1932 ! --- First, get a PBL-top cloud ceiling, if it exists.
1933 ! This value is the first level under the cloud top if
1934 ! the RH is greater than 95%. This should help to identify
1935 ! ceilings that the RUC model doesn't quite catch due to
1936 ! vertical resolution.
1937 
1938 ! - compute relative humidity
1939  do k=1,lm
1940  ll=lm-k+1
1941  tx=t(i,j,ll)-273.15
1942  pol = 0.99999683 + tx*(-0.90826951e-02 + &
1943  tx*(0.78736169e-04 + tx*(-0.61117958e-06 + &
1944  tx*(0.43884187e-08 + tx*(-0.29883885e-10 + &
1945  tx*(0.21874425e-12 + tx*(-0.17892321e-14 + &
1946  tx*(0.11112018e-16 + tx*(-0.30994571e-19)))))))))
1947  esx = 6.1078/pol**8
1948 
1949  es = esx
1950  e = pmid(i,j,ll)/100.*q(i,j,ll)/(0.62197+q(i,j,ll)*0.37803)
1951  rhb(k) = 100.*min(1.,e/es)
1952 !
1953 ! COMPUTE VIRTUAL POTENTIAL TEMPERATURE.
1954 !
1955  enddo
1956 
1957 ! PBL height is computed in INITPOST.f
1958 ! zpbltop is relative to sea level
1959  zsf=zint(i,j,nint(lmh(i,j))+1)
1960  zpbltop = pblh(i,j)+zsf
1961 
1962 ! PBLH(I,J)= zpbltop - FIS(I,J)*GI
1963 ! print *,'I,J,k1,zmid(i,j,lm-k1+1),zmid(i,j,lm-k1),PBLH(I,J)',
1964 ! 1 I,J,k1,zmid(i,j,lm-k1+1),zmid(i,j,lm-k1),PBLH(I,J),RHB(k1)
1965 
1966  loop745: do k2=3,20
1967  if (zpbltop<zmid(i,j,lm-k2+1)) then
1968  if (rhb(k2-1)>95. ) then
1969  zcldbase = zmid(i,j,lm-k2+2)
1970  if (cldz(i,j)<-100.) then
1971  npblcld = npblcld+1
1972  cldz(i,j) = zcldbase
1973  cldp(i,j) = pmid(i,j,lm-k2+2)
1974  exit loop745
1975  end if
1976  if ( zcldbase<cldz(i,j)) then
1977  cldz(i,j) = zcldbase
1978  end if
1979  end if
1980  exit loop745
1981  end if
1982  end do loop745
1983 
1984 !- include convective clouds
1985  ibot=ibotcu(i,j)
1986  if(ibot>0) then
1987  if(cldz(i,j)<-100.) then
1988  cldz(i,j)=zmid(i,j,ibot)
1989  else
1990  if(zmid(i,j,ibot)<cldz(i,j)) then
1991  cldz(i,j)=zmid(i,j,ibot)
1992  endif
1993  endif
1994  endif
1995 
1996  ENDDO !--- End I loop
1997  ENDDO !--- End J loop
1998 
1999  write(6,*)'No. pts with PBL-cloud =',npblcld
2000  write(6,*)'No. pts to eliminate fog =',nfog
2001  do k=2,7
2002  write(6,*)'No. pts with fog below lev',k,' =',nfogn(k)
2003  end do
2004 
2005  nlifr = 0
2006  DO j=jsta,jend
2007  DO i=ista,iend
2008  zcld = cldz(i,j) - fis(i,j)*gi
2009  if (cldz(i,j)>=0..and.zcld<160.) nlifr = nlifr+1
2010  end do
2011  end do
2012  write(6,*)'No. pts w/ LIFR ceiling =',nlifr
2013 
2014 ! Parameter 408: legacy ceiling diagnostic
2015  IF (iget(408)>0) THEN
2016 !!$omp parallel do private(i,j)
2017  DO j=jsta,jend
2018  DO i=ista,iend
2019  grid1(i,j) = cldz(i,j)
2020  ENDDO
2021  ENDDO
2022  if(grib=="grib2" )then
2023  cfld=cfld+1
2024  fld_info(cfld)%ifld=iavblfld(iget(408))
2025  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2026  endif
2027  ENDIF
2028  ENDIF !End of GSD algorithm
2029 
2030 ! BEGIN EXPERIMENTAL GSD CEILING DIAGNOSTICS...
2031 ! J. Kenyon, 4 Feb 2017: this approach uses model-state cloud fractions
2032 ! Parameter 487: experimental ceiling diagnostic #1
2033  IF (iget(487)>0) THEN
2034 ! set some constants for ceiling adjustment in snow (retained from legacy algorithm, also in calvis.f)
2035  rhoice = 970.
2036  coeffp = 10.36
2037  exponfp = 1.
2038  const1 = 3.912
2039 ! set minimum cloud fraction to represent a ceiling
2040  ceiling_thresh_cldfra = 0.5
2041 
2042  DO j=jsta,jend
2043  DO i=ista,iend
2044  ceil(i,j) = spval
2045  zceil = spval
2046  cldfra_max = 0.
2047  do k=1,lm
2048  ll=lm-k+1
2049  cldfra(k) = cfr(i,j,ll)
2050  cldfra_max = max(cldfra_max,cldfra(k)) ! determine the column-maximum cloud fraction
2051  end do
2052 
2053  if (cldfra_max >= ceiling_thresh_cldfra) then ! threshold cloud fraction found in column, get ceiling
2054 
2055 ! threshold cloud fraction (possible ceiling) found somewhere in column, so proceed...
2056 ! first, search for and eliminate fog layers near surface (retained from legacy diagnostic)
2057  do k=2,3 ! Ming, k=3 will never be reached in this logic
2058  if (cldfra(k) < ceiling_thresh_cldfra) then ! these two lines:
2059  if (cldfra(1) > ceiling_thresh_cldfra) then ! ...look for surface-based fog beneath less-cloudy layers
2060  do k1=1,k-1 ! now perform the clearing for k=1 up to k-1
2061  if (cldfra(k1) >= ceiling_thresh_cldfra) then
2062  cldfra(k1)=0.
2063  end if
2064  end do
2065  end if
2066  ! level k=2,3 has no ceiling, and no fog at surface, so skip out of this loop
2067  end if
2068  exit
2069  end do ! k
2070 
2071 ! now search aloft...
2072  loop471:do k=2,lm
2073  k1 = k
2074  if (cldfra(k) >= ceiling_thresh_cldfra) then ! go to 472 ! found ceiling
2075  if (k1 <= 4) then ! within 4 levels of surface, no interpolation
2076  zceil = zmid(i,j,lm-k1+1)
2077  else ! use linear interpolation
2078  zceil = zmid(i,j,lm-k1+1) + (ceiling_thresh_cldfra-cldfra(k1)) &
2079  * (zmid(i,j,lm-k1+2)-zmid(i,j,lm-k1+1)) &
2080  / (cldfra(k1-1) - cldfra(k1))
2081  end if
2082  zceil = max(zceil,fis(i,j)*gi+5.)
2083 
2084 ! consider lowering of ceiling due to falling snow (retained from legacy diagnostic)
2085 ! ...this is extracted from calvis.f (visibility diagnostic)
2086  if (qqs(i,j,lm)>0.) then
2087  tv=t(i,j,lm)*(h1+d608*q(i,j,lm))
2088  rhoair=pmid(i,j,lm)/(rd*tv)
2089  vovermd = (1.+q(i,j,lm))/rhoair + qqs(i,j,lm)/rhoice
2090  concfp = qqs(i,j,lm)/vovermd*1000.
2091  betav = coeffp*concfp**exponfp + 1.e-10
2092  vertvis = 1000.*min(90., const1/betav)
2093  if (vertvis < zceil-fis(i,j)*gi ) then
2094  zceil = fis(i,j)*gi + vertvis
2095  exit loop471
2096  end if
2097  end if
2098 
2099  exit loop471
2100  endif ! cldfra(k) >= ceiling_thresh_cldfra
2101  end do loop471
2102  endif ! cldfra_max >= ceiling_thresh_cldfra
2103  ceil(i,j) = zceil
2104  ENDDO ! i loop
2105  ENDDO ! j loop
2106 
2107 ! Parameter 487: experimental ceiling diagnostic #1
2108  DO j=jsta,jend
2109  DO i=ista,iend
2110  grid1(i,j) = ceil(i,j)
2111  ENDDO
2112  ENDDO
2113  if(grib=="grib2" )then
2114  cfld=cfld+1
2115  fld_info(cfld)%ifld=iavblfld(iget(487))
2116  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2117  endif
2118  ENDIF ! end of parameter-487 conditional code
2119 ! END OF EXPERIMENTAL GSD CEILING DIAGNOSTIC 1
2120 
2121 ! BEGIN EXPERIMENTAL GSD CEILING DIAGNOSTIC 2
2122 ! -- J. Kenyon, 12 Sep 2019
2123 ! Parameter 711 has been developed to eventually replace the GSD
2124 ! legacy ceiling diagnostic, and can be regarded as a ceiling.
2125 ! However, for RAPv5/HRRRv4, paramater 711 will be supplied as
2126 ! the GSD cloud-base height, and parameter 798 will be the
2127 ! corresponding cloud-base pressure. (J. Kenyon, 4 Nov 2019)
2128 ! -- E. James, 15 Dec 2022
2129 ! The above experimental diagnostic, developed for the HRRR with
2130 ! lots of "add-ons" to correct for the HRRR's low bias in cloud
2131 ! cover, needs to be revised for the RRFS with its more extensive
2132 ! cloudiness. For an FAA deliverable due Feb 2023, the diagnostic
2133 ! is being modified to get rid of some of the add ons.
2134 
2135 ! Parameters 711/798: experimental ceiling diagnostic #2 (height and pressure, respectively)
2136  IF ((iget(711)>0) .OR. (iget(798)>0)) THEN
2137  ! set minimum cloud fraction to represent a ceiling
2138 ! ceiling_thresh_cldfra = 0.4
2139  ceiling_thresh_cldfra = 0.5
2140  ! set some constants for ceiling adjustment in snow (retained from legacy algorithm, also in calvis.f)
2141  rhoice = 970.
2142  coeffp = 10.36
2143  exponfp = 1.
2144  const1 = 3.912
2145 
2146  DO j=jsta,jend
2147  DO i=ista,iend
2148  ceil(i,j) = spval
2149  zceil = spval
2150  zceil1 = spval
2151  zceil2 = spval
2152  cldz(i,j) = spval
2153  cldp(i,j) = spval
2154 
2155  !-- Retrieve all cloud fractions in column
2156  do k=1,lm
2157  cldfra(k) = cfr(i,j,lm-k+1)
2158  end do
2159 
2160  !-- Look for surface-based clouds beneath
2161  ! less-cloudy layers. We will regard these
2162  ! instances as surface-based fog, too thin
2163  ! to impose a ceiling.
2164  if (cldfra(1) >= ceiling_thresh_cldfra) then ! possible thin fog; look higher
2165  do k=2,3
2166  if (cldfra(k) < 0.6) then ! confirmed thin fog, extending just below k
2167  cldfra(1:k-1) = 0.0 ! clear fog up to k-1
2168  end if
2169  end do
2170  end if
2171 
2172  !-- Search 1: no summation principle
2173  do k=2,lm
2174  if (cldfra(k) >= ceiling_thresh_cldfra) then ! found ceiling
2175  if (k <= 4) then ! within 4 levels of surface, no interpolation
2176  zceil1 = zmid(i,j,lm-k+1)
2177  else
2178  zceil1 = zmid(i,j,lm-k+1) + (ceiling_thresh_cldfra-cldfra(k)) &
2179  * (zmid(i,j,lm-k+2)-zmid(i,j,lm-k+1)) &
2180  / (cldfra(k-1) - cldfra(k))
2181  end if
2182  exit
2183  end if
2184  end do
2185 
2186  !-- Search 2: apply summation principle; see FAA order
2187  ! JO 7900.5B, Ch. 11 (Sky Condition), available at:
2188  ! https://www.faa.gov/documentLibrary/media/Order/7900_5D.pdf
2189  !
2190  ! and also:
2191  ! http://glossary.ametsoc.org/wiki/Summation_principle
2192  !
2193  ! J. Kenyon, 15 Aug 2019
2194 
2195  ! We seek to identify distinct cloud layers, which is
2196  ! not to be confused with model layers that contain
2197  ! clouds. For instance, a single layer of clouds may be
2198  ! represented across several adjoining model layers.
2199 
2200  cfr_layer_sum(1:lm)=0.0 ! initialize a column of zeros
2201  previous_sum=0.0
2202  do k=4,lm-1
2203  if ( (cldfra(k) >= 0.05 ) .and. & ! criterion 1
2204  (cldfra(k) > cldfra(k-1)) .and. & ! criterion 2
2205  (cldfra(k) >= cldfra(k+1)) ) & ! criterion 3
2206  ! Explanation, by criterion:
2207  ! (1) a reasonably large cloud fraction exists,
2208  ! (2) the cloud fraction is > the adjoining cloud fraction below,
2209  ! (3) the cloud fraction is >= the adjoining cloud fraction above (note that >=
2210  ! is used here, in case k is the lowest of several overcast model layers)
2211  then
2212  ! If all criteria satisfied, then we will consider the local-maximum cldfra(k) as
2213  ! representative of the cloud layer. We then proceed to add this fraction to
2214  ! the accumulated fraction(s) from any lower layer(s).
2215  cfr_layer_sum(k) = min(1.0, previous_sum + cldfra(k))
2216  previous_sum = min(1.0, cfr_layer_sum(k))
2217 
2218  if (cfr_layer_sum(k) >= ceiling_thresh_cldfra) then
2219  zceil2 = zmid(i,j,lm-k+1) + (ceiling_thresh_cldfra-cfr_layer_sum(k)) &
2220  * (zmid(i,j,lm-k+2)-zmid(i,j,lm-k+1)) &
2221  / (cfr_layer_sum(k-1) - cfr_layer_sum(k))
2222  exit ! break from DO K loop
2223  end if
2224 
2225  end if
2226  end do
2227  !-- end of search 2
2228 
2229 ! zceil = min(zceil1,zceil2) ! choose lower of zceil1 and zceil2
2230  zceil = zceil1
2231 
2232  !-- Search for "indefinite ceiling" (vertical visibility) conditions: consider
2233  ! lowering of apparent ceiling due to falling snow (retained from legacy
2234  ! diagnostic); this is extracted from calvis.f (visibility diagnostic)
2235 ! if (QQS(i,j,LM)>1.e-10) then
2236 ! TV=T(I,J,lm)*(H1+D608*Q(I,J,lm))
2237 ! RHOAIR=PMID(I,J,lm)/(RD*TV)
2238 ! vovermd = (1.+Q(i,j,LM))/rhoair + QQS(i,j,LM)/rhoice
2239 ! concfp = QQS(i,j,LM)/vovermd*1000.
2240 ! betav = coeffp*concfp**exponfp + 1.e-10
2241 ! vertvis = 1000.*min(90., const1/betav)
2242 ! if (vertvis < zceil-FIS(I,J)*GI ) then ! if vertvis is more restictive than zceil found above; set zceil to vertvis
2243 ! ! note that FIS is geopotential of the surface (ground), and GI is 1/g
2244 ! zceil = FIS(I,J)*GI + vertvis
2245 ! end if
2246 ! end if
2247 
2248  ceil(i,j) = zceil
2249  ENDDO ! i loop
2250  ENDDO ! j loop
2251 
2252  !-- In order to obtain a somewhat smoothed field of ceiling,
2253  ! we now conduct a horizontal search of neighboring grid
2254  ! boxes and consider each ceiling in AGL. The lowest
2255  ! neighboring AGL value will then be assigned locally.
2256  !
2257  ! In general, the diagnosis of low-AGL ceilings atop hills/peaks
2258  ! will also tend to be "spread" onto the adjacent valleys.
2259  ! However, a neighborhood search using heights in ASL is more
2260  ! problematic, since low ceilings in valleys will tend to be
2261  ! "spread" onto the ajacent hills/peaks as very low ceilings
2262  ! (fog). In actuality, these hills/peaks may exist above the cloud
2263  ! layer.
2264  allocate(full_ceil(im,jm),full_fis(im,jm))
2265  DO j=jsta,jend
2266  DO i=ista,iend
2267  full_ceil(i,j)=ceil(i,j)
2268  full_fis(i,j)=fis(i,j)
2269  ENDDO
2270  ENDDO
2271 ! CALL AllGETHERV(full_ceil)
2272  full_dummy=spval
2273  CALL collect_all(full_ceil(ista:iend,jsta:jend),full_dummy)
2274  full_ceil=full_dummy
2275 ! CALL AllGETHERV(full_fis)
2276  full_dummy=spval
2277  CALL collect_all(full_fis(ista:iend,jsta:jend),full_dummy)
2278  full_fis=full_dummy
2279 
2280  numr = 1
2281  DO j=jsta,jend
2282  DO i=ista,iend
2283  ceil_min = max( ceil(i,j)-fis(i,j)*gi , 5.0) ! ceil_min in AGL
2284  do jc = max(1,j-numr),min(jm,j+numr)
2285  do ic = max(1,i-numr),min(im,i+numr)
2286  ceil_neighbor = max( full_ceil(ic,jc)-full_fis(ic,jc)*gi , 5.0) ! ceil_neighbor in AGL
2287 ! ceil_min = min( ceil_min, ceil_neighbor )
2288  enddo
2289  enddo
2290  cldz(i,j) = ceil_min + fis(i,j)*gi ! convert back to ASL and store
2291  cldz(i,j) = max(min(cldz(i,j), 20000.0),0.0) !set bounds
2292  ! find pressure at CLDZ
2293  do k=2,lm-2
2294  if ( zmid(i,j,lm-k+1) >= cldz(i,j) ) then
2295  cldp(i,j) = pmid(i,j,lm-k+2) + (cldz(i,j)-zmid(i,j,lm-k+2)) &
2296  *(pmid(i,j,lm-k+1)-pmid(i,j,lm-k+2) ) &
2297  /(zmid(i,j,lm-k+1)-zmid(i,j,lm-k+2) )
2298  exit
2299  endif
2300  enddo
2301  ENDDO
2302  ENDDO
2303  if (allocated(full_ceil)) deallocate(full_ceil)
2304  if (allocated(full_fis)) deallocate(full_fis)
2305 
2306  ! Parameters 711/798: experimental ceiling diagnostic #2 (height and pressure, respectively)
2307  IF (iget(711)>0) THEN
2308 !!$omp parallel do private(i,j)
2309  DO j=jsta,jend
2310  DO i=ista,iend
2311  grid1(i,j) = cldz(i,j)
2312  ENDDO
2313  ENDDO
2314  if(grib=="grib2" )then
2315  cfld=cfld+1
2316  fld_info(cfld)%ifld=iavblfld(iget(711))
2317  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2318  endif
2319  ENDIF
2320 
2321  ! Parameters 711/798: experimental ceiling diagnostic #2 (height and pressure, respectively)
2322  IF (iget(798)>0) THEN
2323 !!$omp parallel do private(i,j)
2324  DO j=jsta,jend
2325  DO i=ista,iend
2326  grid1(i,j) = cldp(i,j)
2327  ENDDO
2328  ENDDO
2329  if(grib=="grib2" )then
2330  cfld=cfld+1
2331  fld_info(cfld)%ifld=iavblfld(iget(798))
2332  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2333  endif
2334  ENDIF
2335  ENDIF ! end of parameter-711 and -798 conditional code
2336 
2337 ! END OF EXPERIMENTAL GSD CEILING DIAGNOSTICS
2338 
2339 ! B. ZHOU: CEILING
2340  IF (iget(260)>0) THEN
2341  CALL calceiling(cldz,tcld,ceiling)
2342  DO j=jsta,jend
2343  DO i=ista,iend
2344  grid1(i,j) = ceiling(i,j)
2345  ENDDO
2346  ENDDO
2347  if(grib=="grib2" )then
2348  cfld=cfld+1
2349  fld_info(cfld)%ifld=iavblfld(iget(260))
2350  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2351  endif
2352  ENDIF
2353 ! B. ZHOU: FLIGHT CONDITION RESTRICTION
2354  IF (iget(261) > 0) THEN
2355  CALL calfltcnd(ceiling,grid1(1,jsta))
2356 ! DO J=JSTA,JEND
2357 ! DO I=ISTA,IEND
2358 ! GRID1(I,J) = FLTCND(I,J)
2359 ! ENDDO
2360 ! ENDDO
2361  if(grib=="grib2" )then
2362  cfld=cfld+1
2363  fld_info(cfld)%ifld=iavblfld(iget(261))
2364 !$omp parallel do private(i,j,ii,jj)
2365  do j=1,jend-jsta+1
2366  jj = jsta+j-1
2367  do i=1,iend-ista+1
2368  ii=ista+i-1
2369  datapd(i,j,cfld) = grid1(ii,jj)
2370  enddo
2371  enddo
2372  endif
2373  ENDIF
2374 !
2375 !--- Convective cloud base pressures (deep & shallow; Ferrier, Feb '02)
2376 !
2377  IF (iget(188) > 0) THEN
2378  IF(modelname == 'GFS')THEN
2379 !$omp parallel do private(i,j)
2380  DO j=jsta,jend
2381  DO i=ista,iend
2382  grid1(i,j) = pbot(i,j)
2383  ENDDO
2384  ENDDO
2385  ELSE
2386  DO j=jsta,jend
2387  DO i=ista,iend
2388  ibot=ibotcu(i,j)
2389  IF (ibot>0 .AND. ibot<=nint(lmh(i,j))) THEN
2390  grid1(i,j) = pmid(i,j,ibot)
2391  ELSE
2392  grid1(i,j) = -50000.
2393  ENDIF
2394  ENDDO
2395  ENDDO
2396  END IF
2397  if(grib=="grib2" )then
2398  cfld=cfld+1
2399  fld_info(cfld)%ifld=iavblfld(iget(188))
2400 !$omp parallel do private(i,j,ii,jj)
2401  do j=1,jend-jsta+1
2402  jj = jsta+j-1
2403  do i=1,iend-ista+1
2404  ii=ista+i-1
2405  datapd(i,j,cfld) = grid1(ii,jj)
2406  enddo
2407  enddo
2408  endif
2409  ENDIF
2410 !
2411 !--- Deep convective cloud base pressures (Ferrier, Feb '02)
2412 !
2413  IF (iget(192) > 0) THEN
2414  DO j=jsta,jend
2415  DO i=ista,iend
2416  ibot=ibotdcu(i,j)
2417  IF (ibot>0 .AND. ibot<=nint(lmh(i,j))) THEN
2418  grid1(i,j) = pmid(i,j,ibot)
2419  ELSE
2420  grid1(i,j) = -50000.
2421  ENDIF
2422  ENDDO
2423  ENDDO
2424  if(grib=="grib2" )then
2425  cfld=cfld+1
2426  fld_info(cfld)%ifld=iavblfld(iget(192))
2427  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2428  endif
2429  ENDIF
2430 !--- Shallow convective cloud base pressures (Ferrier, Feb '02)
2431 !
2432  IF (iget(190) > 0) THEN
2433  DO j=jsta,jend
2434  DO i=ista,iend
2435  ibot=ibotscu(i,j)
2436  IF (ibot>0 .AND. ibot<=nint(lmh(i,j))) THEN
2437  grid1(i,j) = pmid(i,j,ibot)
2438  ELSE
2439  grid1(i,j) = -50000.
2440  ENDIF
2441  ENDDO
2442  ENDDO
2443  if(grib=="grib2" )then
2444  cfld=cfld+1
2445  fld_info(cfld)%ifld=iavblfld(iget(190))
2446  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2447  endif
2448  ENDIF
2449 !--- Base of grid-scale cloudiness (Ferrier, Feb '02)
2450 !
2451  IF (iget(194) > 0) THEN
2452  DO j=jsta,jend
2453  DO i=ista,iend
2454  ibot=ibotgr(i,j)
2455  IF (ibot>0 .AND. ibot<=nint(lmh(i,j))) THEN
2456  grid1(i,j) = pmid(i,j,ibot)
2457  ELSE
2458  grid1(i,j) = -50000.
2459  ENDIF
2460  ENDDO
2461  ENDDO
2462  if(grib=="grib2" )then
2463  cfld=cfld+1
2464  fld_info(cfld)%ifld=iavblfld(iget(194))
2465  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2466  endif
2467  ENDIF
2468 
2469  !--- Base of low cloud
2470  !
2471  IF (iget(303) > 0) THEN
2472  DO j=jsta,jend
2473  DO i=ista,iend
2474 ! IF(PBOTL(I,J) > SMALL)THEN
2475  grid1(i,j) = pbotl(i,j)
2476 ! ELSE
2477 ! GRID1(I,J) = SPVAL
2478 ! END IF
2479  ENDDO
2480  ENDDO
2481  id(1:25)=0
2482  itclod = nint(tclod)
2483  IF(itclod /= 0) then
2484  ifincr = mod(ifhr,itclod)
2485  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
2486  ELSE
2487  ifincr = 0
2488  ENDIF
2489  id(19) = ifhr
2490  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2491  id(20) = 3
2492  IF (ifincr==0) THEN
2493  id(18) = ifhr-itclod
2494  ELSE
2495  id(18) = ifhr-ifincr
2496  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2497  ENDIF
2498  IF (id(18)<0) id(18) = 0
2499  if(grib=="grib2" )then
2500  cfld=cfld+1
2501  fld_info(cfld)%ifld=iavblfld(iget(303))
2502  if(itclod==0) then
2503  fld_info(cfld)%ntrange=0
2504  else
2505  fld_info(cfld)%ntrange=1
2506  endif
2507  fld_info(cfld)%tinvstat=ifhr-id(18)
2508 
2509  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2510  endif
2511  ENDIF
2512  !--- Base of middle cloud
2513  !
2514  IF (iget(306) > 0) THEN
2515  DO j=jsta,jend
2516  DO i=ista,iend
2517  IF(pbotm(i,j) > small)THEN
2518  grid1(i,j) = pbotm(i,j)
2519  ELSE
2520  grid1(i,j) = spval
2521  END IF
2522  ENDDO
2523  ENDDO
2524  id(1:25)=0
2525  itclod = nint(tclod)
2526  IF(itclod /= 0) then
2527  ifincr = mod(ifhr,itclod)
2528  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
2529  ELSE
2530  ifincr = 0
2531  ENDIF
2532  id(19) = ifhr
2533  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2534  id(20) = 3
2535  IF (ifincr==0) THEN
2536  id(18) = ifhr-itclod
2537  ELSE
2538  id(18) = ifhr-ifincr
2539  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2540  ENDIF
2541  IF (id(18)<0) id(18) = 0
2542  if(grib=="grib2" )then
2543  cfld=cfld+1
2544  fld_info(cfld)%ifld=iavblfld(iget(306))
2545  if(itclod==0) then
2546  fld_info(cfld)%ntrange=0
2547  else
2548  fld_info(cfld)%ntrange=1
2549  endif
2550  fld_info(cfld)%tinvstat=ifhr-id(18)
2551 
2552  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2553  endif
2554  ENDIF
2555  !--- Base of high cloud
2556  !
2557  IF (iget(309) > 0) THEN
2558  DO j=jsta,jend
2559  DO i=ista,iend
2560  IF(pboth(i,j) > small)THEN
2561  grid1(i,j) = pboth(i,j)
2562  ELSE
2563  grid1(i,j) = spval
2564  END IF
2565  ENDDO
2566  ENDDO
2567  id(1:25)=0
2568  itclod = nint(tclod)
2569  IF(itclod /= 0) then
2570  ifincr = mod(ifhr,itclod)
2571  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
2572  ELSE
2573  ifincr = 0
2574  ENDIF
2575  id(19) = ifhr
2576  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2577  id(20) = 3
2578  IF (ifincr==0) THEN
2579  id(18) = ifhr-itclod
2580  ELSE
2581  id(18) = ifhr-ifincr
2582  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2583  ENDIF
2584  IF (id(18)<0) id(18) = 0
2585  if(grib=="grib2" )then
2586  cfld=cfld+1
2587  fld_info(cfld)%ifld=iavblfld(iget(309))
2588  if(itclod==0) then
2589  fld_info(cfld)%ntrange=0
2590  else
2591  fld_info(cfld)%ntrange=1
2592  endif
2593  fld_info(cfld)%tinvstat=ifhr-id(18)
2594 
2595  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2596  endif
2597  ENDIF
2598 !
2599 !------------------------------------------------
2600 !----------- VARIOUS CLOUD TOP FIELDS ----------
2601 !------------------------------------------------
2602 !
2603 !--- "TOTAL" CLOUD TOP FIELDS (convective + grid-scale; Ferrier, Feb '02)
2604 !
2605  IF ((iget(149)>0) .OR. (iget(179)>0) .OR. &
2606  (iget(168)>0) .OR. (iget(275)>0)) THEN
2607  DO j=jsta,jend
2608  DO i=ista,iend
2609  itop=itopt(i,j)
2610  IF (itop>0 .AND. itop<=nint(lmh(i,j))) THEN
2611  IF(t(i,j,itop)<spval .AND. &
2612  pmid(i,j,itop)<spval .AND. &
2613  zmid(i,j,itop)<spval) THEN
2614  cldp(i,j) = pmid(i,j,itop)
2615  cldz(i,j) = zmid(i,j,itop)
2616  cldt(i,j) = t(i,j,itop)
2617  ELSE
2618  IF(modelname == 'RAPR') then
2619  cldp(i,j) = spval
2620  cldz(i,j) = spval
2621  ELSE
2622  cldp(i,j) = -50000.
2623  cldz(i,j) = -5000.
2624  ENDIF
2625  cldt(i,j) = -500.
2626  ENDIF
2627  ELSE
2628  IF(modelname == 'RAPR') then
2629  cldp(i,j) = spval
2630  cldz(i,j) = spval
2631  ELSE
2632  cldp(i,j) = -50000.
2633  cldz(i,j) = -5000.
2634  ENDIF
2635  cldt(i,j) = -500.
2636  ENDIF !--- End IF (ITOP>0 .AND. ITOP<=LMH(I,J)) ...
2637  ENDDO !--- End DO I loop
2638  ENDDO !--- End DO J loop
2639 !
2640 ! CLOUD TOP PRESSURE
2641 !
2642  IF (iget(149)>0) THEN
2643  DO j=jsta,jend
2644  DO i=ista,iend
2645  grid1(i,j) = cldp(i,j)
2646  ENDDO
2647  ENDDO
2648  if(grib=="grib2" )then
2649  cfld=cfld+1
2650  fld_info(cfld)%ifld=iavblfld(iget(149))
2651  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2652  endif
2653  ENDIF
2654 ! CLOUD TOP HEIGHT
2655 !
2656  IF (iget(179)>0) THEN
2657  DO j=jsta,jend
2658  DO i=ista,iend
2659  grid1(i,j) = cldz(i,j)
2660  ENDDO
2661  ENDDO
2662  if(grib=="grib2" )then
2663  cfld=cfld+1
2664  fld_info(cfld)%ifld=iavblfld(iget(179))
2665  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2666  endif
2667  ENDIF
2668  ENDIF
2669 
2670 ! GSD COULD TOP HEIGHTS AND PRESSURE
2671  IF ((iget(409)>0) .OR. (iget(406)>0)) THEN
2672 
2673  cloud_def_p = 0.0000001
2674 
2675  DO j=jsta,jend
2676  DO i=ista,iend
2677 ! imported from RUC post
2678 ! Cloud top
2679  zcldtop = -5000.
2680  IF(modelname == 'RAPR') zcldtop = spval
2681  do k=1,lm
2682  ll=lm-k+1
2683  watericetotal(k) = qqw(i,j,ll) + qqi(i,j,ll)
2684  enddo
2685 
2686  if (watericetotal(lm)<=cloud_def_p) then
2687  loop373 : do k=lm-1,2,-1
2688  if (watericetotal(k)>cloud_def_p) then
2689  zcldtop = zmid(i,j,lm-k+1) + (cloud_def_p-watericetotal(k)) &
2690  * (zmid(i,j,lm-k)-zmid(i,j,lm-k+1)) &
2691  / (watericetotal(k+1) - watericetotal(k))
2692  exit loop373
2693  end if
2694  end do loop373
2695  else
2696  zcldtop = zmid(i,j,1)
2697  end if
2698 
2699  itop=itopt(i,j)
2700  IF (itop>0 .AND. itop<=nint(lmh(i,j))) THEN
2701  cldp(i,j) = pmid(i,j,itop)
2702  cldt(i,j) = t(i,j,itop)
2703  ELSE
2704  cldp(i,j) = -50000.
2705  IF(modelname == 'RAPR') cldp(i,j) = spval
2706 ! CLDZ(I,J) = -5000.
2707  cldt(i,j) = -500.
2708  ENDIF !--- End IF (ITOP>0 .AND. ITOP<=LMH(I,J)) ...
2709 
2710 !- include convective clouds
2711  itop=itopcu(i,j)
2712  if(itop<lm+1) then
2713 ! print *,'ITOPCu(i,j)',i,j,ITOPCu(i,j)
2714  if(zcldtop <-100.) then
2715 ! print *,'add convective cloud, ITOP,CLDZ(I,J),ZMID(I,J,ITOP)'
2716 ! 1 ,ITOP,zcldtop,ZMID(I,J,ITOP),i,j
2717  zcldtop=zmid(i,j,itop)
2718  else if(zmid(i,j,itop)>zcldtop) then
2719 ! print *,'change cloud top for convective cloud, zcldtop,
2720 ! 1 ZMID(I,J,ITOP),ITOP,i,j'
2721 ! 1 ,zcldtop,ZMID(I,J,ITOP),ITOP,i,j
2722  zcldtop=zmid(i,j,itop)
2723  endif
2724  endif
2725 
2726 ! check consistency of cloud base and cloud top
2727  if(cldz(i,j)>-100. .and. zcldtop<-100.) then
2728  zcldtop = cldz(i,j) + 200.
2729  endif
2730 
2731  cldz(i,j) = zcldtop ! Now CLDZ is cloud top height
2732 
2733  ENDDO !--- End DO I loop
2734  ENDDO !--- End DO J loop
2735 !
2736 ! GSD CLOUD TOP PRESSURE
2737 !
2738  IF (iget(406)>0) THEN
2739  DO j=jsta,jend
2740  DO i=ista,iend
2741  grid1(i,j) = cldp(i,j)
2742  ENDDO
2743  ENDDO
2744  if(grib=="grib2" )then
2745  cfld=cfld+1
2746  fld_info(cfld)%ifld=iavblfld(iget(406))
2747  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2748  endif
2749  ENDIF
2750 ! GSD CLOUD TOP HEIGHT
2751 !
2752  IF (iget(409)>0) THEN
2753  DO j=jsta,jend
2754  DO i=ista,iend
2755  grid1(i,j) = cldz(i,j)
2756  ENDDO
2757  ENDDO
2758  if(grib=="grib2" )then
2759  cfld=cfld+1
2760  fld_info(cfld)%ifld=iavblfld(iget(409))
2761  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2762  endif
2763  ENDIF
2764  ENDIF ! end of GSD algorithm
2765 !
2766 ! CLOUD TOP TEMPS
2767 !
2768  IF (iget(168)>0) THEN
2769  DO j=jsta,jend
2770  DO i=ista,iend
2771  grid1(i,j) = cldt(i,j)
2772  ENDDO
2773  ENDDO
2774  if(grib=="grib2" )then
2775  cfld=cfld+1
2776  fld_info(cfld)%ifld=iavblfld(iget(168))
2777  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2778  endif
2779  ENDIF
2780 !
2781 !huang CLOUD TOP BRIGHTNESS TEMPERATURE
2782  IF (iget(275)>0) THEN
2783  num_thick=0 ! for debug
2784  grid1=spval
2785  DO j=jsta,jend
2786  DO i=ista,iend
2787  opdepth=0.
2788  llmh=nint(lmh(i,j))
2789 !bsf - start
2790 !-- Add subgrid-scale convective clouds for WRF runs
2791  do k=1,llmh
2792  cu_ir(k)=0.
2793  enddo
2794 ! Chuang: GFS specified non-convective grid points as missing
2795  if((hbot(i,j)-spval)>small .and. (htop(i,j)-spval)>small)then
2796  lcbot=nint(hbot(i,j))
2797  lctop=nint(htop(i,j))
2798  if (lcbot-lctop > 1) then
2799  q_conv=cnvcfr(i,j)*qconv
2800  do k=lctop,lcbot
2801  if (t(i,j,k) < trad_ice) then
2802  cu_ir(k)=abscoefi*q_conv
2803  else
2804  cu_ir(k)=abscoef*q_conv
2805  endif
2806  end do !-- do k = lctop,lcbot
2807  endif !-- if (lcbot-lctop > 1) then
2808  end if ! end of check for meaningful hbot and htop
2809  do k=1,llmh
2810 ! if(imp_physics==99 .and. t(i,j,k)<(tfrz-15.))then
2811 ! qqi(i,j,k)=qqw(i,j,k) ! because GFS only uses cloud water
2812 ! qqw(i,j,k)=0.
2813 ! end if
2814  if(pint(i,j,k)<spval.and.qqw(i,j,k)<spval.and. &
2815  qqi(i,j,k)<spval.and.qqs(i,j,k)<spval)then
2816  dp=pint(i,j,k+1)-pint(i,j,k)
2817  opdepth=opdepth+( cu_ir(k) + abscoef*qqw(i,j,k)+ &
2818 !bsf - end
2819  & abscoefi*( qqi(i,j,k)+qqs(i,j,k) ) )*dp
2820  endif
2821  if (opdepth > 1.) exit
2822  enddo
2823  if (opdepth > 1.) num_thick=num_thick+1 ! for debug
2824  k=min(k,llmh)
2825  grid1(i,j)=t(i,j,k)
2826  ENDDO
2827  ENDDO
2828 ! print *,'num_points, num_thick = ',(jend-jsta+1)*im,num_thick
2829 !! k=0
2830 !! 20 opdepthu=opdepthd
2831 !! k=k+1
2832 !!! if(k==1) then
2833 !!! dp=pint(i,j,itop+k)-pmid(i,j,itop)
2834 !!! opdepthd=opdepthu+(abscoef*(0.75*qqw(i,j,itop)+
2835 !!! & 0.25*qqw(i,j,itop+1))+abscoefi*
2836 !!! & (0.75*qqi(i,j,itop)+0.25*qqi(i,j,itop+1)))
2837 !!! & *dp/g
2838 !!! else
2839 !! dp=pint(i,j,k+1)-pint(i,j,k)
2840 !! opdepthd=opdepthu+(abscoef*qqw(i,j,k)+
2841 !! & abscoefi*qqi(i,j,k))*dp
2842 !!! end if
2843 !!
2844 !! lmhh=nint(lmh(i,j))
2845 !! if (opdepthd<1..and. k<lmhh) then
2846 !! goto 20
2847 !! elseif (opdepthd<1..and. k==lmhh) then
2848 !! GRID1(I,J)=T(i,j,lmhh )
2849 !!! prsctt=pmid(i,j,lmhh)
2850 !! else
2851 !!! GRID1(I,J)=T(i,j,k)
2852 !! if(k==1)then
2853 !! GRID1(I,J)=T(i,j,k)
2854 !! else if(k==lmhh)then
2855 !! GRID1(I,J)=T(i,j,k)
2856 !! else
2857 !! fac=(1.-opdepthu)/(opdepthd-opdepthu)
2858 !! GRID1(I,J)=(T(i,j,k)+T(i,j,k-1))/2.0+
2859 !! & (T(i,j,k+1)-T(i,j,k-1))/2.0*fac
2860 !! end if
2861 !!! prsctt=pf(i,j,k-1)+fac*(pf(i,j,k)-pf(i,j,k-1))
2862 !!! prsctt=min(prs(i,j,mkzh),max(prs(i,j,1),prsctt))
2863 !! endif
2864 !!! do 30 k=2,mkzh
2865 !!! if (prsctt>=prs(i,j,k-1).and.prsctt<=prs(i,j,k)) then
2866 !!! fac=(prsctt-prs(i,j,k-1))/(prs(i,j,k)-prs(i,j,k-1))
2867 !!! ctt(i,j)=tmk(i,j,k-1)+
2868 !!! & fac*(tmk(i,j,k)-tmk(i,j,k-1))-celkel
2869 !!! goto 40
2870 !!! endif
2871 !!! 30 continue
2872 !!! 40 continue
2873 !! END DO
2874 !! END DO
2875  if(grib=="grib2" )then
2876  cfld=cfld+1
2877  fld_info(cfld)%ifld=iavblfld(iget(275))
2878  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2879  endif
2880  ENDIF
2881 
2882 !
2883 !--- Convective cloud top pressures (deep & shallow; Ferrier, Feb '02)
2884 !
2885  IF (iget(189) > 0) THEN
2886  IF(modelname == 'GFS')THEN
2887 !$omp parallel do private(i,j)
2888  DO j=jsta,jend
2889  DO i=ista,iend
2890  grid1(i,j) = ptop(i,j)
2891  ENDDO
2892  ENDDO
2893  ELSE
2894  DO j=jsta,jend
2895  DO i=ista,iend
2896  itop=itopcu(i,j)
2897  IF (itop>0 .AND. itop<=nint(lmh(i,j))) THEN
2898  grid1(i,j) = pmid(i,j,itop)
2899  ELSE
2900  grid1(i,j) = -50000.
2901  ENDIF
2902  ENDDO
2903  ENDDO
2904  END IF
2905  if(grib=="grib2" )then
2906  cfld=cfld+1
2907  fld_info(cfld)%ifld=iavblfld(iget(189))
2908 !$omp parallel do private(i,j,ii,jj)
2909  do j=1,jend-jsta+1
2910  jj = jsta+j-1
2911  do i=1,iend-ista+1
2912  ii=ista+i-1
2913  datapd(i,j,cfld) = grid1(ii,jj)
2914  enddo
2915  enddo
2916  endif
2917  END IF
2918 !
2919 !--- Deep convective cloud top pressures (Ferrier, Feb '02)
2920 !
2921  IF (iget(193) > 0) THEN
2922  DO j=jsta,jend
2923  DO i=ista,iend
2924  itop=itopdcu(i,j)
2925  IF (itop>0 .AND. itop<=nint(lmh(i,j))) THEN
2926  grid1(i,j) = pmid(i,j,itop)
2927  ELSE
2928  grid1(i,j) = -50000.
2929  ENDIF
2930  ENDDO
2931  ENDDO
2932  if(grib=="grib2" )then
2933  cfld=cfld+1
2934  fld_info(cfld)%ifld=iavblfld(iget(193))
2935  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2936  endif
2937  END IF
2938 !--- Shallow convective cloud top pressures (Ferrier, Feb '02)
2939 !
2940  IF (iget(191) > 0) THEN
2941  DO j=jsta,jend
2942  DO i=ista,iend
2943  itop=itopscu(i,j)
2944  IF (itop>0 .AND. itop<=nint(lmh(i,j))) THEN
2945  grid1(i,j) = pmid(i,j,itop)
2946  ELSE
2947  grid1(i,j) = -50000.
2948  ENDIF
2949  ENDDO
2950  ENDDO
2951  if(grib=="grib2" )then
2952  cfld=cfld+1
2953  fld_info(cfld)%ifld=iavblfld(iget(191))
2954  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2955  endif
2956  END IF
2957 !
2958 !--- Top of grid-scale cloudiness (Ferrier, Feb '02)
2959 !
2960  IF (iget(195) > 0) THEN
2961  DO j=jsta,jend
2962  DO i=ista,iend
2963  itop=itopgr(i,j)
2964  IF (itop>0 .AND. itop<=nint(lmh(i,j))) THEN
2965  grid1(i,j) = pmid(i,j,itop)
2966  ELSE
2967  grid1(i,j) = -50000.
2968  ENDIF
2969  ENDDO
2970  ENDDO
2971  if(grib=="grib2" )then
2972  cfld=cfld+1
2973  fld_info(cfld)%ifld=iavblfld(iget(195))
2974  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
2975  endif
2976  END IF
2977 
2978  !--- top of low cloud
2979  !
2980  IF (iget(304) > 0) THEN
2981  DO j=jsta,jend
2982  DO i=ista,iend
2983  IF(ptopl(i,j) > small)THEN
2984  grid1(i,j) = ptopl(i,j)
2985  ELSE
2986  grid1(i,j) = spval
2987  END IF
2988  ENDDO
2989  ENDDO
2990  id(1:25)=0
2991  itclod = nint(tclod)
2992  IF(itclod /= 0) then
2993  ifincr = mod(ifhr,itclod)
2994  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
2995  ELSE
2996  ifincr = 0
2997  ENDIF
2998  id(19) = ifhr
2999  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3000  id(20) = 3
3001  IF (ifincr==0) THEN
3002  id(18) = ifhr-itclod
3003  ELSE
3004  id(18) = ifhr-ifincr
3005  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3006  ENDIF
3007  IF (id(18)<0) id(18) = 0
3008  if(grib=="grib2" )then
3009  cfld=cfld+1
3010  fld_info(cfld)%ifld=iavblfld(iget(304))
3011  if(itclod==0) then
3012  fld_info(cfld)%ntrange=0
3013  else
3014  fld_info(cfld)%ntrange=1
3015  endif
3016  fld_info(cfld)%tinvstat=ifhr-id(18)
3017 
3018  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3019  endif
3020  ENDIF
3021  !--- top of middle cloud
3022  !
3023  IF (iget(307) > 0) THEN
3024  DO j=jsta,jend
3025  DO i=ista,iend
3026  grid1(i,j) = ptopm(i,j)
3027  ENDDO
3028  ENDDO
3029  id(1:25)=0
3030  itclod = nint(tclod)
3031  IF(itclod /= 0) then
3032  ifincr = mod(ifhr,itclod)
3033  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
3034  ELSE
3035  ifincr = 0
3036  ENDIF
3037  id(19) = ifhr
3038  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3039  id(20) = 3
3040  IF (ifincr==0) THEN
3041  id(18) = ifhr-itclod
3042  ELSE
3043  id(18) = ifhr-ifincr
3044  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3045  ENDIF
3046  IF (id(18)<0) id(18) = 0
3047  if(grib=="grib2" )then
3048  cfld=cfld+1
3049  fld_info(cfld)%ifld=iavblfld(iget(307))
3050  if(itclod==0) then
3051  fld_info(cfld)%ntrange=0
3052  else
3053  fld_info(cfld)%ntrange=1
3054  endif
3055  fld_info(cfld)%tinvstat=ifhr-id(18)
3056 
3057  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3058  endif
3059  ENDIF
3060  !--- top of high cloud
3061  !
3062  IF (iget(310) > 0) THEN
3063  DO j=jsta,jend
3064  DO i=ista,iend
3065  grid1(i,j) = ptoph(i,j)
3066  ENDDO
3067  ENDDO
3068  id(1:25)=0
3069  itclod = nint(tclod)
3070  IF(itclod /= 0) then
3071  ifincr = mod(ifhr,itclod)
3072  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
3073  ELSE
3074  ifincr = 0
3075  ENDIF
3076  id(19) = ifhr
3077  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3078  id(20) = 3
3079  IF (ifincr==0) THEN
3080  id(18) = ifhr-itclod
3081  ELSE
3082  id(18) = ifhr-ifincr
3083  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3084  ENDIF
3085  IF (id(18)<0) id(18) = 0
3086  if(grib=="grib2" )then
3087  cfld=cfld+1
3088  fld_info(cfld)%ifld=iavblfld(iget(310))
3089  if(itclod==0) then
3090  fld_info(cfld)%ntrange=0
3091  else
3092  fld_info(cfld)%ntrange=1
3093  endif
3094  fld_info(cfld)%tinvstat=ifhr-id(18)
3095 
3096  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3097  endif
3098  ENDIF
3099 
3100  !--- T of low cloud top
3101  !
3102  IF (iget(305) > 0) THEN
3103  DO j=jsta,jend
3104  DO i=ista,iend
3105  grid1(i,j) = ttopl(i,j)
3106  ENDDO
3107  ENDDO
3108  id(1:25)=0
3109  itclod = nint(tclod)
3110  IF(itclod /= 0) then
3111  ifincr = mod(ifhr,itclod)
3112  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
3113  ELSE
3114  ifincr = 0
3115  ENDIF
3116  id(19) = ifhr
3117  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3118  id(20) = 3
3119  IF (ifincr==0) THEN
3120  id(18) = ifhr-itclod
3121  ELSE
3122  id(18) = ifhr-ifincr
3123  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3124  ENDIF
3125  IF (id(18)<0) id(18) = 0
3126  if(grib=="grib2" )then
3127  cfld=cfld+1
3128  fld_info(cfld)%ifld=iavblfld(iget(305))
3129  if(itclod==0) then
3130  fld_info(cfld)%ntrange=0
3131  else
3132  fld_info(cfld)%ntrange=1
3133  endif
3134  fld_info(cfld)%tinvstat=ifhr-id(18)
3135 
3136  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3137  endif
3138  ENDIF
3139  !--- Base of middle cloud
3140  !
3141  IF (iget(308) > 0) THEN
3142  DO j=jsta,jend
3143  DO i=ista,iend
3144  grid1(i,j) = ttopm(i,j)
3145  ENDDO
3146  ENDDO
3147  id(1:25)=0
3148  itclod = nint(tclod)
3149  IF(itclod /= 0) then
3150  ifincr = mod(ifhr,itclod)
3151  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
3152  ELSE
3153  ifincr = 0
3154  ENDIF
3155  id(19) = ifhr
3156  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3157  id(20) = 3
3158  IF (ifincr==0) THEN
3159  id(18) = ifhr-itclod
3160  ELSE
3161  id(18) = ifhr-ifincr
3162  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3163  ENDIF
3164  IF (id(18)<0) id(18) = 0
3165  if(grib=="grib2" )then
3166  cfld=cfld+1
3167  fld_info(cfld)%ifld=iavblfld(iget(308))
3168  if(itclod==0) then
3169  fld_info(cfld)%ntrange=0
3170  else
3171  fld_info(cfld)%ntrange=1
3172  endif
3173  fld_info(cfld)%tinvstat=ifhr-id(18)
3174 
3175  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3176  endif
3177  ENDIF
3178  !--- Base of high cloud
3179  !
3180  IF (iget(311) > 0) THEN
3181  DO j=jsta,jend
3182  DO i=ista,iend
3183  grid1(i,j) = ttoph(i,j)
3184  ENDDO
3185  ENDDO
3186  id(1:25)=0
3187  itclod = nint(tclod)
3188  IF(itclod /= 0) then
3189  ifincr = mod(ifhr,itclod)
3190  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
3191  ELSE
3192  ifincr = 0
3193  ENDIF
3194  id(19) = ifhr
3195  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3196  id(20) = 3
3197  IF (ifincr==0) THEN
3198  id(18) = ifhr-itclod
3199  ELSE
3200  id(18) = ifhr-ifincr
3201  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3202  ENDIF
3203  IF (id(18)<0) id(18) = 0
3204  if(grib=="grib2" )then
3205  cfld=cfld+1
3206  fld_info(cfld)%ifld=iavblfld(iget(311))
3207  if(itclod==0) then
3208  fld_info(cfld)%ntrange=0
3209  else
3210  fld_info(cfld)%ntrange=1
3211  endif
3212  fld_info(cfld)%tinvstat=ifhr-id(18)
3213  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3214  endif
3215  ENDIF
3216 !
3217 !--- Convective cloud fractions from modified Slingo (1987)
3218 !
3219  IF (iget(196) > 0.or.iget(570)>0) THEN
3220  grid1=spval
3221  DO j=jsta,jend
3222  DO i=ista,iend
3223  if(cnvcfr(i,j)/=spval)grid1(i,j)=100.*cnvcfr(i,j) !-- convert to percent
3224  ENDDO
3225  ENDDO
3226  if(iget(196)>0) then
3227  if(grib=="grib2" )then
3228  cfld=cfld+1
3229  fld_info(cfld)%ifld=iavblfld(iget(196))
3230  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3231  endif
3232  elseif(iget(570)>0) then
3233  if(grib=="grib2" )then
3234  cfld=cfld+1
3235  fld_info(cfld)%ifld=iavblfld(iget(570))
3236  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3237  endif
3238  endif
3239  END IF
3240 !
3241 !--- Boundary layer cloud fractions
3242 !
3243  IF (iget(342) > 0) THEN
3244  grid1=spval
3245  DO j=jsta,jend
3246  DO i=ista,iend
3247  if(pblcfr(i,j)/=spval)grid1(i,j)=100.*pblcfr(i,j) !-- convert to percent
3248  ENDDO
3249  ENDDO
3250  id(1:25)=0
3251  itclod = nint(tclod)
3252  IF(itclod /= 0) then
3253  ifincr = mod(ifhr,itclod)
3254  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
3255  ELSE
3256  ifincr = 0
3257  ENDIF
3258  id(19) = ifhr
3259  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3260  id(20) = 3
3261  IF (ifincr==0) THEN
3262  id(18) = ifhr-itclod
3263  ELSE
3264  id(18) = ifhr-ifincr
3265  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3266  ENDIF
3267  IF (id(18)<0) id(18) = 0
3268  if(grib=="grib2" )then
3269  cfld=cfld+1
3270  fld_info(cfld)%ifld=iavblfld(iget(342))
3271  if(itclod==0) then
3272  fld_info(cfld)%ntrange=0
3273  else
3274  fld_info(cfld)%ntrange=1
3275  endif
3276  fld_info(cfld)%tinvstat=ifhr-id(18)
3277 
3278  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3279  endif
3280  END IF
3281 !
3282 !--- Cloud work function
3283 !
3284  IF (iget(313) > 0) THEN
3285  DO j=jsta,jend
3286  DO i=ista,iend
3287  grid1(i,j)=cldwork(i,j)
3288  ENDDO
3289  ENDDO
3290  id(1:25)=0
3291  itclod = nint(tclod)
3292  IF(itclod /= 0) then
3293  ifincr = mod(ifhr,itclod)
3294  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itclod*60)
3295  ELSE
3296  ifincr = 0
3297  ENDIF
3298  id(19) = ifhr
3299  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3300  id(20) = 3
3301  IF (ifincr==0) THEN
3302  id(18) = ifhr-itclod
3303  ELSE
3304  id(18) = ifhr-ifincr
3305  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3306  ENDIF
3307  IF (id(18)<0) id(18) = 0
3308  if(grib=="grib2" )then
3309  cfld=cfld+1
3310  fld_info(cfld)%ifld=iavblfld(iget(313))
3311  if(itclod==0) then
3312  fld_info(cfld)%ntrange=0
3313  else
3314  fld_info(cfld)%ntrange=1
3315  endif
3316  fld_info(cfld)%tinvstat=ifhr-id(18)
3317 
3318  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3319  endif
3320  END IF
3321 !
3322 !*** BLOCK 3. RADIATION FIELDS.
3323 !
3324 !
3325 ! TIME AVERAGED SURFACE SHORT WAVE INCOMING RADIATION.
3326  IF (iget(126)>0) THEN
3327  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3328  grid1=spval
3329  id(1:25)=0
3330  ELSE
3331 ! print*,'ARDSW in CLDRAD=',ARDSW
3332  IF(ardsw>0.)THEN
3333  rrnum=1./ardsw
3334  ELSE
3335  rrnum=0.
3336  ENDIF
3337  DO j=jsta,jend
3338  DO i=ista,iend
3339  IF(aswin(i,j)/=spval)THEN
3340  grid1(i,j) = aswin(i,j)*rrnum
3341  ELSE
3342  grid1(i,j)=aswin(i,j)
3343  END IF
3344  ENDDO
3345  ENDDO
3346  id(1:25)=0
3347  itrdsw = nint(trdsw)
3348  IF(itrdsw /= 0) then
3349  ifincr = mod(ifhr,itrdsw)
3350  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
3351  ELSE
3352  ifincr = 0
3353  endif
3354  id(19) = ifhr
3355  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3356  id(20) = 3
3357  IF (ifincr==0) THEN
3358  id(18) = ifhr-itrdsw
3359  ELSE
3360  id(18) = ifhr-ifincr
3361  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3362  ENDIF
3363  IF (id(18)<0) id(18) = 0
3364  END IF
3365  if(grib=="grib2" )then
3366  cfld=cfld+1
3367  fld_info(cfld)%ifld=iavblfld(iget(126))
3368  if(itrdsw>0) then
3369  fld_info(cfld)%ntrange=1
3370  else
3371  fld_info(cfld)%ntrange=0
3372  endif
3373  fld_info(cfld)%tinvstat=ifhr-id(18)
3374  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3375  endif
3376  ENDIF
3377 !
3378 ! TIME AVERAGED SURFACE UV-B INCOMING RADIATION.
3379  IF (iget(298)>0) THEN
3380  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3381  grid1=spval
3382  id(1:25)=0
3383  ELSE
3384 ! print*,'ARDSW in CLDRAD=',ARDSW
3385  IF(ardsw>0.)THEN
3386  rrnum=1./ardsw
3387  ELSE
3388  rrnum=0.
3389  ENDIF
3390  DO j=jsta,jend
3391  DO i=ista,iend
3392  IF(auvbin(i,j)/=spval)THEN
3393  grid1(i,j) = auvbin(i,j)*rrnum
3394  ELSE
3395  grid1(i,j) = auvbin(i,j)
3396  END IF
3397  ENDDO
3398  ENDDO
3399  id(1:25)=0
3400  id(02)=129
3401  itrdsw = nint(trdsw)
3402  IF(itrdsw /= 0) then
3403  ifincr = mod(ifhr,itrdsw)
3404  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
3405  ELSE
3406  ifincr = 0
3407  endif
3408  id(19) = ifhr
3409  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3410  id(20) = 3
3411  IF (ifincr==0) THEN
3412  id(18) = ifhr-itrdsw
3413  ELSE
3414  id(18) = ifhr-ifincr
3415  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3416  ENDIF
3417  IF (id(18)<0) id(18) = 0
3418  END IF
3419  if(grib=="grib2" )then
3420  cfld=cfld+1
3421  fld_info(cfld)%ifld=iavblfld(iget(298))
3422  if(itrdsw>0) then
3423  fld_info(cfld)%ntrange=1
3424  else
3425  fld_info(cfld)%ntrange=0
3426  endif
3427  fld_info(cfld)%tinvstat=ifhr-id(18)
3428  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3429  endif
3430  ENDIF
3431 !
3432 ! TIME AVERAGED SURFACE UV-B CLEAR SKY INCOMING RADIATION.
3433  IF (iget(297)>0) THEN
3434  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3435  grid1=spval
3436  id(1:25)=0
3437  ELSE
3438 ! print*,'ARDSW in CLDRAD=',ARDSW
3439  IF(ardsw>0.)THEN
3440  rrnum=1./ardsw
3441  ELSE
3442  rrnum=0.
3443  ENDIF
3444  DO j=jsta,jend
3445  DO i=ista,iend
3446  IF(auvbinc(i,j)/=spval)THEN
3447  grid1(i,j) = auvbinc(i,j)*rrnum
3448  ELSE
3449  grid1(i,j) = auvbinc(i,j)
3450  END IF
3451  ENDDO
3452  ENDDO
3453  id(1:25)=0
3454  id(02)=129
3455  itrdsw = nint(trdsw)
3456  IF(itrdsw /= 0) then
3457  ifincr = mod(ifhr,itrdsw)
3458  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
3459  ELSE
3460  ifincr = 0
3461  endif
3462  id(19) = ifhr
3463  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3464  id(20) = 3
3465  IF (ifincr==0) THEN
3466  id(18) = ifhr-itrdsw
3467  ELSE
3468  id(18) = ifhr-ifincr
3469  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3470  ENDIF
3471  IF (id(18)<0) id(18) = 0
3472  END IF
3473  if(grib=="grib2" )then
3474  cfld=cfld+1
3475  fld_info(cfld)%ifld=iavblfld(iget(297))
3476  if(itrdsw>0) then
3477  fld_info(cfld)%ntrange=1
3478  else
3479  fld_info(cfld)%ntrange=0
3480  endif
3481  fld_info(cfld)%tinvstat=ifhr-id(18)
3482  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3483  endif
3484  ENDIF
3485 !
3486 ! TIME AVERAGED SURFACE LONG WAVE INCOMING RADIATION.
3487  IF (iget(127)>0) THEN
3488  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3489  grid1=spval
3490  id(1:25)=0
3491  ELSE
3492  IF(ardlw>0.)THEN
3493  rrnum=1./ardlw
3494  ELSE
3495  rrnum=0.
3496  ENDIF
3497  DO j=jsta,jend
3498  DO i=ista,iend
3499  IF(alwin(i,j)/=spval)THEN
3500  grid1(i,j) = alwin(i,j)*rrnum
3501  ELSE
3502  grid1(i,j)=alwin(i,j)
3503  END IF
3504  ENDDO
3505  ENDDO
3506  id(1:25)=0
3507  itrdlw = nint(trdlw)
3508  IF(itrdlw /= 0) then
3509  ifincr = mod(ifhr,itrdlw)
3510  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdlw*60)
3511  ELSE
3512  ifincr = 0
3513  endif
3514  id(19) = ifhr
3515  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3516  id(20) = 3
3517  IF (ifincr==0) THEN
3518  id(18) = ifhr-itrdlw
3519  ELSE
3520  id(18) = ifhr-ifincr
3521  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3522  ENDIF
3523  IF (id(18)<0) id(18) = 0
3524  END IF
3525  if(grib=="grib2" )then
3526  cfld=cfld+1
3527  fld_info(cfld)%ifld=iavblfld(iget(127))
3528  if(itrdlw>0) then
3529  fld_info(cfld)%ntrange=1
3530  else
3531  fld_info(cfld)%ntrange=0
3532  endif
3533  fld_info(cfld)%tinvstat=ifhr-id(18)
3534  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3535  endif
3536  ENDIF
3537 !
3538 ! TIME AVERAGED SURFACE SHORT WAVE OUTGOING RADIATION.
3539  IF (iget(128)>0) THEN
3540  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3541  grid1=spval
3542  id(1:25)=0
3543  ELSE
3544  IF(ardsw>0.)THEN
3545  rrnum=1./ardsw
3546  ELSE
3547  rrnum=0.
3548  ENDIF
3549  DO j=jsta,jend
3550  DO i=ista,iend
3551  IF(aswout(i,j)/=spval)THEN
3552  grid1(i,j) = -1.0*aswout(i,j)*rrnum
3553  ELSE
3554  grid1(i,j)=aswout(i,j)
3555  END IF
3556  ENDDO
3557  ENDDO
3558  id(1:25)=0
3559  itrdsw = nint(trdsw)
3560  IF(itrdsw /= 0) then
3561  ifincr = mod(ifhr,itrdsw)
3562  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
3563  ELSE
3564  ifincr = 0
3565  endif
3566  id(19) = ifhr
3567  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3568  id(20) = 3
3569  IF (ifincr==0) THEN
3570  id(18) = ifhr-itrdsw
3571  ELSE
3572  id(18) = ifhr-ifincr
3573  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3574  ENDIF
3575  IF (id(18)<0) id(18) = 0
3576  END IF
3577  if(grib=="grib2" )then
3578  cfld=cfld+1
3579  fld_info(cfld)%ifld=iavblfld(iget(128))
3580  if(itrdsw>0) then
3581  fld_info(cfld)%ntrange=1
3582  else
3583  fld_info(cfld)%ntrange=0
3584  endif
3585  fld_info(cfld)%tinvstat=ifhr-id(18)
3586  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3587  endif
3588  ENDIF
3589 !
3590 ! TIME AVERAGED SURFACE LONG WAVE OUTGOING RADIATION.
3591  IF (iget(129)>0) THEN
3592  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3593  grid1=spval
3594  id(1:25)=0
3595  ELSE
3596  IF(ardlw>0.)THEN
3597  rrnum=1./ardlw
3598  ELSE
3599  rrnum=0.
3600  ENDIF
3601  DO j=jsta,jend
3602  DO i=ista,iend
3603  IF(alwout(i,j)/=spval)THEN
3604  grid1(i,j) = -1.0*alwout(i,j)*rrnum
3605  ELSE
3606  grid1(i,j)=alwout(i,j)
3607  END IF
3608  ENDDO
3609  ENDDO
3610  id(1:25)=0
3611  itrdlw = nint(trdlw)
3612  IF(itrdlw /= 0) then
3613  ifincr = mod(ifhr,itrdlw)
3614  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdlw*60)
3615  ELSE
3616  ifincr = 0
3617  endif
3618  id(19) = ifhr
3619  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3620  id(20) = 3
3621  IF (ifincr==0) THEN
3622  id(18) = ifhr-itrdlw
3623  ELSE
3624  id(18) = ifhr-ifincr
3625  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3626  ENDIF
3627  IF (id(18)<0) id(18) = 0
3628  END IF
3629  if(grib=="grib2" )then
3630  cfld=cfld+1
3631  fld_info(cfld)%ifld=iavblfld(iget(129))
3632  if(itrdlw>0) then
3633  fld_info(cfld)%ntrange=1
3634  else
3635  fld_info(cfld)%ntrange=0
3636  endif
3637  fld_info(cfld)%tinvstat=ifhr-id(18)
3638  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3639  endif
3640  ENDIF
3641 !
3642 ! TIME AVERAGED TOP OF THE ATMOSPHERE SHORT WAVE RADIATION.
3643  IF (iget(130)>0) THEN
3644  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3645  grid1=spval
3646  id(1:25)=0
3647  ELSE
3648  IF(ardsw>0.)THEN
3649  rrnum=1./ardsw
3650  ELSE
3651  rrnum=0.
3652  ENDIF
3653  DO j=jsta,jend
3654  DO i=ista,iend
3655  IF(aswtoa(i,j)/=spval)THEN
3656  grid1(i,j) = aswtoa(i,j)*rrnum
3657  ELSE
3658  grid1(i,j)=aswtoa(i,j)
3659  END IF
3660  ENDDO
3661  ENDDO
3662  id(1:25)=0
3663  itrdsw = nint(trdsw)
3664  IF(itrdsw /= 0) then
3665  ifincr = mod(ifhr,itrdsw)
3666  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
3667  ELSE
3668  ifincr = 0
3669  endif
3670  id(19) = ifhr
3671  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3672  id(20) = 3
3673  IF (ifincr==0) THEN
3674  id(18) = ifhr-itrdsw
3675  ELSE
3676  id(18) = ifhr-ifincr
3677  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3678  ENDIF
3679  IF (id(18)<0) id(18) = 0
3680  END IF
3681  if(grib=="grib2" )then
3682  cfld=cfld+1
3683  fld_info(cfld)%ifld=iavblfld(iget(130))
3684  if(itrdsw>0) then
3685  fld_info(cfld)%ntrange=1
3686  else
3687  fld_info(cfld)%ntrange=0
3688  endif
3689  fld_info(cfld)%tinvstat=ifhr-id(18)
3690  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3691  endif
3692  ENDIF
3693 !
3694 ! TIME AVERAGED TOP OF THE ATMOSPHERE LONG WAVE RADIATION.
3695  IF (iget(131)>0) THEN
3696  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3697  grid1=spval
3698  id(1:25)=0
3699  ELSE
3700  IF(ardlw>0.)THEN
3701  rrnum=1./ardlw
3702  ELSE
3703  rrnum=0.
3704  ENDIF
3705  DO j=jsta,jend
3706  DO i=ista,iend
3707  IF(alwtoa(i,j)/=spval)THEN
3708  grid1(i,j) = alwtoa(i,j)*rrnum
3709  ELSE
3710  grid1(i,j)=alwtoa(i,j)
3711  END IF
3712  ENDDO
3713  ENDDO
3714  id(1:25)=0
3715  itrdlw = nint(trdlw)
3716  IF(itrdlw /= 0) then
3717  ifincr = mod(ifhr,itrdlw)
3718  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdlw*60)
3719  ELSE
3720  ifincr = 0
3721  endif
3722  id(19) = ifhr
3723  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3724  id(20) = 3
3725  IF (ifincr==0) THEN
3726  id(18) = ifhr-itrdlw
3727  ELSE
3728  id(18) = ifhr-ifincr
3729  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3730  ENDIF
3731  IF (id(18)<0) id(18) = 0
3732  END IF
3733  if(grib=="grib2" )then
3734  cfld=cfld+1
3735  fld_info(cfld)%ifld=iavblfld(iget(131))
3736  if(itrdlw>0) then
3737  fld_info(cfld)%ntrange=1
3738  else
3739  fld_info(cfld)%ntrange=0
3740  endif
3741  fld_info(cfld)%tinvstat=ifhr-id(18)
3742  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3743  endif
3744  ENDIF
3745 !
3746 ! CURRENT TOP OF THE ATMOSPHERE LONG WAVE RADIATION.
3747  IF (iget(274)>0) THEN
3748  IF(modelname == 'NCAR'.OR.modelname=='RSM')THEN
3749  grid1=spval
3750  ELSE
3751  DO j=jsta,jend
3752  DO i=ista,iend
3753  grid1(i,j) = rlwtoa(i,j)
3754  ENDDO
3755  ENDDO
3756  END IF
3757  if(grib=="grib2" )then
3758  cfld=cfld+1
3759  fld_info(cfld)%ifld=iavblfld(iget(274))
3760  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3761  endif
3762  ENDIF
3763 !
3764 ! CLOUD TOP BRIGHTNESS TEMPERATURE FROM TOA OUTGOING LW.
3765  IF (iget(265)>0) THEN
3766  grid1=spval
3767  IF(modelname == 'NCAR'.OR.modelname=='RSM' .OR. modelname == 'RAPR')THEN
3768  grid1=spval
3769  ELSE
3770  DO j=jsta,jend
3771  DO i=ista,iend
3772  IF(rlwtoa(i,j) < spval) &
3773  & grid1(i,j) = (rlwtoa(i,j)*stbol)**0.25
3774  ENDDO
3775  ENDDO
3776  END IF
3777  if(grib=="grib2" )then
3778  cfld=cfld+1
3779  fld_info(cfld)%ifld=iavblfld(iget(265))
3780  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3781  endif
3782  ENDIF
3783 !
3784 ! CURRENT INCOMING SW RADIATION AT THE SURFACE.
3785  IF (iget(156)>0) THEN
3786  grid1=spval
3787  DO j=jsta,jend
3788  DO i=ista,iend
3789  IF(rswin(i,j)<spval) THEN
3790  IF(czmean(i,j)>1.e-6) THEN
3791  factrs=czen(i,j)/czmean(i,j)
3792  ELSE
3793  factrs=0.0
3794  ENDIF
3795  IF(rswin(i,j)<spval) grid1(i,j)=rswin(i,j)*factrs
3796  ENDIF
3797  ENDDO
3798  ENDDO
3799 !
3800  if(grib=="grib2" )then
3801  cfld=cfld+1
3802  fld_info(cfld)%ifld=iavblfld(iget(156))
3803  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3804  endif
3805  ENDIF
3806 !
3807 ! CURRENT INCOMING LW RADIATION AT THE SURFACE.
3808  IF (iget(157)>0) THEN
3809 ! dong add missing value to DLWRF
3810  grid1 = spval
3811  DO j=jsta,jend
3812  DO i=ista,iend
3813  IF(modelname=='RSM' .OR. modelname == 'RAPR') THEN !add by Binbin: RSM has direct RLWIN output
3814  grid1(i,j)=rlwin(i,j)
3815  ELSE
3816  IF(sigt4(i,j)<spval.and.t(i,j,nint(lmh(i,j)))<spval) THEN
3817  IF(sigt4(i,j)>0.0) THEN
3818  llmh=nint(lmh(i,j))
3819  tlmh=t(i,j,llmh)
3820  factrl=5.67e-8*tlmh*tlmh*tlmh*tlmh/sigt4(i,j)
3821  ELSE
3822  factrl=0.0
3823  ENDIF
3824  IF(rlwin(i,j) < spval) grid1(i,j)=rlwin(i,j)*factrl
3825  ENDIF
3826  ENDIF
3827  ENDDO
3828  ENDDO
3829 !
3830  if(grib=="grib2" )then
3831  cfld=cfld+1
3832  fld_info(cfld)%ifld=iavblfld(iget(157))
3833  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3834  endif
3835  ENDIF
3836 !
3837 ! CURRENT OUTGOING SW RADIATION AT THE SURFACE.
3838  IF (iget(141)>0) THEN
3839  grid1 = spval
3840 !$omp parallel do private(i,j)
3841  DO j=jsta,jend
3842  DO i=ista,iend
3843  IF(rswout(i,j)<spval) THEN
3844  IF(czmean(i,j)>1.e-6) THEN
3845  factrs=czen(i,j)/czmean(i,j)
3846  ELSE
3847  factrs=0.0
3848  ENDIF
3849  IF(rswout(i,j)<spval) grid1(i,j)=rswout(i,j)*factrs
3850  ENDIF
3851  ENDDO
3852  ENDDO
3853 !
3854  if(grib=="grib2" )then
3855  cfld=cfld+1
3856  fld_info(cfld)%ifld=iavblfld(iget(141))
3857  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3858  endif
3859  ENDIF
3860 
3861 ! CURRENT OUTGOING LW RADIATION AT THE SURFACE.
3862  IF (iget(142)>0) THEN
3863 !$omp parallel do private(i,j)
3864  DO j=jsta,jend
3865  DO i=ista,iend
3866  grid1(i,j) = radot(i,j)
3867  ENDDO
3868  ENDDO
3869  if(grib=="grib2" )then
3870  cfld=cfld+1
3871  fld_info(cfld)%ifld=iavblfld(iget(142))
3872  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3873  endif
3874  ENDIF
3875 
3876 ! Instantaneous clear-sky downwelling LW at the surface
3877  IF (iget(764)>0) THEN
3878  DO j=jsta,jend
3879  DO i=ista,iend
3880  grid1(i,j) = lwdnbc(i,j)
3881  ENDDO
3882  ENDDO
3883  if(grib=='grib2') then
3884  cfld=cfld+1
3885  fld_info(cfld)%ifld=iavblfld(iget(764))
3886  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3887  endif
3888  ENDIF
3889 
3890 ! Instantaneous clear-sky upwelling LW at the surface
3891  IF (iget(765)>0) THEN
3892  DO j=jsta,jend
3893  DO i=ista,iend
3894  grid1(i,j) = lwupbc(i,j)
3895  ENDDO
3896  ENDDO
3897  if(grib=='grib2') then
3898  cfld=cfld+1
3899  fld_info(cfld)%ifld=iavblfld(iget(765))
3900  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3901  endif
3902  ENDIF
3903 
3904 ! Instantaneous MEAN_FRP
3905  IF (iget(740)>0) THEN
3906  DO j=jsta,jend
3907  DO i=ista,iend
3908  grid1(i,j) = mean_frp(i,j)
3909  ENDDO
3910  ENDDO
3911  if(grib=='grib2') then
3912  cfld=cfld+1
3913  fld_info(cfld)%ifld=iavblfld(iget(740))
3914  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3915  endif
3916  ENDIF
3917 
3918 ! Biomass burning emissions (EBB)
3919  IF (iget(745)>0) THEN
3920  DO j=jsta,jend
3921  DO i=ista,iend
3922  IF (ebb(i,j)<spval) THEN
3923  grid1(i,j) = ebb(i,j)/(1e9)
3924  ELSE
3925  grid1(i,j) = spval
3926  ENDIF
3927  ENDDO
3928  ENDDO
3929  if(grib=='grib2') then
3930  cfld=cfld+1
3931  fld_info(cfld)%ifld=iavblfld(iget(745))
3932  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3933  endif
3934  ENDIF
3935 
3936 ! Hourly wildfire potential (HWP)
3937  IF (iget(755)>0) THEN
3938  DO j=jsta,jend
3939  DO i=ista,iend
3940  IF (hwp(i,j)<spval) THEN
3941  grid1(i,j) = hwp(i,j)
3942  ELSE
3943  grid1(i,j) = spval
3944  ENDIF
3945  ENDDO
3946  ENDDO
3947  if(grib=='grib2') then
3948  cfld=cfld+1
3949  fld_info(cfld)%ifld=iavblfld(iget(755))
3950  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3951  endif
3952  ENDIF
3953 
3954 ! CURRENT (instantaneous) INCOMING CLEARSKY SW RADIATION AT THE SURFACE.
3955  IF (iget(262)>0) THEN
3956  grid1 = spval
3957 !$omp parallel do private(i,j)
3958  DO j=jsta,jend
3959  DO i=ista,iend
3960  IF(rswinc(i,j)<spval) THEN
3961  IF(czmean(i,j)>1.e-6) THEN
3962  factrs=czen(i,j)/czmean(i,j)
3963  ELSE
3964  factrs=0.0
3965  ENDIF
3966  IF(rswinc(i,j)<spval) grid1(i,j) = rswinc(i,j)*factrs
3967  ENDIF
3968  ENDDO
3969  ENDDO
3970  if(grib=="grib2" )then
3971  cfld=cfld+1
3972  fld_info(cfld)%ifld=iavblfld(iget(262))
3973  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3974  endif
3975  ENDIF
3976 
3977 ! Instantaneous SWDDNI
3978  IF (iget(772)>0)THEN
3979 !$omp parallel do private(i,j)
3980  DO j=jsta,jend
3981  DO i=ista,iend
3982  grid1(i,j) = swddni(i,j)
3983  ENDDO
3984  ENDDO
3985  if(grib=='grib2') then
3986  cfld=cfld+1
3987  fld_info(cfld)%ifld=iavblfld(iget(772))
3988  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3989  endif
3990  ENDIF
3991 
3992 ! Instantaneous clear-sky SWDDNI
3993  IF (iget(796)>0) THEN
3994  DO j=jsta,jend
3995  DO i=ista,iend
3996  grid1(i,j) = swddnic(i,j)
3997  ENDDO
3998  ENDDO
3999  if(grib=='grib2') then
4000  cfld=cfld+1
4001  fld_info(cfld)%ifld=iavblfld(iget(796))
4002  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4003  endif
4004  ENDIF
4005 
4006 ! Instantaneous SWDDIF
4007  IF (iget(773)>0) THEN
4008 !$omp parallel do private(i,j)
4009  DO j=jsta,jend
4010  DO i=ista,iend
4011  grid1(i,j) = swddif(i,j)
4012  ENDDO
4013  ENDDO
4014  if(grib=='grib2') then
4015  cfld=cfld+1
4016  fld_info(cfld)%ifld=iavblfld(iget(773))
4017  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4018  endif
4019  ENDIF
4020 
4021 ! Instantaneous clear-sky SWDDIF
4022  IF (iget(797)>0) THEN
4023  DO j=jsta,jend
4024  DO i=ista,iend
4025  grid1(i,j) = swddifc(i,j)
4026  ENDDO
4027  ENDDO
4028  if(grib=='grib2') then
4029  cfld=cfld+1
4030  fld_info(cfld)%ifld=iavblfld(iget(797))
4031  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4032  endif
4033  ENDIF
4034 
4035 ! TIME AVERAGED INCOMING CLEARSKY SW RADIATION AT THE SURFACE.
4036  IF (iget(383)>0) THEN
4037  DO j=jsta,jend
4038  DO i=ista,iend
4039  grid1(i,j) = aswinc(i,j)
4040  ENDDO
4041  ENDDO
4042  id(1:25)=0
4043  itrdsw = nint(trdsw)
4044  IF(itrdsw /= 0) then
4045  ifincr = mod(ifhr,itrdsw)
4046  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4047  ELSE
4048  ifincr = 0
4049  endif
4050  id(19) = ifhr
4051  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4052  id(20) = 3
4053  IF (ifincr==0) THEN
4054  id(18) = ifhr-itrdsw
4055  ELSE
4056  id(18) = ifhr-ifincr
4057  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4058  ENDIF
4059  IF (id(18)<0) id(18) = 0
4060  if(grib=="grib2" )then
4061  cfld=cfld+1
4062  fld_info(cfld)%ifld=iavblfld(iget(383))
4063  if(itrdsw>0) then
4064  fld_info(cfld)%ntrange=1
4065  else
4066  fld_info(cfld)%ntrange=0
4067  endif
4068  fld_info(cfld)%tinvstat=ifhr-id(18)
4069  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4070  endif
4071  ENDIF
4072 !
4073 ! TIME AVERAGED OUTGOING CLEARSKY SW RADIATION AT THE SURFACE.
4074  IF (iget(386)>0) THEN
4075  DO j=jsta,jend
4076  DO i=ista,iend
4077  grid1(i,j) = aswoutc(i,j)
4078  ENDDO
4079  ENDDO
4080  id(1:25)=0
4081  itrdsw = nint(trdsw)
4082  IF(itrdsw /= 0) then
4083  ifincr = mod(ifhr,itrdsw)
4084  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4085  ELSE
4086  ifincr = 0
4087  endif
4088  id(19) = ifhr
4089  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4090  id(20) = 3
4091  IF (ifincr==0) THEN
4092  id(18) = ifhr-itrdsw
4093  ELSE
4094  id(18) = ifhr-ifincr
4095  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4096  ENDIF
4097  IF (id(18)<0) id(18) = 0
4098  if(grib=="grib2" )then
4099  cfld=cfld+1
4100  fld_info(cfld)%ifld=iavblfld(iget(386))
4101  if(itrdsw>0) then
4102  fld_info(cfld)%ntrange=1
4103  else
4104  fld_info(cfld)%ntrange=0
4105  endif
4106  fld_info(cfld)%tinvstat=ifhr-id(18)
4107  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4108  endif
4109  ENDIF
4110 
4111 ! Instantaneous all-sky outgoing SW flux at the model top
4112  IF (iget(719)>0) THEN
4113  DO j=jsta,jend
4114  DO i=ista,iend
4115  grid1(i,j) = swupt(i,j)
4116  ENDDO
4117  ENDDO
4118  if(grib=='grib2') then
4119  cfld=cfld+1
4120  fld_info(cfld)%ifld=iavblfld(iget(719))
4121  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4122  endif
4123  ENDIF
4124 
4125 ! TIME AVERAGED OUTGOING CLEARSKY SW RADIATION AT THE MODEL TOP
4126  IF (iget(387)>0) THEN
4127  DO j=jsta,jend
4128  DO i=ista,iend
4129  grid1(i,j) = aswtoac(i,j)
4130  ENDDO
4131  ENDDO
4132  id(1:25)=0
4133  itrdsw = nint(trdsw)
4134  IF(itrdsw /= 0) then
4135  ifincr = mod(ifhr,itrdsw)
4136  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4137  ELSE
4138  ifincr = 0
4139  endif
4140  id(19) = ifhr
4141  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4142  id(20) = 3
4143  IF (ifincr==0) THEN
4144  id(18) = ifhr-itrdsw
4145  ELSE
4146  id(18) = ifhr-ifincr
4147  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4148  ENDIF
4149  IF (id(18)<0) id(18) = 0
4150  if(grib=="grib2" )then
4151  cfld=cfld+1
4152  fld_info(cfld)%ifld=iavblfld(iget(387))
4153  if(itrdsw>0) then
4154  fld_info(cfld)%ntrange=1
4155  else
4156  fld_info(cfld)%ntrange=0
4157  endif
4158  fld_info(cfld)%tinvstat=ifhr-id(18)
4159  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4160  endif
4161  ENDIF
4162 !
4163 ! TIME AVERAGED INCOMING SW RADIATION AT THE MODEL TOP
4164  IF (iget(388)>0) THEN
4165  DO j=jsta,jend
4166  DO i=ista,iend
4167  grid1(i,j) = aswintoa(i,j)
4168  ENDDO
4169  ENDDO
4170  id(1:25)=0
4171  itrdsw = nint(trdsw)
4172  IF(itrdsw /= 0) then
4173  ifincr = mod(ifhr,itrdsw)
4174  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4175  ELSE
4176  ifincr = 0
4177  endif
4178  id(19) = ifhr
4179  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4180  id(20) = 3
4181  IF (ifincr==0) THEN
4182  id(18) = ifhr-itrdsw
4183  ELSE
4184  id(18) = ifhr-ifincr
4185  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4186  ENDIF
4187  IF (id(18)<0) id(18) = 0
4188  if(grib=="grib2" )then
4189  cfld=cfld+1
4190  fld_info(cfld)%ifld=iavblfld(iget(388))
4191  if(itrdsw>0) then
4192  fld_info(cfld)%ntrange=1
4193  else
4194  fld_info(cfld)%ntrange=0
4195  endif
4196  fld_info(cfld)%tinvstat=ifhr-id(18)
4197  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4198  endif
4199  ENDIF
4200 !
4201 ! TIME AVERAGED INCOMING CLEARSKY LW RADIATION AT THE SURFACE
4202  IF (iget(382)>0) THEN
4203  DO j=jsta,jend
4204  DO i=ista,iend
4205  grid1(i,j) = alwinc(i,j)
4206  ENDDO
4207  ENDDO
4208  id(1:25)=0
4209  itrdlw = nint(trdlw)
4210  IF(itrdlw /= 0) then
4211  ifincr = mod(ifhr,itrdlw)
4212  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdlw*60)
4213  ELSE
4214  ifincr = 0
4215  endif
4216  id(19) = ifhr
4217  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4218  id(20) = 3
4219  IF (ifincr==0) THEN
4220  id(18) = ifhr-itrdlw
4221  ELSE
4222  id(18) = ifhr-ifincr
4223  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4224  ENDIF
4225  IF (id(18)<0) id(18) = 0
4226  if(grib=="grib2" )then
4227  cfld=cfld+1
4228  fld_info(cfld)%ifld=iavblfld(iget(382))
4229  if(itrdlw>0) then
4230  fld_info(cfld)%ntrange=1
4231  else
4232  fld_info(cfld)%ntrange=0
4233  endif
4234  fld_info(cfld)%tinvstat=ifhr-id(18)
4235  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4236  endif
4237  ENDIF
4238 !
4239 ! TIME AVERAGED OUTGOING CLEARSKY LW RADIATION AT THE SURFACE
4240  IF (iget(384)>0) THEN
4241  DO j=jsta,jend
4242  DO i=ista,iend
4243  grid1(i,j) = alwoutc(i,j)
4244  ENDDO
4245  ENDDO
4246  id(1:25)=0
4247  itrdlw = nint(trdlw)
4248  IF(itrdlw /= 0) then
4249  ifincr = mod(ifhr,itrdlw)
4250  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdlw*60)
4251  ELSE
4252  ifincr = 0
4253  endif
4254  id(19) = ifhr
4255  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4256  id(20) = 3
4257  IF (ifincr==0) THEN
4258  id(18) = ifhr-itrdlw
4259  ELSE
4260  id(18) = ifhr-ifincr
4261  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4262  ENDIF
4263  IF (id(18)<0) id(18) = 0
4264  if(grib=="grib2" )then
4265  cfld=cfld+1
4266  fld_info(cfld)%ifld=iavblfld(iget(384))
4267  if(itrdlw>0) then
4268  fld_info(cfld)%ntrange=1
4269  else
4270  fld_info(cfld)%ntrange=0
4271  endif
4272  fld_info(cfld)%tinvstat=ifhr-id(18)
4273  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4274  endif
4275  ENDIF
4276 !
4277 ! TIME AVERAGED OUTGOING CLEARSKY LW RADIATION AT THE MODEL TOP
4278  IF (iget(385)>0) THEN
4279  DO j=jsta,jend
4280  DO i=ista,iend
4281  grid1(i,j) = alwtoac(i,j)
4282  ENDDO
4283  ENDDO
4284  id(1:25)=0
4285  itrdlw = nint(trdlw)
4286  IF(itrdlw /= 0) then
4287  ifincr = mod(ifhr,itrdlw)
4288  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdlw*60)
4289  ELSE
4290  ifincr = 0
4291  endif
4292  id(19) = ifhr
4293  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4294  id(20) = 3
4295  IF (ifincr==0) THEN
4296  id(18) = ifhr-itrdlw
4297  ELSE
4298  id(18) = ifhr-ifincr
4299  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4300  ENDIF
4301  IF (id(18)<0) id(18) = 0
4302  if(grib=="grib2" )then
4303  cfld=cfld+1
4304  fld_info(cfld)%ifld=iavblfld(iget(385))
4305  if(itrdlw>0) then
4306  fld_info(cfld)%ntrange=1
4307  else
4308  fld_info(cfld)%ntrange=0
4309  endif
4310  fld_info(cfld)%tinvstat=ifhr-id(18)
4311  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4312  endif
4313  ENDIF
4314 !
4315 ! TIME AVERAGED SURFACE VISIBLE BEAM DOWNWARD SOLAR FLUX
4316  IF (iget(401)>0) THEN
4317  DO j=jsta,jend
4318  DO i=ista,iend
4319  grid1(i,j) = avisbeamswin(i,j)
4320  ENDDO
4321  ENDDO
4322  id(1:25)=0
4323  itrdsw = nint(trdsw)
4324  IF(itrdsw /= 0) then
4325  ifincr = mod(ifhr,itrdsw)
4326  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4327  ELSE
4328  ifincr = 0
4329  endif
4330  id(19) = ifhr
4331  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4332  id(20) = 3
4333  IF (ifincr==0) THEN
4334  id(18) = ifhr-itrdsw
4335  ELSE
4336  id(18) = ifhr-ifincr
4337  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4338  ENDIF
4339  IF (id(18)<0) id(18) = 0
4340 ! CFS labels time ave fields as inst in long range forecast
4341  IF(itrdsw < 0)id(1:25)=0
4342  if(grib=="grib2" )then
4343  cfld=cfld+1
4344  fld_info(cfld)%ifld=iavblfld(iget(401))
4345  if(itrdsw>0) then
4346  fld_info(cfld)%ntrange=1
4347  else
4348  fld_info(cfld)%ntrange=0
4349  endif
4350  fld_info(cfld)%tinvstat=ifhr-id(18)
4351  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4352  endif
4353  ENDIF
4354 !
4355 ! TIME AVERAGED SURFACE VISIBLE DIFFUSE DOWNWARD SOLAR FLUX
4356  IF (iget(402)>0) THEN
4357  DO j=jsta,jend
4358  DO i=ista,iend
4359  grid1(i,j) = avisdiffswin(i,j)
4360  ENDDO
4361  ENDDO
4362  id(1:25)=0
4363  itrdsw = nint(trdsw)
4364  IF(itrdsw /= 0) then
4365  ifincr = mod(ifhr,itrdsw)
4366  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4367  ELSE
4368  ifincr = 0
4369  endif
4370  id(19) = ifhr
4371  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4372  id(20) = 3
4373  IF (ifincr==0) THEN
4374  id(18) = ifhr-itrdsw
4375  ELSE
4376  id(18) = ifhr-ifincr
4377  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4378  ENDIF
4379  IF (id(18)<0) id(18) = 0
4380  IF(itrdsw < 0)id(1:25)=0
4381  if(grib=="grib2" )then
4382  cfld=cfld+1
4383  fld_info(cfld)%ifld=iavblfld(iget(402))
4384  if(itrdsw>0) then
4385  fld_info(cfld)%ntrange=1
4386  else
4387  fld_info(cfld)%ntrange=0
4388  endif
4389  fld_info(cfld)%tinvstat=ifhr-id(18)
4390  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4391  endif
4392  ENDIF
4393 !
4394 ! TIME AVERAGED SURFACE VISIBLE BEAM DOWNWARD SOLAR FLUX
4395  IF (iget(403)>0) THEN
4396  DO j=jsta,jend
4397  DO i=ista,iend
4398  grid1(i,j) = airbeamswin(i,j)
4399  ENDDO
4400  ENDDO
4401  id(1:25)=0
4402  itrdsw = nint(trdsw)
4403  IF(itrdsw /= 0) then
4404  ifincr = mod(ifhr,itrdsw)
4405  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4406  ELSE
4407  ifincr = 0
4408  endif
4409  id(19) = ifhr
4410  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4411  id(20) = 3
4412  IF (ifincr==0) THEN
4413  id(18) = ifhr-itrdsw
4414  ELSE
4415  id(18) = ifhr-ifincr
4416  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4417  ENDIF
4418  IF (id(18)<0) id(18) = 0
4419  IF(itrdsw < 0)id(1:25)=0
4420  if(grib=="grib2" )then
4421  cfld=cfld+1
4422  fld_info(cfld)%ifld=iavblfld(iget(403))
4423  if(itrdsw>0) then
4424  fld_info(cfld)%ntrange=1
4425  else
4426  fld_info(cfld)%ntrange=0
4427  endif
4428  fld_info(cfld)%tinvstat=ifhr-id(18)
4429  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4430  endif
4431  ENDIF
4432 !
4433 ! TIME AVERAGED SURFACE VISIBLE DIFFUSE DOWNWARD SOLAR FLUX
4434  IF (iget(404)>0) THEN
4435  DO j=jsta,jend
4436  DO i=ista,iend
4437  grid1(i,j) = airdiffswin(i,j)
4438  ENDDO
4439  ENDDO
4440  id(1:25)=0
4441  itrdsw = nint(trdsw)
4442  IF(itrdsw /= 0) then
4443  ifincr = mod(ifhr,itrdsw)
4444  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itrdsw*60)
4445  ELSE
4446  ifincr = 0
4447  endif
4448  id(19) = ifhr
4449  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
4450  id(20) = 3
4451  IF (ifincr==0) THEN
4452  id(18) = ifhr-itrdsw
4453  ELSE
4454  id(18) = ifhr-ifincr
4455  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
4456  ENDIF
4457  IF (id(18)<0) id(18) = 0
4458  IF(itrdsw < 0)id(1:25)=0
4459  if(grib=="grib2" )then
4460  cfld=cfld+1
4461  fld_info(cfld)%ifld=iavblfld(iget(404))
4462  if(itrdsw>0) then
4463  fld_info(cfld)%ntrange=1
4464  else
4465  fld_info(cfld)%ntrange=0
4466  endif
4467  fld_info(cfld)%tinvstat=ifhr-id(18)
4468  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4469  endif
4470  ENDIF
4471 
4472  !2D AEROSOL OPTICAL DEPTH AT 550 NM
4473  IF(rdaod) then
4474  IF (iget(600).GT.0) THEN
4475  DO j=jsta,jend
4476  DO i=ista,iend
4477  grid1(i,j)=aod550(i,j)
4478  ENDDO
4479  ENDDO
4480  if(grib=="grib2" )then
4481  cfld=cfld+1
4482  fld_info(cfld)%ifld=iavblfld(iget(600))
4483  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4484  endif
4485  ENDIF
4486 
4487  IF (iget(601).GT.0) THEN
4488  DO j=jsta,jend
4489  DO i=ista,iend
4490  grid1(i,j)=du_aod550(i,j)
4491  ENDDO
4492  ENDDO
4493  if(grib=="grib2" )then
4494  cfld=cfld+1
4495  fld_info(cfld)%ifld=iavblfld(iget(601))
4496  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4497  endif
4498  ENDIF
4499 
4500  IF (iget(602).GT.0) THEN
4501  DO j=jsta,jend
4502  DO i=ista,iend
4503  grid1(i,j)=ss_aod550(i,j)
4504  ENDDO
4505  ENDDO
4506  if(grib=="grib2" )then
4507  cfld=cfld+1
4508  fld_info(cfld)%ifld=iavblfld(iget(602))
4509  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4510  endif
4511  ENDIF
4512 
4513  IF (iget(603).GT.0) THEN
4514  DO j=jsta,jend
4515  DO i=ista,iend
4516  grid1(i,j)=su_aod550(i,j)
4517  ENDDO
4518  ENDDO
4519  if(grib=="grib2" )then
4520  cfld=cfld+1
4521  fld_info(cfld)%ifld=iavblfld(iget(603))
4522  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4523  endif
4524  ENDIF
4525 
4526  IF (iget(604).GT.0) THEN
4527  DO j=jsta,jend
4528  DO i=ista,iend
4529  grid1(i,j)=oc_aod550(i,j)
4530  ENDDO
4531  ENDDO
4532  if(grib=="grib2" )then
4533  cfld=cfld+1
4534  fld_info(cfld)%ifld=iavblfld(iget(604))
4535  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4536  endif
4537  ENDIF
4538 
4539 
4540  IF (iget(605).GT.0) THEN
4541  DO j=jsta,jend
4542  DO i=ista,iend
4543  grid1(i,j)=bc_aod550(i,j)
4544  ENDDO
4545  ENDDO
4546  if(grib=="grib2" )then
4547  cfld=cfld+1
4548  fld_info(cfld)%ifld=iavblfld(iget(605))
4549  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4550  endif
4551  ENDIF
4552  END IF !rdaod
4553 
4554  !AQM AEROSOL OPTICAL DEPTH AT 550 NM
4555  IF (aqf_on) THEN
4556  IF (iget(712).GT.0) THEN
4557  DO j=jsta,jend
4558  DO i=ista,iend
4559  grid1(i,j)=aqm_aod550(i,j)
4560  ENDDO
4561  ENDDO
4562  if(grib=="grib2" )then
4563  cfld=cfld+1
4564  fld_info(cfld)%ifld=iavblfld(iget(712))
4565  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4566  endif
4567  ENDIF
4568  END IF !aqf_on
4569 
4570  !2D AEROSOL OPTICAL DEPTH AT 550 NM
4571  IF (iget(715)>0) THEN
4572  DO j=jsta,jend
4573  DO i=ista,iend
4574  grid1(i,j)=taod5502d(i,j)
4575  ENDDO
4576  ENDDO
4577  if(grib=="grib2" )then
4578  cfld=cfld+1
4579  fld_info(cfld)%ifld=iavblfld(iget(715))
4580  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4581  endif
4582  ENDIF
4583 
4584  !AEROSOL ASYMMETRY FACTOR
4585  IF (iget(716)>0) THEN
4586  DO j=jsta,jend
4587  DO i=ista,iend
4588  grid1(i,j)=aerasy2d(i,j)
4589  ENDDO
4590  ENDDO
4591  if(grib=="grib2" )then
4592  cfld=cfld+1
4593  fld_info(cfld)%ifld=iavblfld(iget(716))
4594  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4595  endif
4596  ENDIF
4597 
4598  !AEROSOL SINGLE-SCATTERING ALBEDO
4599  IF (iget(717)>0) THEN
4600  DO j=jsta,jend
4601  DO i=ista,iend
4602  grid1(i,j)=aerssa2d(i,j)
4603  ENDDO
4604  ENDDO
4605  if(grib=="grib2" )then
4606  cfld=cfld+1
4607  fld_info(cfld)%ifld=iavblfld(iget(717))
4608  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
4609  endif
4610  ENDIF
4611 !
4612  if (gocart_on .or. gccpp_on .or. nasa_on) then
4613 !
4614 !*** BLOCK 4. GOCART AEROSOL FIELDS
4615 !
4616 !
4617 !!! ADD AOD AT 550 NM AND OTHER CHANNELS
4618 
4619 !! Aerosol Optical Depth (AOD)
4620 !! CALPW(dust mixing ratio in kg/kg) => column density (kg/m2)
4621 !! CALPW(dust mixing ratio in kg/kg * Qext [aerosol extinction efficiency
4622 !! in m2/g] * 1000. [convert m2/g to m2/kg]) => AOD (no unit)
4623 !!
4624 
4625 !! DETERMINE WHETHER TO COMPUTE AEROSOL OPTICAL PROPERTIES
4626  laeropt = .false.
4627  DO i = 609, 614 ! TOTAL AND SPECIATED AOD AT 550NM
4628  IF ( iget(i)>0 ) laeropt = .true.
4629  ENDDO
4630  DO i = 623, 628 ! AOD AT MULTI-CHANNELS
4631  IF ( iget(i)>0 ) laeropt = .true.
4632  ENDDO
4633  DO i = 648, 656 ! (SSA, ASY AT 340),(SCA AT 550), ANGSTROM
4634  IF ( iget(i)>0 ) laeropt = .true.
4635  ENDDO
4636 
4637 !! DETERMINE WHETHER TO COMPUTE INSTANT SURFACE MASS CONC
4638  laersmass = .false.
4639  DO i = 690, 698 ! TOTAL AND SPECIATED AEROSOL
4640  IF ( iget(i)>0 ) laersmass = .true.
4641  ENDDO
4642 
4643  IF ( laeropt ) THEN
4644  print *, 'COMPUTE AEROSOL OPTICAL PROPERTIES'
4645 
4646 !!! ALLOCATE AEROSOL OPTICAL PROPERTIES
4647  ALLOCATE ( extrhd_du(krhlev,nbin_du,nbdsw))
4648  ALLOCATE ( extrhd_ss(krhlev,nbin_ss,nbdsw))
4649  ALLOCATE ( extrhd_su(krhlev,nbin_su,nbdsw))
4650  ALLOCATE ( extrhd_bc(krhlev,nbin_bc,nbdsw))
4651  ALLOCATE ( extrhd_oc(krhlev,nbin_oc,nbdsw))
4652 
4653  ALLOCATE ( scarhd_du(krhlev,nbin_du,nbdsw))
4654  ALLOCATE ( scarhd_ss(krhlev,nbin_ss,nbdsw))
4655  ALLOCATE ( scarhd_su(krhlev,nbin_su,nbdsw))
4656  ALLOCATE ( scarhd_bc(krhlev,nbin_bc,nbdsw))
4657  ALLOCATE ( scarhd_oc(krhlev,nbin_oc,nbdsw))
4658 
4659  ALLOCATE ( asyrhd_du(krhlev,nbin_du,nbdsw))
4660  ALLOCATE ( asyrhd_ss(krhlev,nbin_ss,nbdsw))
4661  ALLOCATE ( asyrhd_su(krhlev,nbin_su,nbdsw))
4662  ALLOCATE ( asyrhd_bc(krhlev,nbin_bc,nbdsw))
4663  ALLOCATE ( asyrhd_oc(krhlev,nbin_oc,nbdsw))
4664 
4665  ALLOCATE ( ssarhd_du(krhlev,nbin_du,nbdsw))
4666  ALLOCATE ( ssarhd_ss(krhlev,nbin_ss,nbdsw))
4667  ALLOCATE ( ssarhd_su(krhlev,nbin_su,nbdsw))
4668  ALLOCATE ( ssarhd_bc(krhlev,nbin_bc,nbdsw))
4669  ALLOCATE ( ssarhd_oc(krhlev,nbin_oc,nbdsw))
4670 
4671  ALLOCATE ( extrhd_ni(krhlev,nbin_no3,nbdsw))
4672  ALLOCATE ( scarhd_ni(krhlev,nbin_no3,nbdsw))
4673  ALLOCATE ( asyrhd_ni(krhlev,nbin_no3,nbdsw))
4674  ALLOCATE ( ssarhd_ni(krhlev,nbin_no3,nbdsw))
4675 
4676  if (gocart_on .or. gccpp_on) then
4677  naero=kcm1
4678  else if (nasa_on) then
4679  naero=kcm2
4680  endif
4681  print *, 'aft AEROSOL allocate, nbin_du=',nbin_du, &
4682  'nbin_ss=',nbin_ss,'nbin_su=',nbin_su,'nbin_bc=', &
4683  'nbin_oc=',nbin_oc,'nbin_ni=',nbin_no3,'nAero=',naero
4684 
4685 !!! READ AEROSOL LUTS
4686  DO i = 1, naero
4687  CLOSE(unit=noaer)
4688  if (gocart_on) then
4689  aerosol_file='optics_luts_'//aerosolname(i)//'.dat'
4690  else if ( gccpp_on ) then
4691  aerosol_file='optics_luts_'//aerosolname(i)//'_nasa.dat'
4692  else if ( nasa_on ) then
4693  aerosol_file='optics_luts_'//aerosolname(i)//'_nasa.dat'
4694  endif
4695  open(unit=noaer, file=aerosol_file, status='OLD', iostat=ios)
4696  IF (ios > 0) THEN
4697  print *,' ERROR! Non-zero iostat for rd_LUTS ', aerosol_file
4698  stop
4699  ENDIF
4700  if(debugprint)print *,'i=',i,'read aerosol_file=',trim(aerosol_file),'ios=',ios
4701 !
4702  IF (aerosolname(i) == 'DUST') nbin = nbin_du
4703  IF (aerosolname(i) == 'SALT') nbin = nbin_ss
4704  IF (aerosolname(i) == 'SUSO') nbin = nbin_su
4705  IF (aerosolname(i) == 'SOOT') nbin = nbin_bc
4706  IF (aerosolname(i) == 'WASO') nbin = nbin_oc
4707  if (nasa_on) then
4708  IF (aerosolname(i) == 'NITR') nbin = nbin_no3
4709  endif
4710  DO j = 1, nbin
4711  read(noaer,'(2x,a4,1x,i1,1x,a3)')aername_rd,ib, aeropt
4712  IF (aername_rd /= aerosolname(i)) stop
4713  IF (j /= ib ) stop
4714  IF (aeropt /= 'ext' ) stop
4715 
4716  IF (aerosolname(i) == 'DUST') THEN
4717  do ib = 1, nbdsw
4718  read(noaer,'(8f10.5)') (extrhd_du(ii,j,ib), ii=1,krhlev)
4719  enddo
4720  read(noaer,'(2x,a4)') aername_rd
4721  do ib = 1, nbdsw
4722  read(noaer,'(8f10.5)') (scarhd_du(ii,j,ib), ii=1,krhlev)
4723  enddo
4724  read(noaer,'(2x,a4)') aername_rd
4725  do ib = 1, nbdsw
4726  read(noaer,'(8f10.5)') (asyrhd_du(ii,j,ib), ii=1,krhlev)
4727  enddo
4728  read(noaer,'(2x,a4)') aername_rd
4729  do ib = 1, nbdsw
4730  read(noaer,'(8f10.5)') (ssarhd_du(ii,j,ib), ii=1,krhlev)
4731  enddo
4732 
4733  ELSEIF (aerosolname(i) == 'SALT') THEN
4734  do ib = 1, nbdsw
4735  read(noaer,'(8f10.5)') (extrhd_ss(ii,j,ib), ii=1,krhlev)
4736  enddo
4737  read(noaer,'(2x,a4)') aername_rd
4738  do ib = 1, nbdsw
4739  read(noaer,'(8f10.5)') (scarhd_ss(ii,j,ib), ii=1,krhlev)
4740  enddo
4741  read(noaer,'(2x,a4)') aername_rd
4742  do ib = 1, nbdsw
4743  read(noaer,'(8f10.5)') (asyrhd_ss(ii,j,ib), ii=1,krhlev)
4744  enddo
4745  read(noaer,'(2x,a4)') aername_rd
4746  do ib = 1, nbdsw
4747  read(noaer,'(8f10.5)') (ssarhd_ss(ii,j,ib), ii=1,krhlev)
4748  enddo
4749 
4750  ELSEIF (aerosolname(i) == 'SUSO') THEN
4751  do ib = 1, nbdsw
4752  read(noaer,'(8f10.5)') (extrhd_su(ii,j,ib), ii=1,krhlev)
4753  enddo
4754  read(noaer,'(2x,a4)') aername_rd
4755  do ib = 1, nbdsw
4756  read(noaer,'(8f10.5)') (scarhd_su(ii,j,ib), ii=1,krhlev)
4757  enddo
4758  read(noaer,'(2x,a4)') aername_rd
4759  do ib = 1, nbdsw
4760  read(noaer,'(8f10.5)') (asyrhd_su(ii,j,ib), ii=1,krhlev)
4761  enddo
4762  read(noaer,'(2x,a4)') aername_rd
4763  do ib = 1, nbdsw
4764  read(noaer,'(8f10.5)') (ssarhd_su(ii,j,ib), ii=1,krhlev)
4765  enddo
4766 
4767  ELSEIF (aerosolname(i) == 'SOOT') THEN
4768  do ib = 1, nbdsw
4769  read(noaer,'(8f10.5)') (extrhd_bc(ii,j,ib), ii=1,krhlev)
4770  enddo
4771  read(noaer,'(2x,a4)') aername_rd
4772  do ib = 1, nbdsw
4773  read(noaer,'(8f10.5)') (scarhd_bc(ii,j,ib), ii=1,krhlev)
4774  enddo
4775  read(noaer,'(2x,a4)') aername_rd
4776  do ib = 1, nbdsw
4777  read(noaer,'(8f10.5)') (asyrhd_bc(ii,j,ib), ii=1,krhlev)
4778  enddo
4779  read(noaer,'(2x,a4)') aername_rd
4780  do ib = 1, nbdsw
4781  read(noaer,'(8f10.5)') (ssarhd_bc(ii,j,ib), ii=1,krhlev)
4782  enddo
4783 
4784  ELSEIF (aerosolname(i) == 'WASO') THEN
4785  do ib = 1, nbdsw
4786  read(noaer,'(8f10.5)') (extrhd_oc(ii,j,ib), ii=1,krhlev)
4787  enddo
4788  read(noaer,'(2x,a4)') aername_rd
4789  do ib = 1, nbdsw
4790  read(noaer,'(8f10.5)') (scarhd_oc(ii,j,ib), ii=1,krhlev)
4791  enddo
4792  read(noaer,'(2x,a4)') aername_rd
4793  do ib = 1, nbdsw
4794  read(noaer,'(8f10.5)') (asyrhd_oc(ii,j,ib), ii=1,krhlev)
4795  enddo
4796  read(noaer,'(2x,a4)') aername_rd
4797  do ib = 1, nbdsw
4798  read(noaer,'(8f10.5)') (ssarhd_oc(ii,j,ib), ii=1,krhlev)
4799  enddo
4800  ENDIF
4801  if ( nasa_on ) then
4802  IF (aerosolname(i) == 'NITR') THEN
4803  do ib = 1, nbdsw
4804  read(noaer,'(8f10.5)') (extrhd_ni(ii,j,ib), ii=1,krhlev)
4805  enddo
4806  read(noaer,'(2x,a4)') aername_rd
4807  do ib = 1, nbdsw
4808  read(noaer,'(8f10.5)') (scarhd_ni(ii,j,ib), ii=1,krhlev)
4809  enddo
4810  read(noaer,'(2x,a4)') aername_rd
4811  do ib = 1, nbdsw
4812  read(noaer,'(8f10.5)') (asyrhd_ni(ii,j,ib), ii=1,krhlev)
4813  enddo
4814  read(noaer,'(2x,a4)') aername_rd
4815  do ib = 1, nbdsw
4816  read(noaer,'(8f10.5)') (ssarhd_ni(ii,j,ib), ii=1,krhlev)
4817  enddo
4818  ENDIF
4819  endif !nasa_on
4820 
4821  ENDDO ! j-loop for nbin
4822  ENDDO ! i-loop for nAero
4823 ! print *,'finish reading coef'
4824 
4825  CLOSE(unit=noaer)
4826 
4827 !!! COMPUTES RELATIVE HUMIDITY AND RDRH
4828 ! allocate (RH3D(ista:iend,jsta:jend,lm))
4829  allocate (rdrh(ista:iend,jsta:jend,lm))
4830  allocate (ihh(ista:iend,jsta:jend,lm))
4831  DO l=1,lm ! L FROM TOA TO SFC
4832  ll=lm-l+1 ! LL FROM SFC TO TOA
4833 !$omp parallel do private(i,j)
4834  DO j=jsta,jend
4835  DO i=ista,iend
4836  p1d(i,j) = pmid(i,j,ll)
4837  t1d(i,j) = t(i,j,ll)
4838  q1d(i,j) = q(i,j,ll)
4839  ENDDO
4840  ENDDO
4841  !CALL CALRH(P1D,T1D,Q1D,EGRID4)
4842  CALL calrh(p1d,t1d,q1d,egrid4)
4843  DO j=jsta,jend
4844  DO i=ista,iend
4845 ! RH3D(I,J,LL) = EGRID4(I,J)
4846  rh3d = egrid4(i,j)
4847 ! DETERMINE RDRH (wgt for IH2) and IHH (index for IH2)
4848 ! IF ( RH3D(I,J,LL) > RHLEV(KRHLEV) ) THEN
4849  IF ( rh3d > rhlev(krhlev) ) THEN
4850  ih2 = krhlev
4851  ih1 = ih2 - 1
4852  rdrh(i,j,ll) = 1.
4853 ! ELSEIF ( RH3D(I,J,LL) < RHLEV(1)) THEN
4854  ELSEIF ( rh3d < rhlev(1)) THEN
4855  ih2 = 1
4856  ih1 = 1
4857  rdrh(i,j,ll) = 0.
4858  ELSE
4859  ih2 = 1
4860 ! DO WHILE ( RH3D(I,J,LL) > RHLEV(IH2))
4861  DO WHILE ( rh3d > rhlev(ih2))
4862  ih2 = ih2 + 1
4863  IF ( ih2 > krhlev ) EXIT
4864  ENDDO
4865  ih2 = min( krhlev, ih2 )
4866  ih1 = max( 1, ih2-1 )
4867  drh0 = rhlev(ih2) - rhlev(ih1)
4868 ! DRH1 = RH3D(I,J,LL) - RHLEV(IH1)
4869  drh1 = rh3d - rhlev(ih1)
4870  rdrh(i,j,ll) = drh1 / drh0
4871  ENDIF
4872  ihh(i,j,ll) = ih1
4873 !
4874  ENDDO
4875  ENDDO
4876  ENDDO
4877 
4878 !!!
4879 !!! COMPUTE AOD FOR SPECIFIED WAVELENGTHS
4880  DO ib = 1, nbdsw
4881 
4882 ! AOD AT 340 NM
4883  IF (ib == 1 ) indx = 623
4884 ! AOD AT 440 NM
4885  IF (ib == 2 ) indx = 624
4886 ! AOD AT 550 NM
4887  IF (ib == 3 ) indx = 609
4888 ! AOD AT 660 NM
4889  IF (ib == 4 ) indx = 625
4890 ! AOD AT 860 NM
4891  IF (ib == 5 ) indx = 626
4892 ! AOD AT 1630 NM
4893  IF (ib == 6 ) indx = 627
4894 ! AOD AT 11100 NM
4895  IF (ib == 7 ) indx = 628
4896 
4897 ! DETERMINE LEXT AND LSCA (DEFAULT TO F)
4898  lext = .false.
4899  lsca = .false.
4900  lasy = .false.
4901 ! -- CHECK WHETHER TOTAL EXT AOD IS REQUESTED
4902  IF (iget(indx)>0 ) lext =.true.
4903 ! -- CHECK WHETHER SPECIATED AOD AT 550 NM IS REQUESTED
4904  IF ( ib == 3 ) THEN
4905  IF (iget(650)>0 ) lsca =.true. !TOTAL SCA AOD
4906  DO i = 1, naero
4907  IF (iget(indx_ext(i))>0 ) lext = .true.
4908  IF (iget(indx_sca(i))>0 ) lsca = .true.
4909  ENDDO
4910  ENDIF
4911 ! -- CHECK WHETHER ASY AND SSA AT 340NM IS REQUESTED
4912  IF ( ib == 1 ) THEN
4913  IF (iget(648)>0 ) lsca =.true.
4914  IF (iget(649)>0 ) lasy =.true.
4915  ENDIF
4916 ! -- CHECK WHETHER ANGSTROM EXPONENT IS REQUESTED
4917  IF (iget(656)>0 ) THEN
4918  IF ( ib == 2 ) lext = .true.
4919  IF ( ib == 5 ) lext = .true.
4920  ENDIF
4921 ! print *,'LEXT=',LEXT,'LSCA=',LSCA,'LASY=',LASY
4922 ! SKIP IF POST PRODUCT IS NOT REQUESTED
4923  IF ( lext .OR. lsca .OR. lasy ) THEN
4924 ! COMPUTE DUST AOD
4925  aod_du=spval
4926  sca_du=spval
4927  asy_du=spval
4928  ext=0.0
4929  sca=0.0
4930  asy=0.0
4931  DO j=jsta,jend
4932  DO i=ista,iend
4933  DO l=1,lm
4934  DO n=1, nbin_du
4935  ext01 = extrhd_du(1,n,ib)
4936  sca01 = scarhd_du(1,n,ib)
4937  asy01 = asyrhd_du(1,n,ib)
4938  ext(i,j,l) = ext(i,j,l)+1e-9*dust(i,j,l,n) * ext01
4939  sca(i,j,l) = sca(i,j,l)+1e-9*dust(i,j,l,n) * sca01
4940  asy(i,j,l) = asy(i,j,l)+1e-9*dust(i,j,l,n) * sca01*asy01
4941  ENDDO
4942  ext(i,j,l) = ext(i,j,l) * 1000.
4943  sca(i,j,l) = sca(i,j,l) * 1000.
4944  asy(i,j,l) = asy(i,j,l) * 1000.
4945  ENDDO ! L-loop
4946  ENDDO ! I-loop
4947  ENDDO ! J-loop
4948  CALL calpw(aod_du,17)
4949  CALL calpw(sca_du,20)
4950  CALL calpw(asy_du,21)
4951 ! COMPUTE SULFATE AOD
4952  aod_su=spval
4953  sca_su=spval
4954  asy_su=spval
4955  ext=0.0
4956  sca=0.0
4957  asy=0.0
4958  DO j=jsta,jend
4959  DO i=ista,iend
4960  DO l=1,lm
4961  ih1 = ihh(i,j,l)
4962  ih2 = ih1 + 1
4963  DO n = 1, nbin_su
4964  ext01 = extrhd_su(ih1,n,ib) &
4965  & + rdrh(i,j,l)*(extrhd_su(ih2,n,ib)-extrhd_su(ih1,n,ib))
4966  sca01 = scarhd_su(ih1,n,ib) &
4967  & + rdrh(i,j,l)*(scarhd_su(ih2,n,ib)-scarhd_su(ih1,n,ib))
4968  asy01 = asyrhd_su(ih1,n,ib) &
4969  & + rdrh(i,j,l)*(asyrhd_su(ih2,n,ib)-asyrhd_su(ih1,n,ib))
4970  ext(i,j,l) = ext(i,j,l)+1e-9*suso(i,j,l,n) * ext01
4971  sca(i,j,l) = sca(i,j,l)+1e-9*suso(i,j,l,n)*sca01
4972  asy(i,j,l) = asy(i,j,l)+1e-9*suso(i,j,l,n)*sca01*asy01
4973 
4974  ENDDO ! N-loop
4975  ext(i,j,l) = ext(i,j,l) * 1000.
4976  sca(i,j,l) = sca(i,j,l) * 1000.
4977  asy(i,j,l) = asy(i,j,l) * 1000.
4978  ENDDO ! L-loop
4979  ENDDO ! I-loop
4980  ENDDO ! J-loop
4981  CALL calpw(aod_su,17)
4982  CALL calpw(sca_su,20)
4983  CALL calpw(asy_su,21)
4984 
4985 ! COMPUTE SEA SALT AOD
4986  aod_ss=spval
4987  sca_ss=spval
4988  asy_ss=spval
4989  ext=0.0
4990  sca=0.0
4991  asy=0.0
4992  DO j=jsta,jend
4993  DO i=ista,iend
4994  DO l=1,lm
4995  ih1 = ihh(i,j,l)
4996  ih2 = ih1 + 1
4997  DO n = 1, nbin_ss
4998  ext01 = extrhd_ss(ih1,n,ib) &
4999  & + rdrh(i,j,l)*(extrhd_ss(ih2,n,ib)-extrhd_ss(ih1,n,ib))
5000  sca01 = scarhd_ss(ih1,n,ib) &
5001  & + rdrh(i,j,l)*(scarhd_ss(ih2,n,ib)-scarhd_ss(ih1,n,ib))
5002  asy01 = asyrhd_ss(ih1,n,ib) &
5003  & + rdrh(i,j,l)*(asyrhd_ss(ih2,n,ib)-asyrhd_ss(ih1,n,ib))
5004  ext(i,j,l) = ext(i,j,l)+1e-9*salt(i,j,l,n)*ext01
5005  sca(i,j,l) = sca(i,j,l)+1e-9*salt(i,j,l,n)*sca01
5006  asy(i,j,l) = asy(i,j,l)+1e-9*salt(i,j,l,n)*sca01*asy01
5007  ENDDO ! N-loop
5008  ext(i,j,l) = ext(i,j,l) * 1000.
5009  sca(i,j,l) = sca(i,j,l) * 1000.
5010  asy(i,j,l) = asy(i,j,l) * 1000.
5011  ENDDO ! L-loop
5012  ENDDO ! I-loop
5013  ENDDO ! J-loop
5014  CALL calpw(aod_ss,17)
5015  CALL calpw(sca_ss,20)
5016  CALL calpw(asy_ss,21)
5017 
5018 ! COMPUTE BLACK CARBON AOD
5019  aod_bc=spval
5020  sca_bc=spval
5021  asy_bc=spval
5022  ext=0.0
5023  sca=0.0
5024  asy=0.0
5025  DO j=jsta,jend
5026  DO i=ista,iend
5027  DO l=1,lm
5028  ih1 = ihh(i,j,l)
5029  ih2 = ih1 + 1
5030  DO n = 1, nbin_bc
5031  ext01 = extrhd_bc(ih1,n,ib) &
5032  & + rdrh(i,j,l)*(extrhd_bc(ih2,n,ib)-extrhd_bc(ih1,n,ib))
5033  sca01 = scarhd_bc(ih1,n,ib) &
5034  & + rdrh(i,j,l)*(scarhd_bc(ih2,n,ib)-scarhd_bc(ih1,n,ib))
5035  asy01 = asyrhd_bc(ih1,n,ib) &
5036  & + rdrh(i,j,l)*(asyrhd_bc(ih2,n,ib)-asyrhd_bc(ih1,n,ib))
5037  ext(i,j,l) = ext(i,j,l)+1e-9*soot(i,j,l,n)*ext01
5038  sca(i,j,l) = sca(i,j,l)+1e-9*soot(i,j,l,n)*sca01
5039  asy(i,j,l) = asy(i,j,l)+1e-9*soot(i,j,l,n)*sca01*asy01
5040  ENDDO ! N-loop
5041  ext(i,j,l) = ext(i,j,l) * 1000.
5042  sca(i,j,l) = sca(i,j,l) * 1000.
5043  asy(i,j,l) = asy(i,j,l) * 1000.
5044  ENDDO ! L-loop
5045  ENDDO ! I-loop
5046  ENDDO ! J-loop
5047  CALL calpw(aod_bc,17)
5048  CALL calpw(sca_bc,20)
5049  CALL calpw(asy_bc,21)
5050 ! COMPUTE ORGANIC CARBON AOD
5051  aod_oc=spval
5052  sca_oc=spval
5053  asy_oc=spval
5054  ext=0.0
5055  sca=0.0
5056  asy=0.0
5057  DO j=jsta,jend
5058  DO i=ista,iend
5059  DO l=1,lm
5060  ih1 = ihh(i,j,l)
5061  ih2 = ih1 + 1
5062  DO n = 1, nbin_oc
5063  ext01 = extrhd_oc(ih1,n,ib) &
5064  & + rdrh(i,j,l)*(extrhd_oc(ih2,n,ib)-extrhd_oc(ih1,n,ib))
5065  sca01 = scarhd_oc(ih1,n,ib) &
5066  & + rdrh(i,j,l)*(scarhd_oc(ih2,n,ib)-scarhd_oc(ih1,n,ib))
5067  asy01 = asyrhd_oc(ih1,n,ib) &
5068  & + rdrh(i,j,l)*(asyrhd_oc(ih2,n,ib)-asyrhd_oc(ih1,n,ib))
5069  ext(i,j,l) = ext(i,j,l)+1e-9*waso(i,j,l,n)*ext01
5070  sca(i,j,l) = sca(i,j,l)+1e-9*waso(i,j,l,n)*sca01
5071  asy(i,j,l) = asy(i,j,l)+1e-9*waso(i,j,l,n)*sca01*asy01
5072  ENDDO ! N-loop
5073  ext(i,j,l) = ext(i,j,l) * 1000.
5074  sca(i,j,l) = sca(i,j,l) * 1000.
5075  asy(i,j,l) = asy(i,j,l) * 1000.
5076  ENDDO ! L-loop
5077  ENDDO ! I-loop
5078  ENDDO ! J-loop
5079  CALL calpw(aod_oc,17)
5080  CALL calpw(sca_oc,20)
5081  CALL calpw(asy_oc,21)
5082 
5083  if ( nasa_on ) then
5084 ! COMPUTE ORGANIC CARBON AOD
5085  aod_ni=spval
5086  sca_ni=spval
5087  asy_ni=spval
5088  ext=0.0
5089  sca=0.0
5090  asy=0.0
5091  DO j=jsta,jend
5092  DO i=ista,iend
5093  DO l=1,lm
5094  ih1 = ihh(i,j,l)
5095  ih2 = ih1 + 1
5096  DO n = 1, nbin_no3
5097  ext01 = extrhd_ni(ih1,n,ib) &
5098  & + rdrh(i,j,l)*(extrhd_ni(ih2,n,ib)-extrhd_ni(ih1,n,ib))
5099  sca01 = scarhd_ni(ih1,n,ib) &
5100  & + rdrh(i,j,l)*(scarhd_ni(ih2,n,ib)-scarhd_ni(ih1,n,ib))
5101  asy01 = asyrhd_ni(ih1,n,ib) &
5102  & + rdrh(i,j,l)*(asyrhd_ni(ih2,n,ib)-asyrhd_ni(ih1,n,ib))
5103  ext(i,j,l) = ext(i,j,l)+1e-9*no3(i,j,l,n)*ext01
5104  sca(i,j,l) = sca(i,j,l)+1e-9*no3(i,j,l,n)*sca01
5105  asy(i,j,l) = asy(i,j,l)+1e-9*no3(i,j,l,n)*sca01*asy01
5106  ENDDO ! N-loop
5107  ext(i,j,l) = ext(i,j,l) * 1000.
5108  sca(i,j,l) = sca(i,j,l) * 1000.
5109  asy(i,j,l) = asy(i,j,l) * 1000.
5110  ENDDO ! L-loop
5111  ENDDO ! I-loop
5112  ENDDO ! J-loop
5113  CALL calpw(aod_ni,17)
5114  CALL calpw(sca_ni,20)
5115  CALL calpw(asy_ni,21)
5116  endif
5117 
5118 
5119 ! COMPUTE TOTAL AOD
5120  aod=spval
5121  sca=spval
5122  asy=spval
5123  DO j=jsta,jend
5124  DO i=ista,iend
5125  aod_du(i,j) = max(aod_du(i,j), 0.0)
5126  aod_bc(i,j) = max(aod_bc(i,j), 0.0)
5127  aod_oc(i,j) = max(aod_oc(i,j), 0.0)
5128  aod_su(i,j) = max(aod_su(i,j), 0.0)
5129  aod_ss(i,j) = max(aod_ss(i,j), 0.0)
5130 
5131  sca_du(i,j) = max(sca_du(i,j), 0.0)
5132  sca_bc(i,j) = max(sca_bc(i,j), 0.0)
5133  sca_oc(i,j) = max(sca_oc(i,j), 0.0)
5134  sca_su(i,j) = max(sca_su(i,j), 0.0)
5135  sca_ss(i,j) = max(sca_ss(i,j), 0.0)
5136 
5137  asy_du(i,j) = max(asy_du(i,j), 0.0)
5138  asy_bc(i,j) = max(asy_bc(i,j), 0.0)
5139  asy_oc(i,j) = max(asy_oc(i,j), 0.0)
5140  asy_su(i,j) = max(asy_su(i,j), 0.0)
5141  asy_ss(i,j) = max(asy_ss(i,j), 0.0)
5142 
5143  if ( nasa_on ) then
5144  aod_ni(i,j) = max(aod_ni(i,j), 0.0)
5145  sca_ni(i,j) = max(sca_ni(i,j), 0.0)
5146  asy_ni(i,j) = max(asy_ni(i,j), 0.0)
5147 
5148  aod(i,j) = aod_du(i,j) + aod_bc(i,j) + aod_oc(i,j) + &
5149  & aod_su(i,j) + aod_ss(i,j) + aod_ni(i,j)
5150  sca2d(i,j) = sca_du(i,j) + sca_bc(i,j) + sca_oc(i,j) + &
5151  & sca_su(i,j) + sca_ss(i,j) + sca_ni(i,j)
5152  asy2d(i,j) = asy_du(i,j) + asy_bc(i,j) + asy_oc(i,j) + &
5153  & asy_su(i,j) + asy_ss(i,j) + asy_ni(i,j)
5154  endif
5155 
5156  if (gocart_on .or. gccpp_on) then
5157  aod(i,j) = aod_du(i,j) + aod_bc(i,j) + aod_oc(i,j) + &
5158  & aod_su(i,j) + aod_ss(i,j)
5159  sca2d(i,j) = sca_du(i,j) + sca_bc(i,j) + sca_oc(i,j) + &
5160  & sca_su(i,j) + sca_ss(i,j)
5161  asy2d(i,j) = asy_du(i,j) + asy_bc(i,j) + asy_oc(i,j) + &
5162  & asy_su(i,j) + asy_ss(i,j)
5163  endif
5164 
5165  ENDDO ! I-loop
5166  ENDDO ! J-loop
5167 ! FILL UP AOD_440 AND AOD_860, IF ANGSTROM EXP IS REQUESTED
5168  IF ( iget(656) > 0 ) THEN
5169  IF (ib == 2 ) THEN !! AOD AT 440 NM
5170 !$omp parallel do private(i,j)
5171  DO j=jsta,jend
5172  DO i=ista,iend
5173  aod_440(i,j) = aod(i,j)
5174  ENDDO ! I-loop
5175  ENDDO ! J-loop
5176  ENDIF
5177 
5178  IF (ib == 5 ) THEN !! AOD AT 860 NM
5179 !$omp parallel do private(i,j)
5180  DO j=jsta,jend
5181  DO i=ista,iend
5182  aod_860(i,j) = aod(i,j)
5183  ENDDO ! I-loop
5184  ENDDO ! J-loop
5185  ENDIF
5186  ENDIF
5187 
5188 ! WRITE OUT TOTAL AOD
5189  IF ( iget(indx) > 0) THEN
5190 !$omp parallel do private(i,j)
5191  do j=jsta,jend
5192  do i=ista,iend
5193  grid1(i,j) = aod(i,j)
5194  enddo
5195  enddo
5196  CALL bound(grid1,d00,h99999)
5197  if(grib=="grib2" )then
5198  cfld=cfld+1
5199  fld_info(cfld)%ifld=iavblfld(iget(indx))
5200  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5201  endif
5202  ENDIF
5203 !
5204 ! WRITE OUT ASY AND SSA AT 340NM
5205  IF ( ib == 1 ) THEN !!! FOR 340NM ONLY
5206 
5207 ! AER ASYM FACTOR AT 340 NM
5208  IF ( iget(649) > 0 ) THEN
5209  grid1 = spval
5210 !$omp parallel do private(i,j)
5211  DO j=jsta,jend
5212  DO i=ista,iend
5213  IF(sca2d(i,j)<spval.and.asy2d(i,j)<spval) THEN
5214  IF ( sca2d(i,j) > 0.0 ) THEN
5215  asy2d(i,j) = asy2d(i,j) / sca2d(i,j)
5216  ELSE
5217  asy2d(i,j) = 0.
5218  ENDIF
5219  IF(asy2d(i,j)<spval) grid1(i,j)=asy2d(i,j)
5220  ENDIF
5221  ENDDO
5222  ENDDO
5223  CALL bound(grid1,d00,h99999)
5224  if(grib=="grib2" )then
5225  cfld=cfld+1
5226  fld_info(cfld)%ifld=iavblfld(iget(649))
5227  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5228  endif
5229  ENDIF ! IGET(649)
5230 
5231 ! AER SINGLE SCATTER ALB AT 340 NM
5232  IF ( iget(648) > 0 ) THEN
5233  grid1 = spval
5234 !$omp parallel do private(i,j)
5235  DO j=jsta,jend
5236  DO i=ista,iend
5237  IF(aod(i,j)<spval.and.sca2d(i,j)<spval) THEN
5238  IF ( aod(i,j) > 0.0 ) THEN
5239  sca2d(i,j) = sca2d(i,j) / aod(i,j)
5240  ELSE
5241  sca2d(i,j) = 1.0
5242  ENDIF
5243  IF(sca2d(i,j)<spval) grid1(i,j)=sca2d(i,j)
5244  ENDIF
5245  ENDDO
5246  ENDDO
5247  CALL bound(grid1,d00,h99999)
5248  if(grib=="grib2" )then
5249  cfld=cfld+1
5250  fld_info(cfld)%ifld=iavblfld(iget(648))
5251  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5252  endif
5253  ENDIF ! IGET(648)
5254 ! print *,'aft compute sca340'
5255 
5256  ENDIF ! IB IF-BLOCK (340NM)
5257 
5258 
5259 ! WRITE OUT AOD FOR DU, SU, SS, OC, BC for all wavelengths
5260 ! WRITE OUT SPECIATED AEROSOL OPTICAL PROPERTIES
5261  IF ( ib == 3 ) THEN !!! FOR 550NM ONLY
5262 
5263 ! WRITE OUT TOTAL SCATTERING AOD
5264  IF ( iget(650) > 0 ) THEN
5265 !$omp parallel do private(i,j)
5266  DO j=jsta,jend
5267  DO i=ista,iend
5268  grid1(i,j)=sca2d(i,j)
5269  ENDDO
5270  ENDDO
5271  CALL bound(grid1,d00,h99999)
5272  if(grib=="grib2" )then
5273  cfld=cfld+1
5274  fld_info(cfld)%ifld=iavblfld(iget(650))
5275  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5276  endif
5277  ENDIF
5278 ! LOOP THROUGH EACH SPECIES
5279  DO ii = 1, naero
5280 
5281 ! WRITE OUT EXT AOD
5282  jj = indx_ext(ii)
5283  IF ( iget(jj) > 0) THEN ! EXT AOD
5284 !$omp parallel do private(i,j)
5285  DO j=jsta,jend
5286  DO i=ista,iend
5287  IF ( ii == 1 ) grid1(i,j) = aod_du(i,j)
5288  IF ( ii == 2 ) grid1(i,j) = aod_ss(i,j)
5289  IF ( ii == 3 ) grid1(i,j) = aod_su(i,j)
5290  IF ( ii == 4 ) grid1(i,j) = aod_oc(i,j)
5291  IF ( ii == 5 ) grid1(i,j) = aod_bc(i,j)
5292  IF ( ii == 6 ) grid1(i,j) = aod_ni(i,j)
5293  ENDDO
5294  ENDDO
5295  CALL bound(grid1,d00,h99999)
5296  if(grib=="grib2" )then
5297  cfld=cfld+1
5298  fld_info(cfld)%ifld=iavblfld(iget(jj))
5299  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5300  endif
5301  ENDIF
5302 
5303 ! WRITE OUT SCA AOD
5304  jj = indx_sca(ii)
5305  IF ( iget(jj) > 0) THEN ! SCA AOD
5306 !$omp parallel do private(i,j)
5307  DO j=jsta,jend
5308  DO i=ista,iend
5309  IF ( ii == 1 ) grid1(i,j) = sca_du(i,j)
5310  IF ( ii == 2 ) grid1(i,j) = sca_ss(i,j)
5311  IF ( ii == 3 ) grid1(i,j) = sca_su(i,j)
5312  IF ( ii == 4 ) grid1(i,j) = sca_oc(i,j)
5313  IF ( ii == 5 ) grid1(i,j) = sca_bc(i,j)
5314  IF ( ii == 6 ) grid1(i,j) = sca_ni(i,j)
5315  ENDDO
5316  ENDDO
5317  CALL bound(grid1,d00,h99999)
5318  if(grib=="grib2" )then
5319  cfld=cfld+1
5320  fld_info(cfld)%ifld=iavblfld(iget(jj))
5321  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5322  endif
5323  ENDIF
5324 
5325  ENDDO ! II DO-LOOP
5326 
5327  ENDIF ! IB IF-BLOCK (550NM)
5328 
5329  ENDIF ! LEXT IF-BLOCK
5330  ENDDO ! LOOP THROUGH NBDSW CHANNELS
5331 ! COMPUTE AND WRITE OUT ANGSTROM EXPONENT
5332  IF ( iget(656) > 0 ) THEN
5333  angst=spval
5334 ! ANG2 = LOG ( 0.860 / 0.440 )
5335  ang2 = log( 860. / 440. )
5336 !$omp parallel do private(i,j,ang1)
5337  DO j=jsta,jend
5338  DO i=ista,iend
5339  IF (aod_860(i,j) > 0.) THEN
5340  ang1 = log( aod_440(i,j)/aod_860(i,j) )
5341  angst(i,j) = ang1 / ang2
5342  ENDIF
5343  grid1(i,j)=angst(i,j)
5344  ENDDO
5345  ENDDO
5346  if(debugprint)print *,'output angstrom exp,angst=',maxval(angst(ista:iend,jsta:jend)), &
5347  minval(angst(ista:iend,jsta:jend))
5348  CALL bound(grid1,d00,h99999)
5349  if(grib=="grib2" )then
5350  cfld=cfld+1
5351  fld_info(cfld)%ifld=iavblfld(iget(656))
5352  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5353  endif
5354  ENDIF ! ANGSTROM EXPONENT
5355 
5356  ENDIF ! END OF LAEROPT IF-BLOCK
5357 
5358 
5359 !! ADD AEROSOL SURFACE PM25 DUST MASS CONCENTRATION (ug/m3)
5360  IF (iget(686)>0 ) THEN
5361 !$omp parallel do private(i,j)
5362  DO j = jsta,jend
5363  DO i = ista,iend
5364  !GRID1(I,J) = DUSMASS(I,J) * 1.E-6
5365  grid1(i,j) = dustpm(i,j) !ug/m3
5366  END DO
5367  END DO
5368  if(grib=='grib2') then
5369  cfld=cfld+1
5370  fld_info(cfld)%ifld=iavblfld(iget(686))
5371  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5372  endif
5373  ENDIF
5374 
5375  IF (iget(685)>0 ) THEN
5376 !$omp parallel do private(i,j)
5377  DO j = jsta,jend
5378  DO i = ista,iend
5379  grid1(i,j) = dustpm10(i,j) !ug/m3
5380  END DO
5381  END DO
5382  if(grib=='grib2') then
5383  cfld=cfld+1
5384  fld_info(cfld)%ifld=iavblfld(iget(685))
5385  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5386  endif
5387  ENDIF
5388 
5389 
5390 !! ADD AEROSOL SURFACE PM25 SEA SALT MASS CONCENTRATION (ug/m3)
5391  IF (iget(684)>0 ) THEN
5392 !$omp parallel do private(i,j)
5393  DO j = jsta,jend
5394  DO i = ista,iend
5395  !GRID1(I,J) = DUSMASS(I,J) * 1.E-6
5396  grid1(i,j) = sspm(i,j) !ug/m3
5397  END DO
5398  END DO
5399  if(grib=='grib2') then
5400  cfld=cfld+1
5401  fld_info(cfld)%ifld=iavblfld(iget(684))
5402  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5403  endif
5404  ENDIF
5405 !! ADD AEROSOL SURFACE PM10 MASS CONCENTRATION (ug/m3)
5406  IF (iget(619)>0 ) THEN
5407 !$omp parallel do private(i,j)
5408  DO j = jsta,jend
5409  DO i = ista,iend
5410  !GRID1(I,J) = DUSMASS(I,J) * 1.E-6
5411  grid1(i,j) = dusmass(i,j) !ug/m3
5412  END DO
5413  END DO
5414  if(grib=='grib2') then
5415  cfld=cfld+1
5416  fld_info(cfld)%ifld=iavblfld(iget(619))
5417  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5418  endif
5419  ENDIF
5420 
5421 !! ADD AEROSOL SURFACE PM2.5 MASS CONCENTRATION (ug/m3)
5422  IF (iget(620)>0 ) THEN
5423 !$omp parallel do private(i,j)
5424  DO j = jsta,jend
5425  DO i = ista,iend
5426  !GRID1(I,J) = DUSMASS25(I,J) * 1.E-6
5427  grid1(i,j) = dusmass25(i,j) ! ug/m3
5428  END DO
5429  END DO
5430  if(grib=='grib2') then
5431  cfld=cfld+1
5432  fld_info(cfld)%ifld=iavblfld(iget(620))
5433  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5434  endif
5435  ENDIF
5436 !! ADD TOTAL AEROSOL PM10 COLUMN DENSITY (kg/m2) !
5437  IF (iget(621)>0 ) THEN
5438  grid1=spval
5439 !$omp parallel do private(i,j)
5440  DO j = jsta,jend
5441  DO i = ista,iend
5442  !GRID1(I,J) = DUCMASS(I,J) * 1.E-6
5443  IF(ducmass(i,j)<spval) grid1(i,j) = ducmass(i,j) * 1.e-9
5444  END DO
5445  END DO
5446  if(grib=='grib2') then
5447  cfld=cfld+1
5448  fld_info(cfld)%ifld=iavblfld(iget(621))
5449  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5450  endif
5451  ENDIF
5452 
5453 !! ADD TOTAL AEROSOL PM2.5 COLUMN DENSITY (kg/m2)
5454  IF (iget(622)>0 ) THEN
5455  grid1=spval
5456 !$omp parallel do private(i,j)
5457  DO j = jsta,jend
5458  DO i = ista,iend
5459  !GRID1(I,J) = DUCMASS25(I,J) * 1.E-6
5460  IF(ducmass25(i,j)<spval) grid1(i,j) = ducmass25(i,j) * 1.e-9
5461  END DO
5462  END DO
5463  if(grib=='grib2') then
5464  cfld=cfld+1
5465  fld_info(cfld)%ifld=iavblfld(iget(622))
5466  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5467  endif
5468  ENDIF
5469 
5470 !! ADD DUST PM2.5 COLUMN DENSITY (kg/m2)
5471  IF (iget(646)>0 ) THEN
5472  grid1=spval
5473 !$omp parallel do private(i,j)
5474  DO j = jsta,jend
5475  DO i = ista,iend
5476  IF(dustcb(i,j)<spval) grid1(i,j) = dustcb(i,j) * 1.e-9
5477  END DO
5478  END DO
5479  if(grib=='grib2') then
5480  cfld=cfld+1
5481  fld_info(cfld)%ifld=iavblfld(iget(646))
5482  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5483  endif
5484  ENDIF
5485 
5486 !! ADD SEA SALT PM2.5 COLUMN DENSITY (kg/m2)
5487  IF (iget(647)>0 ) THEN
5488  grid1=spval
5489 !$omp parallel do private(i,j)
5490  DO j = jsta,jend
5491  DO i = ista,iend
5492  IF(sscb(i,j)<spval) grid1(i,j) = sscb(i,j) * 1.e-9
5493  END DO
5494  END DO
5495  if(grib=='grib2') then
5496  cfld=cfld+1
5497  fld_info(cfld)%ifld=iavblfld(iget(647))
5498  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5499  endif
5500  ENDIF
5501 !! ADD BC COLUMN DENSITY (kg/m2)
5502  IF (iget(616)>0 ) THEN
5503  grid1=spval
5504 !$omp parallel do private(i,j)
5505  DO j = jsta,jend
5506  DO i = ista,iend
5507  IF(bccb(i,j)<spval) grid1(i,j) = bccb(i,j) * 1.e-9
5508  END DO
5509  END DO
5510  if(grib=='grib2') then
5511  cfld=cfld+1
5512  fld_info(cfld)%ifld=iavblfld(iget(616))
5513  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5514  endif
5515  ENDIF
5516 
5517 !! ADD OC COLUMN DENSITY (kg/m2) !
5518  IF (iget(617)>0 ) THEN
5519  grid1=spval
5520 !$omp parallel do private(i,j)
5521  DO j = jsta,jend
5522  DO i = ista,iend
5523  IF(occb(i,j)<spval) grid1(i,j) = occb(i,j) * 1.e-9
5524  END DO
5525  END DO
5526  if(grib=='grib2') then
5527  cfld=cfld+1
5528  fld_info(cfld)%ifld=iavblfld(iget(617))
5529  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5530  endif
5531  ENDIF
5532 
5533 !! ADD SULF COLUMN DENSITY (kg/m2) !
5534  IF (iget(618)>0 ) THEN
5535  grid1=spval
5536 !$omp parallel do private(i,j)
5537  DO j = jsta,jend
5538  DO i = ista,iend
5539  IF(sulfcb(i,j)<spval) grid1(i,j) = sulfcb(i,j) * 1.e-9
5540  END DO
5541  END DO
5542  if(grib=='grib2') then
5543  cfld=cfld+1
5544  fld_info(cfld)%ifld=iavblfld(iget(618))
5545  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5546  endif
5547  ENDIF
5548 
5549  if (nasa_on) then
5550 !! ADD NO3 COLUMN DENSITY (kg/m2) !
5551  IF (iget(657)>0 ) THEN
5552  grid1=spval
5553 !$omp parallel do private(i,j)
5554  DO j = jsta,jend
5555  DO i = ista,iend
5556  IF(no3cb(i,j)<spval) grid1(i,j) = no3cb(i,j) * 1.e-9
5557  END DO
5558  END DO
5559  if(grib=='grib2') then
5560  cfld=cfld+1
5561  fld_info(cfld)%ifld=iavblfld(iget(657))
5562  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5563  endif
5564  ENDIF
5565 
5566 !! ADD NH4 COLUMN DENSITY (kg/m2) !
5567  IF (iget(658)>0 ) THEN
5568  grid1=spval
5569 !$omp parallel do private(i,j)
5570  DO j = jsta,jend
5571  DO i = ista,iend
5572  IF(nh4cb(i,j)<spval) grid1(i,j) = nh4cb(i,j) * 1.e-9
5573  END DO
5574  END DO
5575  if(grib=='grib2') then
5576  cfld=cfld+1
5577  fld_info(cfld)%ifld=iavblfld(iget(658))
5578  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5579  endif
5580  ENDIF
5581 
5582  endif !nasa_on
5583 
5584  if (gocart_on .or. gccpp_on ) then
5585 !! ADD EMISSION FLUXES,dry depostion, wet/convective depostion (kg/m2/sec)
5586 !! The AER file uses 1.E6 to scale all 2d diagnosis fields
5587 !! Multiply by 1.E-6 to revert these fields back
5588  IF (iget(659)>0) call wrt_aero_diag(659,nbin_du,duem)
5589 ! print *,'aft wrt disg duem'
5590  IF (iget(660)>0) call wrt_aero_diag(660,nbin_du,dusd)
5591  IF (iget(661)>0) call wrt_aero_diag(661,nbin_du,dudp)
5592  IF (iget(662)>0) call wrt_aero_diag(662,nbin_du,duwt)
5593  IF (iget(679)>0) call wrt_aero_diag(679,nbin_du,dusv)
5594 ! print *,'aft wrt disg duwt'
5595 
5596 !! wrt SS diag field
5597  IF (iget(663)>0) call wrt_aero_diag(663,nbin_ss,ssem)
5598  IF (iget(664)>0) call wrt_aero_diag(664,nbin_ss,sssd)
5599  IF (iget(665)>0) call wrt_aero_diag(665,nbin_ss,ssdp)
5600  IF (iget(666)>0) call wrt_aero_diag(666,nbin_ss,sswt)
5601  IF (iget(680)>0) call wrt_aero_diag(680,nbin_ss,sssv)
5602 ! print *,'aft wrt disg sswt'
5603 
5604 !! wrt BC diag field
5605  IF (iget(667)>0) call wrt_aero_diag(667,nbin_bc,bcem)
5606  IF (iget(668)>0) call wrt_aero_diag(668,nbin_bc,bcsd)
5607  IF (iget(669)>0) call wrt_aero_diag(669,nbin_bc,bcdp)
5608  IF (iget(670)>0) call wrt_aero_diag(670,nbin_bc,bcwt)
5609  IF (iget(681)>0) call wrt_aero_diag(681,nbin_bc,bcsv)
5610 ! print *,'aft wrt disg bcwt'
5611 
5612 !! wrt OC diag field
5613  IF (iget(671)>0) call wrt_aero_diag(671,nbin_oc,ocem)
5614  IF (iget(672)>0) call wrt_aero_diag(672,nbin_oc,ocsd)
5615  IF (iget(673)>0) call wrt_aero_diag(673,nbin_oc,ocdp)
5616  IF (iget(674)>0) call wrt_aero_diag(674,nbin_oc,ocwt)
5617  IF (iget(682)>0) call wrt_aero_diag(682,nbin_oc,ocsv)
5618 ! print *,'aft wrt disg ocwt'
5619 !! wrt MIE AOD at 550nm
5620  IF (iget(699).GT.0) call wrt_aero_diag(699,1,maod)
5621  print *,'aft wrt disg maod'
5622  endif !gocart_on
5623 
5624 !! wrt SU diag field
5625 ! IF (IGET(675)>0) call wrt_aero_diag(675,nbin_su,suem)
5626 ! IF (IGET(676)>0) call wrt_aero_diag(676,nbin_su,susd)
5627 ! IF (IGET(677)>0) call wrt_aero_diag(677,nbin_su,sudp)
5628 ! IF (IGET(678)>0) call wrt_aero_diag(678,nbin_su,suwt)
5629 ! print *,'aft wrt disg suwt'
5630  endif ! if gocart_on or nasa_on
5631 
5632 ! CB for WAFS
5633  if(iget(473)>0 .or. iget(474)>0 .or. iget(475)>0) then
5634  ! CB cover is derived from CPRAT (same as #272 in SURFCE.f)
5635  egrid1 = spval
5636  DO j=jsta,jend
5637  DO i=ista,iend
5638  if(avgcprate(i,j) /= spval) then
5639  egrid1(i,j) = avgcprate(i,j)*(1000./dtq2)
5640  end if
5641  END DO
5642  END DO
5643  call cb_cover(egrid1)
5644 
5645  ! CB base(height):derived from convective cloud base (pressure, same as #188 in CLDRAD.f)
5646  ! CB top(height): derived from convective cloud top (pressure, same as #189 in CLDRAD.f)
5647  egrid2 = spval
5648  egrid3 = spval
5649  IF(modelname == 'GFS' .OR. modelname == 'FV3R') then
5650  DO j=jsta,jend
5651  DO i=ista,iend
5652  egrid2(i,j) = pbot(i,j)
5653  egrid3(i,j) = ptop(i,j)
5654  END DO
5655  END DO
5656  END IF
5657 
5658  ! Derive CB base and top, relationship among CB fields
5659  DO j=jsta,jend
5660  DO i=ista,iend
5661  if(egrid1(i,j)<= 0. .or. egrid2(i,j)<= 0. .or. egrid3(i,j) <= 0.) then
5662  egrid1(i,j) = spval
5663  egrid2(i,j) = spval
5664  egrid3(i,j) = spval
5665  end if
5666  END DO
5667  END DO
5668  DO j=jsta,jend
5669  DO i=ista,iend
5670  IF(egrid2(i,j) == spval .or. egrid3(i,j) == spval) cycle
5671  if(egrid3(i,j) < 400.*100. .and. &
5672  (egrid2(i,j)-egrid3(i,j)) > 300.*100) then
5673  ! Convert PBOT to height
5674  if(egrid2(i,j) > pmid(i,j,lm)) then
5675  egrid2(i,j) = 0.
5676  else
5677  do l = lm-1, 1, -1
5678  if(egrid2(i,j) >= pmid(i,j,l)) then
5679  if(egrid2(i,j)-pmid(i,j,l)<0.5) then
5680  egrid2(i,j) = zmid(i,j,l)
5681  else
5682  dp = (log(egrid2(i,j)) - log(pmid(i,j,l)))/ &
5683  max(1.e-6,(log(pmid(i,j,l+1))-log(pmid(i,j,l))))
5684  egrid2(i,j) = zmid(i,j,l)+(zmid(i,j,l+1)-zmid(i,j,l))*dp
5685  end if
5686  exit
5687  end if
5688  end do
5689  end if
5690  ! Convert PTOP to height
5691  if(egrid3(i,j) < pmid(i,j,1)) then
5692  egrid3(i,j) = zmid(i,j,1)
5693  else
5694  do l = 2, lm
5695  if(egrid3(i,j) <= pmid(i,j,l)) then
5696  if(pmid(i,j,l)-egrid3(i,j)<0.5) then
5697  egrid3(i,j) = zmid(i,j,l)
5698  else
5699  dp = (log(egrid3(i,j)) - log(pmid(i,j,l)))/ &
5700  max(1.e-6,(log(pmid(i,j,l))-log(pmid(i,j,l-1))))
5701  egrid3(i,j) = zmid(i,j,l)+(zmid(i,j,l)-zmid(i,j,l-1))*dp
5702  end if
5703  exit
5704  end if
5705  end do
5706  end if
5707  else
5708  egrid1(i,j) = spval
5709  egrid2(i,j) = spval
5710  egrid3(i,j) = spval
5711  end if
5712  END DO
5713  END DO
5714 
5715  IF(iget(473) > 0) THEN
5716 !$omp parallel do private(i,j)
5717  DO j=jsta,jend
5718  DO i=ista,iend
5719  grid1(i,j) = egrid1(i,j)
5720  ENDDO
5721  ENDDO
5722  cfld=cfld+1
5723  fld_info(cfld)%ifld=iavblfld(iget(473))
5724 !$omp parallel do private(i,j,ii,jj)
5725  do j=1,jend-jsta+1
5726  jj = jsta+j-1
5727  do i=1,iend-ista+1
5728  ii=ista+i-1
5729  datapd(i,j,cfld) = grid1(ii,jj)
5730  enddo
5731  enddo
5732  END IF
5733 
5734  IF(iget(474) > 0) THEN
5735 !$omp parallel do private(i,j)
5736  DO j=jsta,jend
5737  DO i=ista,iend
5738  grid1(i,j) = egrid2(i,j)
5739  ENDDO
5740  ENDDO
5741  cfld=cfld+1
5742  fld_info(cfld)%ifld=iavblfld(iget(474))
5743 !$omp parallel do private(i,j,ii,jj)
5744  do j=1,jend-jsta+1
5745  jj = jsta+j-1
5746  do i=1,iend-ista+1
5747  ii=ista+i-1
5748  datapd(i,j,cfld) = grid1(ii,jj)
5749  enddo
5750  enddo
5751  END IF
5752 
5753  IF(iget(475) > 0) THEN
5754 !$omp parallel do private(i,j)
5755  DO j=jsta,jend
5756  DO i=ista,iend
5757  grid1(i,j) = egrid3(i,j)
5758  ENDDO
5759  ENDDO
5760  cfld=cfld+1
5761  fld_info(cfld)%ifld=iavblfld(iget(475))
5762 !$omp parallel do private(i,j,ii,jj)
5763  do j=1,jend-jsta+1
5764  jj = jsta+j-1
5765  do i=1,iend-ista+1
5766  ii=ista+i-1
5767  datapd(i,j,cfld) = grid1(ii,jj)
5768  enddo
5769  enddo
5770  END IF
5771  end if
5772 
5773 !
5774 ! END OF ROUTINE.
5775 !
5776  RETURN
5777  END
5778 
5779  subroutine cb_cover(cbcov)
5780 
5783  use ctlblk_mod, only: spval,jsta,jend,im,ista,iend
5784  implicit none
5785  real, intent(inout) :: cbcov(ista:iend,jsta:jend)
5786 
5787  ! x - convective precipitation [1.0e6*kg/(m2s)]
5788  ! y - cloud cover fraction, between 0 and 1
5789  ! These are original values from Slingo (Table 1):
5790  ! c = -.006 + 0.125*log(p)
5791  ! x = 1.6 3.6 8.1 18.5 39.0 89.0 197.0 440.0 984.0 10000.0
5792  ! y = 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.8
5793  integer, parameter :: np=10
5794  real :: x(np), y(np)
5795 
5796  integer :: i, j, k
5797  real :: val, delta
5798 
5799  x = (/ 1.6,3.6,8.1,18.5,39.0,89.0,197.0,440.0,984.0,10000.0 /)
5800  y = (/ 0.0,0.1,0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.8 /)
5801 
5802  x = log(x)
5803 
5804  do j = jsta, jend
5805  do i = ista, iend
5806  if(cbcov(i,j) == spval) cycle
5807  if(cbcov(i,j) <= 0.) then
5808  cbcov(i,j) = 0.
5809  else
5810  val=log(1.0e6*cbcov(i,j))
5811  if (val <= x(1)) then
5812  cbcov(i,j) = 0.0
5813  else if (val >= x(np)) then
5814  cbcov(i,j) = 0.0
5815  else
5816  do k = 2, np
5817  if (val < x(k)) then
5818  delta = x(k) - x(k-1)
5819  if (delta <= 0.0) then
5820  cbcov(i,j) = y(k-1)
5821  else
5822  cbcov(i,j) = (y(k) * (val-x(k-1)) + &
5823  y(k-1) * (x(k)-val)) / delta
5824  end if
5825  exit
5826  end if
5827  end do
5828  end if
5829  end if
5830  end do
5831  end do
5832  end subroutine cb_cover
5833 
5834  subroutine wrt_aero_diag(igetfld,nbin,data)
5835  use ctlblk_mod, only: jsta, jend, spval, im, jm, grib, &
5836  cfld, datapd, fld_info, jsta_2l, jend_2u,ista_2l,iend_2u,ista,iend
5837  use rqstfld_mod, only: iget, id, lvls, iavblfld
5838  implicit none
5839 !
5840  integer igetfld,nbin
5841  real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u,nbin) :: data
5842 !
5843  integer i,j,k
5844  REAL,dimension(im,jm) :: grid1
5845 !
5846  grid1=spval
5847 !$omp parallel do private(i,j,k)
5848  DO j = jsta,jend
5849  DO i = ista,iend
5850  if(data(i,j,1)<spval) grid1(i,j) = data(i,j,1)
5851  DO k=2,nbin
5852  if(data(i,j,k)<spval)&
5853  grid1(i,j) = grid1(i,j)+ data(i,j,k)
5854  END DO
5855  END DO
5856  END DO
5857  if(grib=='grib2') then
5858  cfld=cfld+1
5859  fld_info(cfld)%ifld=iavblfld(iget(igetfld))
5860  datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
5861  endif
5862 
5863  end subroutine wrt_aero_diag
ctlblk_mod
Definition: CTLBLK.f:1
rqstfld_mod
Definition: RQSTFLD.f:1
calfltcnd
subroutine calfltcnd(CEILING, FLTCND)
Computes Ceiling.
Definition: AVIATION.f:550
masks
Definition: MASKS_mod.f:1
params_mod
Definition: params.F:1
otlift
subroutine otlift(SLINDX)
otlift() computes SFC to 500mb lifted index.
Definition: OTLIFT.f:38
bound
subroutine bound(FLD, FMIN, FMAX)
This routine bounds data in the passed array FLD (im x jm elements long) and clips data values such t...
Definition: BOUND.f:29
calceiling
subroutine calceiling(CLDZ, TCLD, CEILING)
Computes ceiling.
Definition: AVIATION.f:496
upp_physics::calrh
subroutine, public calrh(P1, T1, Q1, RH)
CALRH() computes relative humidity.
Definition: UPP_PHYSICS.f:72
cmassi_mod
Definition: CMASSI.f:1
cb_cover
subroutine cb_cover(cbcov)
Definition: CLDRAD.f:5779
upp_physics::calcape
subroutine, public calcape(ITYPE, DPBND, P1D, T1D, Q1D, L1D, CAPE, CINS, PPARC, ZEQL, THUND)
calcape() computes CAPE and CINS.
Definition: UPP_PHYSICS.f:562
machine_post
Definition: machine.f:1
upp_physics
Definition: UPP_PHYSICS.f:40