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