UPP (develop)
Loading...
Searching...
No Matches
MDLFLD.f
Go to the documentation of this file.
1
2!
99 SUBROUTINE mdlfld
100
101!
102 use vrbls4d, only: dust, salt, suso, waso, soot, no3, nh4, smoke, fv3dust,&
103 coarsepm, ebb
104 use vrbls3d, only: zmid, t, pmid, q, cwm, f_ice, f_rain, f_rimef, qqw, qqi,&
105 qqr, qqs, cfr, cfr_raw, dbz, dbzr, dbzi, dbzc, qqw, nlice, nrain, qqg, qqh, zint,&
106 qqni, qqnr, qqnw, qqnwfa, qqnifa, uh, vh, mcvg, omga, wh, q2, ttnd, rswtt, &
107 rlwtt, train, tcucn, o3, rhomid, dpres, el_pbl, pint, icing_gfip, icing_gfis, &
108 catedr,mwt,gtg,cit, ref_10cm, avgpmtf, avgozcon
109
110 use vrbls2d, only: slp, hbot, htop, cnvcfr, cprate, cnvcfr, sfcshx,sfclhx,ustar,z0,&
111 sr, prec, vis, czen, pblh, pblhgust, u10, v10, avgprec, avgcprate, &
112 ref1km_10cm,ref4km_10cm,refc_10cm,refd_max
113 use masks, only: lmh, gdlat, gdlon,sm,sice,dx,dy
114 use params_mod, only: rd, gi, g, rog, h1, tfrz, d00, dbzmin, d608, small,&
115 h100, h1m12, h99999,pi,erad
116 use pmicrph_mod, only: r1, const1r, qr0, delqr0, const2r, ron, topr, son,&
117 tops, dsnow, drain,const_ng1, const_ng2, gon, topg, dgraupel
118 use ctlblk_mod, only: jsta_2l, jend_2u, lm, jsta, jend, grib, cfld, datapd,&
119 fld_info, modelname, imp_physics, dtq2, spval, icount_calmict,&
120 me, dt, avrain, theat, ifhr, ifmin, avcnvc, lp1, im, jm, &
121 ista, iend, ista_2l, iend_2u, aqf_on, gocart_on, gccpp_on, nasa_on, gtg_on
122 use rqstfld_mod, only: iget, id, lvls, iavblfld, lvlsxml
123 use gridspec_mod, only: gridtype,maptype,dxval
124 use upp_physics, only: calrh, calcape, calvor
125 use upp_math, only: h2u, h2v, u2h, v2h
126
127!
128!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129 implicit none
130!
131 REAL, PARAMETER :: CURATE=24.*1000., ctim1=0., ctim2=24.*3600. &
132 &, raincon=0.8333*1.1787e4, snocon=0.94*1.4594e5 &
133! specify in params now
134!
135!--- 88D reflectivity algorithm, Z = 300.*R**1.4 , R is rain rate in mm/h
136!
137 &, dbzmax=80., zr_a=300., zr_b=1.4
138!
139!--- Modification of Slingo (1987) to enhance convective cloudiness
140!
141 REAL CC(10), PPT(10)
142 DATA cc / 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /
143 DATA ppt/ 0., .14, .31, .70, 1.6, 3.4, 7.7, 17., 38., 85. /
144 INTEGER, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: ICBOT, ICTOP, LPBL
145
146!
147! DECLARE VARIABLES.
148!
149! LOGICAL NORTH,NEED(IM,JM), NMM_GFSmicro
150 LOGICAL NMM_GFSmicro
151 LOGiCAL Model_Radar
152 real, dimension(im,jm) :: GRID1, GRID2
153 real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: EGRID1, EGRID2, EGRID3, EGRID4, EGRID5,&
154 el0, p1d, t1d, q1d, c1d, &
155 fi1d, fr1d, fs1d, qw1, qi1, &
156 qr1, qs1, curefl_s, &
157 curefl, curefl_i, zfrz, dbz1, dbzr1, &
158 dbzi1, dbzc1, egrid6, egrid7, nlice1, &
159 qi, qint, tt, ppp, qv, &
160 qcd, qice1, qrain1, qsno1, refl, &
161 qg1, refl1km, refl4km, rh, gust, nrain1,zm10c, &
162 ustore, vstore
163! T700, TH700
164!
165 REAL, ALLOCATABLE :: EL(:,:,:),RICHNO(:,:,:) ,PBLRI(:,:), PBLREGIME(:,:)
166!
167 integer I,J,L,Lctop,LLMH,IICE,LL,II,JJ,IFINCR,ITHEAT,NC,NMOD,LLL &
168 ,iz1km,iz4km, lcount, hcount, itype, item
169
170 real RDTPHS,CFRdum,PMOD,CC1,CC2,P1,P2,CUPRATE,FACR,RRNUM &
171 ,rainrate,term1,term2,term3,qrold,snorate,dens,delz,fctr,hgt &
172 ,rain,ronv,slor,snow,rhoqs,temp_c,sonv,slos &
173 ,graupel,rhoqg,gonv,slog, alpha, rhod, bb &
174 ,ze_s, ze_r, ze_g, ze_max, ze_nc, ze_conv, ze_sum &
175 ,ze_smax, ze_rmax,ze_gmax, ze_nc_1km, ze_nc_4km, dz &
176 ,lapses, expo,expinv,tsfcnew, gam,gamd,gams, pblhold &
177 ,psfc,tsfc,zsfc,dp,dpbnd,zmin
178
179 real, allocatable :: RH3D(:,:,:)
180
181! for PBL smoothing used in GUST
182 integer ks,nsmooth
183 REAL SDUMMY(IM,2),dxm
184! added to calculate cape and cin for icing
185 real, dimension(ista:iend,jsta:jend) :: dummy, cape, cin
186 integer idummy(ista:iend,jsta:jend)
187
188 real, PARAMETER :: ZSL=0.0, taucr=rd*gi*290.66, const=0.005*g/rd, gord=g/rd
189 logical, parameter :: debugprint = .false.
190
191 gams = 0.0065
192 gamd = 0.0100
193
194 lapses = 0.0065 ! deg K / meter
195 expo = rog*lapses
196 expinv = 1./expo
197
198 zmin=10.**(0.1*dbzmin)
199!
200!
201!*****************************************************************************
202! START SUBROUTINE MDLFLD.
203!
204! ALLOCATE LOCAL ARRAYS
205!
206! Initialize halo regions of USTORE & VSTORE for cases when the halo extends
207! beyond the computational domain boundary.
208!$OMP PARALLEL DO COLLAPSE(2)
209 DO j=jsta_2l,jend_2u
210 DO i=ista_2l,iend_2u
211 ustore(i,j) = 0
212 vstore(i,j) = 0
213 ENDDO
214 ENDDO
215! Set up logical flag to indicate whether model outputs radar directly
216 model_radar = .false.
217! IF (ABS(MAXVAL(REF_10CM)-SPVAL)>SMALL)Model_Radar=.True.
218 check_ref: DO l=1,lm
219 DO j=jsta,jend
220 DO i=ista,iend
221 IF(abs(ref_10cm(i,j,l)-spval)>small) THEN
222 model_radar=.true.
223 exit check_ref
224 ENDIF
225 ENDDO
226 ENDDO
227 ENDDO check_ref
228 if(debugprint .and. me==0)print*,'Did post read in model derived radar ref ',model_radar, &
229 'MODELNAME=',trim(modelname),' imp_physics=',imp_physics
230 ALLOCATE(el(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
231 ALLOCATE(richno(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
232 ALLOCATE(pblri(ista_2l:iend_2u,jsta_2l:jend_2u))
233!
234! SECOND, STANDARD NGM SEA LEVEL PRESSURE.
235 IF (iget(023) > 0 .OR. iget(105) > 0 .OR. iget(445) > 0) THEN
236 CALL ngmslp ! this value is used in some later calculation.
237 ENDIF
238 IF (iget(105) > 0) THEN
239!$omp parallel do private(i,j)
240 DO j=jsta,jend
241 DO i=ista,iend
242 grid1(i,j) = slp(i,j)
243 ENDDO
244 ENDDO
245 if(grib=="grib2") then
246 cfld=cfld+1
247 fld_info(cfld)%ifld=iavblfld(iget(105))
248!$omp parallel do private(i,j,ii,jj)
249 do j=1,jend-jsta+1
250 jj = jsta+j-1
251 do i=1,iend-ista+1
252 ii = ista+i-1
253 datapd(i,j,cfld) = grid1(ii,jj)
254 enddo
255 enddo
256 endif
257
258 ENDIF
259!
260!--- Calculate convective cloud fractions following radiation in
261! NMM; used in subroutine CALRAD_WCLOUD for satellite radiances
262!Both FV3 regional and global output CNVCFR directly
263 IF (modelname=='NMM' .OR. imp_physics==5 .or. &
264 imp_physics==85 .or. imp_physics==95) THEN
265! print*,'DTQ2 in MDLFLD= ',DTQ2
266 rdtphs=24.*3.6e6/dtq2
267 DO j=jsta,jend
268 DO i=ista,iend
269 IF ((hbot(i,j)-htop(i,j)) <= 1.0) THEN
270 icbot(i,j)=0
271 ictop(i,j)=0
272 cnvcfr(i,j)=0.
273 ELSE
274 icbot(i,j)=nint(hbot(i,j))
275 ictop(i,j)=nint(htop(i,j))
276 cfrdum=cc(1)
277 pmod=rdtphs*cprate(i,j) ! mm/day
278 IF (pmod > ppt(1)) THEN
279 DO nc=1,10
280 IF(pmod>ppt(nc)) nmod=nc
281 ENDDO
282 IF (nmod >= 10) THEN
283 cfrdum=cc(10)
284 ELSE
285 cc1=cc(nmod)
286 cc2=cc(nmod+1)
287 p1=ppt(nmod)
288 p2=ppt(nmod+1)
289 cfrdum=cc1+(cc2-cc1)*(pmod-p1)/(p2-p1)
290 ENDIF !--- End IF (NMOD >= 10) ...
291 cfrdum=min(h1, cfrdum)
292 ENDIF !--- End IF (PMOD > PPT(1)) ...
293! CNVCFR(I,J)=100.*CFRdum
294 cnvcfr(i,j)=cfrdum
295 ENDIF !--- End IF (HBOT(I,J)-HTOP(I,J) <= 1.0) ...
296 ENDDO !--- DO I=ista,iend
297 ENDDO !--- DO J=JSTA,JEND
298 ENDIF !-- IF (MODELNAME=='NMM' .OR. imp_physics==5) THEN
299!
300!--- Set
301 IF (modelname=='NMM' .AND. imp_physics==9) THEN
302 nmm_gfsmicro=.true.
303 ELSE
304 nmm_gfsmicro=.false.
305 ENDIF
306!
307! Calculate convective radar reflectivity at the surface (CUREFL_S),
308! and the decrease in reflectivity above the 0C level (CUREFL_I)
309!
310 IF(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95 &
311 .or. nmm_gfsmicro)THEN
312 rdtphs=3.6e6/dtq2
313 DO j=jsta,jend
314 DO i=ista,iend
315 cuprate=rdtphs*cprate(i,j) !--- Cu precip rate, R (mm/h)
316! CUPRATE=CUPPT(I,J)*1000./TRDLW !--- mm/h
317 zfrz(i,j)=zmid(i,j,nint(lmh(i,j))) !-- Initialize to lowest model level
318 DO l=1,nint(lmh(i,j)) !-- Start from the top, work down
319 IF (t(i,j,l) >= tfrz) THEN
320 zfrz(i,j)=zmid(i,j,l) !-- Find highest level where T>0C
321 EXIT
322 ENDIF
323 ENDDO !--- DO L=1,NINT(LMH(I,J))
324! IF (CUPRATE <= 0. .OR. CUPPT(I,J)<=0.) THEN
325 IF (cuprate <= 0. .or. htop(i,j)>=spval) THEN ! bug fix, post doesn not use CUPPT
326 curefl_s(i,j)=0.
327 curefl_i(i,j)=0.
328 ELSE
329 curefl_s(i,j)=zr_a*cuprate**zr_b !--- Use Z=A*R**B
330 lctop=nint(htop(i,j)) !--- Cu cld top level
331!
332!--- Assume convective reflectivity (Z, not dBZ) above 0C level decreases
333! with height by two orders of magnitude (20 dBZ) from the 0C level up
334! to cloud top. If cloud top temperature is above 0C, assume 20 dBZ
335! decrease occurs in the first 1 km above the 0C level.
336!
337 curefl_i(i,j)=-2./max( 1000., zmid(i,j,lctop)-zfrz(i,j) )
338 ENDIF !--- IF (CUPRATE <= 0. .OR. CUPPT(I,J)<=0.) THEN
339 ENDDO !--- End DO I
340 ENDDO
341
342!
343!--- Calculate each hydrometeor category & GRID-SCALE cloud fraction
344! (Jin, Aug-Oct '01; Ferrier, Feb '02)
345!
346
347 if(icount_calmict==0)then !only call calmict once in multiple grid processing
348 DO l=1,lm
349 DO j=jsta,jend
350 DO i=ista,iend
351 p1d(i,j)=pmid(i,j,l)
352 t1d(i,j)=t(i,j,l)
353 q1d(i,j)=q(i,j,l)
354 c1d(i,j)=cwm(i,j,l)
355 fi1d(i,j)=f_ice(i,j,l)
356 fr1d(i,j)=f_rain(i,j,l)
357 fs1d(i,j)=max(h1, f_rimef(i,j,l))
358!
359!--- Estimate radar reflectivity factor at level L
360!
361 curefl(i,j)=0.
362 IF (curefl_s(i,j) > 0.) THEN
363 fctr=0.
364 llmh = nint(lmh(i,j))
365 lctop=nint(htop(i,j)) !--- Cu cld top level
366 IF (l>=lctop .AND. l<=llmh) THEN
367 delz=zmid(i,j,l)-zfrz(i,j)
368 IF (delz <= 0.) THEN
369 fctr=1. !-- Below the highest freezing level
370 ELSE
371 !
372 !--- Reduce convective radar reflectivity above freezing level
373 !
374 fctr=10.**(curefl_i(i,j)*delz)
375 ENDIF !-- End IF (DELZ <= 0.)
376 ENDIF !-- End IF (L>=HTOP(I,J) .OR. L<=LLMH)
377 curefl(i,j)=fctr*curefl_s(i,j)
378 ENDIF !-- End IF (CUREFL_S(I,J) > 0.)
379
380 ENDDO !-- End DO I loop
381 ENDDO !-- End DO J loop
382 IF(imp_physics==5 .or. imp_physics==85 .or. imp_physics==95)THEN
383 fer_mic: IF (imp_physics==5) THEN
384!
385!--- Ferrier-Aligo microphysics in the NMMB
386!
387!--- Determine composition of condensate in terms of cloud water,
388! rain, and ice (cloud ice & precipitation ice) following the
389! *NEWER* the version of the microphysics; radar reflectivity
390! is derived to be consistent with the microphysical assumptions
391!
392 CALL calmict_new(p1d,t1d,q1d,c1d,fi1d,fr1d,fs1d,curefl &
393 & ,qw1,qi1,qr1,qs1,dbz1,dbzr1,dbzi1,dbzc1,nlice1, nrain1)
394 IF(modelname == 'NMM' .and. gridtype=='B')THEN !NMMB
395!
396!--- Use reflectivity from NMMB model output for Ferrier-Aligo (imp_physics=5),
397! add bogused contribution from parameterized convection (CUREFL), and
398! estimate reflectivity from rain (DBZR1) & snow/graupel (DBZI1).
399!
400refl_miss: IF (model_radar) THEN
401 ! - Model output DBZ is present - proceed with calc
402 DO j=jsta,jend
403 DO i=ista,iend
404 IF(p1d(i,j)<spval.and.t1d(i,j)<spval.and.q1d(i,j)<spval)THEN
405 ze_nc=10.**(0.1*ref_10cm(i,j,l))
406 dbz1(i,j)=10.*log10(max(zmin,(ze_nc+curefl(i,j))))
407 dbzr1(i,j)=min(dbzr1(i,j), ref_10cm(i,j,l))
408 dbzi1(i,j)=min(dbzi1(i,j), ref_10cm(i,j,l))
409 ze_max=max(dbzr1(i,j),dbzi1(i,j))
410refl_comp: IF(ref_10cm(i,j,l)>dbzmin .OR. ze_max>dbzmin) THEN
411refl_adj: IF(ref_10cm(i,j,l)<=dbzmin) THEN
412 dbzr1(i,j)=dbzmin
413 dbzi1(i,j)=dbzmin
414 ELSE IF(ze_max<=dbzmin) THEN
415 IF(qr1(i,j)>qs1(i,j)) THEN
416 dbzr1(i,j)=ref_10cm(i,j,l)
417 ELSE IF(qs1(i,j)>qr1(i,j)) THEN
418 dbzi1(i,j)=ref_10cm(i,j,l)
419 ELSE
420 IF(t1d(i,j)>=tfrz) THEN
421 dbzr1(i,j)=ref_10cm(i,j,l)
422 ELSE
423 dbzi1(i,j)=ref_10cm(i,j,l)
424 ENDIF
425 ENDIF
426 ELSE
427 ze_nc=10.**(0.1*ref_10cm(i,j,l))
428 ze_r=10.**(0.1*dbzr1(i,j))
429 ze_s=10.**(0.1*dbzi1(i,j))
430 ze_sum=ze_r+ze_s
431 ze_max=ze_nc/ze_sum
432 ze_r=ze_r*ze_max
433 ze_s=ze_s*ze_max
434 dbzr1(i,j)=10.*log10(ze_r)
435 dbzi1(i,j)=10.*log10(ze_s)
436 ENDIF refl_adj
437 ENDIF refl_comp
438 ELSE
439 dbzr1(i,j)=dbzmin
440 dbzi1(i,j)=dbzmin
441 ENDIF
442 ENDDO
443 ENDDO
444 ELSE
445 ! - Model output dBZ is missing
446 IF (me==0 .AND. l==1) THEN
447 WRITE(6,'(4A,1x,F7.2)') 'WARNING - MDLFLD: REF_10CM NOT ', &
448 'IN NMMB OUTPUT. CHECK ', &
449 'SOLVER_STATE.TXT FILE. USING ', &
450 'REFL OUTPUT FROM CALMICT.'
451 ENDIF
452 ENDIF refl_miss
453 ENDIF
454 ELSE fer_mic
455!
456!--- Determine composition of condensate in terms of cloud water,
457! rain, and ice (cloud ice & precipitation ice) following the
458! *OLDER* the version of the microphysics; radar reflectivity
459! is derived to be consistent with the microphysical assumptions
460!
461 CALL calmict_old(p1d,t1d,q1d,c1d,fi1d,fr1d,fs1d,curefl &
462 & ,qw1,qi1,qr1,qs1,dbz1,dbzr1,dbzi1,dbzc1,nlice1, nrain1)
463 ENDIF fer_mic
464
465 ELSE
466!
467!--- This branch is executed if GFS micro (imp_physics=9) is run in the NMM.
468!
469 DO j=jsta,jend
470 DO i=ista,iend
471 IF(c1d(i,j)<spval.and.fi1d(i,j)<spval)THEN
472 qi1(i,j)=c1d(i,j)*fi1d(i,j)
473 qw1(i,j)=c1d(i,j)-qi1(i,j)
474 ELSE
475 qi1(i,j)=d00
476 qw1(i,j)=d00
477 ENDIF
478 qr1(i,j)=d00
479 qs1(i,j)=d00
480 dbz1(i,j)=dbzmin
481 dbzr1(i,j)=dbzmin
482 dbzi1(i,j)=dbzmin
483 dbzc1(i,j)=dbzmin
484 ENDDO
485 ENDDO
486 ENDIF
487 DO j=jsta,jend
488 DO i=ista,iend
489 llmh = nint(lmh(i,j))
490 IF (l > llmh) THEN
491 qqw(i,j,l) = d00
492 qqi(i,j,l) = d00
493 qqr(i,j,l) = d00
494 qqs(i,j,l) = d00
495 cfr(i,j,l) = d00
496 dbz(i,j,l) = dbzmin
497 dbzr(i,j,l) = dbzmin
498 dbzi(i,j,l) = dbzmin
499 dbzc(i,j,l) = dbzmin
500 ELSE
501 qqw(i,j,l) = max(d00, qw1(i,j))
502 qqi(i,j,l) = max(d00, qi1(i,j))
503 qqr(i,j,l) = max(d00, qr1(i,j))
504 qqs(i,j,l) = max(d00, qs1(i,j))
505 dbz(i,j,l) = max(dbzmin, dbz1(i,j))
506 dbzr(i,j,l) = max(dbzmin, dbzr1(i,j))
507 dbzi(i,j,l) = max(dbzmin, dbzi1(i,j))
508 dbzc(i,j,l) = max(dbzmin, dbzc1(i,j))
509 nlice(i,j,l) = max(d00, nlice1(i,j))
510 nrain(i,j,l) = max(d00, nrain1(i,j))
511 ENDIF !-- End IF (L > LMH(I,J)) ...
512 ENDDO !-- End DO I loop
513 ENDDO !-- End DO J loop
514
515 ENDDO !-- End DO L loop
516 END IF ! end of icount_calmict
517 icount_calmict=icount_calmict+1
518 if(debugprint .and. me==0)print*,'debug calmict:icount_calmict= ',icount_calmict
519
520! Chuang: add the option to compute individual microphysics species
521! for NMMB+Zhao and NMMB+WSM6 which are two of SREF members.
522! Per communication with Ferrier (July 2012), he has set up these
523! 2 runs to output CWM plus fraction arrays instead of individual
524! microphysics species arrays.
525! WRF NMM + non Ferrier still outputs individual microphysics
526! arrays so these 2 if branches are excuted for NMMB only.
527 ELSE IF(modelname == 'NMM' .and. gridtype=='B' .and. imp_physics==99)THEN !NMMB+Zhao
528 DO l=1,lm
529 DO j=jsta,jend
530 DO i=ista,iend
531 llmh = nint(lmh(i,j))
532 IF (l > llmh) THEN
533 qqw(i,j,l) = d00
534 qqi(i,j,l) = d00
535 qqr(i,j,l) = d00
536 qqs(i,j,l) = d00
537 cfr(i,j,l) = d00
538 dbz(i,j,l) = dbzmin
539 dbzr(i,j,l) = dbzmin
540 dbzi(i,j,l) = dbzmin
541 dbzc(i,j,l) = dbzmin
542 ELSE
543 qqi(i,j,l) = max(d00, cwm(i,j,l)*f_ice(i,j,l))
544 qqw(i,j,l) = max(d00, cwm(i,j,l)-qqi(i,j,l))
545 qqr(i,j,l) = d00
546 qqs(i,j,l) = d00
547 dbz(i,j,l) = dbzmin
548 dbzr(i,j,l) = dbzmin
549 dbzi(i,j,l) = dbzmin
550 dbzc(i,j,l) = dbzmin
551 ENDIF !-- End IF (L > LMH(I,J)) ...
552 ENDDO !-- End DO I loop
553 ENDDO ! END DO L LOOP
554 END DO
555 ELSE IF(modelname == 'NMM' .and. gridtype=='B' .and. imp_physics==6)THEN !NMMB+WSM6
556 DO l=1,lm
557 DO j=jsta,jend
558 DO i=ista,iend
559 llmh = nint(lmh(i,j))
560 IF (l > llmh) THEN
561 qqw(i,j,l)=d00
562 qqi(i,j,l)=d00
563 qqr(i,j,l)=d00
564 qqs(i,j,l)=d00
565 cfr(i,j,l)=d00
566 dbz(i,j,l)=dbzmin
567 dbzr(i,j,l)=dbzmin
568 dbzi(i,j,l)=dbzmin
569 dbzc(i,j,l)=dbzmin
570 ELSE
571 qqi(i,j,l)=d00
572 qqw(i,j,l)=max(d00, (1.-f_ice(i,j,l))*cwm(i,j,l)*(1.-f_rain(i,j,l)))
573 qqr(i,j,l)=max(d00,(1.-f_ice(i,j,l))*cwm(i,j,l)*f_rain(i,j,l))
574 qqs(i,j,l)=max(d00, cwm(i,j,l)*f_ice(i,j,l))
575 dens=pmid(i,j,l)/(rd*t(i,j,l)*(q(i,j,l)*d608+1.0)) ! DENSITY
576 dbzr(i,j,l)=((qqr(i,j,l)*dens)**1.75)* &
577 & 3.630803e-9 * 1.e18 ! Z FOR RAIN
578 dbzi(i,j,l)= dbzi(i,j,l)+((qqs(i,j,l)*dens)**1.75)* &
579 & 2.18500e-10 * 1.e18 ! Z FOR SNOW
580 dbz(i,j,l)=dbzr(i,j,l)+dbzi(i,j,l)
581 IF (dbz(i,j,l)>0.) dbz(i,j,l)=10.0*log10(dbz(i,j,l)) ! DBZ
582 IF (dbzr(i,j,l)>0.)dbzr(i,j,l)=10.0*log10(dbzr(i,j,l)) ! DBZ
583 IF (dbzi(i,j,l)>0.) &
584 & dbzi(i,j,l)=10.0*log10(dbzi(i,j,l)) ! DBZ
585 dbz(i,j,l)=max(dbzmin, dbz(i,j,l))
586 dbzr(i,j,l)=max(dbzmin, dbzr(i,j,l))
587 dbzi(i,j,l)=max(dbzmin, dbzi(i,j,l))
588 ENDIF !-- End IF (L > LMH(I,J)) ...
589 ENDDO !-- End DO I loop
590 ENDDO
591 END DO
592
593 ELSE IF(((modelname == 'NMM' .and. gridtype=='B') .OR. modelname == 'FV3R' &
594 .OR. modelname == 'GFS') &
595 .and. (imp_physics==8 .or. imp_physics==17 .or. imp_physics==18))THEN !NMMB or FV3R or GFS +THOMPSON
596 DO l=1,lm
597 DO j=jsta,jend
598 DO i=ista,iend
599 dbz(i,j,l)=ref_10cm(i,j,l)
600 ENDDO
601 ENDDO
602 ENDDO
603 ELSE IF(imp_physics==99 .or. imp_physics==98)THEN ! Zhao MP
604 DO l=1,lm
605 DO j=jsta,jend
606 DO i=ista,iend
607 dbz(i,j,l)=spval
608 ENDDO
609 ENDDO
610 ENDDO
611 ELSE ! compute radar refl for other than NAM/Ferrier or GFS/Zhao microphysics
612 if(debugprint .and. me==0)print*,'calculating radar ref for non-Ferrier/non-Zhao schemes'
613! Determine IICE FLAG
614 IF(imp_physics == 1 .OR. imp_physics == 3)THEN
615 iice = 0
616 ELSE
617 iice = 1
618 END IF
619
620! Chuang: add convective contribution for all MP schemes
621 rdtphs=3.6e6/dtq2
622 DO j=jsta,jend
623 DO i=ista,iend
624 cuprate=rdtphs*cprate(i,j) !--- Cu precip rate, R (mm/h)
625 zfrz(i,j)=zmid(i,j,nint(lmh(i,j))) !-- Initialize to lowest model level
626 DO l=1,nint(lmh(i,j)) !-- Start from the top, work down
627 IF (t(i,j,l) >= tfrz) THEN
628 zfrz(i,j)=zmid(i,j,l) !-- Find highest level where T>0C
629 EXIT
630 ENDIF
631 ENDDO !--- DO L=1,NINT(LMH(I,J))
632! IF (CUPRATE <= 0. .OR. CUPPT(I,J)<=0.) THEN
633 IF (cuprate <= 0. .or. htop(i,j)>=spval) THEN ! bug fix, post doesn not use CUPPT
634 curefl_s(i,j)=0.
635 curefl_i(i,j)=0.
636 ELSE
637 curefl_s(i,j)=zr_a*cuprate**zr_b !--- Use Z=A*R**B
638 lctop=nint(htop(i,j)) !--- Cu cld top level
639!
640!--- Assume convective reflectivity (Z, not dBZ) above 0C level decreases
641! with height by two orders of magnitude (20 dBZ) from the 0C level up
642! to cloud top. If cloud top temperature is above 0C, assume 20 dBZ
643! decrease occurs in the first 1 km above the 0C level.
644!
645 curefl_i(i,j)=-2./max( 1000., zmid(i,j,lctop)-zfrz(i,j) )
646 ENDIF !--- IF (CUPRATE <= 0. .OR. CUPPT(I,J)<=0.) THEN
647 ENDDO !--- End DO I
648 ENDDO
649
650 IF(imp_physics /= 8 .AND. imp_physics /= 9 .and. imp_physics /= 28) THEN
651!tgs - non-Thompson schemes
652!$omp parallel do private(i,j,l,curefl,fctr,dens,llmh,lctop,delz,ze_nc)
653 DO l=1,lm
654 DO j=jsta,jend
655 DO i=ista,iend
656!--- Estimate radar reflectivity factor from convection at level L
657!
658 curefl(i,j)=0.
659 IF (curefl_s(i,j) > 0.) THEN
660 fctr=0.
661 llmh = nint(lmh(i,j))
662 lctop=nint(htop(i,j)) !--- Cu cld top level
663 IF (l>=lctop .AND. l<=llmh) THEN
664 delz=zmid(i,j,l)-zfrz(i,j)
665 IF (delz <= 0.) THEN
666 fctr=1. !-- Below the highest freezing level
667 ELSE
668 !
669 !--- Reduce convective radar reflectivity above freezing level
670 !
671 fctr=10.**(curefl_i(i,j)*delz)
672 ENDIF !-- End IF (DELZ <= 0.)
673 ENDIF !-- End IF (L>=HTOP(I,J) .OR. L<=LLMH)
674 curefl(i,j)=fctr*curefl_s(i,j)
675 dbzc(i,j,l)=curefl(i,j)
676 ENDIF !-- End IF (CUREFL_S(I,J) > 0.)
677
678 IF(t(i,j,l)<spval) THEN
679! IF(T(I,J,L) < 1.0E-3) print*,'ZERO T'
680 IF(t(i,j,l) > 1.0e-3) &
681 & dens = pmid(i,j,l)/(rd*t(i,j,l)*(q(i,j,l)*d608+1.0)) ! DENSITY
682
683! PATCH to se(1.-FI1D(I,J))*C1D(I,J)*FR1D(I,J)t QQR, QQS, AND QQG to
684! zeros if they are negative so that post won't abort
685
686 qqr(i,j,l) = max(qqr(i,j,l),0.0)
687 qqs(i,j,l) = max(qqs(i,j,l),0.0) ! jkw
688 IF (iice == 0) THEN
689 IF (t(i,j,l) >= tfrz) THEN
690 dbz(i,j,l) = ((qqr(i,j,l)*dens)**1.75)* &
691 & 3.630803e-9 * 1.e18 ! Z FOR RAIN
692 dbzr(i,j,l) = dbz(i,j,l)
693 ELSE
694!mptest DBZ(I,J,L) = ((QQR(I,J,L)*DENS)**1.75)* &
695 dbz(i,j,l) = ((qqs(i,j,l)*dens)**1.75)* &
696 & 2.18500e-10 * 1.e18 ! Z FOR SNOW
697 dbzi(i,j,l) = dbz(i,j,l)
698 ENDIF
699 ELSEIF (iice == 1) THEN
700 dbzi(i,j,l) = 0.
701 qqg(i,j,l) = max(qqg(i,j,l),0.0)
702 if(qqr(i,j,l) < spval .and. qqr(i,j,l)> 0.0) then
703 dbzr(i,j,l) = ((qqr(i,j,l)*dens)**1.75) * 3.630803e-9 * 1.e18 ! Z FOR RAIN
704 else
705 dbzr(i,j,l) = 0.
706 endif
707 if(qqs(i,j,l) < spval .and. qqs(i,j,l) > 0.0) then
708 dbzi(i,j,l) = ((qqs(i,j,l)*dens)**1.75) * &
709 & 2.18500e-10 * 1.e18 ! Z FOR SNOW
710 else
711 dbzi(i,j,l) = 0.
712 endif
713 IF (qqg(i,j,l) < spval .and. qqg(i,j,l)> 0.0) then
714 dbzi(i,j,l) = dbzi(i,j,l) + ((qqg(i,j,l)*dens)**1.75) * &
715 & 1.033267e-9 * 1.e18 ! Z FOR GRAUP
716 else
717 dbzi(i,j,l) = dbzi(i,j,l)
718 endif
719 IF (model_radar) THEN
720 ze_nc=10.**(0.1*ref_10cm(i,j,l))
721 dbz(i,j,l) = ze_nc+curefl(i,j)
722 ELSE
723 dbz(i,j,l) = dbzr(i,j,l) + dbzi(i,j,l) + curefl(i,j)
724 END IF
725! IF(L==27.and.QQR(I,J,L)>1.e-4)print*, &
726! 'sample QQR DEN,DBZ= ',QQR(I,J,L),DENS,DBZ(I,J,L)
727 ENDIF
728 IF (dbz(i,j,l) > 0.) dbz(i,j,l) = 10.0*log10(dbz(i,j,l)) ! DBZ
729 IF (dbzr(i,j,l) > 0.) dbzr(i,j,l) = 10.0*log10(dbzr(i,j,l)) ! DBZ
730 IF (dbzi(i,j,l) > 0.) dbzi(i,j,l) = 10.0*log10(dbzi(i,j,l)) ! DBZ
731 IF (dbzc(i,j,l) > 0.) dbzc(i,j,l) = 10.0*log10(dbzc(i,j,l)) ! DBZ
732 llmh = nint(lmh(i,j))
733 IF(l > llmh) THEN
734 dbz(i,j,l) = dbzmin
735 dbzr(i,j,l) = dbzmin
736 dbzi(i,j,l) = dbzmin
737 dbzc(i,j,l) = dbzmin
738 ELSE
739 dbz(i,j,l) = max(dbzmin, dbz(i,j,l))
740 dbzr(i,j,l) = max(dbzmin, dbzr(i,j,l))
741 dbzi(i,j,l) = max(dbzmin, dbzi(i,j,l))
742 dbzc(i,j,l) = max(dbzmin, dbzc(i,j,l))
743 END IF
744 ELSE
745 dbz(i,j,l) = dbzmin
746 dbzr(i,j,l) = dbzmin
747 dbzi(i,j,l) = dbzmin
748 dbzc(i,j,l) = dbzmin
749 ENDIF !(T(I,J,L)<spval)
750 ENDDO
751 ENDDO
752 ENDDO
753!
754!tgs
755 ELSE
756! for Thompson microphisics scheme (option 8), developed at GSD/ESRL
757! 13 January 2009
758 call paramr ! compute constants for reflectivity algorithm
759
760 bb = 0. ! bright band effect - yes or no (0)
761 alpha = 0.224 ! = (1000kg/m^3/917kg/m^3)**2)*(0.176/0.930)
762! 1000kg/m^3 is density of liquid water
763! 917kg/m^3 is density of solid ice
764! 0.176 = dielectric factor of ice
765! 0.930 = dielectric factor of liquid water
766
767 ze_smax = -1.e30
768 ze_rmax = -1.e30
769 ze_gmax = -1.e30
770
771 DO j=jsta,jend
772 DO i=ista,iend
773 refl(i,j) = -10.
774 ze_max = -10.
775
776 iz1km = 0
777 iz4km = 0
778
779 DO l=1,lm
780 ll=lm-l+1
781 IF(t(i,j,ll)<spval)THEN
782 IF(t(i,j,ll) < 1.0e-3)print*,'ZERO T'
783 IF(t(i,j,ll) > 1.0e-3) &
784 rhod=pmid(i,j,ll)/ &
785 (rd*t(i,j,ll)*(q(i,j,ll)*d608+1.0)) ! DENSITY
786 dz=zint(i,j,ll)-zint(i,j,lm+1)
787! Particle size distributions and reflectivity
788! ---------------------------------------------
789! Much of this code borrowed from EXMOISG loop 20 to get particle size
790! distributions
791
792!jmb-- Note that SLOR, SLOS and SLOG are inverse slopes!!!! Also,
793! RONV,SONV,GONV, M-P zero intercept values, normalized by
794! max allowable values.
795
796! Have to set min values of hydrometeors (r1) large enough to avoid
797! underflow problems with log later on.
798
799! -- rain
800 ze_r = 1.e-35
801 if (qqr(i,j,ll) >= 1.e-6) then
802 rain = max(r1,qqr(i,j,ll))
803 ronv = (const1r*tanh((qr0 - rain)/delqr0) + &
804 const2r)/ron
805 slor=(rhod*rain/(topr*ronv))**0.25
806 ze_r = 720.*ronv*ron*slor**7 ! Stoelinga Eq. 2, reflectivity
807 endif
808
809! -- snow
810 ze_s = 1.e-35
811 if (qqs(i,j,ll) >= 1.e-6) then
812 snow = max(r1,qqs(i,j,ll))
813! New SONV formulation based on Fig. 7, curve_3 of Houze et al 1979
814 rhoqs=rhod*snow
815 temp_c = min(-0.001, t(i,j,ll)-273.15)
816 sonv = (min(2.0e8, 2.0e6*exp(-0.12*temp_c)))/son
817 slos=(rhoqs/(tops*sonv))**0.25
818 ze_s = 720.*alpha*sonv*son*slos**7*(dsnow/drain)**2
819! From Stoelinga Eq. 5, reflectivity
820
821! For bright band, increase reflectivity by factor of 5.28,
822! which is ratio of dielectric factors for water/ice (.930/.176)
823 IF (t(i,j,ll) > 273.15) &
824 ze_s = ze_s*(1. + 4.28*bb)
825 endif
826
827! -- graupel
828 ze_g = 1.e-35
829 if (qqg(i,j,ll) >= 1.e-6) then
830 graupel = max(r1,qqg(i,j,ll))
831 rhoqg=rhod*graupel
832 gonv=1.
833 gonv=const_ng1*(rhoqg**const_ng2)
834 gonv = max(1.e4, min(gonv,gon))
835 gonv = gonv/gon
836 slog=(rhoqg/(topg*gonv))**0.25
837 ze_g = 720.*alpha*gonv*gon*slog**7*(dgraupel/drain)**2
838! Stoelinga Eq. 5 applied to graupel
839
840! For bright band
841 IF (t(i,j,ll) > 273.15) &
842 ze_g = ze_g*(1. + 4.28*bb)
843 endif
844
845! -- total grid scale
846 ze_nc = ze_r + ze_s + ze_g
847
848 if (iz1km==0 .and. dz>1000.) then
849 ze_nc_1km = ze_nc
850 iz1km = 1
851 end if
852
853 if (iz4km==0 .and. dz>4000.) then
854 ze_nc_4km = ze_nc
855 iz4km = 1
856 end if
857
858 ze_rmax = max(ze_r,ze_rmax)
859 ze_smax = max(ze_s,ze_smax)
860 ze_gmax = max(ze_g,ze_gmax)
861! Reflectivities are in units of m^6/m^3
862! convert to mm^6/m^3 and take log base 10 to get
863! reflectivities in dbZe (decibels).
864! comp_refl_r(j,k) = 10.*LOG10(ze_r*1.E18)
865! comp_refl_s(j,k) = 10.*LOG10(ze_s*1.E18)
866! comp_refl_g(j,k) = 10.*LOG10(ze_g*1.E18)
867! comp_refl_nc(j,k) = 10.*LOG10(ze_nc*1.E18)
868
869
870! Total composite reflectivity, including convection, in dbZe
871 ze_sum = ze_nc*1.e18 ! + ze_conv
872 ze_max = max(ze_max, ze_sum )
873
874 dbz(i,j,ll) = ze_sum
875 dbzr(i,j,ll) = ze_r*1.e18
876 dbzi(i,j,ll) = (ze_s+ze_g)*1.e18
877 ELSE
878 dbz(i,j,ll) = dbzmin
879 dbzr(i,j,ll) = dbzmin
880 dbzi(i,j,ll) = dbzmin
881 ENDIF !T(I,J,LL)<spval
882 ENDDO
883! parameterized convection
884! -------------------------
885! conv_prate(i,j) is convective pcpn rate, assumed in mm/h
886! ze_conv = 300.*conv_prate**1.4 ! Units: mm^6/m^3
887
888 rdtphs=3.6e6/dt
889 cuprate=rdtphs*cprate(i,j) !--- Cu precip rate, R (mm/h)
890
891! ze_conv= max(0.1,300*(4.*CUPRATE)**1.4)
892! -- switch to time-step conv precip in RR
893 ze_conv= max(0.1,300*(cuprate)**1.4)
894
895! Combine max resolved reflectivity component
896! and sub-grid scale component
897! Total composite reflectivity, including convection, in dbZe
898 ze_sum = ze_max + ze_conv
899 refl(i,j) = 10.*log10(ze_sum)
900 refl1km(i,j) = 10.*log10(ze_nc_1km*1.e18 + ze_conv)
901 refl4km(i,j) = 10.*log10(ze_nc_4km*1.e18 + ze_conv)
902
903 ENDDO
904 ENDDO
905
906 ze_rmax = 10.*log10(ze_rmax*1.e18)
907 ze_smax = 10.*log10(ze_smax*1.e18)
908 ze_gmax = 10.*log10(ze_gmax*1.e18)
909
910 ENDIF !tgs endif for Thompson scheme
911
912 END IF
913!
914! OUTPUT/CALCULATE PRESSURE, OMEGA, POTENTIAL TEMPERATURE,
915! DEWPOINT TEMPERATURE, RELATIVE HUMIDITY, AND
916! ABSOLUTE VORTICITY ON MDL SURFACES.
917!
918!
919 allocate (rh3d(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
920 IF ( (iget(001)>0).OR.(iget(077)>0).OR. &
921 (iget(002)>0).OR.(iget(003)>0).OR. &
922 (iget(004)>0).OR.(iget(005)>0).OR. &
923 (iget(006)>0).OR.(iget(083)>0).OR. &
924 (iget(007)>0).OR.(iget(008)>0).OR. &
925 (iget(009)>0).OR.(iget(010)>0).OR. &
926 (iget(084)>0).OR.(iget(011)>0).OR. &
927 (iget(041)>0).OR.(iget(124)>0).OR. &
928 (iget(078)>0).OR.(iget(079)>0).OR. &
929 (iget(125)>0).OR.(iget(145)>0).OR. &
930 (iget(140)>0).OR.(iget(040)>0).OR. &
931 (iget(181)>0).OR.(iget(182)>0).OR. &
932 (iget(199)>0).OR.(iget(185)>0).OR. &
933 (iget(186)>0).OR.(iget(187)>0).OR. &
934 (iget(250)>0).OR.(iget(252)>0).OR. &
935 (iget(276)>0).OR.(iget(277)>0).OR. &
936 (iget(750)>0).OR.(iget(751)>0).OR. &
937 (iget(752)>0).OR.(iget(754)>0).OR. &
938 (iget(278)>0).OR.(iget(264)>0).OR. &
939 (iget(450)>0).OR.(iget(480)>0).OR. &
940 (iget(479)>0).OR.(iget(481)>0).OR. &
941 (iget(774)>0).OR.(iget(747)>0).OR. &
942 (iget(464)>0).OR.(iget(467)>0).OR. &
943 (iget(470)>0).OR.(iget(476)>0).OR. &
944 (iget(629)>0).OR.(iget(630)>0).OR. &
945 (iget(909)>0).OR.(iget(737)>0).OR. &
946 (iget(742)>0).OR. &
947 (iget(994)>0).OR.(iget(995)>0) ) THEN
948
949 DO 190 l=1,lm
950
951! PRESSURE ON MDL SURFACES.
952 IF (iget(001)>0) THEN
953 IF (lvls(l,iget(001))>0) THEN
954 ll=lm-l+1
955!$omp parallel do private(i,j)
956 DO j=jsta,jend
957 DO i=ista,iend
958 grid1(i,j) = pmid(i,j,ll)
959 ENDDO
960 ENDDO
961 if(grib=="grib2" )then
962 cfld=cfld+1
963 fld_info(cfld)%ifld=iavblfld(iget(001))
964 fld_info(cfld)%lvl=lvlsxml(l,iget(001))
965!$omp parallel do private(i,j,ii,jj)
966 do j=1,jend-jsta+1
967 jj = jsta+j-1
968 do i=1,iend-ista+1
969 ii = ista+i-1
970 datapd(i,j,cfld) = grid1(ii,jj)
971 enddo
972 enddo
973 endif
974 ENDIF
975 ENDIF
976!
977!
978!--- CLOUD WATER on MDL SURFACE (Jin, '01; Ferrier, Feb '02)
979!
980 IF (iget(124) > 0) THEN
981 IF (lvls(l,iget(124)) > 0) THEN
982 ll=lm-l+1
983!$omp parallel do private(i,j)
984 DO j=jsta,jend
985 DO i=ista,iend
986 grid1(i,j) = qqw(i,j,ll)
987 if(grid1(i,j)<1e-20) grid1(i,j) = 0.0
988 ENDDO
989 ENDDO
990 if(grib=="grib2" )then
991 cfld=cfld+1
992 fld_info(cfld)%ifld=iavblfld(iget(124))
993 fld_info(cfld)%lvl=lvlsxml(l,iget(124))
994!$omp parallel do private(i,j,ii,jj)
995 do j=1,jend-jsta+1
996 jj = jsta+j-1
997 do i=1,iend-ista+1
998 ii = ista+i-1
999 datapd(i,j,cfld) = grid1(ii,jj)
1000 enddo
1001 enddo
1002 endif
1003 ENDIF
1004 ENDIF
1005!
1006!--- CLOUD ICE ON MDL SURFACE (Jin, '01; Ferrier, Feb '02)
1007!
1008 IF (iget(125) > 0) THEN
1009 IF (lvls(l,iget(125)) > 0) THEN
1010 ll=lm-l+1
1011!$omp parallel do private(i,j)
1012 DO j=jsta,jend
1013 DO i=ista,iend
1014 grid1(i,j) = qqi(i,j,ll)
1015 if(grid1(i,j)<1e-20) grid1(i,j) = 0.0
1016 ENDDO
1017 ENDDO
1018 if(grib=="grib2" )then
1019 cfld=cfld+1
1020 fld_info(cfld)%ifld=iavblfld(iget(125))
1021 fld_info(cfld)%lvl=lvlsxml(l,iget(125))
1022!$omp parallel do private(i,j,ii,jj)
1023 do j=1,jend-jsta+1
1024 jj = jsta+j-1
1025 do i=1,iend-ista+1
1026 ii = ista+i-1
1027 datapd(i,j,cfld) = grid1(ii,jj)
1028 enddo
1029 enddo
1030 endif
1031 ENDIF
1032 ENDIF
1033!
1034!--- RAIN ON MDL SURFACE (Jin, '01; Ferrier, Feb '02)
1035!
1036 IF (iget(181) > 0) THEN
1037 IF (lvls(l,iget(181)) > 0) THEN
1038 ll=lm-l+1
1039!$omp parallel do private(i,j)
1040 DO j=jsta,jend
1041 DO i=ista,iend
1042 grid1(i,j) = qqr(i,j,ll)
1043 if(grid1(i,j)<1e-20) grid1(i,j) = 0.0
1044 ENDDO
1045 ENDDO
1046 if(grib=="grib2" )then
1047 cfld=cfld+1
1048 fld_info(cfld)%ifld=iavblfld(iget(181))
1049 fld_info(cfld)%lvl=lvlsxml(l,iget(181))
1050!$omp parallel do private(i,j,ii,jj)
1051 do j=1,jend-jsta+1
1052 jj = jsta+j-1
1053 do i=1,iend-ista+1
1054 ii = ista+i-1
1055 datapd(i,j,cfld) = grid1(ii,jj)
1056 enddo
1057 enddo
1058 endif
1059 ENDIF
1060 ENDIF
1061!
1062!--- SNOW ON MDL SURFACE (Jin, '01; Ferrier, Feb '02)
1063!
1064 IF (iget(182) > 0) THEN
1065 IF (lvls(l,iget(182)) > 0)THEN
1066 ll=lm-l+1
1067!$omp parallel do private(i,j)
1068 DO j=jsta,jend
1069 DO i=ista,iend
1070 grid1(i,j) = qqs(i,j,ll)
1071 if(grid1(i,j)<1e-20) grid1(i,j) = 0.0
1072 ENDDO
1073 ENDDO
1074 if(grib=="grib2" )then
1075 cfld=cfld+1
1076 fld_info(cfld)%ifld=iavblfld(iget(182))
1077 fld_info(cfld)%lvl=lvlsxml(l,iget(182))
1078!$omp parallel do private(i,j,ii,jj)
1079 do j=1,jend-jsta+1
1080 jj = jsta+j-1
1081 do i=1,iend-ista+1
1082 ii = ista+i-1
1083 datapd(i,j,cfld) = grid1(ii,jj)
1084 enddo
1085 enddo
1086 endif
1087 ENDIF
1088 ENDIF
1089!
1090!--- GRAUPEL ON MDL SURFACE --tgs
1091!
1092 IF (iget(415) > 0) THEN
1093 IF (lvls(l,iget(415)) > 0)THEN
1094 ll=lm-l+1
1095!$omp parallel do private(i,j)
1096 DO j=jsta,jend
1097 DO i=ista,iend
1098 if(qqg(i,j,ll) < 1.e-12) qqg(i,j,ll) = 0. !tgs
1099 grid1(i,j) = qqg(i,j,ll)
1100 ENDDO
1101 ENDDO
1102 if(grib=="grib2" )then
1103 cfld=cfld+1
1104 fld_info(cfld)%ifld=iavblfld(iget(415))
1105 fld_info(cfld)%lvl=lvlsxml(l,iget(415))
1106!$omp parallel do private(i,j,ii,jj)
1107 do j=1,jend-jsta+1
1108 jj = jsta+j-1
1109 do i=1,iend-ista+1
1110 ii = ista+i-1
1111 datapd(i,j,cfld) = grid1(ii,jj)
1112 enddo
1113 enddo
1114 endif
1115 ENDIF
1116 ENDIF
1117!
1118!--- QNCLOUD ON MDL SURFACE --cra
1119!
1120 IF (iget(747) > 0) THEN
1121 IF (lvls(l,iget(747)) > 0)THEN
1122 ll=lm-l+1
1123!$omp parallel do private(i,j)
1124 DO j=jsta,jend
1125 DO i=ista,iend
1126 if(qqnw(i,j,ll) < 1.e-8) qqnw(i,j,ll) = 0. !tgs
1127 grid1(i,j) = qqnw(i,j,ll)
1128 ENDDO
1129 ENDDO
1130 if(grib=="grib2" )then
1131 cfld=cfld+1
1132 fld_info(cfld)%ifld=iavblfld(iget(747))
1133 fld_info(cfld)%lvl=lvlsxml(l,iget(747))
1134!$omp parallel do private(i,j,ii,jj)
1135 do j=1,jend-jsta+1
1136 jj = jsta+j-1
1137 do i=1,iend-ista+1
1138 ii = ista+i-1
1139 datapd(i,j,cfld) = grid1(ii,jj)
1140 enddo
1141 enddo
1142 endif
1143 ENDIF
1144 ENDIF
1145!
1146!--- QNICE ON MDL SURFACE --tgs
1147!
1148 IF (iget(752) > 0) THEN
1149 IF (lvls(l,iget(752)) > 0)THEN
1150 ll=lm-l+1
1151!$omp parallel do private(i,j)
1152 DO j=jsta,jend
1153 DO i=ista,iend
1154 if(qqni(i,j,ll) < 1.e-8) qqni(i,j,ll) = 0. !tgs
1155 grid1(i,j) = qqni(i,j,ll)
1156 ENDDO
1157 ENDDO
1158 if(grib=="grib2" )then
1159 cfld=cfld+1
1160 fld_info(cfld)%ifld=iavblfld(iget(752))
1161 fld_info(cfld)%lvl=lvlsxml(l,iget(752))
1162!$omp parallel do private(i,j,ii,jj)
1163 do j=1,jend-jsta+1
1164 jj = jsta+j-1
1165 do i=1,iend-ista+1
1166 ii = ista+i-1
1167 datapd(i,j,cfld) = grid1(ii,jj)
1168 enddo
1169 enddo
1170 endif
1171 ENDIF
1172 ENDIF
1173!
1174!--- QNRAIN ON MDL SURFACE --tgs
1175!
1176 IF (iget(754) > 0) THEN
1177 IF (lvls(l,iget(754)) > 0)THEN
1178 ll=lm-l+1
1179!$omp parallel do private(i,j)
1180 DO j=jsta,jend
1181 DO i=ista,iend
1182 if(qqnr(i,j,ll) < 1.e-8) qqnr(i,j,ll) = 0. !tgs
1183 grid1(i,j) = qqnr(i,j,ll)
1184 ENDDO
1185 ENDDO
1186 if(grib=="grib2" )then
1187 cfld=cfld+1
1188 fld_info(cfld)%ifld=iavblfld(iget(754))
1189 fld_info(cfld)%lvl=lvlsxml(l,iget(754))
1190!$omp parallel do private(i,j,ii,jj)
1191 do j=1,jend-jsta+1
1192 jj = jsta+j-1
1193 do i=1,iend-ista+1
1194 ii = ista+i-1
1195 datapd(i,j,cfld) = grid1(ii,jj)
1196 enddo
1197 enddo
1198 endif
1199 ENDIF
1200 ENDIF
1201! QNWFA ON MDL SURFACE --tgs
1202!
1203 IF (iget(766) > 0) THEN
1204 IF (lvls(l,iget(766)) > 0)THEN
1205 ll=lm-l+1
1206 DO j=jsta,jend
1207 DO i=ista,iend
1208 if(qqnwfa(i,j,ll)<1.e-8)qqnwfa(i,j,ll)=0. !tgs
1209 grid1(i,j)=qqnwfa(i,j,ll)
1210 ENDDO
1211 ENDDO
1212 if(grib=="grib2" )then
1213 cfld=cfld+1
1214 fld_info(cfld)%ifld=iavblfld(iget(766))
1215 fld_info(cfld)%lvl=lvlsxml(l,iget(766))
1216 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1217 endif
1218 ENDIF
1219 ENDIF
1220!
1221!--- QNIFA ON MDL SURFACE --tgs
1222!
1223 IF (iget(767) > 0) THEN
1224 IF (lvls(l,iget(767)) > 0)THEN
1225 ll=lm-l+1
1226 DO j=jsta,jend
1227 DO i=ista,iend
1228 if(qqnifa(i,j,ll)<1.e-8)qqnifa(i,j,ll)=0. !tgs
1229 grid1(i,j)=qqnifa(i,j,ll)
1230 ENDDO
1231 ENDDO
1232 if(grib=="grib2" )then
1233 cfld=cfld+1
1234 fld_info(cfld)%ifld=iavblfld(iget(767))
1235 fld_info(cfld)%lvl=lvlsxml(l,iget(767))
1236 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1237 endif
1238 ENDIF
1239 ENDIF
1240!
1241!--- Total cloud fraction on MDL surfaces. (Ferrier, Nov '04)
1242!
1243 IF (iget(145) > 0) THEN
1244 IF (lvls(l,iget(145)) > 0) THEN
1245 ll=lm-l+1
1246!$omp parallel do private(i,j)
1247 DO j=jsta,jend
1248 DO i=ista,iend
1249 IF(abs(cfr(i,j,ll)-spval) > small) THEN
1250 grid1(i,j) = cfr(i,j,ll)*h100
1251 ELSE
1252 grid1(i,j) = spval
1253 ENDIF
1254 ENDDO
1255 ENDDO
1256 CALL bound(grid1,d00,h100)
1257 if(grib=="grib2" )then
1258 cfld=cfld+1
1259 fld_info(cfld)%ifld=iavblfld(iget(145))
1260 fld_info(cfld)%lvl=lvlsxml(l,iget(145))
1261!$omp parallel do private(i,j,ii,jj)
1262 do j=1,jend-jsta+1
1263 jj = jsta+j-1
1264 do i=1,iend-ista+1
1265 ii = ista+i-1
1266 datapd(i,j,cfld) = grid1(ii,jj)
1267 enddo
1268 enddo
1269 endif
1270 ENDIF
1271 ENDIF
1272
1273!--- Model-state cloud fraction (unprocessed) on model surfaces (JSK, 8 Jan 2015)
1274!
1275 IF (iget(774) > 0) THEN
1276 IF (lvls(l,iget(774)) > 0) THEN
1277 ll=lm-l+1
1278!$omp parallel do private(i,j)
1279 DO j=jsta,jend
1280 DO i=ista,iend
1281 IF(modelname == 'RAPR') THEN
1282 grid1(i,j) = cfr(i,j,ll)
1283 ELSE
1284 grid1(i,j) = cfr_raw(i,j,ll)
1285 ENDIF
1286 ENDDO
1287 ENDDO
1288 if(grib=="grib2" )then
1289 cfld=cfld+1
1290 fld_info(cfld)%ifld=iavblfld(iget(774))
1291 fld_info(cfld)%lvl=lvlsxml(l,iget(774))
1292!$omp parallel do private(i,j,ii,jj)
1293 do j=1,jend-jsta+1
1294 jj = jsta+j-1
1295 do i=1,iend-ista+1
1296 ii = ista+i-1
1297 datapd(i,j,cfld) = grid1(ii,jj)
1298 enddo
1299 enddo
1300 endif
1301 ENDIF
1302 ENDIF
1303
1304!--- Equivalent radar reflectivity factor.
1305!
1306 IF (iget(250) > 0) THEN
1307 IF (lvls(l,iget(250)) > 0) THEN
1308 ll=lm-l+1
1309
1310!
1311! CRA Use WRF Thompson reflectivity diagnostic from RAPR model output
1312! Use unipost reflectivity diagnostic otherwise
1313!
1314! Chuang Feb 2015: use Thompson reflectivity direct output for all
1315! models
1316!
1317 IF(imp_physics == 8 .or. imp_physics == 28) THEN
1318!$omp parallel do private(i,j)
1319 DO j=jsta,jend
1320 DO i=ista,iend
1321 grid1(i,j) = ref_10cm(i,j,ll)
1322 ENDDO
1323 ENDDO
1324 ELSE
1325!$omp parallel do private(i,j)
1326 DO j=jsta,jend
1327 DO i=ista,iend
1328 grid1(i,j) = dbz(i,j,ll)
1329 ENDDO
1330 ENDDO
1331 ENDIF
1332! CRA
1333 CALL bound(grid1,dbzmin,dbzmax)
1334 if(grib=="grib2" )then
1335 cfld=cfld+1
1336 fld_info(cfld)%ifld=iavblfld(iget(250))
1337 fld_info(cfld)%lvl=lvlsxml(l,iget(250))
1338!$omp parallel do private(i,j,ii,jj)
1339 do j=1,jend-jsta+1
1340 jj = jsta+j-1
1341 do i=1,iend-ista+1
1342 ii = ista+i-1
1343 datapd(i,j,cfld) = grid1(ii,jj)
1344 enddo
1345 enddo
1346 endif
1347 ENDIF
1348 ENDIF
1349
1350!
1351!--- TOTAL CONDENSATE ON MDL SURFACE (CWM array; Ferrier, Feb '02)
1352!
1353 IF (iget(199)>0) THEN
1354 IF (lvls(l,iget(199))>0) THEN
1355 ll=lm-l+1
1356!$omp parallel do private(i,j)
1357 DO j=jsta,jend
1358 DO i=ista,iend
1359 grid1(i,j) = cwm(i,j,ll)
1360 ENDDO
1361 ENDDO
1362 if(grib=="grib2" )then
1363 cfld=cfld+1
1364 fld_info(cfld)%ifld=iavblfld(iget(199))
1365 fld_info(cfld)%lvl=lvlsxml(l,iget(199))
1366!$omp parallel do private(i,j,ii,jj)
1367 do j=1,jend-jsta+1
1368 jj = jsta+j-1
1369 do i=1,iend-ista+1
1370 ii = ista+i-1
1371 datapd(i,j,cfld) = grid1(ii,jj)
1372 enddo
1373 enddo
1374 endif
1375 ENDIF
1376 ENDIF
1377!
1378!--- F_rain ON MDL SURFACE (Jin, '01; Ferrier, Feb '02)
1379!
1380 IF (iget(185)>0) THEN
1381 IF (lvls(l,iget(185))>0) THEN
1382 ll=lm-l+1
1383!$omp parallel do private(i,j)
1384 DO j=jsta,jend
1385 DO i=ista,iend
1386 grid1(i,j) = f_rain(i,j,ll)
1387 ENDDO
1388 ENDDO
1389 if(grib=="grib2" )then
1390 cfld=cfld+1
1391 fld_info(cfld)%ifld=iavblfld(iget(185))
1392 fld_info(cfld)%lvl=lvlsxml(l,iget(185))
1393!$omp parallel do private(i,j,ii,jj)
1394 do j=1,jend-jsta+1
1395 jj = jsta+j-1
1396 do i=1,iend-ista+1
1397 ii = ista+i-1
1398 datapd(i,j,cfld) = grid1(ii,jj)
1399 enddo
1400 enddo
1401 endif
1402 ENDIF
1403 ENDIF
1404!
1405!--- F_ice ON MDL SURFACE (Jin, '01; Ferrier, Feb '02)
1406!
1407 IF (iget(186)>0) THEN
1408 IF (lvls(l,iget(186))>0) THEN
1409 ll=lm-l+1
1410!$omp parallel do private(i,j)
1411 DO j=jsta,jend
1412 DO i=ista,iend
1413 grid1(i,j) = f_ice(i,j,ll)
1414 ENDDO
1415 ENDDO
1416 if(grib=="grib2" )then
1417 cfld=cfld+1
1418 fld_info(cfld)%ifld=iavblfld(iget(186))
1419 fld_info(cfld)%lvl=lvlsxml(l,iget(186))
1420!$omp parallel do private(i,j,ii,jj)
1421 do j=1,jend-jsta+1
1422 jj = jsta+j-1
1423 do i=1,iend-ista+1
1424 ii = ista+i-1
1425 datapd(i,j,cfld) = grid1(ii,jj)
1426 enddo
1427 enddo
1428 endif
1429 ENDIF
1430 ENDIF
1431!
1432!--- F_RimeF ON MDL SURFACE (Jin, '01; Ferrier, Feb '02)
1433!
1434 IF (iget(187)>0) THEN
1435 IF (lvls(l,iget(187))>0) THEN
1436!--- Filter "rime factor" for non-zero precip rates and % frozen precip
1437 ll=lm-l+1
1438!$omp parallel do private(i,j)
1439 DO j=jsta,jend
1440 DO i=ista,iend
1441 grid1(i,j) = f_rimef(i,j,ll)
1442 ENDDO
1443 ENDDO
1444 if(grib=="grib2" )then
1445 cfld=cfld+1
1446 fld_info(cfld)%ifld=iavblfld(iget(187))
1447 fld_info(cfld)%lvl=lvlsxml(l,iget(187))
1448!$omp parallel do private(i,j,ii,jj)
1449 do j=1,jend-jsta+1
1450 jj = jsta+j-1
1451 do i=1,iend-ista+1
1452 ii = ista+i-1
1453 datapd(i,j,cfld) = grid1(ii,jj)
1454 enddo
1455 enddo
1456 endif
1457 ENDIF
1458
1459 ENDIF
1460!
1461! HEIGHTS ON MDL SURFACES.
1462 IF (iget(077)>0) THEN
1463 IF (lvls(l,iget(077))>0) THEN
1464 ll=lm-l+1
1465!$omp parallel do private(i,j)
1466 DO j=jsta,jend
1467 DO i=ista,iend
1468 grid1(i,j) = zmid(i,j,ll)
1469 ENDDO
1470 ENDDO
1471 if(grib=="grib2" )then
1472 cfld=cfld+1
1473 fld_info(cfld)%ifld=iavblfld(iget(077))
1474 fld_info(cfld)%lvl=lvlsxml(l,iget(077))
1475!$omp parallel do private(i,j,ii,jj)
1476 do j=1,jend-jsta+1
1477 jj = jsta+j-1
1478 do i=1,iend-ista+1
1479 ii = ista+i-1
1480 datapd(i,j,cfld) = grid1(ii,jj)
1481 enddo
1482 enddo
1483 endif
1484 ENDIF
1485
1486 ENDIF
1487!
1488! TEMPERATURE ON MDL SURFACES.
1489 IF (iget(002)>0) THEN
1490 IF (lvls(l,iget(002))>0) THEN
1491 ll=lm-l+1
1492!$omp parallel do private(i,j)
1493 DO j=jsta,jend
1494 DO i=ista,iend
1495 grid1(i,j) = t(i,j,ll)
1496 ENDDO
1497 ENDDO
1498 if(grib=="grib2" )then
1499 cfld=cfld+1
1500 fld_info(cfld)%ifld=iavblfld(iget(002))
1501 fld_info(cfld)%lvl=lvlsxml(l,iget(002))
1502!$omp parallel do private(i,j,ii,jj)
1503 do j=1,jend-jsta+1
1504 jj = jsta+j-1
1505 do i=1,iend-ista+1
1506 ii = ista+i-1
1507 datapd(i,j,cfld) = grid1(ii,jj)
1508 enddo
1509 enddo
1510 endif
1511 ENDIF
1512
1513 ENDIF
1514
1515! VIRTUAL TEMPERATURE ON MDL SURFACES.
1516 IF (iget(909)>0) THEN
1517 IF (lvls(l,iget(909))>0) THEN
1518 ll=lm-l+1
1519!$omp parallel do private(i,j)
1520 DO j=jsta,jend
1521 DO i=ista,iend
1522 IF(t(i,j,ll)<spval.and.q(i,j,ll)<spval)THEN
1523 grid1(i,j)=t(i,j,ll)*(1.+d608*q(i,j,ll))
1524 ELSE
1525 grid1(i,j)=spval
1526 ENDIF
1527 ENDDO
1528 ENDDO
1529 if(grib=="grib2" )then
1530 cfld=cfld+1
1531 fld_info(cfld)%ifld=iavblfld(iget(909))
1532 fld_info(cfld)%lvl=lvlsxml(l,iget(909))
1533 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1534 endif
1535 ENDIF
1536
1537 ENDIF
1538!
1539! POTENTIAL TEMPERATURE ON MDL SURFACES.
1540 IF (iget(003)>0) THEN
1541 IF (lvls(l,iget(003))>0) THEN
1542 ll=lm-l+1
1543!$omp parallel do private(i,j)
1544 DO j=jsta,jend
1545 DO i=ista,iend
1546 p1d(i,j) = pmid(i,j,ll)
1547 t1d(i,j) = t(i,j,ll)
1548 ENDDO
1549 ENDDO
1550 CALL calpot(p1d(ista:iend,jsta:jend),t1d(ista:iend,jsta:jend),egrid3(ista:iend,jsta:jend))
1551
1552!$omp parallel do private(i,j)
1553 DO j=jsta,jend
1554 DO i=ista,iend
1555 grid1(i,j) = egrid3(i,j)
1556 ENDDO
1557 ENDDO
1558 if(grib=="grib2") then
1559 cfld=cfld+1
1560 fld_info(cfld)%ifld=iavblfld(iget(003))
1561 fld_info(cfld)%lvl=lvlsxml(l,iget(003))
1562!$omp parallel do private(i,j,ii,jj)
1563 do j=1,jend-jsta+1
1564 jj = jsta+j-1
1565 do i=1,iend-ista+1
1566 ii = ista+i-1
1567 datapd(i,j,cfld) = grid1(ii,jj)
1568 enddo
1569 enddo
1570 endif
1571!
1572 ENDIF
1573 ENDIF
1574!
1575! VIRTUAL POTENTIAL TEMPERATURE ON MDL SURFACES.
1576 IF (iget(751)>0) THEN
1577 IF (lvls(l,iget(751))>0) THEN
1578 ll=lm-l+1
1579!$omp parallel do private(i,j)
1580 DO j=jsta,jend
1581 DO i=ista,iend
1582 p1d(i,j) = pmid(i,j,ll)
1583 t1d(i,j) = t(i,j,ll)
1584 ENDDO
1585 ENDDO
1586 CALL calpot(p1d(ista:iend,jsta:jend),t1d(ista:iend,jsta:jend),egrid3(ista:iend,jsta:jend))
1587
1588!$omp parallel do private(i,j)
1589 DO j=jsta,jend
1590 DO i=ista,iend
1591 IF(p1d(i,j)<spval.and.t1d(i,j)<spval.and.q(i,j,ll)<spval)THEN
1592 grid1(i,j) = egrid3(i,j) * (1.+d608*q(i,j,ll))
1593 ELSE
1594 grid1(i,j) = spval
1595 ENDIF
1596 ENDDO
1597 ENDDO
1598 if(grib=="grib2") then
1599 cfld=cfld+1
1600 fld_info(cfld)%ifld=iavblfld(iget(751))
1601 fld_info(cfld)%lvl=lvlsxml(l,iget(751))
1602!$omp parallel do private(i,j,ii,jj)
1603 do j=1,jend-jsta+1
1604 jj = jsta+j-1
1605 do i=1,iend-ista+1
1606 ii = ista+i-1
1607 datapd(i,j,cfld) = grid1(ii,jj)
1608 enddo
1609 enddo
1610 endif
1611 ENDIF
1612 ENDIF
1613
1614!
1615! RELATIVE HUMIDITY ON MDL SURFACES.
1616 item = -1
1617 IF (iget(006) > 0) item = lvls(l,iget(006))
1618 IF (item > 0 .OR. iget(450) > 0 .OR. iget(480) > 0 .OR. &
1619 iget(479) > 0 .OR. iget(481) > 0 ) THEN
1620 ll=lm-l+1
1621!$omp parallel do private(i,j)
1622 DO j=jsta,jend
1623 DO i=ista,iend
1624 p1d(i,j) = pmid(i,j,ll)
1625 t1d(i,j) = t(i,j,ll)
1626 q1d(i,j) = q(i,j,ll)
1627 ENDDO
1628 ENDDO
1629
1630 CALL calrh(p1d(ista:iend,jsta:jend),t1d(ista:iend,jsta:jend),q1d(ista:iend,jsta:jend),egrid4(ista:iend,jsta:jend))
1631
1632!$omp parallel do private(i,j)
1633 DO j=jsta,jend
1634 DO i=ista,iend
1635 IF(p1d(i,j)<spval.and.t1d(i,j)<spval.and.q1d(i,j)<spval)THEN
1636 grid1(i,j) = egrid4(i,j)*100.
1637 rh3d(i,j,ll) = grid1(i,j)
1638 egrid2(i,j) = q(i,j,ll)/max(1.e-8,egrid4(i,j)) ! Revert QS to compute cloud cover later
1639 ELSE
1640 grid1(i,j) = spval
1641 rh3d(i,j,ll) = spval
1642 egrid2(i,j) = spval
1643 ENDIF
1644 ENDDO
1645 ENDDO
1646 IF (item > 0) then
1647 if(grib=="grib2") then
1648 cfld=cfld+1
1649 fld_info(cfld)%ifld=iavblfld(iget(006))
1650 fld_info(cfld)%lvl=lvlsxml(l,iget(006))
1651!$omp parallel do private(i,j,ii,jj)
1652 do j=1,jend-jsta+1
1653 jj = jsta+j-1
1654 do i=1,iend-ista+1
1655 ii = ista+i-1
1656 datapd(i,j,cfld) = grid1(ii,jj)
1657 enddo
1658 enddo
1659 endif
1660 ENDIF
1661 ENDIF
1662
1663!
1664! DEWPOINT ON MDL SURFACES.
1665 IF (iget(004)>0) THEN
1666 IF (lvls(l,iget(004))>0) THEN
1667 ll=lm-l+1
1668!$omp parallel do private(i,j)
1669 DO j=jsta,jend
1670 DO i=ista,iend
1671 p1d(i,j) = pmid(i,j,ll)
1672 t1d(i,j) = t(i,j,ll)
1673 q1d(i,j) = q(i,j,ll)
1674 ENDDO
1675 ENDDO
1676 CALL caldwp(p1d(ista:iend,jsta:jend),q1d(ista:iend,jsta:jend),egrid3(ista:iend,jsta:jend),t1d(ista:iend,jsta:jend))
1677!$omp parallel do private(i,j)
1678 DO j=jsta,jend
1679 DO i=ista,iend
1680 IF(p1d(i,j)<spval.and.t1d(i,j)<spval.and.q1d(i,j)<spval)THEN
1681 grid1(i,j) = egrid3(i,j)
1682 ELSE
1683 grid1(i,j) = spval
1684 ENDIF
1685 ENDDO
1686 ENDDO
1687 if(grib=="grib2") then
1688 cfld=cfld+1
1689 fld_info(cfld)%ifld=iavblfld(iget(004))
1690 fld_info(cfld)%lvl=lvlsxml(l,iget(004))
1691!$omp parallel do private(i,j,ii,jj)
1692 do j=1,jend-jsta+1
1693 jj = jsta+j-1
1694 do i=1,iend-ista+1
1695 ii = ista+i-1
1696 datapd(i,j,cfld) = grid1(ii,jj)
1697 enddo
1698 enddo
1699 endif
1700 ENDIF
1701 ENDIF
1702!
1703! SPECIFIC HUMIDITY ON MDL SURFACES.
1704 IF (iget(005)>0) THEN
1705 IF (lvls(l,iget(005))>0) THEN
1706 ll=lm-l+1
1707!$omp parallel do private(i,j)
1708 DO j=jsta,jend
1709 DO i=ista,iend
1710 grid1(i,j) = q(i,j,ll)
1711 ENDDO
1712 ENDDO
1713 CALL bound(grid1,h1m12,h99999)
1714 if(grib=="grib2") then
1715 cfld=cfld+1
1716 fld_info(cfld)%ifld=iavblfld(iget(005))
1717 fld_info(cfld)%lvl=lvlsxml(l,iget(005))
1718!$omp parallel do private(i,j,ii,jj)
1719 do j=1,jend-jsta+1
1720 jj = jsta+j-1
1721 do i=1,iend-ista+1
1722 ii = ista+i-1
1723 datapd(i,j,cfld) = grid1(ii,jj)
1724 enddo
1725 enddo
1726 endif
1727 ENDIF
1728 ENDIF
1729!
1730! WATER VAPOR MIXING RATIO ON MDL SURFACES.
1731 IF (iget(750)>0) THEN
1732 IF (lvls(l,iget(750))>0) THEN
1733 ll=lm-l+1
1734!$omp parallel do private(i,j)
1735 DO j=jsta,jend
1736 DO i=ista,iend
1737 IF(q(i,j,ll)<spval)THEN
1738 grid1(i,j) = q(i,j,ll) / (1.-q(i,j,ll))
1739 ELSE
1740 grid1(i,j) = spval
1741 ENDIF
1742 ENDDO
1743 ENDDO
1744 CALL bound(grid1,h1m12,h99999)
1745 if(grib=="grib2") then
1746 cfld=cfld+1
1747 fld_info(cfld)%ifld=iavblfld(iget(750))
1748 fld_info(cfld)%lvl=lvlsxml(l,iget(750))
1749!$omp parallel do private(i,j,ii,jj)
1750 do j=1,jend-jsta+1
1751 jj = jsta+j-1
1752 do i=1,iend-ista+1
1753 ii = ista+i-1
1754 datapd(i,j,cfld) = grid1(ii,jj)
1755 enddo
1756 enddo
1757 endif
1758 ENDIF
1759 ENDIF
1760!
1761! MOISTURE CONVERGENCE ON MDL SURFACES.
1762 lll = 0
1763 if (iget(083) > 0) lll = lvls(l,iget(083))
1764 IF (iget(083)>0 .OR. iget(295)>0) THEN
1765 IF (lll >0 .OR. iget(295)>0) THEN
1766 ll=lm-l+1
1767!$omp parallel do private(i,j)
1768 DO j=jsta_2l,jend_2u
1769 DO i=ista_2l,iend_2u
1770 q1d(i,j) = q(i,j,ll)
1771 egrid1(i,j) = uh(i,j,ll)
1772 egrid2(i,j) = vh(i,j,ll)
1773 ENDDO
1774 ENDDO
1775 CALL calmcvg(q1d,egrid1,egrid2,egrid3)
1776!$omp parallel do private(i,j)
1777 DO j=jsta,jend
1778 DO i=ista,iend
1779 IF(q1d(i,j)<spval.and.egrid1(i,j)<spval.and.egrid2(i,j)<spval)THEN
1780 grid1(i,j) = egrid3(i,j)
1781 mcvg(i,j,ll) = egrid3(i,j)
1782 ELSE
1783 grid1(i,j) = spval
1784 mcvg(i,j,ll) = spval
1785 ENDIF
1786 ENDDO
1787 ENDDO
1788 IF(iget(083)>0 .AND. lll>0)THEN
1789 if(grib=="grib2") then
1790 cfld=cfld+1
1791 fld_info(cfld)%ifld=iavblfld(iget(083))
1792 fld_info(cfld)%lvl=lvlsxml(l,iget(083))
1793!$omp parallel do private(i,j,ii,jj)
1794 do j=1,jend-jsta+1
1795 jj = jsta+j-1
1796 do i=1,iend-ista+1
1797 ii = ista+i-1
1798 datapd(i,j,cfld) = grid1(ii,jj)
1799 enddo
1800 enddo
1801 endif
1802 ENDIF
1803 ENDIF
1804 ENDIF
1805!
1806! U AND/OR V WIND ON MDL SURFACES.
1807!MEB needs to be modified to do u at u-points and v at v-points
1808 IF (iget(007)>0.OR.iget(008)>0) THEN
1809 IF (lvls(l,iget(007))>0.OR.lvls(l,iget(008))>0 ) THEN
1810 ll=lm-l+1
1811!$omp parallel do private(i,j)
1812 DO j=jsta,jend
1813 DO i=ista,iend
1814 grid1(i,j) = uh(i,j,ll)
1815 grid2(i,j) = vh(i,j,ll)
1816 ENDDO
1817 ENDDO
1818 if(grib=="grib2") then
1819 cfld=cfld+1
1820 fld_info(cfld)%ifld=iavblfld(iget(007))
1821 fld_info(cfld)%lvl=lvlsxml(l,iget(007))
1822!$omp parallel do private(i,j,ii,jj)
1823 do j=1,jend-jsta+1
1824 jj = jsta+j-1
1825 do i=1,iend-ista+1
1826 ii = ista+i-1
1827 datapd(i,j,cfld) = grid1(ii,jj)
1828 enddo
1829 enddo
1830 cfld=cfld+1
1831 fld_info(cfld)%ifld=iavblfld(iget(008))
1832 fld_info(cfld)%lvl=lvlsxml(l,iget(008))
1833!$omp parallel do private(i,j,ii,jj)
1834 do j=1,jend-jsta+1
1835 jj = jsta+j-1
1836 do i=1,iend-ista+1
1837 ii = ista+i-1
1838 datapd(i,j,cfld) = grid2(ii,jj)
1839 enddo
1840 enddo
1841 endif
1842 ENDIF
1843 ENDIF
1844!
1845! OMEGA ON MDL SURFACES.
1846 IF (iget(009)>0) THEN
1847 IF (lvls(l,iget(009))>0) THEN
1848 ll=lm-l+1
1849!$omp parallel do private(i,j)
1850 DO j=jsta,jend
1851 DO i=ista,iend
1852 grid1(i,j) = omga(i,j,ll)
1853 ENDDO
1854 ENDDO
1855 if(grib=="grib2") then
1856 cfld=cfld+1
1857 fld_info(cfld)%ifld=iavblfld(iget(009))
1858 fld_info(cfld)%lvl=lvlsxml(l,iget(009))
1859!$omp parallel do private(i,j,ii,jj)
1860 do j=1,jend-jsta+1
1861 jj = jsta+j-1
1862 do i=1,iend-ista+1
1863 ii = ista+i-1
1864 datapd(i,j,cfld) = grid1(ii,jj)
1865 enddo
1866 enddo
1867 endif
1868 ENDIF
1869 ENDIF
1870!
1871! W ON MDL SURFACES.
1872 IF (iget(264)>0) THEN
1873 IF (lvls(l,iget(264))>0) THEN
1874 ll=lm-l+1
1875!$omp parallel do private(i,j)
1876 DO j=jsta,jend
1877 DO i=ista,iend
1878 grid1(i,j)=wh(i,j,ll)
1879 ENDDO
1880 ENDDO
1881 if(grib=="grib2") then
1882 cfld=cfld+1
1883 fld_info(cfld)%ifld=iavblfld(iget(264))
1884 fld_info(cfld)%lvl=lvlsxml(l,iget(264))
1885!$omp parallel do private(i,j,ii,jj)
1886 do j=1,jend-jsta+1
1887 jj = jsta+j-1
1888 do i=1,iend-ista+1
1889 ii = ista+i-1
1890 datapd(i,j,cfld) = grid1(ii,jj)
1891 enddo
1892 enddo
1893 endif
1894 ENDIF
1895 ENDIF
1896!
1897! ABSOLUTE VORTICITY ON MDL SURFACES.
1898 IF (iget(010)>0) THEN
1899 IF (lvls(l,iget(010))>0) THEN
1900 ll=lm-l+1
1901!$omp parallel do private(i,j)
1902 DO j=jsta_2l,jend_2u
1903 DO i=ista_2l,iend_2u
1904 egrid1(i,j) = uh(i,j,ll)
1905 egrid2(i,j) = vh(i,j,ll)
1906 ENDDO
1907 ENDDO
1908 CALL calvor(egrid1,egrid2,egrid3)
1909!$omp parallel do private(i,j)
1910 DO j=jsta,jend
1911 DO i=ista,iend
1912 IF(egrid3(i,j)<spval)THEN
1913 grid1(i,j) = egrid3(i,j)
1914 ELSE
1915 grid1(i,j) = spval
1916 ENDIF
1917 ENDDO
1918 ENDDO
1919 if(grib=="grib2") then
1920 cfld=cfld+1
1921 fld_info(cfld)%ifld=iavblfld(iget(010))
1922 fld_info(cfld)%lvl=lvlsxml(l,iget(010))
1923!$omp parallel do private(i,j,ii,jj)
1924 do j=1,jend-jsta+1
1925 jj = jsta+j-1
1926 do i=1,iend-ista+1
1927 ii = ista+i-1
1928 datapd(i,j,cfld) = grid1(ii,jj)
1929 enddo
1930 enddo
1931 endif
1932 ENDIF
1933 ENDIF
1934!
1935! GEOSTROPHIC STREAMFUNCTION ON MDL SURFACES.
1936 IF (iget(084)>0) THEN
1937 IF (lvls(l,iget(084))>0) THEN
1938 ll=lm-l+1
1939!$omp parallel do private(i,j)
1940 DO j=jsta,jend
1941 DO i=ista,iend
1942 egrid1(i,j) = zmid(i,j,ll)
1943 ENDDO
1944 ENDDO
1945 CALL calstrm(egrid1(ista:iend,jsta:jend),egrid2(ista:iend,jsta:jend))
1946!$omp parallel do private(i,j)
1947 DO j=jsta,jend
1948 DO i=ista,iend
1949 grid1(i,j) = egrid2(i,j)
1950 ENDDO
1951 ENDDO
1952 if(grib=="grib2") then
1953 cfld=cfld+1
1954 fld_info(cfld)%ifld=iavblfld(iget(084))
1955 fld_info(cfld)%lvl=lvlsxml(l,iget(084))
1956!$omp parallel do private(i,j,ii,jj)
1957 do j=1,jend-jsta+1
1958 jj = jsta+j-1
1959 do i=1,iend-ista+1
1960 ii = ista+i-1
1961 datapd(i,j,cfld) = grid1(ii,jj)
1962 enddo
1963 enddo
1964 endif
1965 ENDIF
1966 ENDIF
1967!
1968! TURBULENT KINETIC ENERGY ON MDL SURFACES.
1969 IF (iget(011)>0) THEN
1970 IF (lvls(l,iget(011))>0) THEN
1971 ll=lm-l+1
1972!$omp parallel do private(i,j)
1973 DO j=jsta,jend
1974 DO i=ista,iend
1975 grid1(i,j) = q2(i,j,ll)
1976 ENDDO
1977 ENDDO
1978 if(grib=="grib2") then
1979 cfld=cfld+1
1980 fld_info(cfld)%ifld=iavblfld(iget(011))
1981 fld_info(cfld)%lvl=lvlsxml(l,iget(011))
1982!$omp parallel do private(i,j,ii,jj)
1983 do j=1,jend-jsta+1
1984 jj = jsta+j-1
1985 do i=1,iend-ista+1
1986 ii = ista+i-1
1987 datapd(i,j,cfld) = grid1(ii,jj)
1988 enddo
1989 enddo
1990 endif
1991 ENDIF
1992 ENDIF
1993!
1994! CLOUD WATER CONTENT
1995!HC IF (IGET(124)>0) THEN
1996!HC IF (LVLS(L,IGET(124))>0) THEN
1997!HC DO J=JSTA,JEND
1998!HC DO I=ista,iend
1999!HC IF(CWM(I,J,L)<0..AND.CWM(I,J,L)>-1.E-10)
2000!HC 1 CWM(I,J,L)=0.
2001!HC GRID1(I,J)=CWM(I,J,L)
2002!HC ENDDO
2003!HC ENDDO
2004!HC ID(1:25) = 0
2005!HC CALL GRIBIT(IGET(124),L,GRIDista,iend,JM)
2006!HC ENDIF
2007!HC ENDIF
2008!
2009! CLOUD ICE CONTENT.
2010!commented out until QICE is brought into post
2011! IF (IGET(125)>0) THEN
2012! IF (LVLS(L,IGET(125))>0) THEN
2013! DO J=JSTA,JEND
2014! DO I=ista,iend
2015! GRID1(I,J)=QICE(I,J,L)
2016! ENDDO
2017! ENDDO
2018! ID(1:25) = 0
2019! CALL GRIBIT(IGET(125),L,GRIDista,iend,JM)
2020! ENDIF
2021! ENDIF
2022!
2023! CLOUD FRACTION
2024!
2025!commented out until CFRC is brought into post
2026! IF (IGET(145)>0) THEN
2027! IF (LVLS(L,IGET(145))>0) THEN
2028! DO J=JSTA,JEND
2029! DO I=ista,iend
2030! GRID1(I,J)=CFRC(I,J,L)
2031! ENDDO
2032! ENDDO
2033! ID(1:25) = 0
2034! CALL GRIBIT(IGET(145),L,GRIDista,iend,JM)
2035! ENDIF
2036! ENDIF
2037!
2038! TEMPERATURE TENDENCY DUE TO RADIATIVE FLUX CONVERGENCE
2039!commented out until TTND is brought into post
2040 IF (iget(140)>0) THEN
2041 IF (lvls(l,iget(140))>0) THEN
2042 ll=lm-l+1
2043!$omp parallel do private(i,j)
2044 DO j=jsta,jend
2045 DO i=ista,iend
2046 grid1(i,j) = ttnd(i,j,ll)
2047 ENDDO
2048 ENDDO
2049 if(grib=="grib2") then
2050 cfld=cfld+1
2051 fld_info(cfld)%ifld=iavblfld(iget(140))
2052 fld_info(cfld)%lvl=lvlsxml(l,iget(140))
2053!$omp parallel do private(i,j,ii,jj)
2054 do j=1,jend-jsta+1
2055 jj = jsta+j-1
2056 do i=1,iend-ista+1
2057 ii = ista+i-1
2058 datapd(i,j,cfld) = grid1(ii,jj)
2059 enddo
2060 enddo
2061 endif
2062 ENDIF
2063 ENDIF
2064!
2065! TEMPERATURE TENDENCY DUE TO SHORT WAVE RADIATION.
2066!commented out until RSWTT is brought into post
2067 IF (iget(040)>0) THEN
2068 IF (lvls(l,iget(040))>0) THEN
2069 ll=lm-l+1
2070!$omp parallel do private(i,j)
2071 DO j=jsta,jend
2072 DO i=ista,iend
2073 grid1(i,j) = rswtt(i,j,ll)
2074 ENDDO
2075 ENDDO
2076 if(grib=="grib2") then
2077 cfld=cfld+1
2078 fld_info(cfld)%ifld=iavblfld(iget(040))
2079 fld_info(cfld)%lvl=lvlsxml(l,iget(040))
2080!$omp parallel do private(i,j,ii,jj)
2081 do j=1,jend-jsta+1
2082 jj = jsta+j-1
2083 do i=1,iend-ista+1
2084 ii = ista+i-1
2085 datapd(i,j,cfld) = grid1(ii,jj)
2086 enddo
2087 enddo
2088 endif
2089 ENDIF
2090 ENDIF
2091!
2092! TEMPERATURE TENDENCY DUE TO LONG WAVE RADIATION.
2093!commented out until RLWTT is brought into post
2094 IF (iget(041)>0) THEN
2095 IF (lvls(l,iget(041))>0) THEN
2096 ll=lm-l+1
2097!$omp parallel do private(i,j)
2098 DO j=jsta,jend
2099 DO i=ista,iend
2100 grid1(i,j) = rlwtt(i,j,ll)
2101 ENDDO
2102 ENDDO
2103 if(grib=="grib2") then
2104 cfld=cfld+1
2105 fld_info(cfld)%ifld=iavblfld(iget(041))
2106 fld_info(cfld)%lvl=lvlsxml(l,iget(041))
2107!$omp parallel do private(i,j,ii,jj)
2108 do j=1,jend-jsta+1
2109 jj = jsta+j-1
2110 do i=1,iend-ista+1
2111 ii = ista+i-1
2112 datapd(i,j,cfld) = grid1(ii,jj)
2113 enddo
2114 enddo
2115 endif
2116 ENDIF
2117 ENDIF
2118!
2119!
2120! PROCESS NEXT MDL LEVEL.
2121!
2122! LATENT HEATING FROM GRID SCALE RAIN/EVAP. (TIME AVE)
2123 IF (iget(078)>0) THEN
2124 IF (lvls(l,iget(078))>0) THEN
2125 ll=lm-l+1
2126 IF(avrain>0.)THEN
2127 rrnum=1./avrain
2128 ELSE
2129 rrnum=0.
2130 ENDIF
2131!$omp parallel do private(i,j)
2132 DO j=jsta,jend
2133 DO i=ista,iend
2134 IF(train(i,j,ll)<spval)THEN
2135 grid1(i,j) = train(i,j,ll)*rrnum
2136 ELSE
2137 grid1(i,j) = spval
2138 ENDIF
2139 ENDDO
2140 ENDDO
2141 id(1:25) = 0
2142 itheat = int(theat)
2143 IF (itheat /= 0) THEN
2144 ifincr = mod(ifhr,itheat)
2145 ELSE
2146 ifincr=0
2147 END IF
2148 id(19) = ifhr
2149 IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2150 id(20) = 3
2151 IF (ifincr==0) THEN
2152 id(18) = ifhr-itheat
2153 ELSE
2154 id(18) = ifhr-ifincr
2155 ENDIF
2156 IF(ifmin >= 1)id(18)=id(18)*60
2157 if(grib=="grib2") then
2158 cfld=cfld+1
2159 fld_info(cfld)%ifld=iavblfld(iget(078))
2160 fld_info(cfld)%lvl=lvlsxml(l,iget(078))
2161 if(itheat==0) then
2162 fld_info(cfld)%ntrange=0
2163 else
2164 fld_info(cfld)%ntrange=1
2165 endif
2166 fld_info(cfld)%tinvstat=ifhr-id(18)
2167!$omp parallel do private(i,j,ii,jj)
2168 do j=1,jend-jsta+1
2169 jj = jsta+j-1
2170 do i=1,iend-ista+1
2171 ii = ista+i-1
2172 datapd(i,j,cfld) = grid1(ii,jj)
2173 enddo
2174 enddo
2175 endif
2176 END IF
2177 ENDIF
2178!
2179! LATENT HEATING FROM CONVECTION. (TIME AVE)
2180 IF (iget(079)>0) THEN
2181 IF (lvls(l,iget(079))>0) THEN
2182 ll=lm-l+1
2183 IF(avcnvc>0.)THEN
2184 rrnum=1./avcnvc
2185 ELSE
2186 rrnum=0.
2187 ENDIF
2188!$omp parallel do private(i,j)
2189 DO j=jsta,jend
2190 DO i=ista,iend
2191 IF(tcucn(i,j,ll)<spval)THEN
2192 grid1(i,j) = tcucn(i,j,ll)*rrnum
2193 ELSE
2194 grid1(i,j) = spval
2195 ENDIF
2196 ENDDO
2197 ENDDO
2198 id(1:25) = 0
2199 itheat = int(theat)
2200 IF (itheat /= 0) THEN
2201 ifincr = mod(ifhr,itheat)
2202 ELSE
2203 ifincr=0
2204 END IF
2205 id(19) = ifhr
2206 IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2207 id(20) = 3
2208 IF (ifincr==0) THEN
2209 id(18) = ifhr-itheat
2210 ELSE
2211 id(18) = ifhr-ifincr
2212 ENDIF
2213 IF(ifmin >= 1)id(18)=id(18)*60
2214 if(grib=="grib2") then
2215 cfld=cfld+1
2216 fld_info(cfld)%ifld=iavblfld(iget(079))
2217 fld_info(cfld)%lvl=lvlsxml(l,iget(079))
2218 if(itheat==0) then
2219 fld_info(cfld)%ntrange=0
2220 else
2221 fld_info(cfld)%ntrange=1
2222 endif
2223 fld_info(cfld)%tinvstat=ifhr-id(18)
2224!$omp parallel do private(i,j,ii,jj)
2225 do j=1,jend-jsta+1
2226 jj = jsta+j-1
2227 do i=1,iend-ista+1
2228 ii = ista+i-1
2229 datapd(i,j,cfld) = grid1(ii,jj)
2230 enddo
2231 enddo
2232 endif
2233 END IF
2234 ENDIF
2235!
2236! OZONE
2237 IF (iget(267)>0) THEN
2238 IF (lvls(l,iget(267))>0) THEN
2239 ll=lm-l+1
2240!$omp parallel do private(i,j)
2241 DO j=jsta,jend
2242 DO i=ista,iend
2243 grid1(i,j) = o3(i,j,ll)
2244 ENDDO
2245 ENDDO
2246 if(grib=="grib2") then
2247 cfld=cfld+1
2248 fld_info(cfld)%ifld=iavblfld(iget(267))
2249 fld_info(cfld)%lvl=lvlsxml(l,iget(267))
2250!$omp parallel do private(i,j,ii,jj)
2251 do j=1,jend-jsta+1
2252 jj = jsta+j-1
2253 do i=1,iend-ista+1
2254 ii = ista+i-1
2255 datapd(i,j,cfld) = grid1(ii,jj)
2256 enddo
2257 enddo
2258 endif
2259 END IF
2260 ENDIF
2261
2262!===============
2263! AQF
2264!===============
2265
2266 if (aqf_on) then
2267
2268 IF (iget(994)>0) THEN
2269 IF (lvls(l,iget(994))>0) THEN
2270 ll=lm-l+1
2271!$omp parallel do private(i,j)
2272 DO j=jsta,jend
2273 DO i=ista,iend
2274 IF(avgozcon(i,j,ll)<spval) THEN
2275 grid1(i,j) = avgozcon(i,j,ll) ! in ppb
2276 ELSE
2277 grid1(i,j) = spval
2278 ENDIF
2279 ENDDO
2280 ENDDO
2281 id(1:25) = 0
2282 itheat = int(theat)
2283 id(19) = ifhr
2284 id(20) = 3
2285 IF (ifhr==0) THEN
2286 id(18) = 0
2287 ELSE
2288 id(18) = ifhr-1
2289 ENDIF
2290 if(grib=="grib2") then
2291 cfld=cfld+1
2292 fld_info(cfld)%ifld=iavblfld(iget(994))
2293 fld_info(cfld)%lvl=lvlsxml(l,iget(994))
2294 if(ifhr==0) then
2295 fld_info(cfld)%ntrange=0
2296 else
2297 fld_info(cfld)%ntrange=1
2298 endif
2299 fld_info(cfld)%tinvstat=ifhr-id(18)
2300!$omp parallel do private(i,j,ii,jj)
2301 do j=1,jend-jsta+1
2302 jj = jsta+j-1
2303 do i=1,iend-ista+1
2304 ii = ista+i-1
2305 datapd(i,j,cfld) = grid1(ii,jj)
2306 enddo
2307 enddo
2308 endif
2309 END IF
2310 ENDIF
2311
2312
2313 !---- PM25 ----
2314
2315 IF (iget(995)>0) THEN
2316 IF (lvls(l,iget(995))>0) THEN
2317 ll=lm-l+1
2318!$omp parallel do private(i,j)
2319 DO j=jsta,jend
2320 DO i=ista,iend
2321 grid1(i,j) = avgpmtf(i,j,ll) !ug/m3
2322 ENDDO
2323 ENDDO
2324 id(1:25) = 0
2325 itheat = int(theat)
2326 id(19) = ifhr
2327 id(20) = 3
2328 IF (ifhr==0) THEN
2329 id(18) = 0
2330 ELSE
2331 id(18) = ifhr-1
2332 ENDIF
2333 if(grib=="grib2") then
2334 cfld=cfld+1
2335 fld_info(cfld)%ifld=iavblfld(iget(995))
2336 fld_info(cfld)%lvl=lvlsxml(l,iget(995))
2337 if(ifhr==0) then
2338 fld_info(cfld)%ntrange=0
2339 else
2340 fld_info(cfld)%ntrange=1
2341 endif
2342 fld_info(cfld)%tinvstat=ifhr-id(18)
2343!$omp parallel do private(i,j,ii,jj)
2344 do j=1,jend-jsta+1
2345 jj = jsta+j-1
2346 do i=1,iend-ista+1
2347 ii = ista+i-1
2348 datapd(i,j,cfld) = grid1(ii,jj)
2349 enddo
2350 enddo
2351 endif
2352 END IF
2353 ENDIF
2354
2355 endif ! -- aqfcmaq_on
2356
2357!===================================
2358
2359!
2360! E. James - 8 Dec 2017: SMOKE from WRF-CHEM
2361! SMOKE
2362 IF (iget(737)>0) THEN
2363 IF (lvls(l,iget(737))>0) THEN
2364 ll=lm-l+1
2365!$omp parallel do private(i,j)
2366 DO j=jsta,jend
2367 DO i=ista,iend
2368 IF(pmid(i,j,ll)<spval.and.t(i,j,ll)<spval.and.smoke(i,j,ll,1)<spval)THEN
2369 grid1(i,j) = (1./rd)*(pmid(i,j,ll)/t(i,j,ll))*smoke(i,j,ll,1)/(1e9)
2370 ELSE
2371 grid1(i,j) = spval
2372 ENDIF
2373 ENDDO
2374 ENDDO
2375 if(grib=="grib2") then
2376 cfld=cfld+1
2377 fld_info(cfld)%ifld=iavblfld(iget(737))
2378 fld_info(cfld)%lvl=lvlsxml(l,iget(737))
2379!$omp parallel do private(i,j,ii,jj)
2380 do j=1,jend-jsta+1
2381 jj = jsta+j-1
2382 do i=1,iend-ista+1
2383 ii = ista+i-1
2384 datapd(i,j,cfld) = grid1(ii,jj)
2385 enddo
2386 enddo
2387 endif
2388 END IF
2389 ENDIF
2390! E. James - 14 Sep 2022: Dust from RRFS
2391! DUST
2392 IF (iget(742)>0) THEN
2393 IF (lvls(l,iget(742))>0) THEN
2394 ll=lm-l+1
2395!$omp parallel do private(i,j)
2396 DO j=jsta,jend
2397 DO i=ista,iend
2398 IF(pmid(i,j,ll)<spval.and.t(i,j,ll)<spval.and.fv3dust(i,j,ll,1)<spval)THEN
2399 grid1(i,j) = (1./rd)*(pmid(i,j,ll)/t(i,j,ll))*fv3dust(i,j,ll,1)/(1e9)
2400 ELSE
2401 grid1(i,j) = spval
2402 ENDIF
2403 ENDDO
2404 ENDDO
2405 if(grib=="grib2") then
2406 cfld=cfld+1
2407 fld_info(cfld)%ifld=iavblfld(iget(742))
2408 fld_info(cfld)%lvl=lvlsxml(l,iget(742))
2409!$omp parallel do private(i,j,ii,jj)
2410 do j=1,jend-jsta+1
2411 jj = jsta+j-1
2412 do i=1,iend-ista+1
2413 ii = ista+i-1
2414 datapd(i,j,cfld) = grid1(ii,jj)
2415 enddo
2416 enddo
2417 endif
2418 END IF
2419 ENDIF
2420! E. James - 23 Feb 2023: COARSEPM from RRFS
2421! DUST
2422 IF (iget(1012)>0) THEN
2423 IF (lvls(l,iget(1012))>0) THEN
2424 ll=lm-l+1
2425!$omp parallel do private(i,j)
2426 DO j=jsta,jend
2427 DO i=ista,iend
2428 IF(pmid(i,j,ll)<spval.and.t(i,j,ll)<spval.and.coarsepm(i,j,ll,1)<spval)THEN
2429 grid1(i,j) = (1./rd)*(pmid(i,j,ll)/t(i,j,ll))*coarsepm(i,j,ll,1)/(1e9)
2430 ELSE
2431 grid1(i,j) = spval
2432 ENDIF
2433 ENDDO
2434 ENDDO
2435 if(grib=="grib2") then
2436 cfld=cfld+1
2437 fld_info(cfld)%ifld=iavblfld(iget(1012))
2438 fld_info(cfld)%lvl=lvlsxml(l,iget(1012))
2439!$omp parallel do private(i,j,ii,jj)
2440 do j=1,jend-jsta+1
2441 jj = jsta+j-1
2442 do i=1,iend-ista+1
2443 ii = ista+i-1
2444 datapd(i,j,cfld) = grid1(ii,jj)
2445 enddo
2446 enddo
2447 endif
2448 END IF
2449 ENDIF
2450!
2451! E. James - 23 Apr 2024: EBB from RRFS
2452!
2453 IF (iget(1015)>0) THEN
2454 IF (lvls(l,iget(1015))>0) THEN
2455 ll=lm-l+1
2456!$omp parallel do private(i,j)
2457 DO j=jsta,jend
2458 DO i=ista,iend
2459 grid1(i,j) = ebb(i,j,ll,1)/(1e9)
2460 ENDDO
2461 ENDDO
2462 if(grib=="grib2") then
2463 cfld=cfld+1
2464 fld_info(cfld)%ifld=iavblfld(iget(1015))
2465 fld_info(cfld)%lvl=lvlsxml(l,iget(1015))
2466!$omp parallel do private(i,j,ii,jj)
2467 do j=1,jend-jsta+1
2468 jj = jsta+j-1
2469 do i=1,iend-ista+1
2470 ii = ista+i-1
2471 datapd(i,j,cfld) = grid1(ii,jj)
2472 enddo
2473 enddo
2474 endif
2475 END IF
2476 ENDIF
2477!
2478 if ( gocart_on .or. gccpp_on .or. nasa_on ) then
2479! DUST 1
2480 IF (iget(629)>0) THEN
2481 IF (lvls(l,iget(629))>0) THEN
2482 ll=lm-l+1
2483!$omp parallel do private(i,j)
2484 DO j=jsta,jend
2485 DO i=ista,iend
2486 IF(dust(i,j,ll,1)<spval.and.rhomid(i,j,ll)<spval)THEN
2487 !GRID1(I,J) = DUST(I,J,LL,1)
2488 grid1(i,j) = dust(i,j,ll,1)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2489 ELSE
2490 grid1(i,j) = spval
2491 ENDIF
2492 ENDDO
2493 ENDDO
2494 if(grib=="grib2") then
2495 cfld=cfld+1
2496 fld_info(cfld)%ifld=iavblfld(iget(629))
2497 fld_info(cfld)%lvl=lvlsxml(l,iget(629))
2498!$omp parallel do private(i,j,ii,jj)
2499 do j=1,jend-jsta+1
2500 jj = jsta+j-1
2501 do i=1,iend-ista+1
2502 ii = ista+i-1
2503 datapd(i,j,cfld) = grid1(ii,jj)
2504 enddo
2505 enddo
2506 endif
2507 END IF
2508 ENDIF
2509
2510! DUST 2
2511 IF (iget(630)>0) THEN
2512 IF (lvls(l,iget(630))>0) THEN
2513 ll=lm-l+1
2514!$omp parallel do private(i,j)
2515 DO j=jsta,jend
2516 DO i=ista,iend
2517 IF(dust(i,j,ll,2)<spval.and.rhomid(i,j,ll)<spval)THEN
2518 !GRID1(I,J) = DUST(I,J,LL,2)
2519 grid1(i,j) = dust(i,j,ll,2)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2520 ELSE
2521 grid1(i,j) = spval
2522 ENDIF
2523 ENDDO
2524 ENDDO
2525 if(grib=="grib2") then
2526 cfld=cfld+1
2527 fld_info(cfld)%ifld=iavblfld(iget(630))
2528 fld_info(cfld)%lvl=lvlsxml(l,iget(630))
2529!$omp parallel do private(i,j,ii,jj)
2530 do j=1,jend-jsta+1
2531 jj = jsta+j-1
2532 do i=1,iend-ista+1
2533 ii = ista+i-1
2534 datapd(i,j,cfld) = grid1(ii,jj)
2535 enddo
2536 enddo
2537 endif
2538 END IF
2539 ENDIF
2540
2541! DUST 3
2542 IF (iget(631)>0) THEN
2543 IF (lvls(l,iget(631))>0) THEN
2544 ll=lm-l+1
2545!$omp parallel do private(i,j)
2546 DO j=jsta,jend
2547 DO i=ista,iend
2548 IF(dust(i,j,ll,3)<spval.and.rhomid(i,j,ll)<spval)THEN
2549 !GRID1(I,J) = DUST(I,J,LL,3)
2550 grid1(i,j) = dust(i,j,ll,3)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2551 ELSE
2552 grid1(i,j) = spval
2553 ENDIF
2554 ENDDO
2555 ENDDO
2556 if(grib=="grib2") then
2557 cfld=cfld+1
2558 fld_info(cfld)%ifld=iavblfld(iget(631))
2559 fld_info(cfld)%lvl=lvlsxml(l,iget(631))
2560!$omp parallel do private(i,j,ii,jj)
2561 do j=1,jend-jsta+1
2562 jj = jsta+j-1
2563 do i=1,iend-ista+1
2564 ii = ista+i-1
2565 datapd(i,j,cfld) = grid1(ii,jj)
2566 enddo
2567 enddo
2568 endif
2569 END IF
2570 ENDIF
2571
2572! DUST 4
2573 IF (iget(632)>0) THEN
2574 IF (lvls(l,iget(632))>0) THEN
2575 ll=lm-l+1
2576!$omp parallel do private(i,j)
2577 DO j=jsta,jend
2578 DO i=ista,iend
2579 IF(dust(i,j,ll,4)<spval.and.rhomid(i,j,ll)<spval)THEN
2580 !GRID1(I,J) = DUST(I,J,LL,4)
2581 grid1(i,j) = dust(i,j,ll,4)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2582 ELSE
2583 grid1(i,j) = spval
2584 ENDIF
2585 ENDDO
2586 ENDDO
2587 if(grib=="grib2") then
2588 cfld=cfld+1
2589 fld_info(cfld)%ifld=iavblfld(iget(632))
2590 fld_info(cfld)%lvl=lvlsxml(l,iget(632))
2591!$omp parallel do private(i,j,ii,jj)
2592 do j=1,jend-jsta+1
2593 jj = jsta+j-1
2594 do i=1,iend-ista+1
2595 ii = ista+i-1
2596 datapd(i,j,cfld) = grid1(ii,jj)
2597 enddo
2598 enddo
2599 endif
2600 END IF
2601 ENDIF
2602
2603! DUST 5
2604 IF (iget(633)>0) THEN
2605 IF (lvls(l,iget(633))>0) THEN
2606 ll=lm-l+1
2607!$omp parallel do private(i,j)
2608 DO j=jsta,jend
2609 DO i=ista,iend
2610 IF(dust(i,j,ll,5)<spval.and.rhomid(i,j,ll)<spval)THEN
2611 !GRID1(I,J) = DUST(I,J,LL,5)
2612 grid1(i,j) = dust(i,j,ll,5)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2613 ELSE
2614 grid1(i,j) = spval
2615 ENDIF
2616 ENDDO
2617 ENDDO
2618 if(grib=="grib2") then
2619 cfld=cfld+1
2620 fld_info(cfld)%ifld=iavblfld(iget(633))
2621 fld_info(cfld)%lvl=lvlsxml(l,iget(633))
2622!$omp parallel do private(i,j,ii,jj)
2623 do j=1,jend-jsta+1
2624 jj = jsta+j-1
2625 do i=1,iend-ista+1
2626 ii = ista+i-1
2627 datapd(i,j,cfld) = grid1(ii,jj)
2628 enddo
2629 enddo
2630 endif
2631 END IF
2632 ENDIF
2633
2634! SEASALT 1
2635 IF (iget(634)>0) THEN
2636 IF (lvls(l,iget(634))>0) THEN
2637 ll=lm-l+1
2638!$omp parallel do private(i,j)
2639 DO j=jsta,jend
2640 DO i=ista,iend
2641 IF(salt(i,j,ll,1)<spval.and.rhomid(i,j,ll)<spval)THEN
2642 grid1(i,j) = salt(i,j,ll,1)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2643 ELSE
2644 grid1(i,j) = spval
2645 ENDIF
2646 ENDDO
2647 ENDDO
2648 if(grib=="grib2") then
2649 cfld=cfld+1
2650 fld_info(cfld)%ifld=iavblfld(iget(634))
2651 fld_info(cfld)%lvl=lvlsxml(l,iget(634))
2652!$omp parallel do private(i,j,ii,jj)
2653 do j=1,jend-jsta+1
2654 jj = jsta+j-1
2655 do i=1,iend-ista+1
2656 ii = ista+i-1
2657 datapd(i,j,cfld) = grid1(ii,jj)
2658 enddo
2659 enddo
2660 endif
2661 END IF
2662 ENDIF
2663
2664! SEASALT 2
2665 IF (iget(635)>0) THEN
2666 IF (lvls(l,iget(635))>0) THEN
2667 ll=lm-l+1
2668!$omp parallel do private(i,j)
2669 DO j=jsta,jend
2670 DO i=ista,iend
2671 IF(salt(i,j,ll,2)<spval.and.rhomid(i,j,ll)<spval)THEN
2672 grid1(i,j) = salt(i,j,ll,2)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2673 ELSE
2674 grid1(i,j) = spval
2675 ENDIF
2676 ENDDO
2677 ENDDO
2678 if(grib=="grib2") then
2679 cfld=cfld+1
2680 fld_info(cfld)%ifld=iavblfld(iget(635))
2681 fld_info(cfld)%lvl=lvlsxml(l,iget(635))
2682!$omp parallel do private(i,j,ii,jj)
2683 do j=1,jend-jsta+1
2684 jj = jsta+j-1
2685 do i=1,iend-ista+1
2686 ii = ista+i-1
2687 datapd(i,j,cfld) = grid1(ii,jj)
2688 enddo
2689 enddo
2690 endif
2691 END IF
2692 ENDIF
2693
2694! SEASALT 3
2695 IF (iget(636)>0) THEN
2696 IF (lvls(l,iget(636))>0) THEN
2697 ll=lm-l+1
2698!$omp parallel do private(i,j)
2699 DO j=jsta,jend
2700 DO i=ista,iend
2701 IF(salt(i,j,ll,3)<spval.and.rhomid(i,j,ll)<spval)THEN
2702 grid1(i,j) = salt(i,j,ll,3)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2703 ELSE
2704 grid1(i,j) = spval
2705 ENDIF
2706 ENDDO
2707 ENDDO
2708 if(grib=="grib2") then
2709 cfld=cfld+1
2710 fld_info(cfld)%ifld=iavblfld(iget(636))
2711 fld_info(cfld)%lvl=lvlsxml(l,iget(636))
2712!$omp parallel do private(i,j,ii,jj)
2713 do j=1,jend-jsta+1
2714 jj = jsta+j-1
2715 do i=1,iend-ista+1
2716 ii = ista+i-1
2717 datapd(i,j,cfld) = grid1(ii,jj)
2718 enddo
2719 enddo
2720 endif
2721 END IF
2722 ENDIF
2723
2724! SEASALT 4
2725 IF (iget(637)>0) THEN
2726 IF (lvls(l,iget(637))>0) THEN
2727 ll=lm-l+1
2728!$omp parallel do private(i,j)
2729 DO j=jsta,jend
2730 DO i=ista,iend
2731 IF(salt(i,j,ll,4)<spval.and.rhomid(i,j,ll)<spval)THEN
2732 grid1(i,j) = salt(i,j,ll,4)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2733 ELSE
2734 grid1(i,j) = spval
2735 ENDIF
2736 ENDDO
2737 ENDDO
2738 if(grib=="grib2") then
2739 cfld=cfld+1
2740 fld_info(cfld)%ifld=iavblfld(iget(637))
2741 fld_info(cfld)%lvl=lvlsxml(l,iget(637))
2742!$omp parallel do private(i,j,ii,jj)
2743 do j=1,jend-jsta+1
2744 jj = jsta+j-1
2745 do i=1,iend-ista+1
2746 ii = ista+i-1
2747 datapd(i,j,cfld) = grid1(ii,jj)
2748 enddo
2749 enddo
2750 endif
2751 END IF
2752 ENDIF
2753
2754! SEASALT 0
2755 IF (iget(638)>0) THEN
2756 IF (lvls(l,iget(638))>0) THEN
2757 ll=lm-l+1
2758!$omp parallel do private(i,j)
2759 DO j=jsta,jend
2760 DO i=ista,iend
2761 IF(salt(i,j,ll,5)<spval.and.rhomid(i,j,ll)<spval)THEN
2762 grid1(i,j) = salt(i,j,ll,5)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2763 ELSE
2764 grid1(i,j) = spval
2765 ENDIF
2766 ENDDO
2767 ENDDO
2768 if(grib=="grib2") then
2769 cfld=cfld+1
2770 fld_info(cfld)%ifld=iavblfld(iget(638))
2771 fld_info(cfld)%lvl=lvlsxml(l,iget(638))
2772!$omp parallel do private(i,j,ii,jj)
2773 do j=1,jend-jsta+1
2774 jj = jsta+j-1
2775 do i=1,iend-ista+1
2776 ii = ista+i-1
2777 datapd(i,j,cfld) = grid1(ii,jj)
2778 enddo
2779 enddo
2780 endif
2781 END IF
2782 ENDIF
2783
2784! SULFATE
2785 IF (iget(639)>0) THEN
2786 IF (lvls(l,iget(639))>0) THEN
2787 ll=lm-l+1
2788!$omp parallel do private(i,j)
2789 DO j=jsta,jend
2790 DO i=ista,iend
2791 IF(suso(i,j,ll,1)<spval.and.rhomid(i,j,ll)<spval)THEN
2792 !GRID1(I,J) = SUSO(I,J,LL,1)
2793 grid1(i,j) = suso(i,j,ll,1)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2794 ELSE
2795 grid1(i,j) = spval
2796 ENDIF
2797 ENDDO
2798 ENDDO
2799 if(grib=="grib2") then
2800 cfld=cfld+1
2801 fld_info(cfld)%ifld=iavblfld(iget(639))
2802 fld_info(cfld)%lvl=lvlsxml(l,iget(639))
2803!$omp parallel do private(i,j,ii,jj)
2804 do j=1,jend-jsta+1
2805 jj = jsta+j-1
2806 do i=1,iend-ista+1
2807 ii = ista+i-1
2808 datapd(i,j,cfld) = grid1(ii,jj)
2809 enddo
2810 enddo
2811 endif
2812 END IF
2813 ENDIF
2814
2815! OC DRY (HYDROPHOBIC ORGANIC CARBON)
2816 IF (iget(640)>0) THEN
2817 IF (lvls(l,iget(640))>0) THEN
2818 ll=lm-l+1
2819!$omp parallel do private(i,j)
2820 DO j=jsta,jend
2821 DO i=ista,iend
2822 IF(waso(i,j,ll,1)<spval.and.rhomid(i,j,ll)<spval)THEN
2823 !GRID1(I,J) = WASO(I,J,LL,1)
2824 grid1(i,j) = waso(i,j,ll,1)*rhomid(i,j,ll) !lzhang
2825 ELSE
2826 grid1(i,j) = spval
2827 ENDIF
2828 ENDDO
2829 ENDDO
2830 if(grib=="grib2") then
2831 cfld=cfld+1
2832 fld_info(cfld)%ifld=iavblfld(iget(640))
2833 fld_info(cfld)%lvl=lvlsxml(l,iget(640))
2834!$omp parallel do private(i,j,ii,jj)
2835 do j=1,jend-jsta+1
2836 jj = jsta+j-1
2837 do i=1,iend-ista+1
2838 ii = ista+i-1
2839 datapd(i,j,cfld) = grid1(ii,jj)
2840 enddo
2841 enddo
2842 endif
2843 END IF
2844 ENDIF
2845
2846! OC WET (HYDROPHILIC ORGANIC CARBON)
2847 IF (iget(641)>0) THEN
2848 IF (lvls(l,iget(641))>0) THEN
2849 ll=lm-l+1
2850!$omp parallel do private(i,j)
2851 DO j=jsta,jend
2852 DO i=ista,iend
2853 IF(waso(i,j,ll,2)<spval.and.rhomid(i,j,ll)<spval)THEN
2854 !GRID1(I,J) = WASO(I,J,LL,2)
2855 grid1(i,j) = waso(i,j,ll,2)*rhomid(i,j,ll) !lzhang
2856 ELSE
2857 grid1(i,j) = spval
2858 ENDIF
2859 ENDDO
2860 ENDDO
2861 if(grib=="grib2") then
2862 cfld=cfld+1
2863 fld_info(cfld)%ifld=iavblfld(iget(641))
2864 fld_info(cfld)%lvl=lvlsxml(l,iget(641))
2865!$omp parallel do private(i,j,ii,jj)
2866 do j=1,jend-jsta+1
2867 jj = jsta+j-1
2868 do i=1,iend-ista+1
2869 ii = ista+i-1
2870 datapd(i,j,cfld) = grid1(ii,jj)
2871 enddo
2872 enddo
2873 endif
2874 END IF
2875 ENDIF
2876
2877! BC DRY (HYDROPHOBIC BLACK CARBON)
2878 IF (iget(642)>0) THEN
2879 IF (lvls(l,iget(642))>0) THEN
2880 ll=lm-l+1
2881!$omp parallel do private(i,j)
2882 DO j=jsta,jend
2883 DO i=ista,iend
2884 IF(soot(i,j,ll,1)<spval.and.rhomid(i,j,ll)<spval)THEN
2885 !GRID1(I,J) = SOOT(I,J,LL,1)
2886 grid1(i,j) = soot(i,j,ll,1)*rhomid(i,j,ll) !lzhang
2887 ELSE
2888 grid1(i,j) = spval
2889 ENDIF
2890 ENDDO
2891 ENDDO
2892 if(grib=="grib2") then
2893 cfld=cfld+1
2894 fld_info(cfld)%ifld=iavblfld(iget(642))
2895 fld_info(cfld)%lvl=lvlsxml(l,iget(642))
2896!$omp parallel do private(i,j,ii,jj)
2897 do j=1,jend-jsta+1
2898 jj = jsta+j-1
2899 do i=1,iend-ista+1
2900 ii = ista+i-1
2901 datapd(i,j,cfld) = grid1(ii,jj)
2902 enddo
2903 enddo
2904 endif
2905 END IF
2906 ENDIF
2907
2908! BC WET (HYDROPHILIC BLACK CARBON)
2909 IF (iget(643)>0) THEN
2910 IF (lvls(l,iget(643))>0) THEN
2911 ll=lm-l+1
2912!$omp parallel do private(i,j)
2913 DO j=jsta,jend
2914 DO i=ista,iend
2915 IF(soot(i,j,ll,2)<spval.and.rhomid(i,j,ll)<spval)THEN
2916 !GRID1(I,J) = SOOT(I,J,LL,2)
2917 grid1(i,j) = soot(i,j,ll,2)*rhomid(i,j,ll) !lzhang
2918 ELSE
2919 grid1(i,j) = spval
2920 ENDIF
2921 ENDDO
2922 ENDDO
2923 if(grib=="grib2") then
2924 cfld=cfld+1
2925 fld_info(cfld)%ifld=iavblfld(iget(643))
2926 fld_info(cfld)%lvl=lvlsxml(l,iget(643))
2927!$omp parallel do private(i,j,ii,jj)
2928 do j=1,jend-jsta+1
2929 jj = jsta+j-1
2930 do i=1,iend-ista+1
2931 ii = ista+i-1
2932 datapd(i,j,cfld) = grid1(ii,jj)
2933 enddo
2934 enddo
2935 endif
2936 END IF
2937 ENDIF
2938
2939
2940 if (nasa_on) then
2941! NITRATE
2942 IF (iget(688)>0) THEN
2943 IF (lvls(l,iget(688))>0) THEN
2944 ll=lm-l+1
2945!$omp parallel do private(i,j)
2946 DO j=jsta,jend
2947 DO i=ista,iend
2948 IF(no3(i,j,ll,1)<spval.and.rhomid(i,j,ll)<spval)THEN
2949 grid1(i,j) = no3(i,j,ll,1)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2950 ELSE
2951 grid1(i,j) = spval
2952 ENDIF
2953 ENDDO
2954 ENDDO
2955 if(grib=="grib2") then
2956 cfld=cfld+1
2957 fld_info(cfld)%ifld=iavblfld(iget(688))
2958 fld_info(cfld)%lvl=lvlsxml(l,iget(688))
2959!$omp parallel do private(i,j,ii,jj)
2960 do j=1,jend-jsta+1
2961 jj = jsta+j-1
2962 do i=1,iend-ista+1
2963 ii = ista+i-1
2964 datapd(i,j,cfld) = grid1(ii,jj)
2965 enddo
2966 enddo
2967 endif
2968 END IF
2969 ENDIF
2970
2971! NH4
2972 IF (iget(689)>0) THEN
2973 IF (lvls(l,iget(689))>0) THEN
2974 ll=lm-l+1
2975!$omp parallel do private(i,j)
2976 DO j=jsta,jend
2977 DO i=ista,iend
2978 IF(nh4(i,j,ll,1)<spval.and.rhomid(i,j,ll)<spval)THEN
2979 grid1(i,j) = nh4(i,j,ll,1)*rhomid(i,j,ll) !lzhang ug/kg-->ug/m3
2980 ELSE
2981 grid1(i,j) = spval
2982 ENDIF
2983 ENDDO
2984 ENDDO
2985 if(grib=="grib2") then
2986 cfld=cfld+1
2987 fld_info(cfld)%ifld=iavblfld(iget(689))
2988 fld_info(cfld)%lvl=lvlsxml(l,iget(689))
2989!$omp parallel do private(i,j,ii,jj)
2990 do j=1,jend-jsta+1
2991 jj = jsta+j-1
2992 do i=1,iend-ista+1
2993 ii = ista+i-1
2994 datapd(i,j,cfld) = grid1(ii,jj)
2995 enddo
2996 enddo
2997 endif
2998 END IF
2999 ENDIF
3000 endif !nasa_on
3001 endif !gocart_on
3002
3003! AIR DENSITY
3004 IF (iget(644)>0) THEN
3005 IF (lvls(l,iget(644))>0) THEN
3006 ll=lm-l+1
3007!$omp parallel do private(i,j)
3008 DO j=jsta,jend
3009 DO i=ista,iend
3010 grid1(i,j) = rhomid(i,j,ll)
3011 ENDDO
3012 ENDDO
3013 if(grib=="grib2") then
3014 cfld=cfld+1
3015 fld_info(cfld)%ifld=iavblfld(iget(644))
3016 fld_info(cfld)%lvl=lvlsxml(l,iget(644))
3017!$omp parallel do private(i,j,ii,jj)
3018 do j=1,jend-jsta+1
3019 jj = jsta+j-1
3020 do i=1,iend-ista+1
3021 ii = ista+i-1
3022 datapd(i,j,cfld) = grid1(ii,jj)
3023 enddo
3024 enddo
3025 endif
3026 END IF
3027 ENDIF
3028
3029! LAYER THICKNESS
3030 IF (iget(645)>0) THEN
3031 IF (lvls(l,iget(645))>0) THEN
3032 ll=lm-l+1
3033!$omp parallel do private(i,j)
3034 DO j=jsta,jend
3035 DO i=ista,iend
3036 grid1(i,j) = dpres(i,j,ll)
3037 ENDDO
3038 ENDDO
3039 if(grib=="grib2") then
3040 cfld=cfld+1
3041 fld_info(cfld)%ifld=iavblfld(iget(645))
3042 fld_info(cfld)%lvl=lvlsxml(l,iget(645))
3043!$omp parallel do private(i,j,ii,jj)
3044 do j=1,jend-jsta+1
3045 jj = jsta+j-1
3046 do i=1,iend-ista+1
3047 ii = ista+i-1
3048 datapd(i,j,cfld) = grid1(ii,jj)
3049 enddo
3050 enddo
3051 endif
3052 END IF
3053 ENDIF
3054
3055! CRA DUST FROM WRF CHEM: Removed ths section because GOCART can output
3056! the same fields above (Chuang 2012-03-07)
3057
3058! CRA
3059!
3060 190 CONTINUE
3061!
3062! END OF MDL SURFACE OUTPUT BLOCK.
3063!
3064 ENDIF
3065! VISIBILITY
3066! IF (IGET(180)>0) THEN
3067!comment out until we get QICE, QSNOW brought into post
3068!MEB RDTPHS= 1./(NPHS*DT)
3069!MEB modifying this Eta-specific code, assuming WRF physics will
3070!MEB explicitly predict vapor/water/ice/rain/snow
3071!MEB comments starting with MEB are lines associated with this
3072!MEB Eta-specific code
3073! NEED TO CALCULATE RAIN WATER AND SNOW MIXING RATIOS
3074! DO J=JSTA,JEND
3075! DO I=ista,iend
3076!MEB IF (PREC(I,J)==0) THEN
3077!MEB QSNO(I,J)=0.
3078!MEB QRAIN(I,J)=0.
3079!MEB ELSE
3080!MEB LLMH=LMH(I,J)
3081!MEB SNORATE=SR(I,J)*PREC(I,J)*RDTPHS
3082!MEB RAINRATE=(1-SR(I,J))*PREC(I,J)*RDTPHS
3083!MEB TERM1=(T(I,J,LM)/PSLP(I,J))**0.4167
3084!MEB TERM2=(T(I,J,LLMH)/PMID(I,J,LMH(I,J)))**0.5833
3085!MEB TERM3=RAINRATE**0.8333
3086!MEB QRAIN(I,J)=RAINCON*TERM1*TERM2*TERM3
3087!MEB TERM4=(T(I,J,LM)/PSLP(I,J))**0.47
3088!MEB TERM5=(T(I,J,LLMH)/PMID(I,J,LMH(I,J)))**0.53
3089!MEB TERM6=SNORATE**0.94
3090!MEB QSNO(I,J)=SNOCON*TERM4*TERM5*TERM6
3091!MEB ENDIF
3092! LLMH=NINT(LMH(I,J))
3093! QRAIN1(I,J)=QRAIN(I,J,LLMH)
3094! QSNO1(I,J)=QSNOW(I,J,LLMH)
3095! TT(I,J)=T(I,J,LLMH)
3096! QV(I,J)=Q(I,J,LLMH)
3097! QCD(I,J)=CWM(I,J,LLMH)
3098! QICE1(I,J)=QICE(I,J,LLMH)
3099! PPP(I,J)=PMID(I,J,LLMH)
3100! ENDDO
3101! ENDDO
3102! CALL CALVIS(QV,QCD,QRAIN1,QICE1,QSNO1,TT,PPP,VIS)
3103! DO J=JSTA,JEND
3104! DO I=ista,iend
3105! GRID1(I,J)=VIS(I,J)
3106! ENDDO
3107! ENDDO
3108! ID(1:25) = 0
3109! CALL GRIBIT(IGET(180),LVLS(1,IGET(180)),
3110! X GRIDista,iend,JM)
3111! ENDIF
3112!
3113! INSTANTANEOUS CONVECTIVE PRECIPITATION RATE.
3114!
3115! IF (IGET(249)>0) THEN
3116! RDTPHS=1000./DTQ2
3117! DO J=JSTA,JEND
3118! DO I=ista,iend
3119! GRID1(I,J)=CPRATE(I,J)*RDTPHS
3120! GRID1(I,J)=SPVAL
3121! ENDDO
3122! ENDDO
3123! ID(1:25) = 0
3124! CALL GRIBIT(IGET(249),LM,GRIDista,iend,JM)
3125! ENDIF
3126!
3127! COMPOSITE RADAR REFLECTIVITY (maximum dBZ in each column)
3128!
3129 IF (iget(252) > 0) THEN
3130 IF(imp_physics /= 8 .and. imp_physics /= 28) THEN
3131!$omp parallel do private(i,j,l)
3132 DO j=jsta,jend
3133 DO i=ista,iend
3134 grid1(i,j) = dbzmin
3135 DO l=1,nint(lmh(i,j))
3136 grid1(i,j) = max( grid1(i,j), dbz(i,j,l) )
3137 ENDDO
3138 ENDDO
3139 ENDDO
3140 ELSE
3141!tgs - for Thompson or Milbrandt scheme
3142!
3143! CRA Use WRF Thompson reflectivity diagnostic from RAPR model output
3144! Use unipost reflectivity diagnostic otherwise
3145!
3146 IF(imp_physics == 8 .or. imp_physics == 28) THEN
3147!NMMB does not have composite radar ref in model output
3148 IF(modelname=='NMM' .and. gridtype=='B' .or. &
3149 modelname=='NCAR'.or. modelname=='FV3R' .or. &
3150 modelname=='GFS' .or. &
3151 modelname=='NMM' .and. gridtype=='E')THEN
3152!$omp parallel do private(i,j,l)
3153 DO j=jsta,jend
3154 DO i=ista,iend
3155 grid1(i,j) = dbzmin
3156 DO l=1,nint(lmh(i,j))
3157 grid1(i,j) = max( grid1(i,j), ref_10cm(i,j,l) )
3158 ENDDO
3159 ENDDO
3160 ENDDO
3161 ELSE
3162!$omp parallel do private(i,j)
3163 DO j=jsta,jend
3164 DO i=ista,iend
3165 grid1(i,j) = refc_10cm(i,j)
3166 ENDDO
3167 ENDDO
3168 END IF
3169 CALL bound(grid1,dbzmin,dbzmax)
3170 ELSE
3171!$omp parallel do private(i,j)
3172 DO j=jsta,jend
3173 DO i=ista,iend
3174 grid1(i,j) = refl(i,j)
3175 ENDDO
3176 ENDDO
3177 ENDIF
3178! CRA
3179 ENDIF
3180 if(grib=="grib2") then
3181 cfld=cfld+1
3182 fld_info(cfld)%ifld=iavblfld(iget(252))
3183!$omp parallel do private(i,j,ii,jj)
3184 do j=1,jend-jsta+1
3185 jj = jsta+j-1
3186 do i=1,iend-ista+1
3187 ii = ista+i-1
3188 datapd(i,j,cfld) = grid1(ii,jj)
3189 enddo
3190 enddo
3191 endif
3192 ENDIF
3193!
3194! COMPUTE VIL (radar derived vertically integrated liquid water in each column)
3195! Per Mei Xu, VIL is radar derived vertically integrated liquid water based
3196! on emprical conversion factors (0.00344)
3197 IF (iget(581)>0) THEN
3198 DO j=jsta,jend
3199 DO i=ista,iend
3200 grid1(i,j)=0.0
3201 DO l=1,nint(lmh(i,j))
3202 if(zint(i,j,l) < spval .and.zint(i,j,l+1)<spval.and.dbz(i,j,l)<spval) then
3203 grid1(i,j)=grid1(i,j)+0.00344* &
3204 (10.**(dbz(i,j,l)/10.))**0.57143*(zint(i,j,l)-zint(i,j,l+1))/1000.
3205 else
3206 grid1(i,j)=spval
3207 endif
3208 ENDDO
3209 ENDDO
3210 ENDDO
3211 if(grib=="grib2") then
3212 cfld=cfld+1
3213 fld_info(cfld)%ifld=iavblfld(iget(581))
3214!$omp parallel do private(i,j,ii,jj)
3215 do j=1,jend-jsta+1
3216 jj = jsta+j-1
3217 do i=1,iend-ista+1
3218 ii = ista+i-1
3219 datapd(i,j,cfld) = grid1(ii,jj)
3220 enddo
3221 enddo
3222 endif
3223 ENDIF
3224!
3225!-- COMPOSITE RADAR REFLECTIVITY FROM RAIN (maximum dBZ in each column due to rain)
3226!
3227 IF (iget(276)>0) THEN
3228 DO j=jsta,jend
3229 DO i=ista,iend
3230 grid1(i,j)=dbzmin
3231 DO l=1,nint(lmh(i,j))
3232 grid1(i,j)=max( grid1(i,j), dbzr(i,j,l) )
3233 ENDDO
3234 ENDDO
3235 ENDDO
3236 if(grib=="grib2") then
3237 cfld=cfld+1
3238 fld_info(cfld)%ifld=iavblfld(iget(276))
3239!$omp parallel do private(i,j,ii,jj)
3240 do j=1,jend-jsta+1
3241 jj = jsta+j-1
3242 do i=1,iend-ista+1
3243 ii = ista+i-1
3244 datapd(i,j,cfld) = grid1(ii,jj)
3245 enddo
3246 enddo
3247 endif
3248 ENDIF
3249!
3250!-- COMPOSITE RADAR REFLECTIVITY FROM ICE
3251! (maximum dBZ in each column due to all ice habits; snow + graupel + etc.)
3252!
3253 IF (iget(277)>0) THEN
3254 DO j=jsta,jend
3255 DO i=ista,iend
3256 grid1(i,j)=dbzmin
3257 DO l=1,nint(lmh(i,j))
3258 grid1(i,j)=max( grid1(i,j), dbzi(i,j,l) )
3259 ENDDO
3260 ENDDO
3261 ENDDO
3262 if(grib=="grib2") then
3263 cfld=cfld+1
3264 fld_info(cfld)%ifld=iavblfld(iget(277))
3265!$omp parallel do private(i,j,ii,jj)
3266 do j=1,jend-jsta+1
3267 jj = jsta+j-1
3268 do i=1,iend-ista+1
3269 ii = ista+i-1
3270 datapd(i,j,cfld) = grid1(ii,jj)
3271 enddo
3272 enddo
3273 endif
3274 ENDIF
3275!
3276!-- COMPOSITE RADAR REFLECTIVITY FROM PARAMETERIZED CONVECTION
3277! (maximum dBZ in each column due to parameterized convection, as bogused into
3278! post assuming a constant reflectivity from the surface to the 0C level,
3279! and decreasing with height at higher levels)
3280!
3281 IF (iget(278)>0) THEN
3282 DO j=jsta,jend
3283 DO i=ista,iend
3284 grid1(i,j)=dbzmin
3285 DO l=1,nint(lmh(i,j))
3286 grid1(i,j)=max( grid1(i,j), dbzc(i,j,l) )
3287 ENDDO
3288 ENDDO
3289 ENDDO
3290 if(grib=="grib2") then
3291 cfld=cfld+1
3292 fld_info(cfld)%ifld=iavblfld(iget(278))
3293!$omp parallel do private(i,j,ii,jj)
3294 do j=1,jend-jsta+1
3295 jj = jsta+j-1
3296 do i=1,iend-ista+1
3297 ii = ista+i-1
3298 datapd(i,j,cfld) = grid1(ii,jj)
3299 enddo
3300 enddo
3301 endif
3302 ENDIF
3303
3304! SRD -- converted to kft
3305! J.Case, ENSCO Inc. (5/26/2008) -- Output Echo Tops (Highest HGT in meters
3306! of the 18-dBZ reflectivity on a model level)
3307
3308 IF (iget(426)>0) THEN
3309 DO j=jsta,jend
3310 DO i=ista,iend
3311 grid1(i,j)=0.0
3312 DO l=1,nint(lmh(i,j))
3313 IF (dbz(i,j,l)>=18.0) THEN
3314 grid1(i,j)=zmid(i,j,l)*3.2808/1000.
3315 EXIT
3316 ENDIF
3317 ENDDO
3318 ENDDO
3319 ENDDO
3320 if(grib=="grib2") then
3321 cfld=cfld+1
3322 fld_info(cfld)%ifld=iavblfld(iget(426))
3323!$omp parallel do private(i,j,ii,jj)
3324 do j=1,jend-jsta+1
3325 jj = jsta+j-1
3326 do i=1,iend-ista+1
3327 ii = ista+i-1
3328 datapd(i,j,cfld) = grid1(ii,jj)
3329 enddo
3330 enddo
3331 endif
3332 ENDIF
3333! J.Case (end mods)
3334! SRD
3335
3336! CRA
3337! NCAR fields
3338! Echo top height (Highest height in meters of 11-dBZ reflectivity
3339! interpolated from adjacent model levels in column containing 18-dBZ)
3340! Use WRF Thompson reflectivity diagnostic from RAPR model output
3341! Use unipost reflectivity diagnostic otherwise
3342!
3343 IF (iget(768) > 0) THEN
3344 IF(modelname == 'RAPR' .AND. (imp_physics == 8 .or. imp_physics == 28)) THEN
3345 DO j=jsta,jend
3346 DO i=ista,iend
3347 grid1(i,j) = -999.
3348 DO l=1,nint(lmh(i,j))
3349 IF (ref_10cm(i,j,l)>=18.0) THEN
3350 grid1(i,j)=zmid(i,j,l)
3351 EXIT
3352 ENDIF
3353 ENDDO
3354 IF(grid1(i,j) >= -900) THEN
3355 DO l=1,nint(lmh(i,j))
3356 IF (ref_10cm(i,j,l) >= 11.0) THEN
3357 IF(l == 1) THEN
3358 grid1(i,j) = zmid(i,j,l)
3359 ELSE IF(ref_10cm(i,j,l-1) == ref_10cm(i,j,l)) THEN
3360 grid1(i,j) = zmid(i,j,l)
3361 ELSE
3362 grid1(i,j) = zmid(i,j,l) + &
3363 (11.0 - ref_10cm(i,j,l)) * &
3364 (zmid(i,j,l-1) - zmid(i,j,l)) / &
3365 (ref_10cm(i,j,l-1) - ref_10cm(i,j,l))
3366 ENDIF
3367 EXIT
3368 ENDIF
3369 ENDDO
3370 ENDIF
3371 ENDDO
3372 ENDDO
3373 ELSE
3374 DO j=jsta,jend
3375 DO i=ista,iend
3376 grid1(i,j) = -999.
3377 DO l=1,nint(lmh(i,j))
3378 IF (dbz(i,j,l) >= 18.0) THEN
3379 grid1(i,j) = zmid(i,j,l)
3380 EXIT
3381 ENDIF
3382 ENDDO
3383 ENDDO
3384 ENDDO
3385 ENDIF
3386 if(grib=="grib2") then
3387 cfld=cfld+1
3388 fld_info(cfld)%ifld=iavblfld(iget(768))
3389!$omp parallel do private(i,j,ii,jj)
3390 do j=1,jend-jsta+1
3391 jj = jsta+j-1
3392 do i=1,iend-ista+1
3393 ii = ista+i-1
3394 datapd(i,j,cfld) = grid1(ii,jj)
3395 enddo
3396 enddo
3397 endif
3398 ENDIF
3399!
3400! Vertically integrated liquid in kg/m^2
3401!
3402 IF (iget(769)>0) THEN
3403 DO j=jsta,jend
3404 DO i=ista,iend
3405 grid1(i,j)=0.0
3406 DO l=1,nint(lmh(i,j))
3407 IF(qqr(i,j,l)<spval.and.qqs(i,j,l)<spval.and.qqg(i,j,l)<spval.and.&
3408 zint(i,j,l)<spval.and.zint(i,j,l+1)<spval.and.&
3409 pmid(i,j,l)<spval.and.t(i,j,l)<spval.and.q(i,j,l)<spval)THEN
3410 IF(qqh(i,j,l)<spval)THEN
3411 grid1(i,j)=grid1(i,j) + (qqr(i,j,l) + qqh(i,j,l) + &
3412 qqs(i,j,l) + qqg(i,j,l))* &
3413 (zint(i,j,l)-zint(i,j,l+1))*pmid(i,j,l)/ &
3414 (rd*t(i,j,l)*(q(i,j,l)*d608+1.0))
3415 ELSE
3416 grid1(i,j)=grid1(i,j) + (qqr(i,j,l) + &
3417 qqs(i,j,l) + qqg(i,j,l))* &
3418 (zint(i,j,l)-zint(i,j,l+1))*pmid(i,j,l)/ &
3419 (rd*t(i,j,l)*(q(i,j,l)*d608+1.0))
3420 ENDIF
3421 ELSE
3422 grid1(i,j)=spval
3423 ENDIF
3424 ENDDO
3425 ENDDO
3426 ENDDO
3427 if(grib=="grib2") then
3428 cfld=cfld+1
3429 fld_info(cfld)%ifld=iavblfld(iget(769))
3430!$omp parallel do private(i,j,ii,jj)
3431 do j=1,jend-jsta+1
3432 jj = jsta+j-1
3433 do i=1,iend-ista+1
3434 ii = ista+i-1
3435 datapd(i,j,cfld) = grid1(ii,jj)
3436 enddo
3437 enddo
3438 endif
3439 ENDIF
3440!
3441! Vertically integrated liquid based on reflectivity factor in kg/m^2
3442! Use WRF Thompson reflectivity diagnostic from RAPR model output
3443! Use unipost reflectivity diagnostic otherwise
3444!
3445 IF (iget(770) > 0) THEN
3446 IF(modelname == 'RAPR' .AND. (imp_physics == 8 .or. imp_physics == 28)) THEN
3447 DO j=jsta,jend
3448 DO i=ista,iend
3449 grid1(i,j) = 0.0
3450 DO l=1,nint(lmh(i,j))
3451 IF (ref_10cm(i,j,l) > -10.0 ) THEN
3452 grid1(i,j) = grid1(i,j) + 0.00344 * &
3453 (10.**(ref_10cm(i,j,l)/10.))**0.57143 * &
3454 (zint(i,j,l)-zint(i,j,l+1))/1000.
3455 ENDIF
3456 ENDDO
3457 ENDDO
3458 ENDDO
3459 ELSE
3460 DO j=jsta,jend
3461 DO i=ista,iend
3462 grid1(i,j) = 0.0
3463 DO l=1,nint(lmh(i,j))
3464 grid1(i,j) = grid1(i,j) + 0.00344 * &
3465 (10.**(dbz(i,j,l)/10.))**0.57143 * &
3466 (zint(i,j,l)-zint(i,j,l+1))/1000.
3467 ENDDO
3468 ENDDO
3469 ENDDO
3470 ENDIF
3471 if(grib=="grib2") then
3472 cfld=cfld+1
3473 fld_info(cfld)%ifld=iavblfld(iget(770))
3474!$omp parallel do private(i,j,ii,jj)
3475 do j=1,jend-jsta+1
3476 jj = jsta+j-1
3477 do i=1,iend-ista+1
3478 ii = ista+i-1
3479 datapd(i,j,cfld) = grid1(ii,jj)
3480 enddo
3481 enddo
3482 endif
3483 ENDIF
3484! CRA
3485
3486!
3487!--- VISIBILITY
3488!
3489 IF (iget(180)>0) THEN
3490 rdtphs=1./dtq2
3491 !
3492 !--- Needed values at 1st level above ground (Jin, '01; Ferrier, Feb '02)
3493 !
3494 DO j=jsta,jend
3495 DO i=ista,iend
3496 llmh=nint(lmh(i,j))
3497 q1d(i,j)=q(i,j,llmh)
3498 if(q1d(i,j)<=0.) q1d(i,j)=0. !tgs
3499 qw1(i,j)=qqw(i,j,llmh)
3500 qr1(i,j)=qqr(i,j,llmh)
3501 qi1(i,j)=qqi(i,j,llmh)
3502 qs1(i,j)=qqs(i,j,llmh)
3503 qg1(i,j)=qqg(i,j,llmh) !tgs
3504 t1d(i,j)=t(i,j,llmh)
3505 p1d(i,j)=pmid(i,j,llmh)
3506
3507!HC July 2012, per communication with Ferrier, modify post to add convective
3508! contribution to visibility for all non GFS models
3509
3510! IF(MODELNAME/='GFS')THEN
3511 IF(imp_physics/=99)THEN
3512 IF (cprate(i,j) > 0. .and. cprate(i,j) < spval &
3513 .and. pmid(i,j,lm) < spval .and. qr1(i,j) < spval) THEN
3514! IF (CUPPT(I,J) > 0.) THEN
3515 rainrate=(1-sr(i,j))*cprate(i,j)*rdtphs
3516! RAINRATE=(1-SR(I,J))*CUPPT(I,J)/(TRDLW*3600.)
3517 term1=(t(i,j,lm)/pmid(i,j,lm))**0.4167
3518 term2=(t1d(i,j)/p1d(i,j))**0.5833
3519 term3=rainrate**0.8333
3520 qrold=1.2*qr1(i,j)
3521 qr1(i,j)=qr1(i,j)+raincon*term1*term2*term3
3522 IF (sr(i,j) > 0. .and. qs1(i,j) < spval) THEN
3523 snorate=sr(i,j)*cprate(i,j)*rdtphs
3524! SNORATE=SR(I,J)*CUPPT(I,J)/(TRDLW*3600.)
3525 term1=(t(i,j,lm)/pmid(i,j,lm))**0.47
3526 term2=(t1d(i,j)/p1d(i,j))**0.53
3527 term3=snorate**0.94
3528 qs1(i,j)=qs1(i,j)+snocon*term1*term2*term3
3529 ENDIF
3530 ENDIF
3531 ELSE !imp_physics is 99
3532! Zhao microphysics option in NMMB is identified as 9
3533! However, microphysics option 9 in WRF is Milbrandt-Yau 2-moment scheme.
3534! 3/14/2013: Ratko comitted NEMS change (r26409) to change mp_physics from 9 to 99 for Zhao
3535! scheme used with NMMB. Post is changing accordingly
3536! IF(imp_physics==99)THEN ! use rain rate for visibility
3537 IF (prec(i,j) < spval .and. prec(i,j) > 0. .and. &
3538 sr(i,j)<spval) THEN
3539! IF (CUPPT(I,J) > 0.) THEN
3540 rainrate=(1-sr(i,j))*prec(i,j)*rdtphs
3541! RAINRATE=(1-SR(I,J))*CUPPT(I,J)/(TRDLW*3600.)
3542 term1=(t(i,j,lm)/pmid(i,j,lm))**0.4167
3543 term2=(t1d(i,j)/p1d(i,j))**0.5833
3544 term3=rainrate**0.8333
3545 qrold=1.2*qr1(i,j)
3546 qr1(i,j)=qr1(i,j)+raincon*term1*term2*term3
3547 IF (sr(i,j) > 0.) THEN
3548 snorate=sr(i,j)*prec(i,j)*rdtphs
3549! SNORATE=SR(I,J)*CUPPT(I,J)/(TRDLW*3600.)
3550 term1=(t(i,j,lm)/pmid(i,j,lm))**0.47
3551 term2=(t1d(i,j)/p1d(i,j))**0.53
3552 term3=snorate**0.94
3553 qs1(i,j)=qs1(i,j)+snocon*term1*term2*term3
3554 ENDIF
3555 ENDIF
3556 END IF
3557
3558 ENDDO
3559 ENDDO
3560!
3561!-- Visibility using Warner-Stoelinga algorithm (Jin, '01)
3562!
3563 ii=(ista+iend)/2
3564 jj=(jsta+jend)/2
3565! print*,'Debug: Visbility ',Q1D(ii,jj),QW1(ii,jj),QR1(ii,jj)
3566! +,QI1(ii,jj) ,QS1(ii,jj),T1D(ii,jj),P1D(ii,jj)
3567
3568 CALL calvis(q1d,qw1,qr1,qi1,qs1,t1d,p1d,vis)
3569
3570! print*,'Debug: Visbility ',Q1D(ii,jj),QW1(ii,jj),QR1(ii,jj),QI1(ii,jj)
3571! +,QS1(ii,jj),T1D(ii,jj),P1D(ii,jj)
3572!
3573
3574 DO j=jsta,jend
3575 DO i=ista,iend
3576 IF(vis(i,j)/=spval.and.abs(vis(i,j))>24135.1)print*,'bad visbility' &
3577 , i,j,q1d(i,j),qw1(i,j),qr1(i,j),qi1(i,j) &
3578 , qs1(i,j),t1d(i,j),p1d(i,j),vis(i,j)
3579
3580 grid1(i,j)=vis(i,j)
3581 END DO
3582 END DO
3583 if(grib=="grib2") then
3584 cfld=cfld+1
3585 fld_info(cfld)%ifld=iavblfld(iget(180))
3586 fld_info(cfld)%lvl=lvlsxml(1,iget(180))
3587 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3588 endif
3589 ENDIF
3590
3591!
3592! --- GSD VISIBILITY
3593!
3594 IF (iget(410)>0) THEN
3595 CALL calvis_gsd(czen,vis)
3596 DO j=jsta,jend
3597 DO i=ista,iend
3598 grid1(i,j)=vis(i,j)
3599 END DO
3600 END DO
3601 if(grib=="grib2") then
3602 cfld=cfld+1
3603 fld_info(cfld)%ifld=iavblfld(iget(410))
3604 fld_info(cfld)%lvl=lvlsxml(1,iget(410))
3605 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3606 endif
3607 ENDIF
3608!
3609! --- RADAR REFLECT - 1km
3610!
3611 IF (iget(748) > 0) THEN
3612!
3613! CRA Use WRF Thompson reflectivity diagnostic from RAPR model output
3614! Use unipost reflectivity diagnostic otherwise
3615!
3616 IF(modelname == 'RAPR' .AND. (imp_physics == 8 .or. imp_physics == 28)) THEN
3617 grid1 = -20.0
3618!$omp parallel do private(i,j)
3619 DO j=jsta,jend
3620 DO i=ista,iend
3621 grid1(i,j) = ref1km_10cm(i,j)
3622 END DO
3623 END DO
3624 CALL bound(grid1,dbzmin,dbzmax)
3625 ELSE
3626!$omp parallel do private(i,j)
3627 DO j=jsta,jend
3628 DO i=ista,iend
3629 grid1(i,j) = refl1km(i,j)
3630 END DO
3631 END DO
3632 ENDIF
3633! CRA
3634! print *,'MAX/MIN radar reflct - 1km ',maxval(grid1),minval(grid1)
3635 if(grib=="grib2") then
3636 cfld=cfld+1
3637 fld_info(cfld)%ifld=iavblfld(iget(748))
3638 fld_info(cfld)%lvl=lvlsxml(1,iget(748))
3639 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3640 endif
3641 ENDIF
3642
3643!
3644! --- RADAR REFLECT - 4km
3645!
3646 IF (iget(757) > 0) THEN
3647!
3648! CRA Use WRF Thompson reflectivity diagnostic from RAPR model output
3649! Use unipost reflectivity diagnostic otherwise
3650!
3651 IF(modelname == 'RAPR' .AND. (imp_physics == 8 .or. imp_physics == 28)) THEN
3652!$omp parallel do private(i,j)
3653 DO j=jsta,jend
3654 DO i=ista,iend
3655 grid1(i,j) = ref4km_10cm(i,j)
3656 END DO
3657 END DO
3658 CALL bound(grid1,dbzmin,dbzmax)
3659 ELSE
3660!$omp parallel do private(i,j)
3661 DO j=jsta,jend
3662 DO i=ista,iend
3663 grid1(i,j) = refl4km(i,j)
3664 END DO
3665 END DO
3666 ENDIF
3667! CRA
3668! print *,'MAX/MIN radar reflct - 4km ',maxval(grid1),minval(grid1)
3669 if(grib=="grib2") then
3670 cfld=cfld+1
3671 fld_info(cfld)%ifld=iavblfld(iget(757))
3672 fld_info(cfld)%lvl=lvlsxml(1,iget(757))
3673 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3674 endif
3675 ENDIF
3676
3677! RADAR REFLECTIVITY AT -10C LEVEL
3678 IF (iget(912)>0) THEN
3679 zm10c=spval
3680 DO j=jsta,jend
3681 DO i=ista,iend
3682! dong handle missing value
3683 if (slp(i,j) < spval) then
3684 zm10c(i,j)=zmid(i,j,nint(lmh(i,j)))
3685 DO l=nint(lmh(i,j)),1,-1
3686 IF (t(i,j,l) <= 263.15) THEN
3687 zm10c(i,j)= l !-- Find lowest level where T<-10C
3688 EXIT
3689 ENDIF
3690 ENDDO
3691 end if ! spval
3692 ENDDO
3693 ENDDO
3694
3695! REFD at -10 C level
3696!
3697! CRA Use WRF Thompson reflectivity diagnostic from RAPR model output
3698! Use unipost reflectivity diagnostic otherwise
3699! Chuang: use Thompson reflectivity direct output for all
3700! models
3701!
3702 IF(imp_physics==8 .or. imp_physics==28) THEN
3703!$omp parallel do private(i,j)
3704 DO j=jsta,jend
3705 DO i=ista,iend
3706 grid1(i,j)=spval
3707! dong handle missing value
3708 if (slp(i,j) < spval) then
3709 grid1(i,j)=ref_10cm(i,j,zm10c(i,j))
3710 end if ! spval
3711 ENDDO
3712 ENDDO
3713 ELSE
3714!$omp parallel do private(i,j)
3715 DO j=jsta,jend
3716 DO i=ista,iend
3717 grid1(i,j)=spval
3718! dong handle missing value
3719 if (slp(i,j) < spval) then
3720 grid1(i,j)=dbz(i,j,zm10c(i,j))
3721 end if ! spval
3722 ENDDO
3723 ENDDO
3724 ENDIF
3725
3726 CALL bound(grid1,dbzmin,dbzmax)
3727
3728 if(grib=="grib2" )then
3729 cfld=cfld+1
3730 fld_info(cfld)%ifld=iavblfld(iget(912))
3731 fld_info(cfld)%lvl=lvlsxml(l,iget(912))
3732 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3733 endif
3734 ENDIF
3735!
3736! ASYMPTOTIC AND FREE ATMOSPHERE MASTER LENGTH SCALE (EL), PLUS
3737! GRADIENT RICHARDSON NUMBER.
3738!
3739 IF ( (iget(111)>0) .OR. (iget(146)>0) .OR. &
3740 (iget(147)>0) ) THEN
3741!
3742! COMPUTE ASYMPTOTIC MASTER LENGTH SCALE.
3743 CALL clmax(el0(1,jsta),egrid2(1,jsta),egrid3(1,jsta),egrid4(1,jsta),egrid5(1,jsta))
3744!
3745! IF REQUESTED, POST ASYMPTOTIC MASTER LENGTH SCALE.
3746 IF (iget(147)>0) THEN
3747!
3748 DO j=jsta,jend
3749 DO i=ista,iend
3750 grid1(i,j) = el0(i,j)
3751 ENDDO
3752 ENDDO
3753 if(grib=="grib2") then
3754 cfld=cfld+1
3755 fld_info(cfld)%ifld=iavblfld(iget(147))
3756 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
3757 endif
3758 ENDIF
3759!
3760! IF REQUESTED, POST FREE ATMOSPHERE MASTER LENGTH SCALE
3761! AND/OR THE GRADIENT RICHARDSON NUMBER.
3762!
3763 IF ( (iget(111)>0) .OR. (iget(146)>0) ) THEN
3764!
3765! COMPUTE FREE ATMOSPHERE MASTER LENGTH SCALE.
3766!$omp parallel do private(i,j,l)
3767 DO l=1,lm
3768 DO j=jsta,jend
3769 DO i=ista,iend
3770 el(i,j,l) = d00
3771 ENDDO
3772 ENDDO
3773 ENDDO
3774
3775 IF(modelname == 'NCAR'.OR.modelname=='RSM'.OR. modelname == 'RAPR')THEN
3776! CALL MIXLEN(EL0,EL)
3777 ELSE IF(modelname == 'NMM')THEN
3778 DO l=1,lm
3779 DO j=jsta,jend
3780 DO i=ista,iend
3781 el(i,j,l)=el_pbl(i,j,l) !NOW EL COMES OUT OF WRF NMM
3782 ENDDO
3783 ENDDO
3784 ENDDO
3785 END IF
3786!
3787! COMPUTE GRADIENT RICHARDSON NUMBER IF REQUESTED.
3788!
3789 IF ( (iget(111)>0) ) CALL calrch(el,richno)
3790!
3791! LOOP OVER MDL LAYERS.
3792 DO 200 l = 1,lm
3793!
3794! POST MIXING LENGTH.
3795!
3796 IF (iget(146)>0) THEN
3797!
3798!
3799 IF (lvls(l,iget(146))>0) THEN
3800 ll=lm-l+1
3801!$omp parallel do private(i,j)
3802 DO j=jsta,jend
3803 DO i=ista,iend
3804 grid1(i,j) = el(i,j,ll)
3805 ENDDO
3806 ENDDO
3807 if(grib=="grib2") then
3808 cfld=cfld+1
3809 fld_info(cfld)%ifld=iavblfld(iget(146))
3810 fld_info(cfld)%lvl=lvlsxml(l,iget(146))
3811!$omp parallel do private(i,j,ii,jj)
3812 do j=1,jend-jsta+1
3813 jj = jsta+j-1
3814 do i=1,iend-ista+1
3815 ii = ista+i-1
3816 datapd(i,j,cfld) = grid1(ii,jj)
3817 enddo
3818 enddo
3819 endif
3820 ENDIF
3821 ENDIF
3822!
3823! POST GRADIENT RICHARDSON NUMBER.
3824!
3825 IF(l < lm)THEN
3826 IF (iget(111)>0) THEN
3827 IF (lvls(l,iget(111))>0) THEN
3828 ll=lm-l+1
3829!$omp parallel do private(i,j)
3830 DO j=jsta,jend
3831 DO i=ista,iend
3832 grid1(i,j) = richno(i,j,ll)
3833 ENDDO
3834 ENDDO
3835 if(grib=="grib2") then
3836 cfld=cfld+1
3837 fld_info(cfld)%ifld=iavblfld(iget(111))
3838 fld_info(cfld)%lvl=lvlsxml(l,iget(111))
3839!$omp parallel do private(i,j,ii,jj)
3840 do j=1,jend-jsta+1
3841 jj = jsta+j-1
3842 do i=1,iend-ista+1
3843 ii = ista+i-1
3844 datapd(i,j,cfld) = grid1(ii,jj)
3845 enddo
3846 enddo
3847 endif
3848 ENDIF
3849 ENDIF
3850 END IF
3851
3852 200 CONTINUE
3853!
3854!
3855 ENDIF
3856 ENDIF
3857!
3858! COMPUTE PBL HEIGHT BASED ON RICHARDSON NUMBER
3859!
3860 IF ( (iget(289)>0) .OR. (iget(389)>0) .OR. (iget(454)>0) &
3861 .OR. (iget(245)>0) .or. iget(464)>0 .or. iget(467)>0 &
3862 .or. iget(470)>0 .or. iget(476)>0) THEN
3863! should only compute pblri if pblh from model is not computed based on Ri
3864! post does not yet read pbl scheme used by model. Will do this soon
3865! For now, compute PBLRI for non GFS models.
3866 IF(modelname == 'GFS')THEN
3867 pblri=pblh
3868 ELSE
3869 CALL calpbl(pblri)
3870 END IF
3871 END IF
3872
3873 IF (iget(289) > 0) THEN
3874!$omp parallel do private(i,j)
3875 DO j=jsta,jend
3876 DO i=ista,iend
3877 grid1(i,j) = pblri(i,j)
3878! PBLH(I,J) = PBLRI(I,J)
3879 ENDDO
3880 ENDDO
3881 if(grib=="grib2") then
3882 cfld=cfld+1
3883 fld_info(cfld)%ifld=iavblfld(iget(289))
3884!$omp parallel do private(i,j,ii,jj)
3885 do j=1,jend-jsta+1
3886 jj = jsta+j-1
3887 do i=1,iend-ista+1
3888 ii = ista+i-1
3889 datapd(i,j,cfld) = grid1(ii,jj)
3890 enddo
3891 enddo
3892 endif
3893 ENDIF
3894! Pyle
3895! COMPUTE TRANSPORT WIND COMPONENTS (AVG WIND OVER MIXED LAYER)
3896!
3897!mp have model layer heights (ZMID, known) so we can average the winds (known) up to the PBLH (known)
3898
3899 IF ( (iget(389) > 0) .OR. (iget(454) > 0) ) THEN
3900!$omp parallel do private(i,j)
3901 DO j=jsta,jend
3902 DO i=ista,iend
3903 IF(pblri(i,j)<spval.and.zint(i,j,lm+1)<spval)THEN
3904 egrid3(i,j) = pblri(i,j) + zint(i,j,lm+1)
3905 ELSE
3906 egrid3(i,j) = spval
3907 ENDIF
3908 END DO
3909 END DO
3910! compute U and V separately because they are on different locations for B grid
3911 CALL h2u(egrid3(ista_2l:iend_2u,jsta_2l:jend_2u),egrid4)
3912!$omp parallel do private(i,j)
3913 DO j=jsta,jend
3914 DO i=ista,iend
3915 egrid1(i,j) = 0.0
3916 egrid2(i,j) = 0.0
3917 END DO
3918 END DO
3919 vert_loopu: DO l=lm,1,-1
3920 CALL h2u(zmid(ista_2l:iend_2u,jsta_2l:jend_2u,l), egrid5)
3921 CALL h2u(pint(ista_2l:iend_2u,jsta_2l:jend_2u,l+1),egrid6)
3922 CALL h2u(pint(ista_2l:iend_2u,jsta_2l:jend_2u,l), egrid7)
3923 hcount=0
3924 DO j=jsta,jend
3925 DO i=ista,iend
3926 if (egrid4(i,j)<spval.and.egrid5(i,j)<spval.and.&
3927 egrid6(i,j)<spval.and.egrid7(i,j)<spval.and.&
3928 uh(i,j,1)<spval)THEN
3929 if (egrid5(i,j) <= egrid4(i,j)) then
3930 hcount = hcount+1
3931 dp = egrid6(i,j) - egrid7(i,j)
3932 egrid1(i,j) = egrid1(i,j) + uh(i,j,l)*dp
3933 egrid2(i,j) = egrid2(i,j) + dp
3934! else
3935! exit vert_loopu
3936 endif
3937 endif
3938 end do
3939 end do
3940 if(hcount < 1 )exit vert_loopu
3941 ENDDO vert_loopu
3942!$omp parallel do private(i,j)
3943 DO j=jsta,jend
3944 DO i=ista,iend
3945 IF(egrid2(i,j) > 0.)THEN
3946 grid1(i,j) = egrid1(i,j)/egrid2(i,j)
3947 ELSE
3948 grid1(i,j) = u10(i,j) ! IF NO MIX LAYER, SPECIFY 10 M WIND, PER DIMEGO,
3949 END IF
3950 ustore(i,j) = grid1(i,j)
3951 END DO
3952 END DO
3953! compute v component now
3954 CALL h2v(egrid3(ista_2l:iend_2u,jsta_2l:jend_2u),egrid4)
3955!$omp parallel do private(i,j)
3956 DO j=jsta,jend
3957 DO i=ista,iend
3958 egrid1(i,j) = 0.
3959 egrid2(i,j) = 0.
3960 egrid5(i,j) = 0.
3961 egrid6(i,j) = 0.
3962 egrid7(i,j) = 0.
3963 END DO
3964 END DO
3965 vert_loopv: DO l=lm,1,-1
3966 CALL h2v(zmid(ista_2l:iend_2u,jsta_2l:jend_2u,l), egrid5)
3967 CALL h2v(pint(ista_2l:iend_2u,jsta_2l:jend_2u,l+1),egrid6)
3968 CALL h2v(pint(ista_2l:iend_2u,jsta_2l:jend_2u,l), egrid7)
3969 hcount=0
3970 DO j=jsta,jend
3971 DO i=ista,iend
3972 if (egrid4(i,j)<spval.and.egrid5(i,j)<spval.and.&
3973 egrid6(i,j)<spval.and.egrid7(i,j)<spval.and.&
3974 vh(i,j,1)<spval)THEN
3975 if (egrid5(i,j) <= egrid4(i,j)) then
3976 hcount=hcount+1
3977 dp = egrid6(i,j) - egrid7(i,j)
3978 egrid1(i,j) = egrid1(i,j) + vh(i,j,l)*dp
3979 egrid2(i,j) = egrid2(i,j) + dp
3980! else
3981! exit vert_loopu
3982 endif
3983 endif
3984 end do
3985 end do
3986 if(hcount<1)exit vert_loopv
3987 ENDDO vert_loopv
3988!$omp parallel do private(i,j)
3989 DO j=jsta,jend
3990 DO i=ista,iend
3991 IF(egrid2(i,j) > 0.)THEN
3992 grid2(i,j) = egrid1(i,j)/egrid2(i,j)
3993 ELSE
3994 grid2(i,j) = v10(i,j) ! IF NO MIX LAYER, SPECIFY 10 M WIND, PER DIMEGO,
3995 END IF
3996 vstore(i,j) = grid2(i,j)
3997 END DO
3998 END DO
3999
4000
4001 CALL u2h(ustore,egrid1)
4002 CALL v2h(vstore,egrid2)
4003!$omp parallel do private(i,j)
4004 DO j=jsta,jend
4005 DO i=ista,iend
4006
4007! EGRID1 is transport wind speed
4008 ! prevent floating overflow if either component is undefined
4009 IF (egrid1(i,j)<spval .and. egrid2(i,j)<spval) THEN
4010 egrid3(i,j) = sqrt((egrid1(i,j)*egrid1(i,j)+egrid2(i,j)*egrid2(i,j)))
4011 ELSe
4012 egrid3(i,j) = spval
4013 END IF
4014 ENDDO
4015 ENDDO
4016
4017 IF(iget(389) > 0)THEN
4018 if(grib=='grib2') then
4019 cfld=cfld+1
4020 fld_info(cfld)%ifld=iavblfld(iget(389))
4021!$omp parallel do private(i,j,ii,jj)
4022 do j=1,jend-jsta+1
4023 jj = jsta+j-1
4024 do i=1,iend-ista+1
4025 ii=ista+i-1
4026 datapd(i,j,cfld) = grid1(ii,jj)
4027 enddo
4028 enddo
4029 cfld=cfld+1
4030 fld_info(cfld)%ifld=iavblfld(iget(390))
4031!$omp parallel do private(i,j,ii,jj)
4032 do j=1,jend-jsta+1
4033 jj = jsta+j-1
4034 do i=1,iend-ista+1
4035 ii=ista+i-1
4036 datapd(i,j,cfld) = grid2(ii,jj)
4037 enddo
4038 enddo
4039 endif
4040 END IF
4041 ENDIF
4042!
4043! COMPUTE VENTILATION RATE (TRANSPORT WIND SPEED * MIXED LAYER HEIGHT)
4044!
4045! OK Mesonet has it in MKS units, so go with it. Ignore South Carolina fire
4046! comments about the winds being in MPH and the mixing height in feet.
4047
4048 IF ( (iget(454) > 0) ) THEN
4049
4050!$omp parallel do private(i,j)
4051 DO j=jsta,jend
4052 DO i=ista,iend
4053
4054 IF (pblri(i,j) /= spval .and. egrid3(i,j)/=spval) then
4055 grid1(i,j) = egrid3(i,j)*pblri(i,j)
4056 else
4057 grid1(i,j) = 0.
4058 ENDIF
4059
4060 ENDDO
4061 ENDDO
4062
4063 if(grib=='grib2') then
4064 cfld=cfld+1
4065 fld_info(cfld)%ifld=iavblfld(iget(454))
4066!$omp parallel do private(i,j,ii,jj)
4067 do j=1,jend-jsta+1
4068 jj = jsta+j-1
4069 do i=1,iend-ista+1
4070 ii = ista+i-1
4071 datapd(i,j,cfld) = grid1(ii,jj)
4072 enddo
4073 enddo
4074 endif
4075
4076
4077 ENDIF
4078!
4079! Calculate 10-m wind gust based on PBL height (as diagnosed from either Ri or theta-v)
4080 IF (iget(245)>0 .or. iget(464)>0 .or. iget(467)>0.or. iget(470)>0 .or. iget(476)>0) THEN
4081 IF (modelname=='RAPR') THEN
4082!tgs - 24may17 - smooth PBLHGUST
4083 if(maptype == 6) then
4084 if(grib=='grib2') then
4085 dxm = (dxval / 360.)*(erad*2.*pi)/1.d6 ! [mm]
4086 endif
4087 else
4088 dxm = dxval
4089 endif
4090 if(grib == 'grib2')then
4091 dxm=dxm/1000.0
4092 endif
4093! if(me==0)print *,'dxm=',dxm
4094 nsmooth = nint(5.*(13500./dxm))
4095 do j = jsta_2l, jend_2u
4096 do i = ista_2l, iend_2u
4097 grid1(i,j)=pblhgust(i,j)
4098 enddo
4099 enddo
4100 call allgetherv(grid1)
4101 do ks=1,nsmooth
4102 CALL smooth(grid1,sdummy,im,jm,0.5)
4103 end do
4104 do j = jsta_2l, jend_2u
4105 do i = ista_2l, iend_2u
4106 pblhgust(i,j)=grid1(i,j)
4107 enddo
4108 enddo
4109 ENDIF
4110
4111 DO j=jsta,jend
4112 DO i=ista,iend
4113 lpbl(i,j)=lm
4114
4115 if(zint(i,j,nint(lmh(i,j))+1) <spval) then
4116
4117 zsfc=zint(i,j,nint(lmh(i,j))+1)
4118 loopl:DO l=nint(lmh(i,j)),1,-1
4119 IF (modelname=='RAPR' .OR. modelname=='FV3R') THEN
4120 hgt=zmid(i,j,l)
4121 pblhold=pblhgust(i,j)
4122 IF(pblhold == spval) THEN
4123 lpbl(i,j) = lm
4124 EXIT loopl
4125 ENDIF
4126 ELSE
4127 hgt=zint(i,j,l)
4128 pblhold=pblri(i,j)
4129 ENDIF
4130 IF(hgt > pblhold+zsfc)THEN
4131 lpbl(i,j)=l+1
4132 IF(lpbl(i,j)>=lp1) lpbl(i,j) = lm
4133 EXIT loopl
4134 END IF
4135 ENDDO loopl
4136
4137 else
4138 lpbl(i,j) = lm
4139 endif
4140 if(lpbl(i,j)<1)print*,'zero lpbl',i,j,pblri(i,j),lpbl(i,j)
4141 ENDDO
4142 ENDDO
4143 IF (modelname=='RAPR' .OR. modelname=='FV3R') THEN
4144 CALL calgust(lpbl,pblhgust,gust)
4145 ELSE
4146 CALL calgust(lpbl,pblri,gust)
4147 END IF
4148 IF (iget(245)>0) THEN
4149!$omp parallel do private(i,j)
4150 DO j=jsta,jend
4151 DO i=ista,iend
4152! if(GUST(I,J) > 200. .and. gust(i,j)<spval) &
4153! print*,'big gust at ',i,j
4154 grid1(i,j) = gust(i,j)
4155 ENDDO
4156 ENDDO
4157 if(grib=='grib2') then
4158 cfld=cfld+1
4159 fld_info(cfld)%ifld=iavblfld(iget(245))
4160!$omp parallel do private(i,j,ii,jj)
4161 do j=1,jend-jsta+1
4162 jj = jsta+j-1
4163 do i=1,iend-ista+1
4164 ii = ista+i-1
4165 datapd(i,j,cfld) = grid1(ii,jj)
4166 enddo
4167 enddo
4168 endif
4169 ENDIF
4170 END IF
4171!
4172! COMPUTE PBL REGIME BASED ON WRF version of BULK RICHARDSON NUMBER
4173!
4174
4175 IF (iget(344)>0) THEN
4176 allocate(pblregime(ista_2l:iend_2u,jsta_2l:jend_2u))
4177 CALL calpblregime(pblregime)
4178!$omp parallel do private(i,j)
4179 DO j=jsta,jend
4180 DO i=ista,iend
4181 grid1(i,j) = pblregime(i,j)
4182 ENDDO
4183 ENDDO
4184 if(grib=="grib2") then
4185 cfld=cfld+1
4186 fld_info(cfld)%ifld=iavblfld(iget(344))
4187!$omp parallel do private(i,j,ii,jj)
4188 do j=1,jend-jsta+1
4189 jj = jsta+j-1
4190 do i=1,iend-ista+1
4191 ii = ista+i-1
4192 datapd(i,j,cfld) = grid1(ii,jj)
4193 enddo
4194 enddo
4195 endif
4196 deallocate(pblregime)
4197 ENDIF
4198!
4199! RADAR ECHO TOP (highest 18.3 dBZ level in each column)
4200!
4201 IF(iget(400)>0)THEN
4202 DO j=jsta,jend
4203 DO i=ista,iend
4204!Initialed as 'undetected'. Nov. 17, 2014, B. ZHOU:
4205!changed from SPVAL to -5000. to distinguish missing grids and undetected
4206! GRID1(I,J) = SPVAL
4207 grid1(i,j) = -5000. !undetected initially
4208 IF(imp_physics == 8.)then ! If Thompson MP
4209 DO l=1,nint(lmh(i,j))
4210 IF(ref_10cm(i,j,l) > 18.3) then
4211 grid1(i,j) = zmid(i,j,l)
4212 EXIT
4213 ENDIF
4214 ENDDO
4215 ELSE ! if other MP than Thompson
4216 DO l=1,nint(lmh(i,j))
4217 IF(dbz(i,j,l) > 18.3) then
4218 grid1(i,j) = zmid(i,j,l)
4219 EXIT
4220 END IF
4221 ENDDO
4222 END IF
4223 201 CONTINUE
4224! if(grid1(i,j)<0.)print*,'bad echo top', &
4225! + i,j,grid1(i,j),dbz(i,j,1:lm)
4226 ENDDO
4227 ENDDO
4228 if(grib=="grib2") then
4229 cfld=cfld+1
4230 fld_info(cfld)%ifld=iavblfld(iget(400))
4231!$omp parallel do private(i,j,ii,jj)
4232 do j=1,jend-jsta+1
4233 jj = jsta+j-1
4234 do i=1,iend-ista+1
4235 ii = ista+i-1
4236 datapd(i,j,cfld) = grid1(ii,jj)
4237 enddo
4238 enddo
4239 endif
4240 ENDIF
4241!
4242!
4243! COMPUTE NCAR GTG turbulence
4244 IF(gtg_on .and. (iget(464) > 0 .or. iget(467) > 0 .or. iget(470) > 0)) then
4245 i=(ista+iend)/2
4246 j=(jsta+jend)/2
4247! if(me == 0) print*,'sending input to GTG i,j,hgt,gust',i,j,ZINT(i,j,LP1),gust(i,j)
4248
4249 ! Use the existing 3D local arrays as cycled variables
4250 richno=spval
4251
4252 call gtg_algo(im,jm,lm,jsta,jend,jsta_2l,jend_2u,&
4253 uh(ista:iend,:,:),vh(ista:iend,:,:),wh(ista:iend,:,:),&
4254 zmid(ista:iend,:,:),pmid(ista:iend,:,:),t(ista:iend,:,:),&
4255 q(ista:iend,:,:),qqw(ista:iend,:,:),qqr(ista:iend,:,:),&
4256 qqs(ista:iend,:,:),qqg(ista:iend,:,:),qqi(ista:iend,:,:),&
4257 q2(ista:iend,:,:),&
4258 zint(ista:iend,:,lp1),pblh(ista:iend,:),sfcshx(ista:iend,:),&
4259 sfclhx(ista:iend,:),ustar(ista:iend,:),&
4260 z0(ista:iend,:),gdlat(ista:iend,:),gdlon(ista:iend,:),&
4261 dx(ista:iend,:),dy(ista:iend,:),u10(ista:iend,:),v10(ista:iend,:),&
4262 gust(ista:iend,:),avgprec(ista:iend,:),sm(ista:iend,:),sice(ista:iend,:),&
4263 catedr(ista:iend,:,:),mwt(ista:iend,:,:),cit(ista:iend,:,:),&
4264 richno(ista:iend,:,:),gtg(ista:iend,:,:),item)
4265
4266 i=iend
4267 j=jend ! 321,541
4268! print*,'GTG output: l,cat,mwt,gtg at',i,j
4269! do l=1,lm
4270! print*,l,catedr(i,j,l),mwt(i,j,l),gtg(i,j,l)
4271! end do
4272 ENDIF
4273
4274 IF (iget(470)>0) THEN
4275 Do l=1,lm
4276 IF (lvls(l,iget(470))>0) THEN
4277 ll=lm-l+1
4278 DO j=jsta,jend
4279 DO i=ista,iend
4280 grid1(i,j)=gtg(i,j,ll)
4281 ENDDO
4282 ENDDO
4283 if(grib=="grib2")then
4284 cfld=cfld+1
4285 fld_info(cfld)%ifld=iavblfld(iget(470))
4286 fld_info(cfld)%lvl=lvlsxml(l,iget(470))
4287!$omp parallel do private(i,j,ii,jj)
4288 do j=1,jend-jsta+1
4289 jj = jsta+j-1
4290 do i=1,iend-ista+1
4291 ii = ista+i-1
4292 datapd(i,j,cfld) = grid1(ii,jj)
4293 enddo
4294 enddo
4295 endif
4296
4297
4298 DO j=jsta,jend
4299 DO i=ista,iend
4300 grid1(i,j)=catedr(i,j,ll)
4301 ENDDO
4302 ENDDO
4303 if(grib=="grib2")then
4304 cfld=cfld+1
4305 fld_info(cfld)%ifld=iavblfld(iget(471))
4306 fld_info(cfld)%lvl=lvlsxml(l,iget(471))
4307!$omp parallel do private(i,j,ii,jj)
4308 do j=1,jend-jsta+1
4309 jj = jsta+j-1
4310 do i=1,iend-ista+1
4311 ii = ista+i-1
4312 datapd(i,j,cfld) = grid1(ii,jj)
4313 enddo
4314 enddo
4315 endif
4316
4317 DO j=jsta,jend
4318 DO i=ista,iend
4319 grid1(i,j)=mwt(i,j,ll)
4320 ENDDO
4321 ENDDO
4322 if(grib=="grib2")then
4323 cfld=cfld+1
4324 fld_info(cfld)%ifld=iavblfld(iget(472))
4325 fld_info(cfld)%lvl=lvlsxml(l,iget(472))
4326!$omp parallel do private(i,j,ii,jj)
4327 do j=1,jend-jsta+1
4328 jj = jsta+j-1
4329 do i=1,iend-ista+1
4330 ii = ista+i-1
4331 datapd(i,j,cfld) = grid1(ii,jj)
4332 enddo
4333 enddo
4334 endif
4335
4336 ENDIF
4337 end do
4338 end IF
4339
4340! COMPUTE NCAR FIP
4341 IF(iget(450)>0 .or. iget(480)>0 .or. iget(479)>0 .or. iget(481)>0)THEN
4342
4343! cape and cin
4344 itype = 1
4345 dpbnd = 300.e2
4346 dummy = 0.
4347 idummy = 0
4348 CALL calcape(itype,dpbnd,dummy,dummy,dummy,idummy,cape,cin, &
4349 dummy,dummy,dummy)
4350
4351 icing_gfip = spval
4352 icing_gfis = spval
4353 DO j=jsta,jend
4354 DO i=ista,iend
4355 if(debugprint .and. i==50 .and. j==jsta .and. me == 0) then
4356 print*,'sending input to FIP ',i,j,lm,gdlat(i,j),gdlon(i,j), &
4357 zint(i,j,lp1),cprate(i,j),prec(i,j),avgcprate(i,j),cape(i,j),cin(i,j)
4358 do l=1,lm
4359 if(debugprint)print*,'l,P,T,RH,CWM,QQW,QQI,QQR,QQS,QQG,OMEG',&
4360 l,pmid(i,j,l),t(i,j,l),rh3d(i,j,l),cwm(i,j,l), &
4361 q(i,j,l),qqw(i,j,l),qqi(i,j,l), &
4362 qqr(i,j,l),qqs(i,j,l),qqg(i,j,l),&
4363 rh3d(i,j,l),zmid(i,j,l),cwm(i,j,l),omga(i,j,l)
4364 end do
4365 end if
4366 CALL icing_algo(i,j,pmid(i,j,1:lm),t(i,j,1:lm),rh3d(i,j,1:lm) &
4367 ,zmid(i,j,1:lm),omga(i,j,1:lm),wh(i,j,1:lm) &
4368 ,q(i,j,1:lm),cwm(i,j,1:lm),qqw(i,j,1:lm),qqi(i,j,1:lm) &
4369 ,qqr(i,j,1:lm),qqs(i,j,1:lm),qqg(i,j,1:lm) &
4370 ,lm,gdlat(i,j),gdlon(i,j),zint(i,j,lp1) &
4371 ,prec(i,j),cprate(i,j),cape(i,j),cin(i,j) &
4372 ,icing_gfip(i,j,1:lm),icing_gfis(i,j,1:lm))
4373! if(gdlon(i,j)>=274. .and. gdlon(i,j)<=277. .and. gdlat(i,j)>=42. &
4374! .and. gdlat(i,j)<=45.)then
4375! print*,'sample FIP profile: l, H, T, RH, CWAT, VV, ICE POT at ' &
4376! , gdlon(i,j),gdlat(i,j)
4377! do l=1,lm
4378! print*,l,zmid(i,j,l),T(i,j,l),rh3d(i,j,l),cwm(i,j,l) &
4379! ,omga(i,j,l),icing_gfip(i,j,l),icing_gfis(i,j,l)
4380! end do
4381! end if
4382 ENDDO
4383 ENDDO
4384! Chuang: Change to output isobaric NCAR icing
4385! do l=1,lm
4386! if(LVLS(L,IGET(450))>0 .or. LVLS(L,IGET(480))>0)then
4387! do j=jsta,jend
4388! do i=ista,iend
4389! grid1(i,j)=icing_gfip(i,j,l)
4390! end do
4391! end do
4392! ID(1:25) = 0
4393! CALL GRIBIT(IGET(450),L,GRIDista,iend,JM)
4394! end if
4395! end do
4396 ENDIF
4397
4398 DEALLOCATE(el, richno, pblri)
4399 if (allocated(rh3d)) deallocate(rh3d)
4400!
4401! END OF ROUTINE.
4402!
4403 RETURN
4404 END
bound
subroutine bound(fld, fmin, fmax)
Clips data in passed array.
Definition BOUND.f:37
caldwp
subroutine caldwp(p1d, q1d, tdwp, t1d)
Computes dewpoint from P, T, and Q.
Definition CALDWP.f:28
calgust
subroutine calgust(lpbl, zpbl, gust)
This routine computes surface wind gust by mixing down momentum from the level at the height of the P...
Definition CALGUST.f:25
calmcvg
subroutine calmcvg(q1d, u1d, v1d, qcnvg)
Subroutine that computes moisture convergence.
Definition CALMCVG.f:44
calmict_new
subroutine calmict_new(p1d, t1d, q1d, c1d, fi1d, fr1d, fs1d, curefl, qw1, qi1, qr1, qs1, dbz1, dbzr1, dbzi1, dbzc1, nlice1, nrain1)
Subroutine that computes hydrometeors.
Definition CALMICT.f:67
calmict_old
subroutine calmict_old(p1d, t1d, q1d, c1d, fi1d, fr1d, fs1d, curefl, qw1, qi1, qr1, qs1, dbz1, dbzr1, dbzi1, dbzc1, nlice1, nrain1)
Definition CALMICT.f:333
calpblregime
subroutine calpblregime(pblregime)
Subroutine that determines the PBL regime.
Definition CALPBLREGIME.f:34
calpbl
subroutine calpbl(pblri)
Subroutine that computes PBL height based on bulk Richardson number.
Definition CALPBL.f:22
calpot
subroutine calpot(p1d, t1d, theta)
Subroutine that computes potential temperature.
Definition CALPOT.f:29
calrch
subroutine calrch(el, richno)
Subroutine that computes GRD RCH number.
Definition CALRCH.f:32
calstrm
subroutine calstrm(z1d, strm)
Subroutine that computes geo streamfunction.
Definition CALSTRM.f:32
mdlfld
subroutine mdlfld
SUBPROGRAM: MDLFLD SLP AND NATIVE LEVEL POSTING PRGRMMR: TREADON ORG: W/NP2 DATE: 92-12-21
Definition MDLFLD.f:100
ngmslp
subroutine ngmslp
SUBPROGRAM: NGMSLP NMC SEA LEVEL PRESSURE REDUCTION PRGRMMR: TREADON ORG: W/NP2 DATE: 93-02-02
Definition NGMSLP.f:91
smooth
subroutine smooth(field, hold, ix, iy, smth)
smooth() smooths a meteorological field using Shapiro smoother.
Definition SMOOTH.f:43