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