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