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