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