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