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