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