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