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