UPP  11.0.0
 All Data Structures Files Functions Variables Pages
MDL2P.f
Go to the documentation of this file.
1 
41  SUBROUTINE mdl2p(iostatusD3D)
42 
43 !
44 !
45  use vrbls4d, only: dust, smoke, fv3dust, coarsepm
46  use vrbls3d, only: pint, o3, pmid, t, q, uh, vh, wh, omga, q2, cwm, &
47  qqw, qqi, qqr, qqs, qqg, dbz, f_rimef, ttnd, cfr, &
48  rlwtt, rswtt, vdifftt, tcucn, tcucns, &
49  train, vdiffmois, dconvmois, sconvmois,nradtt, &
50  o3vdiff, o3prod, o3tndy, mwpv, unknown, vdiffzacce, &
51  zgdrag, cnvctvmmixing, vdiffmacce, mgdrag, &
52  cnvctummixing, ncnvctcfrac, cnvctumflx, cnvctdetmflx, &
53  cnvctzgdrag, cnvctmgdrag, zmid, zint, pmidv, &
54  cnvctdmflx, icing_gfip, icing_gfis,gtg,cat=>catedr,mwt
55  use vrbls2d, only: t500,t700,w_up_max,w_dn_max,w_mean,pslp,fis,z1000,z700,&
56  z500
57  use masks, only: lmh, sm
58  use physcons_post,only: con_fvirt, con_rog, con_eps, con_epsm1
59  use params_mod, only: h1m12, dbzmin, h1, pq0, a2, a3, a4, rhmin, g, &
60  rgamog, rd, d608, gi, erad, pi, small, h100, &
61  h99999, gamma
62  use ctlblk_mod, only: modelname, lp1, me, jsta, jend, lm, spval, spl, &
63  alsl, jend_m, smflag, grib, cfld, fld_info, datapd,&
64  td3d, ifhr, ifmin, im, jm, nbin_du, jsta_2l, &
65  jend_2u, lsm, d3d_on, ioform, nbin_sm, &
66  imp_physics, ista, iend, ista_m, iend_m, ista_2l, &
67  iend_2u, slrutah_on
68  use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml
69  use gridspec_mod, only: gridtype, maptype, dxval
70  use upp_physics, only: fpvsnew, calrh, calvor, calslr_roebber, calslr_uutah
71 
72 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 !
74  implicit none
75 !
76 ! INCLUDE MODEL DIMENSIONS. SET/DERIVE OTHER PARAMETERS.
77 ! GAMMA AND RGAMOG ARE USED IN THE EXTRAPOLATION OF VIRTUAL
78 ! TEMPERATURES BEYOND THE UPPER OF LOWER LIMITS OF DATA.
79 !
80 !
81  real,parameter:: gammam=-1*gamma,zshul=75.,tvshul=290.66
82 !
83 ! DECLARE VARIABLES.
84 !
85  real,PARAMETER :: capa=0.28589641,p1000=1000.e2
86  LOGICAL ioomg,ioall
87  real, dimension(im,jm) :: grid1, grid2
88  real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: fsl, tsl, qsl, osl, usl, vsl &
89  &, Q2SL, WSL, CFRSL, O3SL, TDSL &
90  &, EGRID1, EGRID2 &
91  &, FSL_OLD, USL_OLD, VSL_OLD &
92  &, OSL_OLD, OSL995 &
93  &, ICINGFSL, ICINGVSL &
94  &, GTGSL,CATSL,MWTSL
95  REAL, allocatable :: d3dsl(:,:,:), smokesl(:,:,:), fv3dustsl(:,:,:) &
96  &, COARSEPMSL(:,:,:)
97 !
98  integer,intent(in) :: iostatusd3d
99  INTEGER, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: nl1x, nl1xf
100  real, dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LSM) :: tprs, qprs, fprs
101  real, dimension(ISTA_2L:IEND_2U,JSTA_2L:JEND_2U,LSM) :: rhprs
102 !
103  INTEGER k, nsmooth
104 !
105 !--- Definition of the following 2D (horizontal) dummy variables
106 !
107 ! C1D - total condensate
108 ! QW1 - cloud water mixing ratio
109 ! QI1 - cloud ice mixing ratio
110 ! QR1 - rain mixing ratio
111 ! QS1 - snow mixing ratio
112 ! QG1 - graupel mixing ratio
113 ! DBZ1 - radar reflectivity
114 !
115  REAL, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: c1d, qw1, qi1, qr1, qs1, qg1, dbz1 &
116  , frime, rad, haines
117 
118  REAL sdummy(im,2)
119 
120 ! SAVE RH, U,V, for Icing, CAT, LLWS computation
121  REAL savrh(ista:iend,jsta:jend)
122 !jw
123  integer i,j,l,lp,ll,llmh,jjb,jje,ii,jj,li,ifincr,itd3d,istaa,imois,luhi,la
124  real fact,alpsl,psfc,qblo,pnl1,tblo,tvrl,tvrblo,fac,pslpij, &
125  alpth,ahf,pdv,ql,tvu,tvd,gammas,qsat,rhl,zl,tl,pl,es,part,dum1
126  logical log1
127  real dxm, tem, zero
128 !
129 !******************************************************************************
130 !
131 ! START MDL2P.
132 !
133  if(me==0) print*, 'MDL2P SMFLAG=',smflag
134 
135  if (modelname == 'GFS') then
136  zero = 0.0
137  else
138  zero = h1m12
139  endif
140  if (d3d_on) then
141  if (.not. allocated(d3dsl)) allocate(d3dsl(im,jm,27))
142 !$omp parallel do private(i,j,l)
143  do l=1,27
144  do j=1,jm
145  do i=1,im
146  d3dsl(i,j,l) = spval
147  enddo
148  enddo
149  enddo
150  endif
151  if (.not. allocated(smokesl)) allocate(smokesl(im,jm,nbin_sm))
152 !$omp parallel do private(i,j,l)
153  do l=1,nbin_sm
154  do j=1,jm
155  do i=1,im
156  smokesl(i,j,l) = spval
157  enddo
158  enddo
159  enddo
160  if (.not. allocated(fv3dustsl)) allocate(fv3dustsl(im,jm,nbin_sm))
161 !$omp parallel do private(i,j,l)
162  do l=1,nbin_sm
163  do j=1,jm
164  do i=1,im
165  fv3dustsl(i,j,l) = spval
166  enddo
167  enddo
168  enddo
169  if (.not. allocated(coarsepmsl)) allocate(coarsepmsl(im,jm,nbin_sm))
170 !$omp parallel do private(i,j,l)
171  do l=1,nbin_sm
172  do j=1,jm
173  do i=1,im
174  coarsepmsl(i,j,l) = spval
175  enddo
176  enddo
177  enddo
178 !
179 ! SET TOTAL NUMBER OF POINTS ON OUTPUT GRID.
180 !
181 !---------------------------------------------------------------
182 !
183 ! *** PART I ***
184 !
185 ! VERTICAL INTERPOLATION OF EVERYTHING ELSE. EXECUTE ONLY
186 ! IF THERE'S SOMETHING WE WANT.
187 !
188  IF((iget(012) > 0) .OR. (iget(013) > 0) .OR. &
189  (iget(014) > 0) .OR. (iget(015) > 0) .OR. &
190  (iget(016) > 0) .OR. (iget(017) > 0) .OR. &
191  (iget(018) > 0) .OR. (iget(019) > 0) .OR. &
192  (iget(020) > 0) .OR. (iget(030) > 0) .OR. &
193  (iget(021) > 0) .OR. (iget(022) > 0) .OR. &
194  (iget(023) > 0) .OR. (iget(085) > 0) .OR. &
195  (iget(086) > 0) .OR. (iget(284) > 0) .OR. &
196  (iget(153) > 0) .OR. (iget(166) > 0) .OR. &
197  (iget(183) > 0) .OR. (iget(184) > 0) .OR. &
198  (iget(198) > 0) .OR. (iget(251) > 0) .OR. &
199  (iget(257) > 0) .OR. (iget(258) > 0) .OR. &
200  (iget(294) > 0) .OR. (iget(268) > 0) .OR. &
201  (iget(331) > 0) .OR. (iget(326) > 0) .OR. &
202 ! add D3D fields
203  (iget(354) > 0) .OR. (iget(355) > 0) .OR. &
204  (iget(356) > 0) .OR. (iget(357) > 0) .OR. &
205  (iget(358) > 0) .OR. (iget(359) > 0) .OR. &
206  (iget(360) > 0) .OR. (iget(361) > 0) .OR. &
207  (iget(362) > 0) .OR. (iget(363) > 0) .OR. &
208  (iget(364) > 0) .OR. (iget(365) > 0) .OR. &
209  (iget(366) > 0) .OR. (iget(367) > 0) .OR. &
210  (iget(368) > 0) .OR. (iget(369) > 0) .OR. &
211  (iget(370) > 0) .OR. (iget(371) > 0) .OR. &
212  (iget(372) > 0) .OR. (iget(373) > 0) .OR. &
213  (iget(374) > 0) .OR. (iget(375) > 0) .OR. &
214  (iget(391) > 0) .OR. (iget(392) > 0) .OR. &
215  (iget(393) > 0) .OR. (iget(394) > 0) .OR. &
216  (iget(395) > 0) .OR. (iget(379) > 0) .OR. &
217 ! ADD DUST FIELDS
218  (iget(455) > 0) .OR. &
219 ! Add WAFS hazard fields: Icing and GTG turbulence
220  (iget(464) > 0) .OR. (iget(465) > 0) .OR. &
221  (iget(466) > 0) .OR. (iget(450) > 0) .OR. &
222  (iget(480) > 0) .OR. &
223 ! ADD SMOKE FIELDS
224  (iget(738) > 0) .OR. (iget(743) > 0) .OR. &
225  (modelname == 'RAPR') .OR.&
226 ! LIFTED INDEX needs 500 mb T
227  (iget(030)>0) .OR. (iget(031)>0) .OR. (iget(075)>0)) THEN
228 !
229 !---------------------------------------------------------------------
230 !***
231 !*** BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND, DO ALL
232 !*** INTERPOLATION ABOVE GROUND NOW.
233 !***
234 !
235 ! print*,'LSM= ',lsm
236 
237  if(gridtype == 'B' .or. gridtype == 'E') &
238  call exch(pint(ista_2l:iend_2u,jsta_2l:jend_2u,lp1))
239 
240  DO lp=1,lsm
241 
242 ! if(me == 0) print *,'in LP loop me=',me,'UH=',UH(1:10,JSTA,LP), &
243 ! 'JSTA_2L=',JSTA_2L,'JEND_2U=',JEND_2U,'JSTA=',JSTA,JEND, &
244 ! 'PMID(1,1,L)=',(PMID(1,1,LI),LI=1,LM),'SPL(LP)=',SPL(LP)
245 
246 ! if(me ==0) print *,'in mdl2p,LP loop o3=',maxval(o3(1:im,jsta:jend,lm))
247 !
248 !$omp parallel do private(i,j,l)
249  DO j=jsta_2l,jend_2u
250  DO i=ista_2l,iend_2u
251  tsl(i,j) = spval
252  qsl(i,j) = spval
253  fsl(i,j) = spval
254  osl(i,j) = spval
255  usl(i,j) = spval
256  vsl(i,j) = spval
257 ! dong initialize wsl
258  wsl(i,j) = spval
259  q2sl(i,j) = spval
260  c1d(i,j) = spval ! Total condensate
261  qw1(i,j) = spval ! Cloud water
262  qi1(i,j) = spval ! Cloud ice
263  qr1(i,j) = spval ! Rain
264  qs1(i,j) = spval ! Snow (precip ice)
265  qg1(i,j) = spval ! Graupel
266  dbz1(i,j) = spval
267  frime(i,j) = spval
268  rad(i,j) = spval
269  o3sl(i,j) = spval
270  cfrsl(i,j) = spval
271  icingfsl(i,j) = spval
272  icingvsl(i,j) = spval
273  gtgsl(i,j) = spval
274  catsl(i,j) = spval
275  mwtsl(i,j) = spval
276 !
277 !*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER JUST BELOW
278 !*** THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
279 !
280  nl1x(i,j) = lp1
281  DO l=2,lm
282  IF(nl1x(i,j) == lp1 .AND. pmid(i,j,l) > spl(lp)) THEN
283  nl1x(i,j) = l
284  ENDIF
285  ENDDO
286 !
287 ! IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
288 ! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
289 ! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
290 ! WILL EXTRAPOLATE TO THAT POINT
291 !
292  IF(nl1x(i,j) == lp1 .AND. pint(i,j,lp1) > spl(lp)) THEN
293  nl1x(i,j) = lm
294  ENDIF
295 
296  nl1xf(i,j) = lp1 + 1
297  DO l=2,lp1
298  IF(nl1xf(i,j) == (lp1+1) .AND. pint(i,j,l) > spl(lp)) THEN
299  nl1xf(i,j) = l
300  ENDIF
301  ENDDO
302 ! if(NL1X(I,J) == LMP1)print*,'Debug: NL1X=LMP1 AT ' 1 ,i,j,lp
303  ENDDO
304  ENDDO ! end of j loop
305 
306 !
307 !mptest IF(NHOLD == 0)GO TO 310
308 !
309 !!$omp parallel do
310 !!$omp& private(nn,i,j,ll,fact,qsat,rhl)
311 !hc DO 220 NN=1,NHOLD
312 !hc I=IHOLD(NN)
313 !hc J=JHOLD(NN)
314 ! DO 220 J=JSTA,JEND
315 
316  ii = (ista+iend)/2
317  jj = (jsta+jend)/2
318 
319 !$omp parallel do private(i,j,k,l,ll,llmh,la,tvd,tvu,fact,fac,ahf,rhl,tl,pl,ql,zl,es,qsat,part,tvrl,tvrblo,tblo,qblo,gammas,pnl1)
320  DO j=jsta,jend
321  DO i=ista,iend
322 !---------------------------------------------------------------------
323 !*** VERTICAL INTERPOLATION OF GEOPOTENTIAL, TEMPERATURE, SPECIFIC
324 !*** HUMIDITY, CLOUD WATER/ICE, OMEGA, WINDS, AND TKE.
325 !---------------------------------------------------------------------
326 !
327  ll = nl1x(i,j)
328  llmh = nint(lmh(i,j))
329 
330 !HC IF(NL1X(I,J)<=LM)THEN
331 
332  IF(spl(lp) < pint(i,j,2)) THEN ! Above second interface
333  IF(t(i,j,1) < spval) tsl(i,j) = t(i,j,1)
334  IF(q(i,j,1) < spval) qsl(i,j) = q(i,j,1)
335 
336  IF(gridtype == 'A')THEN
337  IF(uh(i,j,1) < spval) usl(i,j) = uh(i,j,1)
338  IF(vh(i,j,1) < spval) vsl(i,j) = vh(i,j,1)
339  END IF
340 
341 ! if ( J == JSTA.and. I == 1.and.me == 0) &
342 ! print *,'1 USL=',USL(I,J),UH(I,J,1)
343 
344  IF(wh(i,j,1) < spval) wsl(i,j) = wh(i,j,1)
345  IF(omga(i,j,1) < spval) osl(i,j) = omga(i,j,1)
346  IF(q2(i,j,1) < spval) q2sl(i,j) = q2(i,j,1)
347  IF(cwm(i,j,1) < spval) c1d(i,j) = cwm(i,j,1)
348  c1d(i,j) = max(c1d(i,j),zero) ! Total condensate
349  IF(qqw(i,j,1) < spval) qw1(i,j) = qqw(i,j,1)
350  qw1(i,j) = max(qw1(i,j),zero) ! Cloud water
351  IF(qqi(i,j,1) < spval) qi1(i,j) = qqi(i,j,1)
352  qi1(i,j) = max(qi1(i,j),zero) ! Cloud ice
353  IF(qqr(i,j,1) < spval) qr1(i,j) = qqr(i,j,1)
354  qr1(i,j) = max(qr1(i,j),zero) ! Rain
355  IF(qqs(i,j,1) < spval) qs1(i,j) = qqs(i,j,1)
356  qs1(i,j) = max(qs1(i,j),zero) ! Snow (precip ice)
357  IF(qqg(i,j,1) < spval) qg1(i,j) = qqg(i,j,1)
358  qg1(i,j) = max(qg1(i,j),zero) ! Graupel (precip ice)
359  IF(dbz(i,j,1) < spval) dbz1(i,j) = dbz(i,j,1)
360  dbz1(i,j) = max(dbz1(i,j),dbzmin)
361  IF(f_rimef(i,j,1) < spval) frime(i,j) = f_rimef(i,j,1)
362  frime(i,j) = max(frime(i,j),h1)
363  IF(ttnd(i,j,1) < spval) rad(i,j) = ttnd(i,j,1)
364  IF(o3(i,j,1) < spval) o3sl(i,j) = o3(i,j,1)
365  IF(cfr(i,j,1) < spval) cfrsl(i,j) = cfr(i,j,1)
366 !GFIP
367  IF(icing_gfip(i,j,1) < spval) icingfsl(i,j) = icing_gfip(i,j,1)
368  IF(icing_gfis(i,j,1) < spval) icingvsl(i,j) = icing_gfis(i,j,1)
369 !GTG
370  IF(gtg(i,j,1) < spval) gtgsl(i,j) = gtg(i,j,1)
371  IF(cat(i,j,1) < spval) catsl(i,j) = cat(i,j,1)
372  IF(mwt(i,j,1) < spval) mwtsl(i,j) = mwt(i,j,1)
373  DO k = 1, nbin_sm
374  IF(smoke(i,j,1,k) < spval) smokesl(i,j,k)=smoke(i,j,1,k)
375  IF(fv3dust(i,j,1,k) < spval) fv3dustsl(i,j,k)=fv3dust(i,j,1,k)
376  IF(coarsepm(i,j,1,k) < spval) coarsepmsl(i,j,k)=coarsepm(i,j,1,k)
377  ENDDO
378 
379 ! only interpolate GFS d3d fields when reqested
380 ! if(iostatusD3D ==0 .and. d3d_on)then
381  if (d3d_on) then
382  IF((iget(354) > 0) .OR. (iget(355) > 0) .OR. &
383  (iget(356) > 0) .OR. (iget(357) > 0) .OR. &
384  (iget(358) > 0) .OR. (iget(359) > 0) .OR. &
385  (iget(360) > 0) .OR. (iget(361) > 0) .OR. &
386  (iget(362) > 0) .OR. (iget(363) > 0) .OR. &
387  (iget(364) > 0) .OR. (iget(365) > 0) .OR. &
388  (iget(366) > 0) .OR. (iget(367) > 0) .OR. &
389  (iget(368) > 0) .OR. (iget(369) > 0) .OR. &
390  (iget(370) > 0) .OR. (iget(371) > 0) .OR. &
391  (iget(372) > 0) .OR. (iget(373) > 0) .OR. &
392  (iget(374) > 0) .OR. (iget(375) > 0) .OR. &
393  (iget(391) > 0) .OR. (iget(392) > 0) .OR. &
394  (iget(393) > 0) .OR. (iget(394) > 0) .OR. &
395  (iget(395) > 0) .OR. (iget(379) > 0)) THEN
396  d3dsl(i,j,1) = rlwtt(i,j,1)
397  d3dsl(i,j,2) = rswtt(i,j,1)
398  d3dsl(i,j,3) = vdifftt(i,j,1)
399  d3dsl(i,j,4) = tcucn(i,j,1)
400  d3dsl(i,j,5) = tcucns(i,j,1)
401  d3dsl(i,j,6) = train(i,j,1)
402  d3dsl(i,j,7) = vdiffmois(i,j,1)
403  d3dsl(i,j,8) = dconvmois(i,j,1)
404  d3dsl(i,j,9) = sconvmois(i,j,1)
405  d3dsl(i,j,10) = nradtt(i,j,1)
406  d3dsl(i,j,11) = o3vdiff(i,j,1)
407  d3dsl(i,j,12) = o3prod(i,j,1)
408  d3dsl(i,j,13) = o3tndy(i,j,1)
409  d3dsl(i,j,14) = mwpv(i,j,1)
410  d3dsl(i,j,15) = unknown(i,j,1)
411  d3dsl(i,j,16) = vdiffzacce(i,j,1)
412  d3dsl(i,j,17) = zgdrag(i,j,1)
413  d3dsl(i,j,18) = cnvctummixing(i,j,1)
414  d3dsl(i,j,19) = vdiffmacce(i,j,1)
415  d3dsl(i,j,20) = mgdrag(i,j,1)
416  d3dsl(i,j,21) = cnvctvmmixing(i,j,1)
417  d3dsl(i,j,22) = ncnvctcfrac(i,j,1)
418  d3dsl(i,j,23) = cnvctumflx(i,j,1)
419  d3dsl(i,j,24) = cnvctdmflx(i,j,1)
420  d3dsl(i,j,25) = cnvctdetmflx(i,j,1)
421  d3dsl(i,j,26) = cnvctzgdrag(i,j,1)
422  d3dsl(i,j,27) = cnvctmgdrag(i,j,1)
423  end if
424  end if
425 
426  ELSE IF(ll <= llmh)THEN
427 !
428 !---------------------------------------------------------------------
429 ! INTERPOLATE LINEARLY IN LOG(P)
430 !*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
431 !*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
432 !*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
433 !---------------------------------------------------------------------
434 !
435 !KRF: Need ncar and nmm wrf core checks as well?
436  IF (modelname == 'RAPR' .OR. modelname == 'NCAR' .OR. modelname == 'NMM') THEN
437  fact = (alsl(lp)-log(pmid(i,j,ll)))/ &
438  max(1.e-6,(log(pmid(i,j,ll))-log(pmid(i,j,ll-1))))
439  fact = max(-10.0,min(fact, 10.0))
440  ELSEIF (modelname == 'GFS' .OR. modelname == 'FV3R') THEN
441  fact = (alsl(lp)-log(pmid(i,j,ll)))/ &
442  max(1.e-6,(log(pmid(i,j,ll))-log(pmid(i,j,ll-1))))
443  fact = max(-10.0,min(fact, 10.0))
444  IF ( abs(pmid(i,j,ll)-pmid(i,j,ll-1)) < 0.5 ) THEN
445  fact = -1.
446  ENDIF
447  ELSE
448  fact = (alsl(lp)-log(pmid(i,j,ll)))/ &
449  (log(pmid(i,j,ll))-log(pmid(i,j,ll-1)))
450  ENDIF
451  IF(t(i,j,ll) < spval .AND. t(i,j,ll-1) < spval) &
452  tsl(i,j) = t(i,j,ll)+(t(i,j,ll)-t(i,j,ll-1))*fact
453  IF(q(i,j,ll) < spval .AND. q(i,j,ll-1) < spval) &
454  qsl(i,j) = q(i,j,ll)+(q(i,j,ll)-q(i,j,ll-1))*fact
455 
456  IF(gridtype=='A')THEN
457  IF(uh(i,j,ll) < spval .AND. uh(i,j,ll-1) < spval) &
458  usl(i,j) = uh(i,j,ll)+(uh(i,j,ll)-uh(i,j,ll-1))*fact
459  IF(vh(i,j,ll) < spval .AND. vh(i,j,ll-1) < spval) &
460  vsl(i,j) = vh(i,j,ll)+(vh(i,j,ll)-vh(i,j,ll-1))*fact
461  END IF
462 
463  IF(wh(i,j,ll) < spval .AND. wh(i,j,ll-1) < spval) &
464  wsl(i,j) = wh(i,j,ll)+(wh(i,j,ll)-wh(i,j,ll-1))*fact
465  IF(omga(i,j,ll) < spval .AND. omga(i,j,ll-1) < spval) &
466  osl(i,j) = omga(i,j,ll)+(omga(i,j,ll)-omga(i,j,ll-1))*fact
467  IF(q2(i,j,ll) < spval .AND. q2(i,j,ll-1) < spval) &
468  q2sl(i,j) = q2(i,j,ll)+(q2(i,j,ll)-q2(i,j,ll-1))*fact
469 
470 ! IF(ZMID(I,J,LL) < SPVAL .AND. ZMID(I,J,LL-1) < SPVAL) &
471 ! & FSL(I,J) = ZMID(I,J,LL)+(ZMID(I,J,LL)-ZMID(I,J,LL-1))*FACT
472 ! FSL(I,J) = FSL(I,J)*G
473 
474  if (modelname == 'GFS') then
475  es = min(fpvsnew(tsl(i,j)), spl(lp))
476  qsat = con_eps*es/(spl(lp)+con_epsm1*es)
477  else
478  qsat = pq0/spl(lp)*exp(a2*(tsl(i,j)-a3)/(tsl(i,j)-a4))
479  endif
480 !
481  rhl = max(rhmin, min(1.0, qsl(i,j)/qsat))
482  qsl(i,j) = rhl*qsat
483 
484 ! if(tsl(i,j) > 330. .or. tsl(i,j) < 100.)print*, &
485 ! 'bad isobaric T Q',i,j,lp,tsl(i,j),qsl(i,j) &
486 ! ,T(I,J,LL),T(I,J,LL-1),Q(I,J,LL),Q(I,J,LL-1)
487 
488  IF(q2sl(i,j) < 0.0) q2sl(i,j) = 0.0
489 !
490 !HC ADD FERRIER'S HYDROMETEOR
491  IF(cwm(i,j,ll) < spval .AND. cwm(i,j,ll-1) < spval) &
492  c1d(i,j) = cwm(i,j,ll) + (cwm(i,j,ll)-cwm(i,j,ll-1))*fact
493  c1d(i,j) = max(c1d(i,j),zero) ! Total condensate
494 
495  IF(qqw(i,j,ll) < spval .AND. qqw(i,j,ll-1) < spval) &
496  qw1(i,j) = qqw(i,j,ll) + (qqw(i,j,ll)-qqw(i,j,ll-1))*fact
497  qw1(i,j) = max(qw1(i,j),zero) ! Cloud water
498 
499  IF(qqi(i,j,ll) < spval .AND. qqi(i,j,ll-1) < spval) &
500  qi1(i,j) = qqi(i,j,ll) + (qqi(i,j,ll)-qqi(i,j,ll-1))*fact
501  qi1(i,j) = max(qi1(i,j),zero) ! Cloud ice
502 
503  IF(qqr(i,j,ll) < spval .AND. qqr(i,j,ll-1) < spval) &
504  qr1(i,j) = qqr(i,j,ll) + (qqr(i,j,ll)-qqr(i,j,ll-1))*fact
505  qr1(i,j) = max(qr1(i,j),zero) ! Rain
506 
507  IF(qqs(i,j,ll) < spval .AND. qqs(i,j,ll-1) < spval) &
508  qs1(i,j) = qqs(i,j,ll) + (qqs(i,j,ll)-qqs(i,j,ll-1))*fact
509  qs1(i,j) = max(qs1(i,j),zero) ! Snow (precip ice)
510 
511  IF(qqg(i,j,ll) < spval .AND. qqg(i,j,ll-1) < spval) &
512  qg1(i,j) = qqg(i,j,ll) + (qqg(i,j,ll)-qqg(i,j,ll-1))*fact
513  qg1(i,j) = max(qg1(i,j),zero) ! GRAUPEL (precip ice)
514 
515  IF(dbz(i,j,ll) < spval .AND. dbz(i,j,ll-1) < spval) &
516  dbz1(i,j) = dbz(i,j,ll) + (dbz(i,j,ll)-dbz(i,j,ll-1))*fact
517  dbz1(i,j) = max(dbz1(i,j),dbzmin)
518 
519  IF(f_rimef(i,j,ll) < spval .AND. f_rimef(i,j,ll-1) < spval) &
520  frime(i,j) = f_rimef(i,j,ll) + (f_rimef(i,j,ll) - f_rimef(i,j,ll-1))*fact
521  frime(i,j)=max(frime(i,j),h1)
522 
523  IF(ttnd(i,j,ll) < spval .AND. ttnd(i,j,ll-1) < spval) &
524  rad(i,j) = ttnd(i,j,ll) + (ttnd(i,j,ll)-ttnd(i,j,ll-1))*fact
525 
526  IF(o3(i,j,ll) < spval .AND. o3(i,j,ll-1) < spval) &
527  o3sl(i,j) = o3(i,j,ll) + (o3(i,j,ll)-o3(i,j,ll-1))*fact
528 
529  IF(cfr(i,j,ll) < spval .AND. cfr(i,j,ll-1) < spval) &
530  cfrsl(i,j) = cfr(i,j,ll) + (cfr(i,j,ll)-cfr(i,j,ll-1))*fact
531 !GFIP
532  IF(icing_gfip(i,j,ll) < spval .AND. icing_gfip(i,j,ll-1) < spval) &
533  icingfsl(i,j) = icing_gfip(i,j,ll) + (icing_gfip(i,j,ll)-icing_gfip(i,j,ll-1))*fact
534  icingfsl(i,j) = max(0.0, icingfsl(i,j))
535  icingfsl(i,j) = min(1.0, icingfsl(i,j))
536  IF(icing_gfis(i,j,ll) < spval .AND. icing_gfis(i,j,ll-1) < spval) &
537  icingvsl(i,j) = icing_gfis(i,j,ll) + (icing_gfis(i,j,ll)-icing_gfis(i,j,ll-1))*fact
538 ! Icing severity categories
539 ! 0 = none (0, 0.08)
540 ! 1 = trace [0.08, 0.21]
541 ! 2 = light (0.21, 0.37]
542 ! 3 = moderate (0.37, 0.67]
543 ! 4 = severe (0.67, 1]
544 ! https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-228.shtml
545  if (icingvsl(i,j) < 0.08) then
546  icingvsl(i,j) = 0.0
547  elseif (icingvsl(i,j) <= 0.21) then
548  icingvsl(i,j) = 1.
549  else if(icingvsl(i,j) <= 0.37) then
550  icingvsl(i,j) = 2.0
551  else if(icingvsl(i,j) <= 0.67) then
552  icingvsl(i,j) = 3.0
553  else
554  icingvsl(i,j) = 4.0
555  endif
556  if(icingfsl(i,j)< 0.001) icingvsl(i,j) = 0.
557 ! GTG
558  IF(gtg(i,j,ll) < spval .AND. gtg(i,j,ll-1) < spval) THEN
559  gtgsl(i,j) = gtg(i,j,ll) + (gtg(i,j,ll)-gtg(i,j,ll-1))*fact
560  gtgsl(i,j) = max(0.0, gtgsl(i,j))
561  gtgsl(i,j) = min(1.0, gtgsl(i,j))
562  ENDIF
563  IF(cat(i,j,ll) < spval .AND. cat(i,j,ll-1) < spval) THEN
564  catsl(i,j) = cat(i,j,ll) + (cat(i,j,ll)-cat(i,j,ll-1))*fact
565  catsl(i,j) = max(0.0, catsl(i,j))
566  catsl(i,j) = min(1.0, catsl(i,j))
567  ENDIF
568  IF(mwt(i,j,ll) < spval .AND. mwt(i,j,ll-1) < spval) THEN
569  mwtsl(i,j) = mwt(i,j,ll) + (mwt(i,j,ll)-mwt(i,j,ll-1))*fact
570  mwtsl(i,j) = max(0.0, mwtsl(i,j))
571  mwtsl(i,j) = min(1.0, mwtsl(i,j))
572  ENDIF
573  DO k = 1, nbin_sm
574  IF(smoke(i,j,ll,k) < spval .AND. smoke(i,j,ll-1,k) < spval) &
575  smokesl(i,j,k)=smoke(i,j,ll,k)+(smoke(i,j,ll,k)-smoke(i,j,ll-1,k))*fact
576  IF(fv3dust(i,j,ll,k) < spval .AND. fv3dust(i,j,ll-1,k) < spval) &
577  fv3dustsl(i,j,k)=fv3dust(i,j,ll,k)+(fv3dust(i,j,ll,k)-fv3dust(i,j,ll-1,k))*fact
578  IF(coarsepm(i,j,ll,k) < spval .AND. coarsepm(i,j,ll-1,k) < spval) &
579  coarsepmsl(i,j,k)=coarsepm(i,j,ll,k)+(coarsepm(i,j,ll,k)-coarsepm(i,j,ll-1,k))*fact
580  ENDDO
581 
582 ! only interpolate GFS d3d fields when == ested
583 ! if(iostatusD3D==0)then
584  if (d3d_on) then
585  IF((iget(354) > 0) .OR. (iget(355) > 0) .OR. &
586  (iget(356) > 0) .OR. (iget(357) > 0) .OR. &
587  (iget(358) > 0) .OR. (iget(359) > 0) .OR. &
588  (iget(360) > 0) .OR. (iget(361) > 0) .OR. &
589  (iget(362) > 0) .OR. (iget(363) > 0) .OR. &
590  (iget(364) > 0) .OR. (iget(365) > 0) .OR. &
591  (iget(366) > 0) .OR. (iget(367) > 0) .OR. &
592  (iget(368) > 0) .OR. (iget(369) > 0) .OR. &
593  (iget(370) > 0) .OR. (iget(371) > 0) .OR. &
594  (iget(372) > 0) .OR. (iget(373) > 0) .OR. &
595  (iget(374) > 0) .OR. (iget(375) > 0) .OR. &
596  (iget(391) > 0) .OR. (iget(392) > 0) .OR. &
597  (iget(393) > 0) .OR. (iget(394) > 0) .OR. &
598  (iget(395) > 0) .OR. (iget(379) > 0))THEN
599  d3dsl(i,j,1) = rlwtt(i,j,ll)+(rlwtt(i,j,ll) &
600  - rlwtt(i,j,ll-1))*fact
601  d3dsl(i,j,2) = rswtt(i,j,ll)+(rswtt(i,j,ll) &
602  - rswtt(i,j,ll-1))*fact
603  d3dsl(i,j,3) = vdifftt(i,j,ll)+(vdifftt(i,j,ll) &
604  - vdifftt(i,j,ll-1))*fact
605  d3dsl(i,j,4) = tcucn(i,j,ll)+(tcucn(i,j,ll) &
606  - tcucn(i,j,ll-1))*fact
607  d3dsl(i,j,5) = tcucns(i,j,ll)+(tcucns(i,j,ll) &
608  - tcucns(i,j,ll-1))*fact
609  d3dsl(i,j,6) = train(i,j,ll)+(train(i,j,ll) &
610  - train(i,j,ll-1))*fact
611  d3dsl(i,j,7) = vdiffmois(i,j,ll)+ &
612  (vdiffmois(i,j,ll)-vdiffmois(i,j,ll-1))*fact
613  d3dsl(i,j,8) = dconvmois(i,j,ll)+ &
614  (dconvmois(i,j,ll)-dconvmois(i,j,ll-1))*fact
615  d3dsl(i,j,9) = sconvmois(i,j,ll)+ &
616  (sconvmois(i,j,ll)-sconvmois(i,j,ll-1))*fact
617  d3dsl(i,j,10) = nradtt(i,j,ll)+ &
618  (nradtt(i,j,ll)-nradtt(i,j,ll-1))*fact
619  d3dsl(i,j,11) = o3vdiff(i,j,ll)+ &
620  (o3vdiff(i,j,ll)-o3vdiff(i,j,ll-1))*fact
621  d3dsl(i,j,12) = o3prod(i,j,ll)+ &
622  (o3prod(i,j,ll)-o3prod(i,j,ll-1))*fact
623  d3dsl(i,j,13) = o3tndy(i,j,ll)+ &
624  (o3tndy(i,j,ll)-o3tndy(i,j,ll-1))*fact
625  d3dsl(i,j,14) = mwpv(i,j,ll)+ &
626  (mwpv(i,j,ll)-mwpv(i,j,ll-1))*fact
627  d3dsl(i,j,15) = unknown(i,j,ll)+ &
628  (unknown(i,j,ll)-unknown(i,j,ll-1))*fact
629  d3dsl(i,j,16) = vdiffzacce(i,j,ll)+ &
630  (vdiffzacce(i,j,ll)-vdiffzacce(i,j,ll-1))*fact
631  d3dsl(i,j,17) = zgdrag(i,j,ll)+ &
632  (zgdrag(i,j,ll)-zgdrag(i,j,ll-1))*fact
633  d3dsl(i,j,18) = cnvctummixing(i,j,ll)+ &
634  (cnvctummixing(i,j,ll)-cnvctummixing(i,j,ll-1))*fact
635  d3dsl(i,j,19) = vdiffmacce(i,j,ll)+ &
636  (vdiffmacce(i,j,ll)-vdiffmacce(i,j,ll-1))*fact
637  d3dsl(i,j,20) = mgdrag(i,j,ll)+ &
638  (mgdrag(i,j,ll)-mgdrag(i,j,ll-1))*fact
639  d3dsl(i,j,21) = cnvctvmmixing(i,j,ll)+ &
640  (cnvctvmmixing(i,j,ll)-cnvctvmmixing(i,j,ll-1))*fact
641  d3dsl(i,j,22) = ncnvctcfrac(i,j,ll)+ &
642  (ncnvctcfrac(i,j,ll)-ncnvctcfrac(i,j,ll-1))*fact
643  d3dsl(i,j,23) = cnvctumflx(i,j,ll)+ &
644  (cnvctumflx(i,j,ll)-cnvctumflx(i,j,ll-1))*fact
645  d3dsl(i,j,24) = cnvctdmflx(i,j,ll)+ &
646  (cnvctdmflx(i,j,ll)-cnvctdmflx(i,j,ll-1))*fact
647  d3dsl(i,j,25) = cnvctdetmflx(i,j,ll)+ &
648  (cnvctdetmflx(i,j,ll)-cnvctdetmflx(i,j,ll-1))*fact
649  d3dsl(i,j,26) = cnvctzgdrag(i,j,ll)+ &
650  (cnvctzgdrag(i,j,ll)-cnvctzgdrag(i,j,ll-1))*fact
651  d3dsl(i,j,27) = cnvctmgdrag(i,j,ll)+ &
652  (cnvctmgdrag(i,j,ll)-cnvctmgdrag(i,j,ll-1))*fact
653  end if
654  end if ! if d3d_on test
655 
656 ! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE
657 ! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
658 ! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
659 ! GOUND
660  ELSE ! underground
661  IF(modelname == 'GFS')THEN ! GFS deduce T and H using Shuell
662  tvu = t(i,j,lm) * (1.+con_fvirt*q(i,j,lm))
663  if(zmid(i,j,lm) > zshul) then
664  tvd = tvu + gamma*zmid(i,j,lm)
665  if(tvd > tvshul) then
666  if(tvu > tvshul) then
667  tvd = tvshul - 5.e-3*(tvu-tvshul)*(tvu-tvshul)
668  else
669  tvd = tvshul
670  endif
671  endif
672  gammas = (tvu-tvd)/zmid(i,j,lm)
673  else
674  gammas = 0.
675  endif
676  part = con_rog*(alsl(lp)-log(pmid(i,j,lm)))
677  fsl(i,j) = zmid(i,j,lm) - tvu*part/(1.+0.5*gammas*part)
678 ! tp(k) = t(1)+gammas*(hp(k)-h(1))
679  tsl(i,j) = t(i,j,lm) - gamma*(fsl(i,j)-zmid(i,j,lm))
680  fsl(i,j) = fsl(i,j)*g ! just use NAM G for now since FSL will be divided by GI later
681 !
682 ! Compute RH at lowest model layer because Iredell and Chuang decided to compute
683 ! underground GFS Q to maintain RH
684  es = min(fpvsnew(t(i,j,lm)), pmid(i,j,lm))
685  qsat = con_eps*es/(pmid(i,j,lm)+con_epsm1*es)
686  rhl = q(i,j,lm)/qsat
687 ! compute saturation water vapor at isobaric level
688  es = min(fpvsnew(tsl(i,j)), spl(lp))
689  qsat = con_eps*es/(spl(lp)+con_epsm1*es)
690 ! Q at isobaric level is computed by maintaining constant RH
691  qsl(i,j) = rhl*qsat
692 
693  ELSE
694  pl = pint(i,j,lm-1)
695  zl = zint(i,j,lm-1)
696  tl = 0.5*(t(i,j,lm-2)+t(i,j,lm-1))
697  ql = 0.5*(q(i,j,lm-2)+q(i,j,lm-1))
698 ! TMT0=TL-A0
699 ! TMT15=MIN(TMT0,-15.)
700 ! AI=0.008855
701 ! BI=1.
702 ! IF(TMT0 < -20.)THEN
703 ! AI=0.007225
704 ! BI=0.9674
705 ! ENDIF
706 
707  qsat = pq0/pl*exp(a2*(tl-a3)/(tl-a4))
708  rhl = ql/qsat
709 !
710  IF(rhl > 1.)THEN
711  rhl = 1.
712  ql = rhl*qsat
713  ENDIF
714 !
715  IF(rhl < rhmin)THEN
716  rhl = rhmin
717  ql = rhl*qsat
718  ENDIF
719 !
720  tvrl = tl*(1.+0.608*ql)
721  tvrblo = tvrl*(spl(lp)/pl)**rgamog
722  tblo = tvrblo/(1.+0.608*ql)
723 !
724 ! TMT0=TBLO-A3
725 ! TMT15=MIN(TMT0,-15.)
726 ! AI=0.008855
727 ! BI=1.
728 ! IF(TMT0 < -20.)THEN
729 ! AI=0.007225
730 ! BI=0.9674
731 ! ENDIF
732 
733  qsat = pq0/spl(lp)*exp(a2*(tblo-a3)/(tblo-a4))
734  tsl(i,j) = tblo
735  qblo = rhl*qsat
736  qsl(i,j) = max(1.e-12,qblo)
737  END IF ! endif loop for deducing T and H differently for GFS
738 
739 ! if(tsl(i,j) > 330. .or. tsl(i,j) < 100.)print*, &
740 ! 'bad isobaric T Q',i,j,lp,tsl(i,j),qsl(i,j),tl,ql,pl
741 
742  IF(gridtype == 'A')THEN
743  usl(i,j) = uh(i,j,llmh)
744  vsl(i,j) = vh(i,j,llmh)
745  END IF
746 ! if ( J == JSTA.and. I == 1.and.me == 0) &
747 ! & print *,'3 USL=',USL(I,J),UH(I,J,LLMH),LLMH
748  wsl(i,j) = wh(i,j,llmh)
749  osl(i,j) = omga(i,j,llmh)
750  q2sl(i,j) = max(0.0,0.5*(q2(i,j,llmh-1)+q2(i,j,llmh)))
751  pnl1 = pint(i,j,ll)
752  fac = 0.0
753  ahf = 0.0
754 
755 ! FSL(I,J)=(PNL1-SPL(LP))/(SPL(LP)+PNL1)
756 ! 1 *(TSL(I,J))*(1.+0.608*QSL(I,J))
757 ! 2 *RD*2.+ZINT(I,J,NL1X(I,J))*G
758 
759 ! FSL(I,J)=FPRS(I,J,LP-1)-RD*(TPRS(I,J,LP-1)
760 ! 1 *(H1+D608*QPRS(I,J,LP-1))
761 ! 2 +TSL(I,J)*(H1+D608*QSL(I,J)))
762 ! 3 *LOG(SPL(LP)/SPL(LP-1))/2.0
763 
764 ! if(abs(SPL(LP)-97500.0) < 0.01)then
765 ! if(gdlat(i,j) > 35.0.and.gdlat(i,j)<=37.0 .and. &
766 ! gdlon(i,j) > -100.0 .and. gdlon(i,j) < -96.0)print*, &
767 ! 'Debug:I,J,FPRS(LP-1),TPRS(LP-1),TSL,SPL(LP),SPL(LP-1)=' &
768 ! ,i,j,FPRS(I,J,LP-1),TPRS(I,J,LP-1),TSL(I,J),SPL(LP) &
769 ! ,SPL(LP-1)
770 ! if(gdlat(i,j) > 35.0.and.gdlat(i,j)<=37.0 .and.
771 ! 1 gdlon(i,j) > -100.0 .and. gdlon(i,j) < -96.0)print*,
772 ! 2 'Debug:I,J,PNL1,TSL,NL1X,ZINT,FSL= ',I,J,PNL1,TSL(I,J)
773 ! 3 ,NL1X(I,J),ZINT(I,J,NL1X(I,J)),FSL(I,J)/G
774 ! end if
775 ! if(lp == lsm)print*,'Debug:undergound T,Q,U,V,FSL='
776 ! 1,TSL(I,J),QSL(I,J),USL(I,J),VSL(I,J),FSL(I,J)
777 !
778 !--- Set hydrometeor fields to zero below ground
779  c1d(i,j) = 0.
780  qw1(i,j) = 0.
781  qi1(i,j) = 0.
782  qr1(i,j) = 0.
783  qs1(i,j) = 0.
784  qg1(i,j) = 0.
785  dbz1(i,j) = dbzmin
786  frime(i,j) = 1.
787  rad(i,j) = 0.
788  o3sl(i,j) = o3(i,j,llmh)
789  IF(cfr(i,j,1)<spval)cfrsl(i,j) = 0.
790  END IF
791 ! Compute heights by interpolating from heights on interface for NAM but
792 ! hydrostaticJ integration for GFS
793 
794  IF(modelname == 'GFS') then
795  l=ll
796  IF(spl(lp) < pmid(i,j,1)) THEN ! above model domain
797  tvd = t(i,j,1)*(1+con_fvirt*q(i,j,1))
798  fsl(i,j) = zmid(i,j,1)-con_rog*tvd *(alsl(lp)-log(pmid(i,j,1)))
799  fsl(i,j) = fsl(i,j)*g
800  ELSE IF(l <= llmh)THEN
801  tvd = t(i,j,l)*(1+con_fvirt*q(i,j,l))
802  tvu = tsl(i,j)*(1+con_fvirt*qsl(i,j))
803  fsl(i,j) = zmid(i,j,l)-con_rog*0.5*(tvd+tvu) &
804  * (alsl(lp)-log(pmid(i,j,l)))
805  fsl(i,j) = fsl(i,j)*g
806  END IF
807  ELSE
808  la = nl1xf(i,j)
809  IF(nl1xf(i,j)<=(llmh+1)) THEN
810  fact = (alsl(lp)-log(pint(i,j,la)))/ &
811  (log(pint(i,j,la))-log(pint(i,j,la-1)))
812  IF(zint(i,j,la) < spval .AND. zint(i,j,la-1) < spval) &
813  fsl(i,j) = zint(i,j,la)+(zint(i,j,la)-zint(i,j,la-1))*fact
814  fsl(i,j) = fsl(i,j)*g
815  ELSE
816  fsl(i,j) = fprs(i,j,lp-1)-rd*(tprs(i,j,lp-1) &
817  * (h1+d608*qprs(i,j,lp-1)) &
818  + tsl(i,j)*(h1+d608*qsl(i,j))) &
819  * log(spl(lp)/spl(lp-1))/2.0
820  END IF
821  END IF
822 
823  enddo ! End of i loop
824  enddo ! End of J loop
825 
826 !
827 !*** FILL THE 3-D-IN-PRESSURE ARRAYS FOR THE MEMBRANE SLP REDUCTION
828 !
829 !$omp parallel do private(i,j)
830  DO j=jsta,jend
831  DO i=ista,iend
832  tprs(i,j,lp) = tsl(i,j)
833  qprs(i,j,lp) = qsl(i,j)
834  fprs(i,j,lp) = fsl(i,j)
835  ENDDO
836  ENDDO
837 !
838 ! VERTICAL INTERPOLATION FOR WIND FOR E GRID
839 !
840  IF(gridtype == 'E')THEN
841  DO j=jsta,jend
842 ! DO I=2,IM-MOD(J,2)
843  DO i=ista_m,iend-mod(j,2)
844 ! IF(i == im/2 .and. j == (jsta+jend)/2)then
845 ! do l=1,lm
846 ! print*,'PMIDV=',PMIDV(i,j,l)
847 ! end do
848 ! end if
849 !
850 !*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER FOR V POINT JUST BELOW
851 !*** THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
852 !
853  nl1x(i,j) = lp1
854  DO l=2,lm
855 
856 ! IF(J == 1 .AND. I < IM)THEN !SOUTHERN BC
857 ! PDV=0.5*(PMID(I,J,L)+PMID(I+1,J,L))
858 ! ELSE IF(J == JM .AND. I < IM)THEN !NORTHERN BC
859 ! PDV=0.5*(PMID(I,J,L)+PMID(I+1,J,L))
860 ! ELSE IF(I == 1 .AND. MOD(J,2) == 0) THEN !WESTERN EVEN BC
861 ! PDV=0.5*(PMID(I,J-1,L)+PMID(I,J+1,L))
862 ! ELSE IF(I == IM .AND. MOD(J,2) == 0) THEN !EASTERN EVEN BC
863 ! PDV=0.5*(PMID(I,J-1,L)+PMID(I,J+1,L))
864 ! ELSE IF (MOD(J,2) < 1) THEN
865 ! PDV=0.25*(PMID(I,J,L)+PMID(I-1,J,L)
866 ! & +PMID(I,J+1,L)+PMID(I,J-1,L))
867 ! ELSE
868 ! PDV=0.25*(PMID(I,J,L)+PMID(I+1,J,L)
869 ! & +PMID(I,J+1,L)+PMID(I,J-1,L))
870 ! END IF
871 ! JJB=JSTA
872 ! IF(MOD(JSTA,2) == 0)JJB=JSTA+1
873 ! JJE=JEND
874 ! IF(MOD(JEND,2) == 0)JJE=JEND-1
875 ! DO J=JJB,JJE,2 !chc
876 ! PDV(IM,J)=PDV(IM-1,J)
877 ! END DO
878 
879  IF(nl1x(i,j) == lp1.AND.pmidv(i,j,l) > spl(lp))THEN
880  nl1x(i,j) = l
881 ! IF(i == im/2 .and. j == jm/2)print*, &
882 ! 'Wind Debug:LP,NL1X',LP,NL1X(I,J)
883  ENDIF
884  ENDDO
885 !
886 ! IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
887 ! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
888 ! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
889 ! WILL EXTRAPOLATE TO THAT POINT
890 !
891 ! IF(NL1X(I,J) == LMP1.AND.PINT(I,J,LMP1) > SPL(LP))THEN
892  IF(nl1x(i,j) == lp1)THEN
893  IF(j == jsta .AND. i < iend)THEN !SOUTHERN BC
894  pdv = 0.5*(pint(i,j,lp1)+pint(i+1,j,lp1))
895  ELSE IF(j == jend .AND. i < iend)THEN !NORTHERN BC
896  pdv = 0.5*(pint(i,j,lp1)+pint(i+1,j,lp1))
897  ELSE IF(i == ista .AND. mod(j,2) == 0) THEN !WESTERN EVEN BC
898  pdv = 0.5*(pint(i,j-1,lp1)+pint(i,j+1,lp1))
899  ELSE IF(i == iend .AND. mod(j,2) == 0) THEN !EASTERN EVEN BC
900  pdv = 0.5*(pint(i,j-1,lp1)+pint(i,j+1,lp1))
901  ELSE IF (mod(j,2) < 1) THEN
902  pdv = 0.25*(pint(i,j,lp1)+pint(i-1,j,lp1) &
903  + pint(i,j+1,lp1)+pint(i,j-1,lp1))
904  ELSE
905  pdv = 0.25*(pint(i,j,lp1)+pint(i+1,j,lp1) &
906  + pint(i,j+1,lp1)+pint(i,j-1,lp1))
907  END IF
908  IF(pdv > spl(lp))THEN
909  nl1x(i,j) = lm
910  END IF
911  ENDIF
912 !
913  ENDDO
914  ENDDO
915 !
916  DO j=jsta,jend
917 ! DO I=1,IM-MOD(j,2)
918  DO i=ista,iend-mod(j,2)
919  ll = nl1x(i,j)
920 !---------------------------------------------------------------------
921 !*** VERTICAL INTERPOLATION OF WINDS FOR A-E GRID
922 !---------------------------------------------------------------------
923 !
924 !HC IF(NL1X(I,J)<=LM)THEN
925  llmh = nint(lmh(i,j))
926 
927  IF(spl(lp) < pint(i,j,2))THEN ! Above second interface
928  IF(uh(i,j,1) < spval) usl(i,j) = uh(i,j,1)
929  IF(vh(i,j,1) < spval) vsl(i,j) = vh(i,j,1)
930 
931  ELSE IF(nl1x(i,j)<=llmh)THEN
932 !
933 !---------------------------------------------------------------------
934 ! INTERPOLATE LINEARLY IN LOG(P)
935 !*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
936 !*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
937 !*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
938 !---------------------------------------------------------------------
939 !
940 
941  fact = (alsl(lp)-log(pmidv(i,j,ll)))/ &
942  (log(pmidv(i,j,ll))-log(pmidv(i,j,ll-1)))
943  IF(uh(i,j,ll) < spval .AND. uh(i,j,ll-1) < spval) &
944  usl(i,j) = uh(i,j,ll)+(uh(i,j,ll)-uh(i,j,ll-1))*fact
945  IF(vh(i,j,ll) < spval .AND. vh(i,j,ll-1) < spval) &
946  vsl(i,j) = vh(i,j,ll)+(vh(i,j,ll)-vh(i,j,ll-1))*fact
947 ! IF(i == im/2 .and. j == jm/2)print*, &
948 ! 'Wind Debug:LP,NL1X,FACT=',LP,NL1X(I,J),FACT
949 !
950 ! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE
951 ! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
952 ! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
953 ! GOUND
954  ELSE
955  IF(uh(i,j,llmh) < spval) usl(i,j) = uh(i,j,llmh)
956  IF(vh(i,j,llmh) < spval) vsl(i,j) = vh(i,j,llmh)
957  END IF
958  ENDDO ! end of i loop
959  ENDDO ! end of j loop
960 
961 ! if(me == 0) print *,'after 230 me=',me,'USL=',USL(1:10,JSTA)
962  jjb = jsta
963  IF(mod(jsta,2) == 0) jjb = jsta+1
964  jje = jend
965  IF(mod(jend,2) == 0) jje = jend-1
966  DO j=jjb,jje,2 !chc
967  usl(iend,j) = usl(iend-1,j)
968  vsl(iend,j) = vsl(iend-1,j)
969  END DO
970  ELSE IF(gridtype=='B')THEN ! B grid wind interpolation
971  DO j=jsta,jend_m
972 ! DO I=1,IM-1
973  DO i=ista,iend_m
974 !*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER FOR V POINT JUST BELOW
975 !*** THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
976 !
977  nl1x(i,j)=lp1
978  DO l=2,lm
979  IF(nl1x(i,j) == lp1.AND.pmidv(i,j,l) > spl(lp))THEN
980  nl1x(i,j) = l
981 ! IF(i == im/2 .and. j == jm/2)print*, &
982 ! 'Wind Debug for B grid:LP,NL1X',LP,NL1X(I,J)
983  ENDIF
984  ENDDO
985 !
986 ! IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
987 ! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
988 ! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
989 ! WILL EXTRAPOLATE TO THAT POINT
990 !
991  IF(nl1x(i,j)==lp1)THEN
992  pdv = 0.25*(pint(i,j,lp1)+pint(i+1,j,lp1) &
993  + pint(i,j+1,lp1)+pint(i+1,j+1,lp1))
994  IF(pdv > spl(lp))THEN
995  nl1x(i,j)=lm
996  END IF
997  ENDIF
998 !
999  ENDDO
1000  ENDDO
1001 !
1002  DO j=jsta,jend_m
1003 ! DO I=1,IM-1
1004  DO i=ista,iend_m
1005  ll = nl1x(i,j)
1006 !---------------------------------------------------------------------
1007 !*** VERTICAL INTERPOLATION OF WINDS FOR A-E GRID
1008 !---------------------------------------------------------------------
1009 !
1010 !HC IF(NL1X(I,J)<=LM)THEN
1011  llmh = nint(lmh(i,j))
1012 
1013  IF(spl(lp) < pint(i,j,2))THEN ! Above second interface
1014  IF(uh(i,j,1) < spval) usl(i,j) = uh(i,j,1)
1015  IF(vh(i,j,1) < spval) vsl(i,j) = vh(i,j,1)
1016 
1017  ELSE IF(nl1x(i,j)<=llmh)THEN
1018 !
1019 !---------------------------------------------------------------------
1020 ! INTERPOLATE LINEARLY IN LOG(P)
1021 !*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
1022 !*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
1023 !*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
1024 !---------------------------------------------------------------------
1025 !
1026 
1027  fact = (alsl(lp)-log(pmidv(i,j,ll)))/ &
1028  (log(pmidv(i,j,ll))-log(pmidv(i,j,ll-1)))
1029  IF(uh(i,j,ll) < spval .AND. uh(i,j,ll-1) < spval) &
1030  usl(i,j)=uh(i,j,ll)+(uh(i,j,ll)-uh(i,j,ll-1))*fact
1031  IF(vh(i,j,ll) < spval .AND. vh(i,j,ll-1) < spval) &
1032  vsl(i,j)=vh(i,j,ll)+(vh(i,j,ll)-vh(i,j,ll-1))*fact
1033 ! IF(i == im/2 .and. j == jm/2)print*, &
1034 ! 'Wind Debug:LP,NL1X,FACT=',LP,NL1X(I,J),FACT
1035 !
1036 ! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE
1037 ! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
1038 ! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
1039 ! GOUND
1040  ELSE
1041  IF(uh(i,j,llmh) < spval)usl(i,j)=uh(i,j,llmh)
1042  IF(vh(i,j,llmh) < spval)vsl(i,j)=vh(i,j,llmh)
1043  END IF
1044  enddo
1045  enddo
1046  END IF ! END OF WIND INTERPOLATION FOR NMM
1047 ! if(me == 0) print *,'after 230 if me=',me,'USL=',USL(1:10,JSTA)
1048 
1049 
1050 !
1051 !---------------------------------------------------------------------
1052 ! LOAD GEOPOTENTIAL AND TEMPERATURE INTO STANDARD LEVEL
1053 ! ARRAYS FOR THE NEXT PASS.
1054 !---------------------------------------------------------------------
1055 !
1056 !*** SAVE 500MB TEMPERATURE FOR LIFTED INDEX.
1057 !
1058  IF(nint(spl(lp)) == 50000)THEN
1059 !$omp parallel do private(i,j)
1060  DO j=jsta,jend
1061  DO i=ista,iend
1062  t500(i,j) = tsl(i,j)
1063  z500(i,j) = fsl(i,j)*gi
1064  ENDDO
1065  ENDDO
1066  ENDIF
1067 
1068 !
1069 !*** SAVE 700MB TEMPERATURE FOR LIFTED INDEX.
1070 !
1071  IF(nint(spl(lp)) == 70000)THEN
1072 !$omp parallel do private(i,j)
1073  DO j=jsta,jend
1074  DO i=ista,iend
1075  t700(i,j) = tsl(i,j)
1076  z700(i,j) = fsl(i,j)*gi
1077  ENDDO
1078  ENDDO
1079  ENDIF
1080 
1081 !
1082 !---------------------------------------------------------------------
1083 !*** CALCULATE 1000MB GEOPOTENTIALS CONSISTENT WITH SLP OBTAINED
1084 !*** FROM THE MESINGER OR NWS SHUELL SLP REDUCTION.
1085 !---------------------------------------------------------------------
1086 !
1087 !*** FROM MESINGER SLP
1088 !
1089 !HC MOVE THIS PART TO THE END OF THIS SUBROUTINE AFTER PSLP IS COMPUTED
1090 !HC IF(IGET(023) > 0.AND.NINT(SPL(LP)) == 100000)THEN
1091 !HC ALPTH=LOG(1.E5)
1092 !HC!$omp parallel do private(i,j)
1093 !HC DO J=JSTA,JEND
1094 !HC DO I=1,IM
1095 !HC IF(FSL(I,J) < SPVAL) THEN
1096 !HC PSLPIJ=PSLP(I,J)
1097 !HC ALPSL=LOG(PSLPIJ)
1098 !HC PSFC=PINT(I,J,NINT(LMH(I,J))+1)
1099 !HC IF(ABS(PSLPIJ-PSFC) < 5.E2) THEN
1100 !HC FSL(I,J)=R*TSL(I,J)*(ALPSL-ALPTH)
1101 !HC ELSE
1102 !HC FSL(I,J)=FIS(I,J)/(ALPSL-LOG(PSFC))*
1103 !HC 1 (ALPSL-ALPTH)
1104 !HC ENDIF
1105 !HC Z1000(I,J)=FSL(I,J)*GI
1106 !HC ELSE
1107 !HC Z1000(I,J)=SPVAL
1108 !HC ENDIF
1109 !HC ENDDO
1110 !HC ENDDO
1111 !
1112 !*** FROM NWS SHUELL SLP. NGMSLP2 COMPUTES 1000MB GEOPOTENTIAL.
1113 !
1114 !HC ELSEIF(IGET(023)<=0.AND.LP == LSM)THEN
1115 !HC IF(IGET(023)<=0.AND.LP == LSM)THEN
1116 !!$omp parallel do private(i,j)
1117 !HC DO J=JSTA,JEND
1118 !HC DO I=1,IM
1119 !HC IF(Z1000(I,J) < SPVAL) THEN
1120 !HC FSL(I,J)=Z1000(I,J)*G
1121 !HC ELSE
1122 !HC FSL(I,J)=SPVAL
1123 !HC ENDIF
1124 !HC ENDDO
1125 !HC ENDDO
1126 !HC ENDIF
1127 !---------------------------------------------------------------------
1128 !---------------------------------------------------------------------
1129 ! *** PART II ***
1130 !---------------------------------------------------------------------
1131 !---------------------------------------------------------------------
1132 !
1133 ! INTERPOLATE/OUTPUT SELECTED FIELDS.
1134 !
1135 !---------------------------------------------------------------------
1136 !
1137 !*** OUTPUT GEOPOTENTIAL (SCALE BY GI)
1138 !
1139  IF(iget(012) > 0)THEN
1140  IF(lvls(lp,iget(012)) > 0)THEN
1141  IF((iget(023) > 0 .OR. iget(445) > 0) .AND. nint(spl(lp)) == 100000) THEN
1142  ! GO TO 222
1143  ELSE
1144 !$omp parallel do private(i,j)
1145  DO j=jsta,jend
1146  DO i=ista,iend
1147  IF(fsl(i,j) < spval) THEN
1148  grid1(i,j) = fsl(i,j)*gi
1149  ELSE
1150  grid1(i,j) = spval
1151  ENDIF
1152  ENDDO
1153  ENDDO
1154 
1155  IF (smflag) THEN
1156 !tgs - smoothing of geopotential heights
1157  if(maptype == 6) then
1158  if(grib=='grib2') then
1159  dxm = (dxval / 360.)*(erad*2.*pi)/1.d6 ! [mm]
1160  endif
1161  else
1162  dxm = dxval
1163  endif
1164  if(grib == 'grib2')then
1165  dxm=dxm/1000.0
1166  endif
1167 ! print *,'dxm=',dxm
1168  nsmooth = nint(5.*(13500./dxm))
1169  call allgetherv(grid1)
1170  do k=1,nsmooth
1171  CALL smooth(grid1,sdummy,im,jm,0.5)
1172  end do
1173  ENDIF
1174  if(grib == 'grib2')then
1175  cfld = cfld + 1
1176  fld_info(cfld)%ifld=iavblfld(iget(012))
1177  fld_info(cfld)%lvl=lvlsxml(lp,iget(012))
1178 !$omp parallel do private(i,j,ii,jj)
1179  do j=1,jend-jsta+1
1180  jj = jsta+j-1
1181  do i=1,iend-ista+1
1182  ii=ista+i-1
1183  datapd(i,j,cfld) = grid1(ii,jj)
1184  enddo
1185  enddo
1186  endif
1187  END IF
1188  ENDIF
1189  ENDIF
1190 !222 CONTINUE
1191 !
1192 !*** TEMPERATURE
1193 !
1194  IF(iget(013) > 0) THEN
1195  IF(lvls(lp,iget(013)) > 0)THEN
1196 !$omp parallel do private(i,j)
1197  DO j=jsta,jend
1198  DO i=ista,iend
1199  grid1(i,j) = tsl(i,j)
1200  ENDDO
1201  ENDDO
1202 
1203  IF (smflag) THEN
1204  nsmooth = nint(3.*(13500./dxm))
1205  call allgetherv(grid1)
1206  do k=1,nsmooth
1207  CALL smooth(grid1,sdummy,im,jm,0.5)
1208  end do
1209  ENDIF
1210 
1211  if(grib == 'grib2')then
1212  cfld = cfld + 1
1213  fld_info(cfld)%ifld = iavblfld(iget(013))
1214  fld_info(cfld)%lvl = lvlsxml(lp,iget(013))
1215 !$omp parallel do private(i,j,ii,jj)
1216  do j=1,jend-jsta+1
1217  jj = jsta+j-1
1218  do i=1,iend-ista+1
1219  ii=ista+i-1
1220  datapd(i,j,cfld) = grid1(ii,jj)
1221  enddo
1222  enddo
1223  endif
1224  ENDIF
1225  ENDIF
1226 
1227 !*** virtual TEMPERATURE
1228 !
1229  IF(iget(910)>0) THEN
1230  IF(lvls(lp,iget(910))>0)THEN
1231 !$omp parallel do private(i,j)
1232  DO j=jsta,jend
1233  DO i=ista,iend
1234  IF(tsl(i,j) < spval .AND. qsl(i,j) < spval) THEN
1235  grid1(i,j) = tsl(i,j)*(1.+0.608*qsl(i,j))
1236  ELSE
1237  grid1(i,j) = spval
1238  ENDIF
1239  ENDDO
1240  ENDDO
1241 
1242  IF (smflag) THEN
1243  nsmooth = nint(3.*(13500./dxm))
1244  call allgetherv(grid1)
1245  do k=1,nsmooth
1246  CALL smooth(grid1,sdummy,im,jm,0.5)
1247  end do
1248  ENDIF
1249 
1250  if(grib=='grib2')then
1251  cfld=cfld+1
1252  fld_info(cfld)%ifld = iavblfld(iget(910))
1253  fld_info(cfld)%lvl = lvlsxml(lp,iget(910))
1254 !$omp parallel do private(i,j,ii,jj)
1255  do j=1,jend-jsta+1
1256  jj = jsta+j-1
1257  do i=1,iend-ista+1
1258  ii=ista+i-1
1259  datapd(i,j,cfld) = grid1(ii,jj)
1260  enddo
1261  enddo
1262  endif
1263  ENDIF
1264  ENDIF
1265 
1266 !
1267 !*** POTENTIAL TEMPERATURE.
1268 !
1269  IF(iget(014) > 0)THEN
1270  IF(lvls(lp,iget(014)) > 0)THEN
1271 
1272  tem = (p1000/spl(lp)) ** capa
1273 !$omp parallel do private(i,j)
1274  DO j=jsta,jend
1275  DO i=ista,iend
1276  IF(tsl(i,j) < spval) THEN
1277  grid1(i,j) = tsl(i,j) * tem
1278  ELSE
1279  grid1(i,j) = spval
1280  ENDIF
1281  ENDDO
1282  ENDDO
1283 !!$omp parallel do private(i,j)
1284 ! DO J=JSTA,JEND
1285 ! DO I=1,IM
1286 ! EGRID2(I,J) = SPL(LP)
1287 ! ENDDO
1288 ! ENDDO
1289 !
1290 ! CALL CALPOT(EGRID2,TSL,EGRID1)
1291 !!$omp parallel do private(i,j)
1292 ! DO J=JSTA,JEND
1293 ! DO I=1,IM
1294 ! GRID1(I,J) = EGRID1(I,J)
1295 ! ENDDO
1296 ! ENDDO
1297 
1298  if(grib == 'grib2')then
1299  cfld = cfld + 1
1300  fld_info(cfld)%ifld=iavblfld(iget(014))
1301  fld_info(cfld)%lvl=lvlsxml(lp,iget(014))
1302 !$omp parallel do private(i,j,ii,jj)
1303  do j=1,jend-jsta+1
1304  jj = jsta+j-1
1305  do i=1,iend-ista+1
1306  ii=ista+i-1
1307  datapd(i,j,cfld) = grid1(ii,jj)
1308  enddo
1309  enddo
1310  endif
1311  ENDIF
1312  ENDIF
1313 !
1314 !*** RELATIVE HUMIDITY.
1315 !
1316 
1317  IF(iget(017) > 0 .OR. iget(257) > 0 .OR. iget(1006) > 0)THEN
1318 ! if ( me == 0) print *,'IGET(17)=',IGET(017),'LP=',LP,IGET(257), &
1319 ! 'LVLS=',LVLS(1,4)
1320  log1=.false.
1321  IF(iget(017) > 0.) then
1322  if(lvls(lp,iget(017)) > 0 ) log1=.true.
1323  endif
1324  IF(iget(257) > 0) then
1325  if(lvls(lp,iget(257)) > 0 ) log1=.true.
1326  endif
1327 
1328 !$omp parallel do private(i,j)
1329  DO j=jsta,jend
1330  DO i=ista,iend
1331  egrid2(i,j) = spl(lp)
1332  ENDDO
1333  ENDDO
1334 !
1335  CALL calrh(egrid2(ista:iend,jsta:jend),tsl(ista:iend,jsta:jend),qsl(ista:iend,jsta:jend),egrid1(ista:iend,jsta:jend))
1336 
1337 !$omp parallel do private(i,j)
1338  DO j=jsta,jend
1339  DO i=ista,iend
1340  IF(egrid1(i,j) < spval) THEN
1341  grid1(i,j) = egrid1(i,j)*100.
1342  ELSE
1343  grid1(i,j) = egrid1(i,j)
1344  ENDIF
1345  ENDDO
1346  ENDDO
1347 
1348  IF (smflag) THEN
1349  nsmooth=nint(2.*(13500./dxm))
1350  call allgetherv(grid1)
1351  do k=1,nsmooth
1352  CALL smooth(grid1,sdummy,im,jm,0.5)
1353  end do
1354  ENDIF
1355 
1356  if ( log1 ) then
1357  if(grib == 'grib2')then
1358  cfld = cfld + 1
1359  fld_info(cfld)%ifld=iavblfld(iget(017))
1360  fld_info(cfld)%lvl=lvlsxml(lp,iget(017))
1361 !$omp parallel do private(i,j,ii,jj)
1362  do j=1,jend-jsta+1
1363  jj = jsta+j-1
1364  do i=1,iend-ista+1
1365  ii=ista+i-1
1366  datapd(i,j,cfld) = grid1(ii,jj)
1367  enddo
1368  enddo
1369  endif
1370 
1371 !$omp parallel do private(i,j)
1372  DO j=jsta,jend
1373  DO i=ista,iend
1374  savrh(i,j) = grid1(i,j)
1375  ENDDO
1376  ENDDO
1377  ENDIF !if (log1 )
1378 
1379 !$omp parallel do private(i,j)
1380  DO j=jsta,jend
1381  DO i=ista,iend
1382  rhprs(i,j,lp) = grid1(i,j)
1383  ENDDO
1384  ENDDO
1385  ENDIF
1386 !
1387 !*** CLOUD FRACTION.
1388 !
1389  IF(iget(331) > 0)THEN
1390  IF(lvls(lp,iget(331)) > 0)THEN
1391 !$omp parallel do private(i,j)
1392  DO j=jsta,jend
1393  DO i=ista,iend
1394  grid1(i,j) = spval
1395  IF(abs(cfrsl(i,j)-spval) > small) THEN
1396  cfrsl(i,j) = min(max(0.0,cfrsl(i,j)),1.0)
1397  grid1(i,j) = cfrsl(i,j)*h100
1398  ENDIF
1399  ENDDO
1400  ENDDO
1401  if(grib == 'grib2')then
1402  cfld = cfld + 1
1403  fld_info(cfld)%ifld = iavblfld(iget(331))
1404  fld_info(cfld)%lvl = lvlsxml(lp,iget(331))
1405 !$omp parallel do private(i,j,ii,jj)
1406  do j=1,jend-jsta+1
1407  jj = jsta+j-1
1408  do i=1,iend-ista+1
1409  ii=ista+i-1
1410  datapd(i,j,cfld) = grid1(ii,jj)
1411  enddo
1412  enddo
1413  endif
1414  ENDIF
1415  ENDIF
1416 !
1417 !*** DEWPOINT TEMPERATURE.
1418 !
1419  IF(iget(015) > 0)THEN
1420  IF(lvls(lp,iget(015)) > 0)THEN
1421 !$omp parallel do private(i,j)
1422  DO j=jsta,jend
1423  DO i=ista,iend
1424  egrid2(i,j) = spl(lp)
1425  ENDDO
1426  ENDDO
1427 !
1428  CALL caldwp(egrid2(ista:iend,jsta:jend),qsl(ista:iend,jsta:jend),egrid1(ista:iend,jsta:jend),tsl(ista:iend,jsta:jend))
1429 !$omp parallel do private(i,j)
1430  DO j=jsta,jend
1431  DO i=ista,iend
1432  IF(tsl(i,j) < spval) THEN
1433  grid1(i,j) = egrid1(i,j)
1434  ELSE
1435  grid1(i,j) = spval
1436  ENDIF
1437  ENDDO
1438  ENDDO
1439  if(grib == 'grib2')then
1440  cfld = cfld + 1
1441  fld_info(cfld)%ifld=iavblfld(iget(015))
1442  fld_info(cfld)%lvl=lvlsxml(lp,iget(015))
1443 !$omp parallel do private(i,j,ii,jj)
1444  do j=1,jend-jsta+1
1445  jj = jsta+j-1
1446  do i=1,iend-ista+1
1447  ii=ista+i-1
1448  datapd(i,j,cfld) = grid1(ii,jj)
1449  enddo
1450  enddo
1451  endif
1452  ENDIF
1453  ENDIF
1454 !
1455 !*** SPECIFIC HUMIDITY.
1456 !
1457  IF(iget(016) > 0)THEN
1458  IF(lvls(lp,iget(016)) > 0)THEN
1459 !$omp parallel do private(i,j)
1460  DO j=jsta,jend
1461  DO i=ista,iend
1462  grid1(i,j) = qsl(i,j)
1463  ENDDO
1464  ENDDO
1465  CALL bound(grid1,zero,h99999)
1466  if(grib == 'grib2')then
1467  cfld = cfld + 1
1468  fld_info(cfld)%ifld=iavblfld(iget(016))
1469  fld_info(cfld)%lvl=lvlsxml(lp,iget(016))
1470 !$omp parallel do private(i,j,ii,jj)
1471  do j=1,jend-jsta+1
1472  jj = jsta+j-1
1473  do i=1,iend-ista+1
1474  ii=ista+i-1
1475  datapd(i,j,cfld) = grid1(ii,jj)
1476  enddo
1477  enddo
1478  endif
1479  ENDIF
1480  ENDIF
1481 !
1482 !*** OMEGA
1483 !
1484  IF(iget(020) > 0)THEN
1485  IF(lvls(lp,iget(020)) > 0)THEN
1486 !$omp parallel do private(i,j)
1487  DO j=jsta,jend
1488  DO i=ista,iend
1489  grid1(i,j) = osl(i,j)
1490  ENDDO
1491  ENDDO
1492 
1493  IF (smflag .or. ioform == 'binarympiio' ) THEN
1494  call allgetherv(grid1)
1495  if (ioform == 'binarympiio') then
1496 ! nsmooth = max(2, min(30,nint(jm/94.0)))
1497 ! do k=1,5
1498  CALL smoothc(grid1,sdummy,im,jm,0.5)
1499  CALL smoothc(grid1,sdummy,im,jm,-0.5)
1500 ! enddo
1501  else
1502  nsmooth = nint(3.*(13500./dxm))
1503 ! endif
1504  do k=1,nsmooth
1505  CALL smooth(grid1,sdummy,im,jm,0.5)
1506  end do
1507  endif
1508  ENDIF
1509 
1510  if(grib == 'grib2')then
1511  cfld = cfld + 1
1512  fld_info(cfld)%ifld=iavblfld(iget(020))
1513  fld_info(cfld)%lvl=lvlsxml(lp,iget(020))
1514 !$omp parallel do private(i,j,ii,jj)
1515  do j=1,jend-jsta+1
1516  jj = jsta+j-1
1517  do i=1,iend-ista+1
1518  ii=ista+i-1
1519  datapd(i,j,cfld) = grid1(ii,jj)
1520  enddo
1521  enddo
1522  endif
1523  ENDIF
1524  ENDIF
1525 !
1526 !*** W
1527 !
1528  IF(iget(284) > 0)THEN
1529  IF(lvls(lp,iget(284)) > 0)THEN
1530 !$omp parallel do private(i,j)
1531  DO j=jsta,jend
1532  DO i=ista,iend
1533  grid1(i,j) = wsl(i,j)
1534  ENDDO
1535  ENDDO
1536  if(grib == 'grib2')then
1537  cfld = cfld + 1
1538  fld_info(cfld)%ifld=iavblfld(iget(284))
1539  fld_info(cfld)%lvl=lvlsxml(lp,iget(284))
1540 !$omp parallel do private(i,j,ii,jj)
1541  do j=1,jend-jsta+1
1542  jj = jsta+j-1
1543  do i=1,iend-ista+1
1544  ii=ista+i-1
1545  datapd(i,j,cfld) = grid1(ii,jj)
1546  enddo
1547  enddo
1548  endif
1549  ENDIF
1550  ENDIF
1551 !
1552 !*** MOISTURE CONVERGENCE
1553 !
1554  IF(iget(085) > 0)THEN
1555  IF(lvls(lp,iget(085)) > 0)THEN
1556  CALL calmcvg(qsl(ista_2l,jsta_2l),usl(ista_2l,jsta_2l),vsl(ista_2l,jsta_2l),egrid1(ista_2l,jsta_2l))
1557 ! if(me == 0) print *,'after calmcvgme=',me,'USL=',USL(1:10,JSTA)
1558 !$omp parallel do private(i,j)
1559  DO j=jsta,jend
1560  DO i=ista,iend
1561  grid1(i,j) = egrid1(i,j)
1562  ENDDO
1563  ENDDO
1564 !MEB NOT SURE IF I STILL NEED THIS
1565 ! CONVERT TO DIVERGENCE FOR GRIB UNITS
1566 !
1567 ! CALL SCLFLD(GRID1(ista:iend,jsta:jend),-1.0,IM,JM)
1568 !MEB NOT SURE IF I STILL NEED THIS
1569  if(grib == 'grib2')then
1570  cfld = cfld + 1
1571  fld_info(cfld)%ifld=iavblfld(iget(085))
1572  fld_info(cfld)%lvl=lvlsxml(lp,iget(085))
1573 !$omp parallel do private(i,j,ii,jj)
1574  do j=1,jend-jsta+1
1575  jj = jsta+j-1
1576  do i=1,iend-ista+1
1577  ii=ista+i-1
1578  datapd(i,j,cfld) = grid1(ii,jj)
1579  enddo
1580  enddo
1581 ! if(me==0) print *,'in mdl2p,mconv, lp=',fld_info(cfld)%lvl,'lp=',lp
1582  endif
1583  ENDIF
1584  ENDIF
1585 !
1586 !*** U AND/OR V WIND
1587 !
1588  IF(iget(018) > 0.OR.iget(019) > 0)THEN
1589  log1=.false.
1590  IF(iget(018) > 0.) then
1591  if(lvls(lp,iget(018)) > 0 ) log1=.true.
1592  endif
1593  IF(iget(019) > 0) then
1594  if(lvls(lp,iget(019)) > 0 ) log1=.true.
1595  endif
1596  if ( log1 ) then
1597 !$omp parallel do private(i,j)
1598  DO j=jsta,jend
1599  DO i=ista,iend
1600  grid1(i,j) = usl(i,j)
1601  grid2(i,j) = vsl(i,j)
1602  ENDDO
1603  ENDDO
1604 
1605  IF (smflag) THEN
1606  nsmooth=nint(5.*(13500./dxm))
1607  call allgetherv(grid1)
1608  do k=1,nsmooth
1609  CALL smooth(grid1,sdummy,im,jm,0.5)
1610  end do
1611  nsmooth=nint(5.*(13500./dxm))
1612  call allgetherv(grid2)
1613  do k=1,nsmooth
1614  CALL smooth(grid2,sdummy,im,jm,0.5)
1615  end do
1616  ENDIF
1617 
1618  if(grib == 'grib2')then
1619  cfld = cfld + 1
1620  fld_info(cfld)%ifld=iavblfld(iget(018))
1621  fld_info(cfld)%lvl=lvlsxml(lp,iget(018))
1622 !$omp parallel do private(i,j,ii,jj)
1623  do j=1,jend-jsta+1
1624  jj = jsta+j-1
1625  do i=1,iend-ista+1
1626  ii=ista+i-1
1627  datapd(i,j,cfld) = grid1(ii,jj)
1628  enddo
1629  enddo
1630 
1631  cfld = cfld + 1
1632  fld_info(cfld)%ifld=iavblfld(iget(019))
1633  fld_info(cfld)%lvl=lvlsxml(lp,iget(019))
1634 !$omp parallel do private(i,j,ii,jj)
1635  do j=1,jend-jsta+1
1636  jj = jsta+j-1
1637  do i=1,iend-ista+1
1638  ii=ista+i-1
1639  datapd(i,j,cfld) = grid2(ii,jj)
1640  enddo
1641  enddo
1642  endif
1643  ENDIF
1644  ENDIF
1645 !
1646 !*** ABSOLUTE VORTICITY
1647 !
1648  IF (iget(021) > 0) THEN
1649  IF (lvls(lp,iget(021)) > 0) THEN
1650  CALL calvor(usl,vsl,egrid1)
1651 ! print *,'me=',me,'EGRID1=',EGRID1(1:10,JSTA)
1652 !$omp parallel do private(i,j)
1653  DO j=jsta,jend
1654  DO i=ista,iend
1655  grid1(i,j) = egrid1(i,j)
1656  ENDDO
1657  ENDDO
1658 
1659  IF (smflag .or. ioform == 'binarympiio' ) THEN
1660  call allgetherv(grid1)
1661  if (ioform == 'binarympiio') then
1662 ! nsmooth = max(2, min(30,nint(jm/94.0)))
1663 ! do k=1,5
1664  CALL smoothc(grid1,sdummy,im,jm,0.5)
1665  CALL smoothc(grid1,sdummy,im,jm,-0.5)
1666 ! enddo
1667  else
1668  nsmooth = nint(4.*(13500./dxm))
1669 ! endif
1670  do k=1,nsmooth
1671  CALL smooth(grid1,sdummy,im,jm,0.5)
1672  end do
1673  endif
1674  ENDIF
1675 
1676  if(grib == 'grib2')then
1677  cfld = cfld + 1
1678  fld_info(cfld)%ifld=iavblfld(iget(021))
1679  fld_info(cfld)%lvl=lvlsxml(lp,iget(021))
1680 !$omp parallel do private(i,j,ii,jj)
1681  do j=1,jend-jsta+1
1682  jj = jsta+j-1
1683  do i=1,iend-ista+1
1684  ii=ista+i-1
1685  datapd(i,j,cfld) = grid1(ii,jj)
1686  enddo
1687  enddo
1688  endif
1689  ENDIF
1690  ENDIF
1691 !
1692 ! GEOSTROPHIC STREAMFUNCTION.
1693  IF (iget(086) > 0) THEN
1694  IF (lvls(lp,iget(086)) > 0) THEN
1695 !$omp parallel do private(i,j)
1696  DO j=jsta,jend
1697  DO i=ista,iend
1698  IF(fsl(i,j)<spval)THEN
1699  egrid2(i,j) = fsl(i,j)*gi
1700  ENDIF
1701  ENDDO
1702  ENDDO
1703  CALL calstrm(egrid2(ista:iend,jsta:jend),egrid1(ista:iend,jsta:jend))
1704 !$omp parallel do private(i,j)
1705  DO j=jsta,jend
1706  DO i=ista,iend
1707  IF(fsl(i,j) < spval) THEN
1708  grid1(i,j) = egrid1(i,j)
1709  ELSE
1710  grid1(i,j) = spval
1711  ENDIF
1712  ENDDO
1713  ENDDO
1714  if(grib == 'grib2')then
1715  cfld = cfld + 1
1716  fld_info(cfld)%ifld=iavblfld(iget(086))
1717  fld_info(cfld)%lvl=lvlsxml(lp,iget(086))
1718 !$omp parallel do private(i,j,ii,jj)
1719  do j=1,jend-jsta+1
1720  jj = jsta+j-1
1721  do i=1,iend-ista+1
1722  ii=ista+i-1
1723  datapd(i,j,cfld) = grid1(ii,jj)
1724  enddo
1725  enddo
1726  endif
1727  ENDIF
1728  ENDIF
1729 !
1730 !*** TURBULENT KINETIC ENERGY
1731 !
1732  IF (iget(022) > 0) THEN
1733  IF (lvls(lp,iget(022)) > 0) THEN
1734 !$omp parallel do private(i,j)
1735  DO j=jsta,jend
1736  DO i=ista,iend
1737  grid1(i,j) = q2sl(i,j)
1738  ENDDO
1739  ENDDO
1740  if(grib == 'grib2')then
1741  cfld = cfld + 1
1742  fld_info(cfld)%ifld=iavblfld(iget(022))
1743  fld_info(cfld)%lvl=lvlsxml(lp,iget(022))
1744 !$omp parallel do private(i,j,ii,jj)
1745  do j=1,jend-jsta+1
1746  jj = jsta+j-1
1747  do i=1,iend-ista+1
1748  ii=ista+i-1
1749  datapd(i,j,cfld) = grid1(ii,jj)
1750  enddo
1751  enddo
1752  endif
1753  ENDIF
1754  ENDIF
1755 !
1756 !*** CLOUD WATER
1757 !
1758  IF (iget(153) > 0) THEN
1759  IF (lvls(lp,iget(153)) > 0) THEN
1760  IF(imp_physics==99 .or. imp_physics==98)then
1761 ! GFS does not seperate cloud water from ice, hoping to do that in Feb 08 implementation
1762 !$omp parallel do private(i,j)
1763  DO j=jsta,jend
1764  DO i=ista,iend
1765  IF(qw1(i,j) < spval .AND. qi1(i,j) < spval) THEN
1766  grid1(i,j) = qw1(i,j) + qi1(i,j)
1767  qi1(i,j) = spval
1768  ELSE
1769  grid1(i,j) = spval
1770  ENDIF
1771  ENDDO
1772  ENDDO
1773  ELSE
1774 !$omp parallel do private(i,j)
1775  DO j=jsta,jend
1776  DO i=ista,iend
1777  grid1(i,j) = qw1(i,j)
1778  ENDDO
1779  ENDDO
1780  END IF
1781  if(grib == 'grib2')then
1782  cfld = cfld + 1
1783  fld_info(cfld)%ifld=iavblfld(iget(153))
1784  fld_info(cfld)%lvl=lvlsxml(lp,iget(153))
1785 !$omp parallel do private(i,j,ii,jj)
1786  do j=1,jend-jsta+1
1787  jj = jsta+j-1
1788  do i=1,iend-ista+1
1789  ii=ista+i-1
1790  datapd(i,j,cfld) = grid1(ii,jj)
1791  enddo
1792  enddo
1793  endif
1794  ENDIF
1795  ENDIF
1796 !
1797 !*** CLOUD ICE
1798 !
1799  IF (iget(166) > 0) THEN
1800  IF (lvls(lp,iget(166)) > 0) THEN
1801 !$omp parallel do private(i,j)
1802  DO j=jsta,jend
1803  DO i=ista,iend
1804  grid1(i,j) = qi1(i,j)
1805  ENDDO
1806  ENDDO
1807  if(grib == 'grib2')then
1808  cfld = cfld + 1
1809  fld_info(cfld)%ifld=iavblfld(iget(166))
1810  fld_info(cfld)%lvl=lvlsxml(lp,iget(166))
1811 !$omp parallel do private(i,j,ii,jj)
1812  do j=1,jend-jsta+1
1813  jj = jsta+j-1
1814  do i=1,iend-ista+1
1815  ii=ista+i-1
1816  datapd(i,j,cfld) = grid1(ii,jj)
1817  enddo
1818  enddo
1819  endif
1820  ENDIF
1821  ENDIF
1822 !
1823 !--- RAIN
1824  IF (iget(183) > 0) THEN
1825  IF (lvls(lp,iget(183)) > 0) THEN
1826 !$omp parallel do private(i,j)
1827  DO j=jsta,jend
1828  DO i=ista,iend
1829  grid1(i,j) = qr1(i,j)
1830  ENDDO
1831  ENDDO
1832  if(grib == 'grib2')then
1833  cfld = cfld + 1
1834  fld_info(cfld)%ifld=iavblfld(iget(183))
1835  fld_info(cfld)%lvl=lvlsxml(lp,iget(183))
1836 !$omp parallel do private(i,j,ii,jj)
1837  do j=1,jend-jsta+1
1838  jj = jsta+j-1
1839  do i=1,iend-ista+1
1840  ii=ista+i-1
1841  datapd(i,j,cfld) = grid1(ii,jj)
1842  enddo
1843  enddo
1844  endif
1845  ENDIF
1846  ENDIF
1847 !
1848 !--- SNOW
1849  IF (iget(184) > 0) THEN
1850  IF (lvls(lp,iget(184)) > 0) THEN
1851 !$omp parallel do private(i,j)
1852  DO j=jsta,jend
1853  DO i=ista,iend
1854  grid1(i,j) = qs1(i,j)
1855  ENDDO
1856  ENDDO
1857  if(grib == 'grib2')then
1858  cfld = cfld + 1
1859  fld_info(cfld)%ifld=iavblfld(iget(184))
1860  fld_info(cfld)%lvl=lvlsxml(lp,iget(184))
1861 !$omp parallel do private(i,j,ii,jj)
1862  do j=1,jend-jsta+1
1863  jj = jsta+j-1
1864  do i=1,iend-ista+1
1865  ii=ista+i-1
1866  datapd(i,j,cfld) = grid1(ii,jj)
1867  enddo
1868  enddo
1869  endif
1870  ENDIF
1871  ENDIF
1872 !
1873 !--- GRAUPEL
1874  IF (iget(416) > 0) THEN
1875  IF (lvls(lp,iget(416)) > 0) THEN
1876 !$omp parallel do private(i,j)
1877  DO j=jsta,jend
1878  DO i=ista,iend
1879  grid1(i,j) = qg1(i,j)
1880  ENDDO
1881  ENDDO
1882  if(grib == 'grib2')then
1883  cfld = cfld + 1
1884  fld_info(cfld)%ifld=iavblfld(iget(416))
1885  fld_info(cfld)%lvl=lvlsxml(lp,iget(416))
1886 !$omp parallel do private(i,j,ii,jj)
1887  do j=1,jend-jsta+1
1888  jj = jsta+j-1
1889  do i=1,iend-ista+1
1890  ii=ista+i-1
1891  datapd(i,j,cfld) = grid1(ii,jj)
1892  enddo
1893  enddo
1894  endif
1895  ENDIF
1896  ENDIF
1897 
1898 !
1899 !--- TOTAL CONDENSATE
1900  IF (iget(198) > 0) THEN
1901  IF (lvls(lp,iget(198)) > 0) THEN
1902 !$omp parallel do private(i,j)
1903  DO j=jsta,jend
1904  DO i=ista,iend
1905  grid1(i,j) = c1d(i,j)
1906  ENDDO
1907  ENDDO
1908  if(grib == 'grib2')then
1909  cfld = cfld + 1
1910  fld_info(cfld)%ifld=iavblfld(iget(198))
1911  fld_info(cfld)%lvl=lvlsxml(lp,iget(198))
1912 !$omp parallel do private(i,j,ii,jj)
1913  do j=1,jend-jsta+1
1914  jj = jsta+j-1
1915  do i=1,iend-ista+1
1916  ii=ista+i-1
1917  datapd(i,j,cfld) = grid1(ii,jj)
1918  enddo
1919  enddo
1920  endif
1921  ENDIF
1922  ENDIF
1923 !
1924 !--- RIME FACTOR
1925  IF (iget(263) > 0) THEN
1926  IF (lvls(lp,iget(263)) > 0) THEN
1927 !$omp parallel do private(i,j)
1928  DO j=jsta,jend
1929  DO i=ista,iend
1930  grid1(i,j) = frime(i,j)
1931  ENDDO
1932  ENDDO
1933  if(grib == 'grib2')then
1934  cfld = cfld + 1
1935  fld_info(cfld)%ifld=iavblfld(iget(263))
1936  fld_info(cfld)%lvl=lvlsxml(lp,iget(263))
1937 !$omp parallel do private(i,j,ii,jj)
1938  do j=1,jend-jsta+1
1939  jj = jsta+j-1
1940  do i=1,iend-ista+1
1941  ii=ista+i-1
1942  datapd(i,j,cfld) = grid1(ii,jj)
1943  enddo
1944  enddo
1945  endif
1946  ENDIF
1947  ENDIF
1948 !
1949 !--- Temperature tendency by all radiation: == ested by AFWA
1950  IF (iget(294) > 0) THEN
1951  IF (lvls(lp,iget(294)) > 0) THEN
1952 !$omp parallel do private(i,j)
1953  DO j=jsta,jend
1954  DO i=ista,iend
1955  grid1(i,j) = rad(i,j)
1956  ENDDO
1957  ENDDO
1958  if(grib == 'grib2')then
1959  cfld = cfld + 1
1960  fld_info(cfld)%ifld=iavblfld(iget(294))
1961  fld_info(cfld)%lvl=lvlsxml(lp,iget(294))
1962 !$omp parallel do private(i,j,ii,jj)
1963  do j=1,jend-jsta+1
1964  jj = jsta+j-1
1965  do i=1,iend-ista+1
1966  ii=ista+i-1
1967  datapd(i,j,cfld) = grid1(ii,jj)
1968  enddo
1969  enddo
1970  endif
1971  ENDIF
1972  ENDIF
1973 !
1974 !--- Radar Reflectivity
1975  IF (iget(251) > 0) THEN
1976  IF (lvls(lp,iget(251)) > 0) THEN
1977 !$omp parallel do private(i,j)
1978  DO j=jsta,jend
1979  DO i=ista,iend
1980  grid1(i,j) = dbz1(i,j)
1981  ENDDO
1982  ENDDO
1983  if(grib == 'grib2')then
1984  cfld = cfld + 1
1985  fld_info(cfld)%ifld=iavblfld(iget(251))
1986  fld_info(cfld)%lvl=lvlsxml(lp,iget(251))
1987 !$omp parallel do private(i,j,ii,jj)
1988  do j=1,jend-jsta+1
1989  jj = jsta+j-1
1990  do i=1,iend-ista+1
1991  ii=ista+i-1
1992  datapd(i,j,cfld) = grid1(ii,jj)
1993  enddo
1994  enddo
1995  endif
1996  ENDIF
1997  ENDIF
1998 !
1999 !--- IN-FLIGHT ICING CONDITION: ADD BY B. ZHOU
2000  IF(iget(257) > 0)THEN
2001  IF(lvls(lp,iget(257)) > 0)THEN
2002  CALL calicing(tsl(ista:iend,jsta:jend), savrh, osl(ista:iend,jsta:jend), egrid1(ista:iend,jsta:jend))
2003 
2004 !$omp parallel do private(i,j)
2005  DO j=jsta,jend
2006  DO i=ista,iend
2007  grid1(i,j) = egrid1(i,j)
2008  ENDDO
2009  ENDDO
2010  if(grib == 'grib2')then
2011  cfld = cfld + 1
2012  fld_info(cfld)%ifld=iavblfld(iget(257))
2013  fld_info(cfld)%lvl=lvlsxml(lp,iget(257))
2014 !$omp parallel do private(i,j,ii,jj)
2015  do j=1,jend-jsta+1
2016  jj = jsta+j-1
2017  do i=1,iend-ista+1
2018  ii=ista+i-1
2019  datapd(i,j,cfld) = grid1(ii,jj)
2020  enddo
2021  enddo
2022  endif
2023  ENDIF
2024  ENDIF
2025 
2026 !
2027 
2028 !--- CLEAR AIR TURBULENCE (CAT): ADD BY B. ZHOU
2029  IF (lp > 1) THEN
2030  IF(iget(258) > 0)THEN
2031  IF(lvls(lp,iget(258)) > 0)THEN
2032 !$omp parallel do private(i,j)
2033  DO j=jsta,jend
2034  DO i=ista,iend
2035  IF(fsl(i,j)<spval)THEN
2036  grid1(i,j) = fsl(i,j)*gi
2037  ELSE
2038  grid1(i,j) = spval
2039  ENDIF
2040  egrid1(i,j) = spval
2041  ENDDO
2042  ENDDO
2043  CALL calcat(usl(ista_2l,jsta_2l),vsl(ista_2l,jsta_2l),grid1(ista_2l,jsta_2l) &
2044  ,usl_old(ista_2l,jsta_2l),vsl_old(ista_2l,jsta_2l) &
2045  ,fsl_old(ista_2l,jsta_2l),egrid1(ista_2l,jsta_2l))
2046 !$omp parallel do private(i,j)
2047  DO j=jsta,jend
2048  DO i=ista,iend
2049  grid1(i,j) = egrid1(i,j)
2050 ! IF(GRID1(I,J) > 3. .OR. GRID1(I,J) < 0.)
2051 ! + print*,'bad CAT',i,j,GRID1(I,J)
2052  ENDDO
2053  ENDDO
2054  if(grib == 'grib2')then
2055  cfld = cfld + 1
2056  fld_info(cfld)%ifld=iavblfld(iget(258))
2057  fld_info(cfld)%lvl=lvlsxml(lp,iget(258))
2058 !$omp parallel do private(i,j,ii,jj)
2059  do j=1,jend-jsta+1
2060  jj = jsta+j-1
2061  do i=1,iend-ista+1
2062  ii=ista+i-1
2063  datapd(i,j,cfld) = grid1(ii,jj)
2064  enddo
2065  enddo
2066  endif
2067  end if
2068  end if
2069  end if
2070 !
2071 
2072 !--- GFIP IN-FLIGHT ICING POTENTIAL: ADDED BY H CHUANG
2073  IF(iget(450) > 0)THEN
2074  IF(lvls(lp,iget(450)) > 0)THEN
2075 !$omp parallel do private(i,j)
2076  DO j=jsta,jend
2077  DO i=ista,iend
2078  grid1(i,j) = icingfsl(i,j)
2079  ENDDO
2080  ENDDO
2081  if(grib == 'grib2') then
2082  cfld = cfld + 1
2083  fld_info(cfld)%ifld=iavblfld(iget(450))
2084  fld_info(cfld)%lvl=lvlsxml(lp,iget(450))
2085 !$omp parallel do private(i,j,jj)
2086  do j=1,jend-jsta+1
2087  jj = jsta+j-1
2088  do i=1,iend-ista+1
2089  ii=ista+i-1
2090  datapd(i,j,cfld) = grid1(ii,jj)
2091  enddo
2092  enddo
2093  endif
2094  ENDIF
2095  ENDIF
2096 !--- GFIP IN-FLIGHT ICING SEVERITY: ADDED BY Y MAO
2097  IF(iget(480) > 0) THEN
2098  IF(lvls(lp,iget(480)) > 0) THEN
2099 !$omp parallel do private(i,j)
2100  DO j=jsta,jend
2101  DO i=ista,iend
2102  grid1(i,j) = icingvsl(i,j)
2103  ENDDO
2104  ENDDO
2105  if(grib == 'grib2') then
2106  cfld = cfld+1
2107  fld_info(cfld)%ifld=iavblfld(iget(480))
2108  fld_info(cfld)%lvl=lvlsxml(lp,iget(480))
2109 !$omp parallel do private(i,j,jj)
2110  do j=1,jend-jsta+1
2111  jj = jsta+j-1
2112  do i=1,iend-ista+1
2113  ii=ista+i-1
2114  datapd(i,j,cfld) = grid1(ii,jj)
2115  enddo
2116  enddo
2117  endif
2118  ENDIF
2119  ENDIF
2120 !--- GTG EDR turbulence: ADDED BY Y. MAO
2121  IF(iget(464) > 0) THEN
2122  IF(lvls(lp,iget(464)) > 0) THEN
2123 !$omp parallel do private(i,j)
2124  DO j=jsta,jend
2125  DO i=ista,iend
2126  grid1(i,j) = gtgsl(i,j)
2127  ENDDO
2128  ENDDO
2129  if(grib == 'grib2') then
2130  cfld = cfld+1
2131  fld_info(cfld)%ifld=iavblfld(iget(464))
2132  fld_info(cfld)%lvl=lvlsxml(lp,iget(464))
2133 !$omp parallel do private(i,j,jj)
2134  do j=1,jend-jsta+1
2135  jj = jsta+j-1
2136  do i=1,iend-ista+1
2137  ii=ista+i-1
2138  datapd(i,j,cfld) = grid1(ii,jj)
2139  enddo
2140  enddo
2141  endif
2142  ENDIF
2143  ENDIF
2144 !--- GTG CAT turbulence: ADDED BY Y. MAO
2145  IF(iget(465) > 0) THEN
2146  IF(lvls(lp,iget(465)) > 0) THEN
2147 !$omp parallel do private(i,j)
2148  DO j=jsta,jend
2149  DO i=ista,iend
2150  grid1(i,j) = catsl(i,j)
2151  ENDDO
2152  ENDDO
2153  if(grib == 'grib2') then
2154  cfld = cfld+1
2155  fld_info(cfld)%ifld=iavblfld(iget(465))
2156  fld_info(cfld)%lvl=lvlsxml(lp,iget(465))
2157 !$omp parallel do private(i,j,jj)
2158  do j=1,jend-jsta+1
2159  jj = jsta+j-1
2160  do i=1,iend-ista+1
2161  ii=ista+i-1
2162  datapd(i,j,cfld) = grid1(ii,jj)
2163  enddo
2164  enddo
2165  endif
2166  ENDIF
2167  ENDIF
2168 !--- GTG MWT turbulence: ADDED BY Y. MAO
2169  IF(iget(466) > 0) THEN
2170  IF(lvls(lp,iget(466)) > 0) THEN
2171 !$omp parallel do private(i,j)
2172  DO j=jsta,jend
2173  DO i=ista,iend
2174  grid1(i,j) = mwtsl(i,j)
2175  ENDDO
2176  ENDDO
2177  if(grib == 'grib2') then
2178  cfld = cfld+1
2179  fld_info(cfld)%ifld=iavblfld(iget(466))
2180  fld_info(cfld)%lvl=lvlsxml(lp,iget(466))
2181 !$omp parallel do private(i,j,jj)
2182  do j=1,jend-jsta+1
2183  jj = jsta+j-1
2184  do i=1,iend-ista+1
2185  ii=ista+i-1
2186  datapd(i,j,cfld) = grid1(ii,jj)
2187  enddo
2188  enddo
2189  endif
2190  ENDIF
2191  ENDIF
2192 
2193 !$omp parallel do private(i,j)
2194  DO j=jsta_2l,jend_2u
2195  DO i=ista_2l,iend_2u
2196  usl_old(i,j) = usl(i,j)
2197  vsl_old(i,j) = vsl(i,j)
2198  IF(fsl(i,j)<spval)THEN
2199  fsl_old(i,j) = fsl(i,j)*gi
2200  ELSE
2201  fsl_old(i,j) = spval
2202  ENDIF
2203  ENDDO
2204  ENDDO
2205 !
2206 !--- OZONE
2207  IF (iget(268) > 0) THEN
2208  IF (lvls(lp,iget(268)) > 0) THEN
2209 !$omp parallel do private(i,j)
2210  DO j=jsta,jend
2211  DO i=ista,iend
2212  grid1(i,j) = o3sl(i,j)
2213  ENDDO
2214  ENDDO
2215 ! print *,'in mdl2p,o3sl=',minval(o3sl(1:im,jsta:jend)), &
2216 ! minval(o3sl(1:im,jsta:jend))
2217  if(grib == 'grib2')then
2218  cfld = cfld + 1
2219  fld_info(cfld)%ifld=iavblfld(iget(268))
2220  fld_info(cfld)%lvl=lvlsxml(lp,iget(268))
2221 !$omp parallel do private(i,j,ii,jj)
2222  do j=1,jend-jsta+1
2223  jj = jsta+j-1
2224  do i=1,iend-ista+1
2225  ii=ista+i-1
2226  datapd(i,j,cfld) = grid1(ii,jj)
2227  enddo
2228  enddo
2229  endif
2230  ENDIF
2231  ENDIF
2232 ! E. James - 8 Dec 2017: SMOKE from WRF-CHEM
2233 !--- SMOKE
2234  IF (iget(738) > 0) THEN
2235  IF (lvls(lp,iget(738)) > 0) THEN
2236 !$omp parallel do private(i,j)
2237  DO j=jsta,jend
2238  DO i=ista,iend
2239  IF(smokesl(i,j,1)<spval.and.spl(lp)<spval.and.tsl(i,j)<spval)THEN
2240  grid1(i,j) = (1./rd)*smokesl(i,j,1)*(spl(lp)/(tsl(i,j)*(1e9)))
2241  ELSE
2242  grid1(i,j) = spval
2243  ENDIF
2244  ENDDO
2245  ENDDO
2246  if(grib == 'grib2')then
2247  cfld = cfld + 1
2248  fld_info(cfld)%ifld=iavblfld(iget(738))
2249  fld_info(cfld)%lvl=lvlsxml(lp,iget(738))
2250 !$omp parallel do private(i,j,ii,jj)
2251  do j=1,jend-jsta+1
2252  jj = jsta+j-1
2253  do i=1,iend-ista+1
2254  ii=ista+i-1
2255  datapd(i,j,cfld) = grid1(ii,jj)
2256  enddo
2257  enddo
2258  endif
2259  ENDIF
2260  ENDIF
2261 ! E. James - 14 Sep 2022: DUST from RRFS
2262  IF (iget(743) > 0) THEN
2263  IF (lvls(lp,iget(743)) > 0) THEN
2264 !$omp parallel do private(i,j)
2265  DO j=jsta,jend
2266  DO i=ista,iend
2267  IF(fv3dustsl(i,j,1)<spval.and.spl(lp)<spval.and.tsl(i,j)<spval)THEN
2268  grid1(i,j) = (1./rd)*fv3dustsl(i,j,1)*(spl(lp)/(tsl(i,j)*(1e9)))
2269  ELSE
2270  grid1(i,j) = spval
2271  ENDIF
2272  ENDDO
2273  ENDDO
2274  if(grib == 'grib2')then
2275  cfld = cfld + 1
2276  fld_info(cfld)%ifld=iavblfld(iget(743))
2277  fld_info(cfld)%lvl=lvlsxml(lp,iget(743))
2278 !$omp parallel do private(i,j,ii,jj)
2279  do j=1,jend-jsta+1
2280  jj = jsta+j-1
2281  do i=1,iend-ista+1
2282  ii=ista+i-1
2283  datapd(i,j,cfld) = grid1(ii,jj)
2284  enddo
2285  enddo
2286  endif
2287  ENDIF
2288  ENDIF
2289 ! E. James - 23 Feb 2023: COARSEPM from RRFS
2290  IF (iget(1013) > 0) THEN
2291  IF (lvls(lp,iget(1013)) > 0) THEN
2292 !$omp parallel do private(i,j)
2293  DO j=jsta,jend
2294  DO i=ista,iend
2295  IF(coarsepmsl(i,j,1)<spval.and.spl(lp)<spval.and.tsl(i,j)<spval)THEN
2296  grid1(i,j) = (1./rd)*coarsepmsl(i,j,1)*(spl(lp)/(tsl(i,j)*(1e9)))
2297  ELSE
2298  grid1(i,j) = spval
2299  ENDIF
2300  ENDDO
2301  ENDDO
2302  if(grib == 'grib2')then
2303  cfld = cfld + 1
2304  fld_info(cfld)%ifld=iavblfld(iget(1013))
2305  fld_info(cfld)%lvl=lvlsxml(lp,iget(1013))
2306 !$omp parallel do private(i,j,ii,jj)
2307  do j=1,jend-jsta+1
2308  jj = jsta+j-1
2309  do i=1,iend-ista+1
2310  ii=ista+i-1
2311  datapd(i,j,cfld) = grid1(ii,jj)
2312  enddo
2313  enddo
2314  endif
2315  ENDIF
2316  ENDIF
2317 
2318  if(iostatusd3d==0 .and. d3d_on) then
2319 !--- longwave tendency
2320  IF (iget(355) > 0) THEN
2321  IF (lvls(lp,iget(355)) > 0) THEN
2322 !$omp parallel do private(i,j)
2323  DO j=jsta,jend
2324  DO i=ista,iend
2325  grid1(i,j) = d3dsl(i,j,1)
2326  ENDDO
2327  ENDDO
2328  id(1:25)=0
2329  itd3d = nint(td3d)
2330  if (itd3d /= 0) then
2331  ifincr = mod(ifhr,itd3d)
2332  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2333  else
2334  ifincr = 0
2335  endif
2336  id(18) = 0
2337  id(19) = ifhr
2338  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2339  id(20) = 3
2340  IF (ifincr == 0) THEN
2341  id(18) = ifhr-itd3d
2342  ELSE
2343  id(18) = ifhr-ifincr
2344  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2345  ENDIF
2346  if(grib == 'grib2')then
2347  cfld = cfld + 1
2348  fld_info(cfld)%ifld=iavblfld(iget(355))
2349  fld_info(cfld)%lvl=lvlsxml(lp,iget(355))
2350  if(itd3d==0) then
2351  fld_info(cfld)%ntrange=0
2352  else
2353  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2354  endif
2355  fld_info(cfld)%tinvstat=itd3d
2356 !$omp parallel do private(i,j,ii,jj)
2357  do j=1,jend-jsta+1
2358  jj = jsta+j-1
2359  do i=1,iend-ista+1
2360  ii=ista+i-1
2361  datapd(i,j,cfld) = grid1(ii,jj)
2362  enddo
2363  enddo
2364  endif
2365  ENDIF
2366  ENDIF
2367 !--- longwave tendency
2368  IF (iget(354) > 0) THEN
2369  IF (lvls(lp,iget(354)) > 0) THEN
2370 !$omp parallel do private(i,j)
2371  DO j=jsta,jend
2372  DO i=ista,iend
2373  grid1(i,j) = d3dsl(i,j,2)
2374  ENDDO
2375  ENDDO
2376  id(1:25)=0
2377  itd3d = nint(td3d)
2378  if (itd3d /= 0) then
2379  ifincr = mod(ifhr,itd3d)
2380  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2381  else
2382  ifincr = 0
2383  endif
2384  id(18) = 0
2385  id(19) = ifhr
2386  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2387  id(20) = 3
2388  IF (ifincr == 0) THEN
2389  id(18) = ifhr-itd3d
2390  ELSE
2391  id(18) = ifhr-ifincr
2392  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2393  ENDIF
2394  if(grib == 'grib2')then
2395  cfld = cfld + 1
2396  fld_info(cfld)%ifld=iavblfld(iget(354))
2397  fld_info(cfld)%lvl=lvlsxml(lp,iget(354))
2398  if(itd3d==0) then
2399  fld_info(cfld)%ntrange=0
2400  else
2401  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2402  endif
2403  fld_info(cfld)%tinvstat=itd3d
2404 !$omp parallel do private(i,j,ii,jj)
2405  do j=1,jend-jsta+1
2406  jj = jsta+j-1
2407  do i=1,iend-ista+1
2408  ii=ista+i-1
2409  datapd(i,j,cfld) = grid1(ii,jj)
2410  enddo
2411  enddo
2412  endif
2413  ENDIF
2414  ENDIF
2415 !--- longwave tendency
2416  IF (iget(356) > 0) THEN
2417  IF (lvls(lp,iget(356)) > 0) THEN
2418 !$omp parallel do private(i,j)
2419  DO j=jsta,jend
2420  DO i=ista,iend
2421  grid1(i,j) = d3dsl(i,j,3)
2422  ENDDO
2423  ENDDO
2424  id(1:25)=0
2425  itd3d = nint(td3d)
2426  if (itd3d /= 0) then
2427  ifincr = mod(ifhr,itd3d)
2428  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2429  else
2430  ifincr = 0
2431  endif
2432  id(18) = 0
2433  id(19) = ifhr
2434  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2435  id(20) = 3
2436  IF (ifincr == 0) THEN
2437  id(18) = ifhr-itd3d
2438  ELSE
2439  id(18) = ifhr-ifincr
2440  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2441  ENDIF
2442  if(grib == 'grib2')then
2443  cfld = cfld + 1
2444  fld_info(cfld)%ifld=iavblfld(iget(356))
2445  fld_info(cfld)%lvl=lvlsxml(lp,iget(356))
2446  if(itd3d==0) then
2447  fld_info(cfld)%ntrange=0
2448  else
2449  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2450  endif
2451  fld_info(cfld)%tinvstat=itd3d
2452 !$omp parallel do private(i,j,ii,jj)
2453  do j=1,jend-jsta+1
2454  jj = jsta+j-1
2455  do i=1,iend-ista+1
2456  ii=ista+i-1
2457  datapd(i,j,cfld) = grid1(ii,jj)
2458  enddo
2459  enddo
2460  endif
2461  ENDIF
2462  ENDIF
2463 !--- longwave tendency
2464  IF (iget(357) > 0) THEN
2465  IF (lvls(lp,iget(357)) > 0) THEN
2466 !$omp parallel do private(i,j)
2467  DO j=jsta,jend
2468  DO i=ista,iend
2469  grid1(i,j) = d3dsl(i,j,4)
2470  ENDDO
2471  ENDDO
2472  id(1:25)=0
2473  itd3d = nint(td3d)
2474  if (itd3d /= 0) then
2475  ifincr = mod(ifhr,itd3d)
2476  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2477  else
2478  ifincr = 0
2479  endif
2480  id(18) = 0
2481  id(19) = ifhr
2482  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2483  id(20) = 3
2484  IF (ifincr == 0) THEN
2485  id(18) = ifhr-itd3d
2486  ELSE
2487  id(18) = ifhr-ifincr
2488  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2489  ENDIF
2490  if(grib == 'grib2')then
2491  cfld = cfld + 1
2492  fld_info(cfld)%ifld=iavblfld(iget(357))
2493  fld_info(cfld)%lvl=lvlsxml(lp,iget(357))
2494  if(itd3d==0) then
2495  fld_info(cfld)%ntrange=0
2496  else
2497  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2498  endif
2499  fld_info(cfld)%tinvstat=itd3d
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  ENDIF
2510  ENDIF
2511 !--- longwave tendency
2512  IF (iget(358) > 0) THEN
2513  IF (lvls(lp,iget(358)) > 0) THEN
2514 !$omp parallel do private(i,j)
2515  DO j=jsta,jend
2516  DO i=ista,iend
2517  grid1(i,j) = d3dsl(i,j,5)
2518  ENDDO
2519  ENDDO
2520  id(1:25)=0
2521  itd3d = nint(td3d)
2522  if (itd3d /= 0) then
2523  ifincr = mod(ifhr,itd3d)
2524  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2525  else
2526  ifincr = 0
2527  endif
2528  id(18) = 0
2529  id(19) = ifhr
2530  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2531  id(20) = 3
2532  IF (ifincr == 0) THEN
2533  id(18) = ifhr-itd3d
2534  ELSE
2535  id(18) = ifhr-ifincr
2536  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2537  ENDIF
2538  if(grib == 'grib2')then
2539  cfld = cfld + 1
2540  fld_info(cfld)%ifld=iavblfld(iget(358))
2541  fld_info(cfld)%lvl=lvlsxml(lp,iget(358))
2542  if(itd3d==0) then
2543  fld_info(cfld)%ntrange=0
2544  else
2545  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2546  endif
2547  fld_info(cfld)%tinvstat=itd3d
2548 !$omp parallel do private(i,j,ii,jj)
2549  do j=1,jend-jsta+1
2550  jj = jsta+j-1
2551  do i=1,iend-ista+1
2552  ii=ista+i-1
2553  datapd(i,j,cfld) = grid1(ii,jj)
2554  enddo
2555  enddo
2556  endif
2557  ENDIF
2558  ENDIF
2559 !--- longwave tendency
2560  IF (iget(359) > 0) THEN
2561  IF (lvls(lp,iget(359)) > 0) THEN
2562 !$omp parallel do private(i,j)
2563  DO j=jsta,jend
2564  DO i=ista,iend
2565  grid1(i,j) = d3dsl(i,j,6)
2566  ENDDO
2567  ENDDO
2568  id(1:25)=0
2569  itd3d = nint(td3d)
2570  if (itd3d /= 0) then
2571  ifincr = mod(ifhr,itd3d)
2572  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2573  else
2574  ifincr = 0
2575  endif
2576  id(18) = 0
2577  id(19) = ifhr
2578  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2579  id(20) = 3
2580  IF (ifincr == 0) THEN
2581  id(18) = ifhr-itd3d
2582  ELSE
2583  id(18) = ifhr-ifincr
2584  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2585  ENDIF
2586  if(grib == 'grib2')then
2587  cfld = cfld + 1
2588  fld_info(cfld)%ifld=iavblfld(iget(359))
2589  fld_info(cfld)%lvl=lvlsxml(lp,iget(359))
2590  if(itd3d==0) then
2591  fld_info(cfld)%ntrange=0
2592  else
2593  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2594  endif
2595  fld_info(cfld)%tinvstat=itd3d
2596 !$omp parallel do private(i,j,ii,jj)
2597  do j=1,jend-jsta+1
2598  jj = jsta+j-1
2599  do i=1,iend-ista+1
2600  ii=ista+i-1
2601  datapd(i,j,cfld) = grid1(ii,jj)
2602  enddo
2603  enddo
2604  endif
2605  ENDIF
2606  ENDIF
2607 !--- longwave tendency
2608  IF (iget(360) > 0) THEN
2609  IF (lvls(lp,iget(360)) > 0) THEN
2610 !$omp parallel do private(i,j)
2611  DO j=jsta,jend
2612  DO i=ista,iend
2613  grid1(i,j) = d3dsl(i,j,7)
2614  ENDDO
2615  ENDDO
2616  id(1:25)=0
2617  itd3d = nint(td3d)
2618  if (itd3d /= 0) then
2619  ifincr = mod(ifhr,itd3d)
2620  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2621  else
2622  ifincr = 0
2623  endif
2624  id(18) = 0
2625  id(19) = ifhr
2626  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2627  id(20) = 3
2628  IF (ifincr == 0) THEN
2629  id(18) = ifhr-itd3d
2630  ELSE
2631  id(18) = ifhr-ifincr
2632  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2633  ENDIF
2634  if(grib == 'grib2')then
2635  cfld = cfld + 1
2636  fld_info(cfld)%ifld=iavblfld(iget(360))
2637  fld_info(cfld)%lvl=lvlsxml(lp,iget(360))
2638  if(itd3d==0) then
2639  fld_info(cfld)%ntrange=0
2640  else
2641  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2642  endif
2643  fld_info(cfld)%tinvstat=itd3d
2644 !$omp parallel do private(i,j,ii,jj)
2645  do j=1,jend-jsta+1
2646  jj = jsta+j-1
2647  do i=1,iend-ista+1
2648  ii=ista+i-1
2649  datapd(i,j,cfld) = grid1(ii,jj)
2650  enddo
2651  enddo
2652  endif
2653  ENDIF
2654  ENDIF
2655 !--- longwave tendency
2656  IF (iget(361) > 0) THEN
2657  IF (lvls(lp,iget(361)) > 0) THEN
2658 !$omp parallel do private(i,j)
2659  DO j=jsta,jend
2660  DO i=ista,iend
2661  grid1(i,j) = d3dsl(i,j,8)
2662  ENDDO
2663  ENDDO
2664  id(1:25)=0
2665  itd3d = nint(td3d)
2666  if (itd3d /= 0) then
2667  ifincr = mod(ifhr,itd3d)
2668  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2669  else
2670  ifincr = 0
2671  endif
2672  id(18) = 0
2673  id(19) = ifhr
2674  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2675  id(20) = 3
2676  IF (ifincr == 0) THEN
2677  id(18) = ifhr-itd3d
2678  ELSE
2679  id(18) = ifhr-ifincr
2680  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2681  ENDIF
2682  if(grib == 'grib2')then
2683  cfld = cfld + 1
2684  fld_info(cfld)%ifld=iavblfld(iget(361))
2685  fld_info(cfld)%lvl=lvlsxml(lp,iget(361))
2686  if(itd3d==0) then
2687  fld_info(cfld)%ntrange=0
2688  else
2689  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2690  endif
2691  fld_info(cfld)%tinvstat=itd3d
2692 !$omp parallel do private(i,j,ii,jj)
2693  do j=1,jend-jsta+1
2694  jj = jsta+j-1
2695  do i=1,iend-ista+1
2696  ii=ista+i-1
2697  datapd(i,j,cfld) = grid1(ii,jj)
2698  enddo
2699  enddo
2700  endif
2701  ENDIF
2702  ENDIF
2703 !--- longwave tendency
2704  IF (iget(362) > 0) THEN
2705  IF (lvls(lp,iget(362)) > 0) THEN
2706 !$omp parallel do private(i,j)
2707  DO j=jsta,jend
2708  DO i=ista,iend
2709  grid1(i,j) = d3dsl(i,j,9)
2710  ENDDO
2711  ENDDO
2712  id(1:25)=0
2713  itd3d = nint(td3d)
2714  if (itd3d /= 0) then
2715  ifincr = mod(ifhr,itd3d)
2716  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2717  else
2718  ifincr = 0
2719  endif
2720  id(18) = 0
2721  id(19) = ifhr
2722  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2723  id(20) = 3
2724  IF (ifincr == 0) THEN
2725  id(18) = ifhr-itd3d
2726  ELSE
2727  id(18) = ifhr-ifincr
2728  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2729  ENDIF
2730  if(grib == 'grib2')then
2731  cfld = cfld + 1
2732  fld_info(cfld)%ifld=iavblfld(iget(362))
2733  fld_info(cfld)%lvl=lvlsxml(lp,iget(362))
2734  if(itd3d==0) then
2735  fld_info(cfld)%ntrange=0
2736  else
2737  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2738  endif
2739  fld_info(cfld)%tinvstat=itd3d
2740 !$omp parallel do private(i,j,ii,jj)
2741  do j=1,jend-jsta+1
2742  jj = jsta+j-1
2743  do i=1,iend-ista+1
2744  ii=ista+i-1
2745  datapd(i,j,cfld) = grid1(ii,jj)
2746  enddo
2747  enddo
2748  endif
2749  ENDIF
2750  ENDIF
2751 !--- longwave tendency
2752  IF (iget(363) > 0) THEN
2753  IF (lvls(lp,iget(363)) > 0) THEN
2754 !$omp parallel do private(i,j)
2755  DO j=jsta,jend
2756  DO i=ista,iend
2757  grid1(i,j) = d3dsl(i,j,10)
2758  ENDDO
2759  ENDDO
2760  id(1:25)=0
2761  itd3d = nint(td3d)
2762  if (itd3d /= 0) then
2763  ifincr = mod(ifhr,itd3d)
2764  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2765  else
2766  ifincr = 0
2767  endif
2768  id(02)=133 ! Table 133
2769  id(18) = 0
2770  id(19) = ifhr
2771  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2772  id(20) = 3
2773  IF (ifincr == 0) THEN
2774  id(18) = ifhr-itd3d
2775  ELSE
2776  id(18) = ifhr-ifincr
2777  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2778  ENDIF
2779  if(grib == 'grib2')then
2780  cfld = cfld + 1
2781  fld_info(cfld)%ifld=iavblfld(iget(363))
2782  fld_info(cfld)%lvl=lvlsxml(lp,iget(363))
2783  if(itd3d==0) then
2784  fld_info(cfld)%ntrange=0
2785  else
2786  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2787  endif
2788  fld_info(cfld)%tinvstat=itd3d
2789 !$omp parallel do private(i,j,ii,jj)
2790  do j=1,jend-jsta+1
2791  jj = jsta+j-1
2792  do i=1,iend-ista+1
2793  ii=ista+i-1
2794  datapd(i,j,cfld) = grid1(ii,jj)
2795  enddo
2796  enddo
2797  endif
2798  ENDIF
2799  ENDIF
2800 !--- longwave tendency
2801  IF (iget(364) > 0) THEN
2802  IF (lvls(lp,iget(364)) > 0) THEN
2803 !$omp parallel do private(i,j)
2804  DO j=jsta,jend
2805  DO i=ista,iend
2806  grid1(i,j) = d3dsl(i,j,11)
2807  ENDDO
2808  ENDDO
2809  id(1:25)=0
2810  itd3d = nint(td3d)
2811  if (itd3d /= 0) then
2812  ifincr = mod(ifhr,itd3d)
2813  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2814  else
2815  ifincr = 0
2816  endif
2817  id(02)=133 ! Table 133
2818  id(18) = 0
2819  id(19) = ifhr
2820  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2821  id(20) = 3
2822  IF (ifincr == 0) THEN
2823  id(18) = ifhr-itd3d
2824  ELSE
2825  id(18) = ifhr-ifincr
2826  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2827  ENDIF
2828  if(grib == 'grib2')then
2829  cfld = cfld + 1
2830  fld_info(cfld)%ifld=iavblfld(iget(364))
2831  fld_info(cfld)%lvl=lvlsxml(lp,iget(364))
2832  if(itd3d==0) then
2833  fld_info(cfld)%ntrange=0
2834  else
2835  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2836  endif
2837  fld_info(cfld)%tinvstat=itd3d
2838 !$omp parallel do private(i,j,ii,jj)
2839  do j=1,jend-jsta+1
2840  jj = jsta+j-1
2841  do i=1,iend-ista+1
2842  ii=ista+i-1
2843  datapd(i,j,cfld) = grid1(ii,jj)
2844  enddo
2845  enddo
2846  endif
2847  ENDIF
2848  ENDIF
2849 !--- longwave tendency
2850  IF (iget(365) > 0) THEN
2851  IF (lvls(lp,iget(365)) > 0) THEN
2852 !$omp parallel do private(i,j)
2853  DO j=jsta,jend
2854  DO i=ista,iend
2855  grid1(i,j) = d3dsl(i,j,12)
2856  ENDDO
2857  ENDDO
2858  id(1:25)=0
2859  itd3d = nint(td3d)
2860  if (itd3d /= 0) then
2861  ifincr = mod(ifhr,itd3d)
2862  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2863  else
2864  ifincr = 0
2865  endif
2866  id(02)=133 ! Table 133
2867  id(18) = 0
2868  id(19) = ifhr
2869  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2870  id(20) = 3
2871  IF (ifincr == 0) THEN
2872  id(18) = ifhr-itd3d
2873  ELSE
2874  id(18) = ifhr-ifincr
2875  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2876  ENDIF
2877  if(grib == 'grib2')then
2878  cfld = cfld + 1
2879  fld_info(cfld)%ifld=iavblfld(iget(365))
2880  fld_info(cfld)%lvl=lvlsxml(lp,iget(365))
2881  if(itd3d==0) then
2882  fld_info(cfld)%ntrange=0
2883  else
2884  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2885  endif
2886  fld_info(cfld)%tinvstat=itd3d
2887 !$omp parallel do private(i,j,ii,jj)
2888  do j=1,jend-jsta+1
2889  jj = jsta+j-1
2890  do i=1,iend-ista+1
2891  ii=ista+i-1
2892  datapd(i,j,cfld) = grid1(ii,jj)
2893  enddo
2894  enddo
2895  endif
2896  ENDIF
2897  ENDIF
2898 !--- longwave tendency
2899  IF (iget(366) > 0) THEN
2900  IF (lvls(lp,iget(366)) > 0) THEN
2901 !$omp parallel do private(i,j)
2902  DO j=jsta,jend
2903  DO i=ista,iend
2904  grid1(i,j) = d3dsl(i,j,13)
2905  ENDDO
2906  ENDDO
2907  id(1:25)=0
2908  itd3d = nint(td3d)
2909  if (itd3d /= 0) then
2910  ifincr = mod(ifhr,itd3d)
2911  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2912  else
2913  ifincr = 0
2914  endif
2915  id(02)=133 ! Table 133
2916  id(18) = 0
2917  id(19) = ifhr
2918  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2919  id(20) = 3
2920  IF (ifincr == 0) THEN
2921  id(18) = ifhr-itd3d
2922  ELSE
2923  id(18) = ifhr-ifincr
2924  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2925  ENDIF
2926  if(grib == 'grib2')then
2927  cfld = cfld + 1
2928  fld_info(cfld)%ifld=iavblfld(iget(366))
2929  fld_info(cfld)%lvl=lvlsxml(lp,iget(366))
2930  if(itd3d==0) then
2931  fld_info(cfld)%ntrange=0
2932  else
2933  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2934  endif
2935  fld_info(cfld)%tinvstat=itd3d
2936 !$omp parallel do private(i,j,ii,jj)
2937  do j=1,jend-jsta+1
2938  jj = jsta+j-1
2939  do i=1,iend-ista+1
2940  ii=ista+i-1
2941  datapd(i,j,cfld) = grid1(ii,jj)
2942  enddo
2943  enddo
2944  endif
2945  ENDIF
2946  ENDIF
2947 !--- longwave tendency
2948  IF (iget(367) > 0) THEN
2949  IF (lvls(lp,iget(367)) > 0) THEN
2950 !$omp parallel do private(i,j)
2951  DO j=jsta,jend
2952  DO i=ista,iend
2953  grid1(i,j) = d3dsl(i,j,14)
2954  ENDDO
2955  ENDDO
2956  id(1:25)=0
2957  itd3d = nint(td3d)
2958  if (itd3d /= 0) then
2959  ifincr = mod(ifhr,itd3d)
2960  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
2961  else
2962  ifincr = 0
2963  endif
2964  id(02)=133 ! Table 133
2965  id(18) = 0
2966  id(19) = ifhr
2967  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
2968  id(20) = 3
2969  IF (ifincr == 0) THEN
2970  id(18) = ifhr-itd3d
2971  ELSE
2972  id(18) = ifhr-ifincr
2973  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
2974  ENDIF
2975  if(grib == 'grib2')then
2976  cfld = cfld + 1
2977  fld_info(cfld)%ifld=iavblfld(iget(367))
2978  fld_info(cfld)%lvl=lvlsxml(lp,iget(367))
2979  if(itd3d==0) then
2980  fld_info(cfld)%ntrange=0
2981  else
2982  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
2983  endif
2984  fld_info(cfld)%tinvstat=itd3d
2985 !$omp parallel do private(i,j,ii,jj)
2986  do j=1,jend-jsta+1
2987  jj = jsta+j-1
2988  do i=1,iend-ista+1
2989  ii=ista+i-1
2990  datapd(i,j,cfld) = grid1(ii,jj)
2991  enddo
2992  enddo
2993  endif
2994  ENDIF
2995  ENDIF
2996 !--- longwave tendency
2997  IF (iget(368) > 0) THEN
2998  IF (lvls(lp,iget(368)) > 0) THEN
2999 !$omp parallel do private(i,j)
3000  DO j=jsta,jend
3001  DO i=ista,iend
3002  grid1(i,j) = d3dsl(i,j,15)
3003  ENDDO
3004  ENDDO
3005  id(1:25)=0
3006  itd3d = nint(td3d)
3007  if (itd3d /= 0) then
3008  ifincr = mod(ifhr,itd3d)
3009  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3010  else
3011  ifincr = 0
3012  endif
3013  id(02)=133 ! Table 133
3014  id(18) = 0
3015  id(19) = ifhr
3016  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3017  id(20) = 3
3018  IF (ifincr == 0) THEN
3019  id(18) = ifhr-itd3d
3020  ELSE
3021  id(18) = ifhr-ifincr
3022  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3023  ENDIF
3024  if(grib == 'grib2')then
3025  cfld = cfld + 1
3026  fld_info(cfld)%ifld=iavblfld(iget(368))
3027  fld_info(cfld)%lvl=lvlsxml(lp,iget(368))
3028  if(itd3d==0) then
3029  fld_info(cfld)%ntrange=0
3030  else
3031  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3032  endif
3033  fld_info(cfld)%tinvstat=itd3d
3034 !$omp parallel do private(i,j,ii,jj)
3035  do j=1,jend-jsta+1
3036  jj = jsta+j-1
3037  do i=1,iend-ista+1
3038  ii=ista+i-1
3039  datapd(i,j,cfld) = grid1(ii,jj)
3040  enddo
3041  enddo
3042  endif
3043  ENDIF
3044  ENDIF
3045 !--- longwave tendency
3046  IF (iget(369) > 0) THEN
3047  IF (lvls(lp,iget(369)) > 0) THEN
3048 !$omp parallel do private(i,j)
3049  DO j=jsta,jend
3050  DO i=ista,iend
3051  grid1(i,j) = d3dsl(i,j,16)
3052  ENDDO
3053  ENDDO
3054  id(1:25)=0
3055  itd3d = nint(td3d)
3056  if (itd3d /= 0) then
3057  ifincr = mod(ifhr,itd3d)
3058  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3059  else
3060  ifincr = 0
3061  endif
3062  id(18) = 0
3063  id(19) = ifhr
3064  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3065  id(20) = 3
3066  IF (ifincr == 0) THEN
3067  id(18) = ifhr-itd3d
3068  ELSE
3069  id(18) = ifhr-ifincr
3070  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3071  ENDIF
3072  if(grib == 'grib2')then
3073  cfld = cfld + 1
3074  fld_info(cfld)%ifld=iavblfld(iget(369))
3075  fld_info(cfld)%lvl=lvlsxml(lp,iget(369))
3076  if(itd3d==0) then
3077  fld_info(cfld)%ntrange=0
3078  else
3079  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3080  endif
3081  fld_info(cfld)%tinvstat=itd3d
3082 !$omp parallel do private(i,j,ii,jj)
3083  do j=1,jend-jsta+1
3084  jj = jsta+j-1
3085  do i=1,iend-ista+1
3086  ii=ista+i-1
3087  datapd(i,j,cfld) = grid1(ii,jj)
3088  enddo
3089  enddo
3090  endif
3091  ENDIF
3092  ENDIF
3093 !--- longwave tendency
3094  IF (iget(370) > 0) THEN
3095  IF (lvls(lp,iget(370)) > 0) THEN
3096 !$omp parallel do private(i,j)
3097  DO j=jsta,jend
3098  DO i=ista,iend
3099  grid1(i,j) = d3dsl(i,j,17)
3100  ENDDO
3101  ENDDO
3102  id(1:25)=0
3103  itd3d = nint(td3d)
3104  if (itd3d /= 0) then
3105  ifincr = mod(ifhr,itd3d)
3106  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3107  else
3108  ifincr = 0
3109  endif
3110  id(02)=133 ! Table 133
3111  id(18) = 0
3112  id(19) = ifhr
3113  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3114  id(20) = 3
3115  IF (ifincr == 0) THEN
3116  id(18) = ifhr-itd3d
3117  ELSE
3118  id(18) = ifhr-ifincr
3119  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3120  ENDIF
3121  if(grib == 'grib2')then
3122  cfld = cfld + 1
3123  fld_info(cfld)%ifld=iavblfld(iget(370))
3124  fld_info(cfld)%lvl=lvlsxml(lp,iget(370))
3125  if(itd3d==0) then
3126  fld_info(cfld)%ntrange=0
3127  else
3128  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3129  endif
3130  fld_info(cfld)%tinvstat=itd3d
3131 !$omp parallel do private(i,j,ii,jj)
3132  do j=1,jend-jsta+1
3133  jj = jsta+j-1
3134  do i=1,iend-ista+1
3135  ii=ista+i-1
3136  datapd(i,j,cfld) = grid1(ii,jj)
3137  enddo
3138  enddo
3139  endif
3140  ENDIF
3141  ENDIF
3142 !--- longwave tendency
3143  IF (iget(371) > 0) THEN
3144  IF (lvls(lp,iget(371)) > 0) THEN
3145 !$omp parallel do private(i,j)
3146  DO j=jsta,jend
3147  DO i=ista,iend
3148  grid1(i,j) = d3dsl(i,j,18)
3149  ENDDO
3150  ENDDO
3151  id(1:25)=0
3152  itd3d = nint(td3d)
3153  if (itd3d /= 0) then
3154  ifincr = mod(ifhr,itd3d)
3155  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3156  else
3157  ifincr = 0
3158  endif
3159  id(02)=133 ! Table 133
3160  id(18) = 0
3161  id(19) = ifhr
3162  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3163  id(20) = 3
3164  IF (ifincr == 0) THEN
3165  id(18) = ifhr-itd3d
3166  ELSE
3167  id(18) = ifhr-ifincr
3168  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3169  ENDIF
3170  if(grib == 'grib2')then
3171  cfld = cfld + 1
3172  fld_info(cfld)%ifld=iavblfld(iget(371))
3173  fld_info(cfld)%lvl=lvlsxml(lp,iget(371))
3174  if(itd3d==0) then
3175  fld_info(cfld)%ntrange=0
3176  else
3177  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3178  endif
3179  fld_info(cfld)%tinvstat=itd3d
3180 !$omp parallel do private(i,j,ii,jj)
3181  do j=1,jend-jsta+1
3182  jj = jsta+j-1
3183  do i=1,iend-ista+1
3184  ii=ista+i-1
3185  datapd(i,j,cfld) = grid1(ii,jj)
3186  enddo
3187  enddo
3188  endif
3189  ENDIF
3190  ENDIF
3191 !--- longwave tendency
3192  IF (iget(372) > 0) THEN
3193  IF (lvls(lp,iget(372)) > 0) THEN
3194 !$omp parallel do private(i,j)
3195  DO j=jsta,jend
3196  DO i=ista,iend
3197  grid1(i,j) = d3dsl(i,j,19)
3198  ENDDO
3199  ENDDO
3200  id(1:25)=0
3201  itd3d = nint(td3d)
3202  if (itd3d /= 0) then
3203  ifincr = mod(ifhr,itd3d)
3204  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3205  else
3206  ifincr = 0
3207  endif
3208  id(18) = 0
3209  id(19) = ifhr
3210  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3211  id(20) = 3
3212  IF (ifincr == 0) THEN
3213  id(18) = ifhr-itd3d
3214  ELSE
3215  id(18) = ifhr-ifincr
3216  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3217  ENDIF
3218  if(grib == 'grib2')then
3219  cfld = cfld + 1
3220  fld_info(cfld)%ifld=iavblfld(iget(372))
3221  fld_info(cfld)%lvl=lvlsxml(lp,iget(372))
3222  if(itd3d==0) then
3223  fld_info(cfld)%ntrange=0
3224  else
3225  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3226  endif
3227  fld_info(cfld)%tinvstat=itd3d
3228 !$omp parallel do private(i,j,ii,jj)
3229  do j=1,jend-jsta+1
3230  jj = jsta+j-1
3231  do i=1,iend-ista+1
3232  ii=ista+i-1
3233  datapd(i,j,cfld) = grid1(ii,jj)
3234  enddo
3235  enddo
3236  endif
3237  ENDIF
3238  ENDIF
3239 !--- longwave tendency
3240  IF (iget(373) > 0) THEN
3241  IF (lvls(lp,iget(373)) > 0) THEN
3242 !$omp parallel do private(i,j)
3243  DO j=jsta,jend
3244  DO i=ista,iend
3245  grid1(i,j) = d3dsl(i,j,20)
3246  ENDDO
3247  ENDDO
3248  id(1:25)=0
3249  itd3d = nint(td3d)
3250  if (itd3d /= 0) then
3251  ifincr = mod(ifhr,itd3d)
3252  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3253  else
3254  ifincr = 0
3255  endif
3256  id(02)=133 ! Table 133
3257  id(18) = 0
3258  id(19) = ifhr
3259  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3260  id(20) = 3
3261  IF (ifincr == 0) THEN
3262  id(18) = ifhr-itd3d
3263  ELSE
3264  id(18) = ifhr-ifincr
3265  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3266  ENDIF
3267  if(grib == 'grib2')then
3268  cfld = cfld + 1
3269  fld_info(cfld)%ifld=iavblfld(iget(373))
3270  fld_info(cfld)%lvl=lvlsxml(lp,iget(373))
3271  if(itd3d==0) then
3272  fld_info(cfld)%ntrange=0
3273  else
3274  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3275  endif
3276  fld_info(cfld)%tinvstat=itd3d
3277 !$omp parallel do private(i,j,ii,jj)
3278  do j=1,jend-jsta+1
3279  jj = jsta+j-1
3280  do i=1,iend-ista+1
3281  ii=ista+i-1
3282  datapd(i,j,cfld) = grid1(ii,jj)
3283  enddo
3284  enddo
3285  endif
3286  ENDIF
3287  ENDIF
3288 !--- longwave tendency
3289  IF (iget(374) > 0) THEN
3290  IF (lvls(lp,iget(374)) > 0) THEN
3291 !$omp parallel do private(i,j)
3292  DO j=jsta,jend
3293  DO i=ista,iend
3294  grid1(i,j) = d3dsl(i,j,21)
3295  ENDDO
3296  ENDDO
3297  id(1:25)=0
3298  itd3d = nint(td3d)
3299  if (itd3d /= 0) then
3300  ifincr = mod(ifhr,itd3d)
3301  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3302  else
3303  ifincr = 0
3304  endif
3305  id(02)=133 ! Table 133
3306  id(18) = 0
3307  id(19) = ifhr
3308  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3309  id(20) = 3
3310  IF (ifincr == 0) THEN
3311  id(18) = ifhr-itd3d
3312  ELSE
3313  id(18) = ifhr-ifincr
3314  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3315  ENDIF
3316  if(grib == 'grib2')then
3317  cfld = cfld + 1
3318  fld_info(cfld)%ifld=iavblfld(iget(374))
3319  fld_info(cfld)%lvl=lvlsxml(lp,iget(374))
3320  if(itd3d==0) then
3321  fld_info(cfld)%ntrange=0
3322  else
3323  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3324  endif
3325  fld_info(cfld)%tinvstat=itd3d
3326 !$omp parallel do private(i,j,ii,jj)
3327  do j=1,jend-jsta+1
3328  jj = jsta+j-1
3329  do i=1,iend-ista+1
3330  ii=ista+i-1
3331  datapd(i,j,cfld) = grid1(ii,jj)
3332  enddo
3333  enddo
3334  endif
3335  ENDIF
3336  ENDIF
3337 !--- longwave tendency
3338  IF (iget(375) > 0) THEN
3339  IF (lvls(lp,iget(375)) > 0) THEN
3340 !$omp parallel do private(i,j)
3341  DO j=jsta,jend
3342  DO i=ista,iend
3343  grid1(i,j) = d3dsl(i,j,22)
3344  ENDDO
3345  ENDDO
3346  id(1:25)=0
3347  itd3d = nint(td3d)
3348  if (itd3d /= 0) then
3349  ifincr = mod(ifhr,itd3d)
3350  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3351  else
3352  ifincr = 0
3353  endif
3354  id(18) = 0
3355  id(19) = ifhr
3356  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3357  id(20) = 3
3358  IF (ifincr == 0) THEN
3359  id(18) = ifhr-itd3d
3360  ELSE
3361  id(18) = ifhr-ifincr
3362  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3363  ENDIF
3364  if(grib == 'grib2') then
3365  cfld = cfld + 1
3366  fld_info(cfld)%ifld=iavblfld(iget(375))
3367  fld_info(cfld)%lvl=lvlsxml(lp,iget(375))
3368  if(itd3d==0) then
3369  fld_info(cfld)%ntrange=0
3370  else
3371  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3372  endif
3373  fld_info(cfld)%tinvstat=itd3d
3374 !$omp parallel do private(i,j,ii,jj)
3375  do j=1,jend-jsta+1
3376  jj = jsta+j-1
3377  do i=1,iend-ista+1
3378  ii=ista+i-1
3379  datapd(i,j,cfld) = grid1(ii,jj)
3380  enddo
3381  enddo
3382  endif
3383  ENDIF
3384  ENDIF
3385 !--- total diabatic heating
3386  IF (iget(379) > 0) THEN
3387  IF (lvls(lp,iget(379)) > 0) THEN
3388 !$omp parallel do private(i,j)
3389  DO j=jsta,jend
3390  DO i=ista,iend
3391  IF(d3dsl(i,j,1)/=spval)THEN
3392  grid1(i,j) = d3dsl(i,j,1) + d3dsl(i,j,2) &
3393  + d3dsl(i,j,3) + d3dsl(i,j,4) &
3394  + d3dsl(i,j,5) + d3dsl(i,j,6)
3395  ELSE
3396  grid1(i,j) = spval
3397  END IF
3398  ENDDO
3399  ENDDO
3400  id(1:25)=0
3401  itd3d = nint(td3d)
3402  if (itd3d /= 0) then
3403  ifincr = mod(ifhr,itd3d)
3404  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3405  else
3406  ifincr = 0
3407  endif
3408  id(18) = 0
3409  id(19) = ifhr
3410  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3411  id(20) = 3
3412  IF (ifincr == 0) THEN
3413  id(18) = ifhr-itd3d
3414  ELSE
3415  id(18) = ifhr-ifincr
3416  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3417  ENDIF
3418  if(grib == 'grib2') then
3419  cfld = cfld + 1
3420  fld_info(cfld)%ifld=iavblfld(iget(379))
3421  fld_info(cfld)%lvl=lvlsxml(lp,iget(379))
3422  if(itd3d==0) then
3423  fld_info(cfld)%ntrange=0
3424  else
3425  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3426  endif
3427  fld_info(cfld)%tinvstat=itd3d
3428 !$omp parallel do private(i,j,ii,jj)
3429  do j=1,jend-jsta+1
3430  jj = jsta+j-1
3431  do i=1,iend-ista+1
3432  ii=ista+i-1
3433  datapd(i,j,cfld) = grid1(ii,jj)
3434  enddo
3435  enddo
3436  endif
3437  ENDIF
3438  ENDIF
3439 !--- convective updraft
3440  IF (iget(391) > 0) THEN
3441  IF (lvls(lp,iget(391)) > 0) THEN
3442 !$omp parallel do private(i,j)
3443  DO j=jsta,jend
3444  DO i=ista,iend
3445  grid1(i,j) = d3dsl(i,j,23)
3446  ENDDO
3447  ENDDO
3448  id(1:25)=0
3449  itd3d = nint(td3d)
3450  if (itd3d /= 0) then
3451  ifincr = mod(ifhr,itd3d)
3452  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3453  else
3454  ifincr = 0
3455  endif
3456  id(02)=133 ! Table 133
3457  id(18) = 0
3458  id(19) = ifhr
3459  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3460  id(20) = 3
3461  IF (ifincr == 0) THEN
3462  id(18) = ifhr-itd3d
3463  ELSE
3464  id(18) = ifhr-ifincr
3465  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3466  ENDIF
3467  if(grib == 'grib2') then
3468  cfld = cfld + 1
3469  fld_info(cfld)%ifld=iavblfld(iget(391))
3470  fld_info(cfld)%lvl=lvlsxml(lp,iget(391))
3471  if(itd3d==0) then
3472  fld_info(cfld)%ntrange=0
3473  else
3474  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3475  endif
3476  fld_info(cfld)%tinvstat=itd3d
3477 !$omp parallel do private(i,j,ii,jj)
3478  do j=1,jend-jsta+1
3479  jj = jsta+j-1
3480  do i=1,iend-ista+1
3481  ii=ista+i-1
3482  datapd(i,j,cfld) = grid1(ii,jj)
3483  enddo
3484  enddo
3485  endif
3486  ENDIF
3487  ENDIF
3488 !--- convective downdraft
3489  IF (iget(392) > 0) THEN
3490  IF (lvls(lp,iget(392)) > 0) THEN
3491 !$omp parallel do private(i,j)
3492  DO j=jsta,jend
3493  DO i=ista,iend
3494  grid1(i,j) = d3dsl(i,j,24)
3495  ENDDO
3496  ENDDO
3497  id(1:25)=0
3498  itd3d = nint(td3d)
3499  if (itd3d /= 0) then
3500  ifincr = mod(ifhr,itd3d)
3501  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3502  else
3503  ifincr = 0
3504  endif
3505  id(02)=133 ! Table 133
3506  id(18) = 0
3507  id(19) = ifhr
3508  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3509  id(20) = 3
3510  IF (ifincr == 0) THEN
3511  id(18) = ifhr-itd3d
3512  ELSE
3513  id(18) = ifhr-ifincr
3514  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3515  ENDIF
3516  if(grib == 'grib2') then
3517  cfld = cfld + 1
3518  fld_info(cfld)%ifld=iavblfld(iget(392))
3519  fld_info(cfld)%lvl=lvlsxml(lp,iget(392))
3520  if(itd3d==0) then
3521  fld_info(cfld)%ntrange=0
3522  else
3523  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3524  endif
3525  fld_info(cfld)%tinvstat=itd3d
3526 !$omp parallel do private(i,j,ii,jj)
3527  do j=1,jend-jsta+1
3528  jj = jsta+j-1
3529  do i=1,iend-ista+1
3530  ii=ista+i-1
3531  datapd(i,j,cfld) = grid1(ii,jj)
3532  enddo
3533  enddo
3534  endif
3535  ENDIF
3536  ENDIF
3537 !--- convective detraintment
3538  IF (iget(393) > 0) THEN
3539  IF (lvls(lp,iget(393)) > 0) THEN
3540 !$omp parallel do private(i,j)
3541  DO j=jsta,jend
3542  DO i=ista,iend
3543  grid1(i,j) = d3dsl(i,j,25)
3544  ENDDO
3545  ENDDO
3546  id(1:25)=0
3547  itd3d = nint(td3d)
3548  if (itd3d /= 0) then
3549  ifincr = mod(ifhr,itd3d)
3550  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3551  else
3552  ifincr = 0
3553  endif
3554  id(02)=133 ! Table 133
3555  id(18) = 0
3556  id(19) = ifhr
3557  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3558  id(20) = 3
3559  IF (ifincr == 0) THEN
3560  id(18) = ifhr-itd3d
3561  ELSE
3562  id(18) = ifhr-ifincr
3563  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3564  ENDIF
3565  if(grib == 'grib2') then
3566  cfld = cfld + 1
3567  fld_info(cfld)%ifld=iavblfld(iget(393))
3568  fld_info(cfld)%lvl=lvlsxml(lp,iget(393))
3569  if(itd3d==0) then
3570  fld_info(cfld)%ntrange=0
3571  else
3572  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3573  endif
3574  fld_info(cfld)%tinvstat=itd3d
3575 !$omp parallel do private(i,j,ii,jj)
3576  do j=1,jend-jsta+1
3577  jj = jsta+j-1
3578  do i=1,iend-ista+1
3579  ii=ista+i-1
3580  datapd(i,j,cfld) = grid1(ii,jj)
3581  enddo
3582  enddo
3583  endif
3584  ENDIF
3585  ENDIF
3586 !--- convective gravity drag zonal acce
3587  IF (iget(394) > 0) THEN
3588  IF (lvls(lp,iget(394)) > 0) THEN
3589 !$omp parallel do private(i,j)
3590  DO j=jsta,jend
3591  DO i=ista,iend
3592  grid1(i,j) = d3dsl(i,j,26)
3593  ENDDO
3594  ENDDO
3595  id(1:25)=0
3596  itd3d = nint(td3d)
3597  if (itd3d /= 0) then
3598  ifincr = mod(ifhr,itd3d)
3599  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3600  else
3601  ifincr = 0
3602  endif
3603  id(02)=133 ! Table 133
3604  id(18) = 0
3605  id(19) = ifhr
3606  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3607  id(20) = 3
3608  IF (ifincr == 0) THEN
3609  id(18) = ifhr-itd3d
3610  ELSE
3611  id(18) = ifhr-ifincr
3612  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3613  ENDIF
3614  if(grib == 'grib2') then
3615  cfld = cfld + 1
3616  fld_info(cfld)%ifld=iavblfld(iget(394))
3617  fld_info(cfld)%lvl=lvlsxml(lp,iget(394))
3618  if(itd3d==0) then
3619  fld_info(cfld)%ntrange=0
3620  else
3621  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3622  endif
3623  fld_info(cfld)%tinvstat=itd3d
3624 !$omp parallel do private(i,j,ii,jj)
3625  do j=1,jend-jsta+1
3626  jj = jsta+j-1
3627  do i=1,iend-ista+1
3628  ii=ista+i-1
3629  datapd(i,j,cfld) = grid1(ii,jj)
3630  enddo
3631  enddo
3632  endif
3633  ENDIF
3634  ENDIF
3635 !--- convective gravity drag meridional acce
3636  IF (iget(395) > 0) THEN
3637  IF (lvls(lp,iget(395)) > 0) THEN
3638 !$omp parallel do private(i,j)
3639  DO j=jsta,jend
3640  DO i=ista,iend
3641  grid1(i,j) = d3dsl(i,j,27)
3642  ENDDO
3643  ENDDO
3644  id(1:25)=0
3645  itd3d = nint(td3d)
3646  if (itd3d /= 0) then
3647  ifincr = mod(ifhr,itd3d)
3648  IF(ifmin >= 1)ifincr= mod(ifhr*60+ifmin,itd3d*60)
3649  else
3650  ifincr = 0
3651  endif
3652  id(02)=133 ! Table 133
3653  id(18) = 0
3654  id(19) = ifhr
3655  IF(ifmin >= 1)id(19)=ifhr*60+ifmin
3656  id(20) = 3
3657  IF (ifincr == 0) THEN
3658  id(18) = ifhr-itd3d
3659  ELSE
3660  id(18) = ifhr-ifincr
3661  IF(ifmin >= 1)id(18)=ifhr*60+ifmin-ifincr
3662  ENDIF
3663  if(grib == 'grib2') then
3664  cfld = cfld + 1
3665  fld_info(cfld)%ifld=iavblfld(iget(395))
3666  fld_info(cfld)%lvl=lvlsxml(lp,iget(395))
3667  if(itd3d==0) then
3668  fld_info(cfld)%ntrange=0
3669  else
3670  fld_info(cfld)%ntrange=(ifhr-id(18))/itd3d
3671  endif
3672  fld_info(cfld)%tinvstat=itd3d
3673 !$omp parallel do private(i,j,ii,jj)
3674  do j=1,jend-jsta+1
3675  jj = jsta+j-1
3676  do i=1,iend-ista+1
3677  ii=ista+i-1
3678  datapd(i,j,cfld) = grid1(ii,jj)
3679  enddo
3680  enddo
3681  endif
3682  ENDIF
3683  ENDIF
3684  END IF ! end of d3d output
3685 
3686 ! CHUANG: COMPUTE HAINES INDEX
3687  IF (iget(455) > 0) THEN
3688  ii=(ista+iend)/2+100
3689  jj=(jsta+jend)/2-100
3690  IF(abs(spl(lp)-50000.)<small) luhi=lp
3691  IF(abs(spl(lp)-70000.)<small) THEN ! high evevation
3692 ! HAINES=SPVAL
3693 ! print*,'computing dew point for Haine Index at ',SPL(LP)
3694 !$omp parallel do private(i,j)
3695  DO j=jsta,jend
3696  DO i=ista,iend
3697  haines(i,j) = spval
3698  egrid2(i,j) = spl(lp)
3699  ENDDO
3700  ENDDO
3701  CALL caldwp(egrid2(ista:iend,jsta:jend),qsl(ista:iend,jsta:jend),tdsl(ista:iend,jsta:jend),tsl(ista:iend,jsta:jend))
3702 
3703 !$omp parallel do private(i,j,dum1,istaa,imois)
3704  DO j=jsta,jend
3705  DO i=ista,iend
3706  IF(sm(i,j) < 1.0 .AND. zint(i,j,lm+1) < fsl(i,j)*gi) THEN
3707  dum1 = tsl(i,j)-tprs(i,j,luhi)
3708  IF(dum1 <= 17.)THEN
3709  istaa = 1
3710  ELSE IF(dum1 > 17. .AND. dum1 <= 21.) THEN
3711  istaa = 2
3712  ELSE
3713  istaa = 3
3714  END IF
3715  dum1 = tsl(i,j)-tdsl(i,j)
3716  IF(dum1 <= 14.) THEN
3717  imois = 1
3718  ELSE IF(dum1>14. .AND. dum1<=20.) THEN
3719  imois = 2
3720  ELSE
3721  imois = 3
3722  END IF
3723  IF(tsl(i,j)<spval.and.tprs(i,j,luhi)<spval.and.tdsl(i,j)<spval)THEN
3724  haines(i,j) = istaa + imois
3725  ELSE
3726  haines(i,j) = spval
3727  ENDIF
3728 ! if(i==570 .and. j==574)print*,'high hainesindex:',i,j,luhi,tsl(i,j) &
3729 ! ,tprs(i,j,luhi),tdsl(i,j),ista,imois,spl(luhi),spl(lp),haines(i,j)
3730  END IF
3731  END DO
3732  END DO
3733 
3734  luhi = lp
3735  ENDIF
3736 
3737  IF(abs(spl(lp)-85000.)<small)THEN ! mid evevation
3738 ! print*,'computing dew point for Haine Index at ',SPL(LP)
3739 !$omp parallel do private(i,j)
3740  DO j=jsta,jend
3741  DO i=ista,iend
3742  egrid2(i,j) = spl(lp)
3743  ENDDO
3744  ENDDO
3745  CALL caldwp(egrid2(ista:iend,jsta:jend),qsl(ista:iend,jsta:jend),tdsl(ista:iend,jsta:jend),tsl(ista:iend,jsta:jend))
3746 
3747 !$omp parallel do private(i,j,dum1,istaa,imois)
3748  DO j=jsta,jend
3749  DO i=ista,iend
3750  IF(sm(i,j) < 1.0 .AND. zint(i,j,lm+1) < fsl(i,j)*gi) THEN
3751  dum1 = tsl(i,j)-tprs(i,j,luhi)
3752  IF(dum1 <=5. ) THEN
3753  istaa = 1
3754  ELSE IF(dum1 > 5. .AND. dum1 <= 10.) THEN
3755  istaa = 2
3756  ELSE
3757  istaa = 3
3758  END IF
3759  dum1 = tsl(i,j)-tdsl(i,j)
3760  IF(dum1 <= 5.) THEN
3761  imois = 1
3762  ELSE IF(dum1 > 5. .AND. dum1 <= 12.) THEN
3763  imois = 2
3764  ELSE
3765  imois = 3
3766  END IF
3767 ! if(i==570 .and. j==574)print*,'mid haines index:',i,j,luhi,tsl(i,j) &
3768 ! ,tprs(i,j,luhi),tdsl(i,j),ista,imois,spl(luhi),spl(lp),haines(i,j)
3769  IF(tsl(i,j)<spval.and.tprs(i,j,luhi)<spval.and.tdsl(i,j)<spval)THEN
3770  haines(i,j) = istaa + imois
3771  ELSE
3772  haines(i,j) = spval
3773  ENDIF
3774  END IF
3775  END DO
3776  END DO
3777 
3778  luhi = lp
3779  ENDIF
3780 
3781  IF(abs(spl(lp)-95000.)<small)THEN ! LOW evevation
3782 ! print*,'computing dew point for Haine Index at ',SPL(LP)
3783 !$omp parallel do private(i,j)
3784  DO j=jsta,jend
3785  DO i=ista,iend
3786  egrid2(i,j)=spl(lp)
3787  ENDDO
3788  ENDDO
3789  CALL caldwp(egrid2(ista:iend,jsta:jend),qsl(ista:iend,jsta:jend),tdsl(ista:iend,jsta:jend),tsl(ista:iend,jsta:jend))
3790 
3791 !$omp parallel do private(i,j,dum1,istaa,imois)
3792  DO j=jsta,jend
3793  DO i=ista,iend
3794  IF(sm(i,j) < 1.0 .AND. zint(i,j,lm+1) < fsl(i,j)*gi) THEN
3795  dum1 = tsl(i,j)-tprs(i,j,luhi)
3796  IF(dum1 <= 3.)THEN
3797  istaa = 1
3798  ELSE IF(dum1 > 3. .AND. dum1 <=7. ) THEN
3799  istaa = 2
3800  ELSE
3801  istaa = 3
3802  END IF
3803  dum1 = tsl(i,j)-tdsl(i,j)
3804  IF(dum1 <=5. ) THEN
3805  imois = 1
3806  ELSE IF(dum1 > 5. .AND. dum1 <= 9.)THEN
3807  imois = 2
3808  ELSE
3809  imois = 3
3810  END IF
3811 ! if(i==570 .and. j==574)print*,'low haines index:',i,j,luhi,tsl(i,j) &
3812 ! ,tprs(i,j,luhi),tdsl(i,j),ista,imois,spl(luhi),spl(lp),haines(i,j)
3813  IF(tsl(i,j)<spval.and.tprs(i,j,luhi)<spval.and.tdsl(i,j)<spval)THEN
3814  haines(i,j) = istaa + imois
3815  ELSE
3816  haines(i,j) = spval
3817  ENDIF
3818  END IF
3819  END DO
3820  END DO
3821 
3822  if(grib == 'grib2') then
3823  cfld = cfld + 1
3824  fld_info(cfld)%ifld=iavblfld(iget(455))
3825 !$omp parallel do private(i,j,ii,jj)
3826  do j=1,jend-jsta+1
3827  jj = jsta+j-1
3828  do i=1,iend-ista+1
3829  ii=ista+i-1
3830  datapd(i,j,cfld) = haines(ii,jj)
3831  enddo
3832  enddo
3833  endif
3834 
3835  ENDIF
3836 
3837  ENDIF
3838 !
3839  ENDDO !*** END OF MAIN VERTICAL LOOP DO LP=1,LSM
3840 !*** ENDIF FOR IF TEST SEEING IF WE WANT ANY OTHER VARIABLES
3841  ENDIF
3842 ! SRD
3843 !
3844 ! MAX VERTICAL VELOCITY UPDRAFT
3845 !
3846  IF (iget(423) > 0) THEN
3847 ! LP=22 ! 400 MB
3848  lp=46 ! 1000 MB
3849 !$omp parallel do private(i,j)
3850  DO j=jsta,jend
3851  DO i=ista,iend
3852  grid1(i,j) = w_up_max(i,j)
3853 ! print *,' writing w_up_max, i,j, = ', w_up_max(i,j)
3854  ENDDO
3855  ENDDO
3856  if(grib == 'grib2')then
3857  cfld = cfld + 1
3858  fld_info(cfld)%ifld = iavblfld(iget(423))
3859  fld_info(cfld)%lvl = lvlsxml(lp,iget(423))
3860  if (ifhr > 0) then
3861  fld_info(cfld)%tinvstat=1
3862  else
3863  fld_info(cfld)%tinvstat=0
3864  endif
3865  fld_info(cfld)%ntrange=1
3866 !$omp parallel do private(i,j,ii,jj)
3867  do j=1,jend-jsta+1
3868  jj = jsta+j-1
3869  do i=1,iend-ista+1
3870  ii=ista+i-1
3871  datapd(i,j,cfld) = grid1(ii,jj)
3872  enddo
3873  enddo
3874  endif
3875  ENDIF
3876 !
3877 ! MAX VERTICAL VELOCITY DOWNDRAFT
3878 !
3879  IF (iget(424) > 0) THEN
3880  lp = 46 ! 1000 MB
3881 !$omp parallel do private(i,j)
3882  DO j=jsta,jend
3883  DO i=ista,iend
3884  grid1(i,j) = w_dn_max(i,j)
3885  ENDDO
3886  ENDDO
3887  if(grib == 'grib2')then
3888  cfld = cfld + 1
3889  fld_info(cfld)%ifld=iavblfld(iget(424))
3890  fld_info(cfld)%lvl=lvlsxml(lp,iget(424))
3891  if (ifhr > 0) then
3892  fld_info(cfld)%tinvstat=1
3893  else
3894  fld_info(cfld)%tinvstat=0
3895  endif
3896  fld_info(cfld)%ntrange=1
3897 !$omp parallel do private(i,j,ii,jj)
3898  do j=1,jend-jsta+1
3899  jj = jsta+j-1
3900  do i=1,iend-ista+1
3901  ii=ista+i-1
3902  datapd(i,j,cfld) = grid1(ii,jj)
3903  enddo
3904  enddo
3905  endif
3906  ENDIF
3907 !
3908 ! MEAN VERTICAL VELOCITY
3909 !
3910 ! This hourly mean field is, in fact, based on the column
3911 ! mean bounded by sigma levels, but I chose to keep the
3912 ! output here with the other diagnostic vertical
3913 ! velocity fields
3914 !
3915  IF (iget(425) > 0) THEN
3916  lp = 46 ! 1000 MB
3917 !$omp parallel do private(i,j)
3918  DO j=jsta,jend
3919  DO i=ista,iend
3920  grid1(i,j) = w_mean(i,j)
3921  ENDDO
3922  ENDDO
3923  if(grib == 'grib2')then
3924  cfld = cfld + 1
3925  fld_info(cfld)%ifld = iavblfld(iget(425))
3926  fld_info(cfld)%lvl = lvlsxml(lp,iget(425))
3927  if (ifhr == 0) then
3928  fld_info(cfld)%tinvstat = 0
3929  else
3930  fld_info(cfld)%tinvstat = 1
3931  endif
3932  fld_info(cfld)%ntrange = 1
3933 !$omp parallel do private(i,j,ii,jj)
3934  do j=1,jend-jsta+1
3935  jj = jsta+j-1
3936  do i=1,iend-ista+1
3937  ii=ista+i-1
3938  datapd(i,j,cfld) = grid1(ii,jj)
3939  enddo
3940  enddo
3941  endif
3942  ENDIF
3943 ! SRD
3944 
3945 !
3946 ! CALL MEMBRANE SLP REDUCTION IF REQESTED IN CONTROL FILE
3947 !
3948 ! OUTPUT MEMBRANCE SLP
3949  IF(iget(023) > 0)THEN
3950  IF(gridtype == 'A'.OR. gridtype == 'B') then
3951  if(me==0)print*,'CALLING MEMSLP for A or B grid'
3952  CALL memslp(tprs,qprs,fprs)
3953  if(me==0)print*,'aft CALLING MEMSLP for A or B grid,pslp=', &
3954  maxval(pslp(ista:iend,jsta:jend)),minval(pslp(ista:iend,jsta:jend)),pslp((ista+iend)/2,(jsta+jend)/2)
3955  ELSE IF (gridtype == 'E')THEN
3956  if(me==0)print*,'CALLING MEMSLP_NMM for E grid'
3957 ! CALL MEMSLP_NMM(TPRS,QPRS,FPRS)
3958  ELSE
3959  print*,'unknow grid type-> WONT DERIVE MESINGER SLP'
3960  END IF
3961 !$omp parallel do private(i,j)
3962  DO j=jsta,jend
3963  DO i=ista,iend
3964  grid1(i,j) = pslp(i,j)
3965  ENDDO
3966  ENDDO
3967 ! print *,'inmdl2p,pslp=',maxval(pslp(1:im,jsta:jend)),minval(pslp(1:im,jsta:jend))
3968 ! print *,'inmdl2p,point pslp=',pslp(im/2,(jsta+jend)/2),pslp(1,jsta),'cfld=',cfld
3969  if(grib == 'grib2')then
3970  cfld = cfld + 1
3971  fld_info(cfld)%ifld = iavblfld(iget(023))
3972 !$omp parallel do private(i,j,ii,jj)
3973  do j=1,jend-jsta+1
3974  jj = jsta+j-1
3975  do i=1,iend-ista+1
3976  ii=ista+i-1
3977  datapd(i,j,cfld) = grid1(ii,jj)
3978  enddo
3979  enddo
3980  endif
3981  ENDIF
3982 
3983 ! OUTPUT of MAPS SLP
3984  IF(iget(445) > 0)THEN
3985  if(me==0)print*,'CALLING MAPS SLP'
3986  CALL mapsslp(tprs)
3987 !$omp parallel do private(i,j)
3988  DO j=jsta,jend
3989  DO i=ista,iend
3990  grid1(i,j) = pslp(i,j)
3991  ENDDO
3992  ENDDO
3993  if(grib == 'grib2') then
3994  cfld = cfld + 1
3995  fld_info(cfld)%ifld = iavblfld(iget(445))
3996 !$omp parallel do private(i,j,ii,jj)
3997  do j=1,jend-jsta+1
3998  jj = jsta+j-1
3999  do i=1,iend-ista+1
4000  ii=ista+i-1
4001  datapd(i,j,cfld) = grid1(ii,jj)
4002  enddo
4003  enddo
4004  endif
4005  ENDIF
4006 
4007 
4008 ! ADJUST 1000 MB HEIGHT TO MEMBEANCE SLP
4009  IF(iget(023) > 0.OR.iget(445) > 0)THEN
4010  IF(iget(012) > 0)THEN
4011 ! dong add missing value to 1000 mb hgt
4012  grid1= spval
4013  DO lp=lsm,1,-1
4014  IF(abs(spl(lp)-1.0e5) <= 1.0e-5)THEN
4015  IF(lvls(lp,iget(012)) > 0)THEN
4016  alpth = log(1.e5)
4017  IF(modelname == 'GFS')THEN
4018 ! GFS does not want to adjust 1000 mb H to membrane SLP
4019 ! because MOS can't adjust to the much lower H
4020 !$omp parallel do private(i,j)
4021  DO j=jsta,jend
4022  DO i=ista,iend
4023  IF(fsl(i,j)<spval)THEN
4024  grid1(i,j) = fsl(i,j)*gi
4025  ELSE
4026  grid1(i,j) = spval
4027  ENDIF
4028  ENDDO
4029  ENDDO
4030  ELSE
4031 !$omp parallel do private(i,j,PSLPIJ,ALPSL,PSFC)
4032  DO j=jsta,jend
4033  DO i=ista,iend
4034  IF(pslp(i,j) < spval) THEN
4035  pslpij = pslp(i,j)
4036  alpsl = log(pslpij)
4037  psfc = pint(i,j,nint(lmh(i,j))+1)
4038  IF(abs(pslpij-psfc) < 5.e2) THEN
4039  grid1(i,j) = rd*tprs(i,j,lp)*(alpsl-alpth)
4040  ELSE
4041  grid1(i,j) = fis(i,j)/(alpsl-log(psfc))*(alpsl-alpth)
4042  ENDIF
4043  z1000(i,j) = grid1(i,j)*gi
4044  grid1(i,j) = z1000(i,j)
4045  ELSE
4046  grid1(i,j) = spval
4047  END IF
4048  ENDDO
4049  ENDDO
4050  END IF
4051 
4052  IF (smflag) THEN
4053 !tgs - smoothing of geopotential heights
4054  nsmooth = nint(5.*(13500./dxm))
4055  call allgetherv(grid1)
4056  do k=1,nsmooth
4057  CALL smooth(grid1,sdummy,im,jm,0.5)
4058  end do
4059  ENDIF
4060 
4061  if(grib == 'grib2') then
4062  cfld = cfld + 1
4063  fld_info(cfld)%ifld = iavblfld(iget(012))
4064  fld_info(cfld)%lvl = lvlsxml(lp,iget(012))
4065 !$omp parallel do private(i,j,ii,jj)
4066  do j=1,jend-jsta+1
4067  jj = jsta+j-1
4068  do i=1,iend-ista+1
4069  ii=ista+i-1
4070  datapd(i,j,cfld) = grid1(ii,jj)
4071  enddo
4072  enddo
4073  endif
4074  exit
4075  ENDIF
4076  ENDIF
4077  END DO
4078  ENDIF
4079  ENDIF
4080 
4081 ! SNOW DESITY SOLID-LIQUID-RATION SLR
4082  IF ( iget(1006)>0 ) THEN
4083  if(me==0)print*,'CALLING SLR'
4084  egrid1=spval
4085  if(slrutah_on) then
4086  call calslr_uutah(egrid1)
4087  else
4088  call calslr_roebber(tprs,rhprs,egrid1)
4089  endif
4090 !$omp parallel do private(i,j)
4091  do j=jsta,jend
4092  do i=ista,iend
4093  grid1(i,j)=spval
4094  if(egrid1(i,j) < spval) then
4095  if(egrid1(i,j)>=1.) then
4096  grid1(i,j)=1000./egrid1(i,j)
4097  else
4098  grid1(i,j)=spval
4099  endif
4100  endif
4101  enddo
4102  enddo
4103  if(grib=='grib2') then
4104  cfld=cfld+1
4105  fld_info(cfld)%ifld=iavblfld(iget(1006))
4106 !$omp parallel do private(i,j,ii,jj)
4107  do j=1,jend-jsta+1
4108  jj = jsta+j-1
4109  do i=1,iend-ista+1
4110  ii=ista+i-1
4111  datapd(i,j,cfld) = grid1(ii,jj)
4112  enddo
4113  enddo
4114  endif
4115  ENDIF
4116 !
4117 if(allocated(d3dsl)) deallocate(d3dsl)
4118 if(allocated(smokesl)) deallocate(smokesl)
4119 if(allocated(fv3dustsl)) deallocate(fv3dustsl)
4120 if(allocated(coarsepmsl)) deallocate(coarsepmsl)
4121  if(me==0)print*,'MDL2P completed'
4122 ! END OF ROUTINE.
4123 !
4124  RETURN
4125  END
ctlblk_mod
Definition: CTLBLK.f:1
calstrm
subroutine calstrm(Z1D, STRM)
Subroutine that computes geo streamfunction.
Definition: CALSTRM.f:31
caldwp
subroutine caldwp(P1D, Q1D, TDWP, T1D)
Computes dewpoint from P, T, and Q.
Definition: CALDWP.f:21
rqstfld_mod
Definition: RQSTFLD.f:1
upp_physics::calslr_roebber
subroutine, public calslr_roebber(tprs, rhprs, slr)
calslr_roebber() computes snow solid-liquid-ratio slr using the Roebber algorithm.
Definition: UPP_PHYSICS.f:2705
upp_physics::calslr_uutah
subroutine, public calslr_uutah(SLR)
calslr_uutah() computes snow solid-liquid-ratio slr using the Steenburgh algorithm.
Definition: UPP_PHYSICS.f:4374
masks
Definition: MASKS_mod.f:1
smooth
subroutine smooth(FIELD, HOLD, IX, IY, SMTH)
smooth() smooths a meteorological field using Shapiro smoother.
Definition: SMOOTH.f:42
physcons_post
Definition: physcons.f:1
params_mod
Definition: params.F:1
upp_physics::calvor
subroutine, public calvor(UWND, VWND, ABSV)
CALVOR() computes absolute vorticity.
Definition: UPP_PHYSICS.f:1744
bound
subroutine bound(FLD, FMIN, FMAX)
This routine bounds data in the passed array FLD (im x jm elements long) and clips data values such t...
Definition: BOUND.f:29
calcat
subroutine calcat(U, V, H, U_OLD, V_OLD, H_OLD, CAT)
Computes Clear Air Turbulence Index.
Definition: AVIATION.f:264
calmcvg
subroutine calmcvg(Q1D, U1D, V1D, QCNVG)
Subroutine that computes moisture convergence.
Definition: CALMCVG.f:43
upp_physics::calrh
subroutine, public calrh(P1, T1, Q1, RH)
CALRH() computes relative humidity.
Definition: UPP_PHYSICS.f:72
smoothc
subroutine smoothc(FIELD, HOLD, IX, IY, SMTH)
smoothc() smooths a meteorological field using Shapiro smoother.
Definition: SMOOTH.f:134
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:378
calicing
subroutine calicing(T1, RH, OMGA, ICING)
Computes In-Flight Icing.
Definition: AVIATION.f:160
upp_physics
Definition: UPP_PHYSICS.f:40