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