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