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