UPP (develop)
Loading...
Searching...
No Matches
MDL2SIGMA.f
Go to the documentation of this file.
1
2!
54 SUBROUTINE mdl2sigma
55
56!
57!
58 use vrbls3d, only: pint, t, q, zint, alpint, pmid, exch_h, uh, &
59 vh, omga, q2, cwm, qqw, qqi, qqr, qqs, cfr, &
60 f_rimef, pmidv
61! use vrbls2d, only:
62 use masks, only: lmh
63 use params_mod, only: d50 , pq0, a2, a3, a4, h1, d01, d608, rgamog,&
64 h1m12, d00, h2, rd, g, gi, h99999
65 use ctlblk_mod, only: jsta_2l, jend_2u, spval, lp1, jsta, jend, lm, &
66 grib, cfld, datapd, fld_info, me, jend_m, im, &
67 jm, im_jm, ista, iend, ista_2l, iend_2u, ista_m, iend_m
68 use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml
69 use gridspec_mod, only :gridtype
70!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 implicit none
72!
73! INCLUDE MODEL DIMENSIONS. SET/DERIVE OTHER PARAMETERS.
74! GAMMA AND RGAMOG ARE USED IN THE EXTRAPOLATION OF VIRTUAL
75! TEMPERATURES BEYOND THE UPPER OF LOWER LIMITS OF DATA.
76!
77 integer,PARAMETER :: LSIG=22
78 real,PARAMETER :: PTSIGO=1.0e4
79!
80! DECLARE VARIABLES.
81!
82 LOGICAL READTHK
83 LOGICAL IOOMG,IOALL
84 LOGICAL DONEFSL1,TSLDONE
85 real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: FSL, TSL, QSL, OSL, USL, VSL, Q2SL, &
86 fsl1, cfrsig, egrid1, egrid2
87 REAL GRID1(IM,JM)
88 real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: grid2
89
90 REAL SIGO(LSIG+1),DSIGO(LSIG),ASIGO(LSIG)
91!
92 INTEGER,dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: NL1X,NL1XF
93!
94!
95!--- Definition of the following 2D (horizontal) dummy variables
96!
97! C1D - total condensate
98! QW1 - cloud water mixing ratio
99! QI1 - cloud ice mixing ratio
100! QR1 - rain mixing ratio
101! QS1 - snow mixing ratio
102!
103 real, dimension(ista_2l:iend_2u,jsta_2l:jend_2u) :: C1D, QW1, QI1, QR1, QS1, QG1, AKH
104!
105 integer I,J,L,LL,LP,LLMH,II,JJ,JJB,JJE,NHOLD
106 real PFSIGO,APFSIGO,PSIGO,APSIGO,PNL1,PU,ZU,TU,QU,QSAT, &
107 rhu,tvru,tvrabv,tabv,qabv,b,ahf,fac,pl,zl,tl,ql, &
108 rhl,tmt0,ai,bi,tvrl,tvrblo,tblo,qblo,fact, &
109 px,bf,facf,ahff,dpsig,tv,pdv,denom,denomf,pnl1f,dum
110!
111!
112!******************************************************************************
113!
114! START MDL2P.
115!
116! SET TOTAL NUMBER OF POINTS ON OUTPUT GRID.
117!
118!---------------------------------------------------------------
119!
120! *** PART I ***
121!
122! VERTICAL INTERPOLATION OF EVERYTHING ELSE. EXECUTE ONLY
123! IF THERE'S SOMETHING WE WANT.
124!
125 IF((iget(205)>0).OR.(iget(206)>0).OR. &
126 (iget(207)>0).OR.(iget(208)>0).OR. &
127 (iget(209)>0).OR.(iget(210)>0).OR. &
128 (iget(216)>0).OR.(iget(217)>0).OR. &
129 (iget(211)>0).OR.(iget(212)>0).OR. &
130 (iget(213)>0).OR.(iget(214)>0).OR. &
131 (iget(215)>0).OR.(iget(222)>0).OR. &
132 (iget(243)>0) ) THEN !!Air Quality (Plee Oct2003)
133!
134!---------------------------------------------------------------------
135!
136!--- VERTICAL INTERPOLATION OF GEOPOTENTIAL, SPECIFIC HUMIDITY, TEMPERATURE,
137! OMEGA, TKE, & CLOUD FIELDS. START AT THE UPPERMOST TARGET SIGMA LEVEL.
138!
139 readthk=.false.
140 IF(readthk)THEN ! EITHER READ DSG THICKNESS
141 READ(41)dsigo !DSIGO FROM TOP TO BOTTOM
142!
143 sigo(1)=0.0
144 DO l=2,lsig+1
145 sigo(l)=sigo(l-1)+dsigo(lsig-l+2)
146 END DO
147 sigo(lsig+1)=1.0
148 DO l=1,lsig
149 asigo(l)=0.5*(sigo(l)+sigo(l+1))
150 END DO
151 ELSE ! SPECIFY SIGO
152 asigo( 1)= 0.0530
153 asigo( 2)= 0.1580
154 asigo( 3)= 0.2605
155 asigo( 4)= 0.3595
156 asigo( 5)= 0.4550
157 asigo( 6)= 0.5470
158 asigo( 7)= 0.6180
159 asigo( 8)= 0.6690
160 asigo( 9)= 0.7185
161 asigo(10)= 0.7585
162 asigo(11)= 0.7890
163 asigo(12)= 0.8190
164 asigo(13)= 0.8480
165 asigo(14)= 0.8755
166 asigo(15)= 0.9015
167 asigo(16)= 0.9260
168 asigo(17)= 0.9490
169 asigo(18)= 0.9650
170 asigo(19)= 0.9745
171 asigo(20)= 0.9835
172 asigo(21)= 0.9915
173 asigo(22)= 0.9975
174!
175 sigo( 1)= 0.0
176 sigo( 2)= 0.1060
177 sigo( 3)= 0.2100
178 sigo( 4)= 0.3110
179 sigo( 5)= 0.4080
180 sigo( 6)= 0.5020
181 sigo( 7)= 0.5920
182 sigo( 8)= 0.6440
183 sigo( 9)= 0.6940
184 sigo(10)= 0.7430
185 sigo(11)= 0.7740
186 sigo(12)= 0.8040
187 sigo(13)= 0.8340
188 sigo(14)= 0.8620
189 sigo(15)= 0.8890
190 sigo(16)= 0.9140
191 sigo(17)= 0.9380
192 sigo(18)= 0.9600
193 sigo(19)= 0.9700
194 sigo(20)= 0.9790
195 sigo(21)= 0.9880
196 sigo(22)= 0.9950
197 sigo(23)= 1.0
198 END IF
199! OBTAIN GEOPOTENTIAL AT 1ST LEVEL
200 DO j=jsta_2l,jend_2u
201 DO i=ista_2l,iend_2u
202 fsl(i,j)=spval
203 akh(i,j)=spval
204 nl1xf(i,j)=lp1
205 DO l=1,lp1
206 IF(nl1xf(i,j)==lp1.AND.pint(i,j,l)>ptsigo)THEN
207 nl1xf(i,j)=l
208 ENDIF
209 ENDDO
210 END DO
211 END DO
212 DO 167 j=jsta,jend
213 DO 167 i=ista_2l,iend_2u
214 donefsl1=.false.
215 pfsigo=ptsigo
216 apfsigo=log(pfsigo)
217 pnl1=pint(i,j,nl1xf(i,j))
218 ll=nl1xf(i,j)
219 llmh = nint(lmh(i,j))
220 IF(nl1xf(i,j)==1 .AND. t(i,j,1)<spval &
221 .AND. t(i,j,2)<spval .AND. q(i,j,1)<spval &
222 .AND. q(i,j,2)<spval)THEN
223 pu=pint(i,j,2)
224 zu=zint(i,j,2)
225 tu=d50*(t(i,j,1)+t(i,j,2))
226 qu=d50*(q(i,j,1)+q(i,j,2))
227 qsat=pq0/pu*exp(a2*(tu-a3)/(tu-a4))
228 rhu =qu/qsat
229 IF(rhu>h1)THEN
230 rhu=h1
231 qu =rhu*qsat
232 ENDIF
233 IF(rhu<d01)THEN
234 rhu=d01
235 qu =rhu*qsat
236 ENDIF
237!
238 tvru=tu*(h1+d608*qu)
239 tvrabv=tvru*(pfsigo/pu)**rgamog
240 tabv=tvrabv/(h1+d608*qu)
241 qsat=pq0/pfsigo*exp(a2*(tabv-a3)/(tabv-a4))
242 qabv =rhu*qsat
243 qabv =max(h1m12,qabv)
244!== B =TABV
245 b =tvrabv !Marina Tsidulko Dec22, 2003
246 ahf =d00
247 fac =d00
248 donefsl1=.true.
249 ELSEIF(nl1xf(i,j)==lp1 .AND. t(i,j,lm-1)<spval &
250 .AND. t(i,j,lm-2)<spval .AND. q(i,j,lm-1)<spval &
251 .AND. q(i,j,lm-2)<spval)THEN
252 pl=pint(i,j,lm-1)
253 zl=zint(i,j,lm-1)
254 tl=d50*(t(i,j,lm-2)+t(i,j,lm-1))
255 ql=d50*(q(i,j,lm-2)+q(i,j,lm-1))
256!
257!--- RH w/r/t water for all conditions (Ferrier, 25 Jan 02)
258!
259 qsat=pq0/pl*exp(a2*(tl-a3)/(tl-a4))
260 rhl=ql/qsat
261 IF(rhl>h1)THEN
262 rhl=h1
263 ql =rhl*qsat
264 ENDIF
265 IF(rhl<d01)THEN
266 rhl=d01
267 ql =rhl*qsat
268 ENDIF
269!
270 tvrl =tl*(h1+d608*ql)
271 tvrblo=tvrl*(pfsigo/pl)**rgamog
272 tblo =tvrblo/(h1+d608*ql)
273 qsat=pq0/pfsigo*exp(a2*(tblo-a3)/(tblo-a4))
274 qblo =rhl*qsat
275 qblo =max(h1m12,qblo)
276!== B =TBLO
277 b =tvrblo !Marina Tsidulko Dec22, 2003
278 ahf =d00
279 fac =d00
280 donefsl1=.true.
281 ELSEIF(nl1xf(i,j)<lp1) THEN ! must check the "<LP1" first; do not combine these into a single "if"
282 IF(t(i,j,nl1xf(i,j))<spval &
283 & .AND. q(i,j,nl1xf(i,j))<spval)THEN
284!== B =T(I,J,NL1XF(I,J)) !Marina Tsidulko Dec22, 2003
285 b =t(i,j,nl1xf(i,j))*(h1+d608*q(i,j,nl1xf(i,j)))
286 denom=(alpint(i,j,nl1xf(i,j)+1)-alpint(i,j,nl1xf(i,j)-1))
287!== AHF =(B-T(I,J,NL1XF(I,J)-1))/DENOM
288 ahf =(b-t(i,j,nl1xf(i,j)-1)*(h1+d608*q(i,j,nl1xf(i,j)-1))) &
289 /denom !Marina Tsidulko Dec22, 2003
290 fac =h2*log(pmid(i,j,nl1xf(i,j)))
291 donefsl1=.true.
292 END IF
293 END IF
294!
295 if(donefsl1)fsl1(i,j)=(pnl1-pfsigo)/(pfsigo+pnl1) &
296 *((apfsigo+alpint(i,j,nl1xf(i,j))-fac)*ahf+b)*rd*h2 &
297 +zint(i,j,nl1xf(i,j))*g
298! COMPUTE EXCHANGE COEFFICIENT ON FIRST INTERFACTE
299 IF(nl1xf(i,j)<=2 .OR. nl1xf(i,j)>(llmh+1))THEN
300 akh(i,j)=0.0
301 ELSE
302 fact=(apfsigo-log(pint(i,j,ll)))/ &
303 & (log(pint(i,j,ll))-log(pint(i,j,ll-1)))
304! EXCH_H is on the bottom of model interfaces
305 IF(exch_h(i,j,ll-2)<spval .AND. exch_h(i,j,ll-1)<spval) &
306 & akh(i,j)=exch_h(i,j,ll-1)+(exch_h(i,j,ll-1) &
307 & -exch_h(i,j,ll-2))*fact
308 END IF
309 167 CONTINUE
310! OUTPUT FIRST LAYER GEOPOTENTIAL
311! GEOPOTENTIAL (SCALE BY GI)
312 IF (iget(205)>0) THEN
313 IF (lvls(1,iget(205))>0) THEN
314!$omp parallel do
315 DO j=jsta,jend
316 DO i=ista,iend
317 IF(fsl1(i,j)<spval) THEN
318 grid1(i,j)=fsl1(i,j)*gi
319 ELSE
320 grid1(i,j)=spval
321 ENDIF
322 ENDDO
323 ENDDO
324 if(grib=='grib2')then
325 cfld=cfld+1
326 fld_info(cfld)%ifld=iavblfld(iget(205))
327 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
328 endif
329 ENDIF
330 ENDIF
331
332! OUTPUT FIRST INTERFACE KH Heat Diffusivity
333 IF (iget(243)>0) THEN !!Air Quality (Plee Oct2003) ^^^^^
334 IF (lvls(1,iget(243))>0) THEN
335!$omp parallel do
336 DO j=jsta,jend
337 DO i=ista,iend
338 grid1(i,j)=akh(i,j)
339 ENDDO
340 ENDDO
341 if(grib=="grib2" )then
342 cfld=cfld+1
343 fld_info(cfld)%ifld=iavblfld(iget(243))
344 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
345 endif
346
347 if(me==0)print*,'output Heat Diffusivity'
348 ENDIF
349 ENDIF
350
351!***
352!*** BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND, DO ALL
353!*** INTERPOLATION ABOVE GROUND NOW.
354!***
355!
356 DO 310 lp=1,lsig
357 nhold=0
358!
359 DO j=jsta_2l,jend_2u
360 DO i=ista_2l,iend_2u
361
362!
363 tsl(i,j)=spval
364 qsl(i,j)=spval
365 fsl(i,j)=spval
366 osl(i,j)=spval
367 usl(i,j)=spval
368 vsl(i,j)=spval
369 q2sl(i,j)=spval
370 c1d(i,j)=spval ! Total condensate
371 qw1(i,j)=spval ! Cloud water
372 qi1(i,j)=spval ! Cloud ice
373 qr1(i,j)=spval ! Rain
374 qs1(i,j)=spval ! Snow (precip ice)
375 qg1(i,j)=spval
376 cfrsig(i,j)=spval
377!
378!*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER JUST BELOW
379!*** THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
380!
381 nl1x(i,j)=lp1
382 DO l=2,lm
383 llmh = nint(lmh(i,j))
384 psigo=ptsigo+asigo(lp)*(pint(i,j,llmh+1)-ptsigo)
385 IF(nl1x(i,j)==lp1.AND.pmid(i,j,l)>psigo)THEN
386 nl1x(i,j)=l
387 ENDIF
388 ENDDO
389!
390! IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
391! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
392! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
393! WILL EXTRAPOLATE TO THAT POINT
394!
395 IF(nl1x(i,j)==lp1.AND.pint(i,j,llmh+1)>=psigo)THEN
396 nl1x(i,j)=lm
397 ENDIF
398!
399! if(NL1X(I,J)==LP1)print*,'Debug: NL1X=LP1 AT '
400! 1 ,i,j,lp
401 ENDDO
402 ENDDO
403!
404!mptest IF(NHOLD==0)GO TO 310
405!
406!$omp parallel do private(i,j,ll,llmh,psigo,apsigo,fact,dum,pl, &
407!$omp & zl,tl,ql,ai,bi,qsat,rhl,tvrl,tvrblo,tblo,tmt0, &
408!$omp & qblo,pnl1,fac,ahf)
409!hc DO 220 NN=1,NHOLD
410!hc I=IHOLD(NN)
411!hc J=JHOLD(NN)
412 DO 220 j=jsta,jend ! Moorthi on Nov 26 2014
413! DO 220 J=JSTA_2L,JEND_2U
414 DO 220 i=ista,iend
415 ll=nl1x(i,j)
416!---------------------------------------------------------------------
417!*** VERTICAL INTERPOLATION OF GEOPOTENTIAL, TEMPERATURE, SPECIFIC
418!*** HUMIDITY, CLOUD WATER/ICE, OMEGA, WINDS, AND TKE.
419!---------------------------------------------------------------------
420!
421!HC IF(NL1X(I,J)<=LM)THEN
422 llmh = nint(lmh(i,j))
423 psigo=ptsigo+asigo(lp)*(pint(i,j,llmh+1)-ptsigo)
424 apsigo=log(psigo)
425 IF(nl1x(i,j)<=llmh)THEN
426!
427!---------------------------------------------------------------------
428! INTERPOLATE LINEARLY IN LOG(P)
429!*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
430!*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
431!*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
432!---------------------------------------------------------------------
433!
434
435 fact=(apsigo-log(pmid(i,j,ll)))/ &
436 & (log(pmid(i,j,ll))-log(pmid(i,j,ll-1)))
437 tsl(i,j)=t(i,j,ll)+(t(i,j,ll)-t(i,j,ll-1))*fact
438 IF(q(i,j,ll)<spval .AND. q(i,j,ll-1)<spval) &
439 & qsl(i,j)=q(i,j,ll)+(q(i,j,ll)-q(i,j,ll-1))*fact
440 IF(gridtype=='A')THEN
441 IF(uh(i,j,ll)<spval .AND. uh(i,j,ll-1)<spval) &
442 & usl(i,j)=uh(i,j,ll)+(uh(i,j,ll)-uh(i,j,ll-1))*fact
443 IF(vh(i,j,ll)<spval .AND. vh(i,j,ll-1)<spval) &
444 & vsl(i,j)=vh(i,j,ll)+(vh(i,j,ll)-vh(i,j,ll-1))*fact
445 END IF
446 IF(omga(i,j,ll)<spval .AND. omga(i,j,ll-1)<spval) &
447 & osl(i,j)=omga(i,j,ll)+(omga(i,j,ll)-omga(i,j,ll-1))*fact
448 IF(q2(i,j,ll)<spval .AND. q2(i,j,ll-1)<spval) &
449 & q2sl(i,j)=q2(i,j,ll)+(q2(i,j,ll)-q2(i,j,ll-1))*fact
450! FSL(I,J)=ZMID(I,J,LL)+(ZMID(I,J,LL)-ZMID(I,J,LL-1))*FACT
451! FSL(I,J)=FSL(I,J)*G
452!hc QSAT=PQ0/PSIGO
453!hc 1 *EXP(A2*(TSL(I,J)-A3)/(TSL(I,J)-A4))
454!
455!hc RHL=QSL(I,J)/QSAT
456!
457!hc IF(RHL>1.) QSL(I,J)=QSAT
458!hc IF(RHL<0.01) QSL(I,J)=0.01*QSAT
459 IF(q2sl(i,j)<0.0) q2sl(i,j)=0.0
460!
461!HC ADD FERRIER'S HYDROMETEOR
462 IF(cwm(i,j,ll)<spval .AND. cwm(i,j,ll-1)<spval) &
463 & c1d(i,j)=cwm(i,j,ll)+(cwm(i,j,ll)-cwm(i,j,ll-1))*fact
464 c1d(i,j)=max(c1d(i,j),h1m12) ! Total condensate
465 IF(qqw(i,j,ll)<spval .AND. qqw(i,j,ll-1)<spval) &
466 & qw1(i,j)=qqw(i,j,ll)+(qqw(i,j,ll)-qqw(i,j,ll-1))*fact
467 qw1(i,j)=max(qw1(i,j),h1m12) ! Cloud water
468 IF(qqi(i,j,ll)<spval .AND. qqi(i,j,ll-1)<spval) &
469 & qi1(i,j)=qqi(i,j,ll)+(qqi(i,j,ll)-qqi(i,j,ll-1))*fact
470 qi1(i,j)=max(qi1(i,j),h1m12) ! Cloud ice
471 IF(qqr(i,j,ll)<spval .AND. qqr(i,j,ll-1)<spval) &
472 & qr1(i,j)=qqr(i,j,ll)+(qqr(i,j,ll)-qqr(i,j,ll-1))*fact
473 qr1(i,j)=max(qr1(i,j),h1m12) ! Rain
474 IF(qqs(i,j,ll)<spval .AND. qqs(i,j,ll-1)<spval) &
475 & qs1(i,j)=qqs(i,j,ll)+(qqs(i,j,ll)-qqs(i,j,ll-1))*fact
476 qs1(i,j)=max(qs1(i,j),h1m12) ! Snow (precip ice)
477 IF(cfr(i,j,ll)<spval .AND. cfr(i,j,ll-1)<spval) &
478 & cfrsig(i,j)=cfr(i,j,ll)+(cfr(i,j,ll)-cfr(i,j,ll-1))*fact
479 cfrsig(i,j)=max(cfrsig(i,j),h1m12)
480 IF(qqs(i,j,ll)<spval .AND. qqs(i,j,ll-1)<spval)THEN
481 dum=f_rimef(i,j,ll)+(f_rimef(i,j,ll)-f_rimef(i,j,ll-1))*fact
482 IF(dum <= 5.0)THEN
483 qg1(i,j)=h1m12
484 ELSE
485 qg1(i,j)=qs1(i,j)
486 qs1(i,j)=h1m12
487 END IF
488 END IF
489! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE
490! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
491! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
492! GOUND
493 ELSE
494! ii=91
495! jj=13
496! if(i==ii.and.j==jj)print*,'Debug: underg extra at i,j,lp' &
497! &, i,j,lp
498 pl=pint(i,j,lm-1)
499 zl=zint(i,j,lm-1)
500 tl=0.5*(t(i,j,lm-2)+t(i,j,lm-1))
501 ql=0.5*(q(i,j,lm-2)+q(i,j,lm-1))
502 ai=0.008855
503 bi=1.
504 tmt0=tl-a3
505 IF(tmt0<-20.)THEN
506 ai=0.007225
507 bi=0.9674
508 ENDIF
509 qsat=pq0/pl*exp(a2*(tl-a3)/(tl-a4))
510!
511 rhl=ql/qsat
512!
513 IF(rhl>1.)THEN
514 rhl=1.
515 ql =rhl*qsat
516 ENDIF
517!
518 IF(rhl<0.01)THEN
519 rhl=0.01
520 ql =rhl*qsat
521 ENDIF
522!
523 tvrl =tl*(1.+0.608*ql)
524 tvrblo=tvrl*(psigo/pl)**rgamog
525 tblo =tvrblo/(1.+0.608*ql)
526!
527 tmt0=tblo-a3
528 ai=0.008855
529 bi=1.
530 IF(tmt0<-20.)THEN
531 ai=0.007225
532 bi=0.9674
533 ENDIF
534 qsat=pq0/psigo*exp(a2*(tblo-a3)/(tblo-a4))
535!
536! TSL(I,J)=TBLO
537 qblo = rhl*qsat
538 qsl(i,j) = max(1.e-12,qblo)
539 IF(gridtype=='A')THEN
540 usl(i,j) = uh(i,j,llmh)
541 vsl(i,j) = vh(i,j,llmh)
542 END IF
543 osl(i,j) = omga(i,j,llmh)
544 q2sl(i,j) = max(0.0,0.5*(q2(i,j,llmh-1)+q2(i,j,llmh)))
545 pnl1 = pint(i,j,nl1x(i,j))
546 fac = 0.
547 ahf = 0.0
548!
549!--- Set hydrometeor fields to zero below ground
550 c1d(i,j)=0.
551 qw1(i,j)=0.
552 qi1(i,j)=0.
553 qr1(i,j)=0.
554 qs1(i,j)=0.
555 qg1(i,j)=0.
556 cfrsig(i,j)=0.
557 END IF
558 220 CONTINUE
559!
560! OBTAIN GEOPOTENTIAL AND KH ON INTERFACES
561 DO j=jsta_2l,jend_2u
562 DO i=ista_2l,iend_2u
563 fsl(i,j)=spval
564 akh(i,j)=spval
565 nl1xf(i,j)=lp1
566 llmh = nint(lmh(i,j))
567 psigo=ptsigo+sigo(lp+1)*(pint(i,j,llmh+1)-ptsigo)
568 DO l=1,lp1
569 IF(nl1xf(i,j)==lp1.AND.pint(i,j,l)>psigo)THEN
570 nl1xf(i,j)=l
571 ENDIF
572 ENDDO
573 END DO
574 END DO
575!
576! DO J=JSTA_2L,JEND_2U
577 DO j=jsta,jend ! Moorthi on 26 Nov 2014
578 DO i=ista,iend
579 donefsl1=.false.
580 tsldone=.false.
581 llmh = nint(lmh(i,j))
582 pfsigo=ptsigo+sigo(lp+1)*(pint(i,j,llmh+1)-ptsigo)
583 psigo=ptsigo+asigo(lp)*(pint(i,j,llmh+1)-ptsigo)
584 apfsigo=log(pfsigo)
585 pnl1f=pint(i,j,nl1xf(i,j))
586 ll=nl1xf(i,j)
587 IF(nl1xf(i,j)==1 .AND. t(i,j,1)<spval &
588 & .AND. t(i,j,2)<spval .AND. q(i,j,1)<spval &
589 & .AND. q(i,j,2)<spval)THEN
590 pu=pint(i,j,2)
591 zu=zint(i,j,2)
592 tu=d50*(t(i,j,1)+t(i,j,2))
593 qu=d50*(q(i,j,1)+q(i,j,2))
594!
595!--- RH w/r/t water for all conditions
596!
597 qsat=pq0/pu*exp(a2*(tu-a3)/(tu-a4))
598 rhu =qu/qsat
599 IF(rhu>h1)THEN
600 rhu=h1
601 qu =rhu*qsat
602 ENDIF
603 IF(rhu<d01)THEN
604 rhu=d01
605 qu =rhu*qsat
606 ENDIF
607!
608 tvru=tu*(h1+d608*qu)
609! TVRABV=TVRU*(SPL(L)/PU)**RGAMOG
610!== TVRABV=TVRU*(PFSIGO/PU)**RGAMOG
611 px=(pfsigo+pnl1f)*0.5 !Marina Tsidulko Dec22, 2003
612 tvrabv=tvru*(px/pu)**rgamog!Marina Tsidulko Dec22, 2003
613 tabv=tvrabv/(h1+d608*qu)
614! ADD ADDITIONAL BLOCK FOR INTERPOLATING HEIGHT ONTO INTERFACES
615!== BF =TABV
616 bf =tvrabv !Marina Tsidulko Dec22, 2003
617 facf =d00
618 ahff =d00
619 donefsl1=.true.
620!
621!
622 ELSEIF(nl1xf(i,j)==lp1 .AND. t(i,j,lm-1)<spval &
623 & .AND. t(i,j,lm-2)<spval .AND. q(i,j,lm-1)<spval &
624 & .AND. q(i,j,lm-2)<spval)THEN
625!
626! EXTRAPOLATION AT LOWER BOUND. THE LOWER BOUND IS
627! LM IF OLDRD=.FALSE. IF OLDRD=.TRUE. THE LOWER
628! BOUND IS THE FIRST ATMOSPHERIC ETA LAYER.
629!
630 pl=pint(i,j,lm-1)
631 zl=zint(i,j,lm-1)
632 tl=d50*(t(i,j,lm-2)+t(i,j,lm-1))
633 ql=d50*(q(i,j,lm-2)+q(i,j,lm-1))
634!
635!--- RH w/r/t water for all conditions (Ferrier, 25 Jan 02)
636!
637 qsat=pq0/pl*exp(a2*(tl-a3)/(tl-a4))
638 rhl=ql/qsat
639 IF(rhl>h1)THEN
640 rhl=h1
641 ql =rhl*qsat
642 ENDIF
643 IF(rhl<d01)THEN
644 rhl=d01
645 ql =rhl*qsat
646 ENDIF
647!
648 tvrl =tl*(h1+d608*ql)
649!== TVRBLO=TVRL*(PFSIGO/PL)**RGAMOG
650 px=(pfsigo+pnl1f)*0.5 !Marina Tsidulko Dec22, 2003
651 tvrblo=tvrl*(px/pl)**rgamog !Marina Tsidulko Dec22, 2003
652 tblo =tvrblo/(h1+d608*ql)
653! ADD ADDITIONAL BLOCK FOR INTERPOLATING HEIGHT ONTO INTERFACES
654!== BF =TBLO
655 bf =tvrblo !Marina Tsidulko Dec22, 2003
656 facf =d00
657 ahff =d00
658 donefsl1=.true.
659 tsldone=.true.
660!
661!
662 ELSEIF(nl1xf(i,j)<lp1) THEN ! must check the "<LP1" first; do not combine these into a single "if"
663 IF(t(i,j,nl1xf(i,j))<spval &
664 & .AND. q(i,j,nl1xf(i,j))<spval)THEN
665!
666! INTERPOLATION BETWEEN LOWER AND UPPER BOUNDS.
667!
668! ADD ADDITIONAL BLOCK FOR INTERPOLATING HEIGHT ONTO INTERFACES
669! if(NL1XF(I,J)==LMp1)
670! + print*,'Debug: using bad temp at',i,j
671!== BF =T(I,J,NL1XF(I,J)) !Marina Tsidulko Dec22, 2003
672 bf =t(i,j,nl1xf(i,j))*(h1+d608*q(i,j,nl1xf(i,j)))
673!HC FACF =H2*LOG(PT+PDSL(I,J)*AETA(NL1XF(I,J)))
674 facf =h2*log(pmid(i,j,nl1xf(i,j)))
675 denomf=(alpint(i,j,nl1xf(i,j)+1)-alpint(i,j,nl1xf(i,j)-1))
676!== AHFF=(BF-T(I,J,NL1XF(I,J)-1))/DENOMF
677 ahff=(bf-t(i,j,nl1xf(i,j)-1)*(h1+d608*q(i,j,nl1xf(i,j)-1))) &
678 /denomf !Marina Tsidulko Dec22, 2003
679!
680 donefsl1=.true.
681 ENDIF
682 ENDIF
683!
684 IF(donefsl1)then
685 fsl(i,j)=(pnl1f-pfsigo)/(pfsigo+pnl1f) &
686 *((apfsigo+alpint(i,j,nl1xf(i,j))-facf)*ahff+bf)*rd*h2 &
687 +zint(i,j,nl1xf(i,j))*g
688! OBTAIN TEMPERATURE AT MID-SIGMA LAYER BASED ON HYDROSTATIC
689 dpsig=(sigo(lp+1)-sigo(lp))*(pint(i,j,llmh+1)-ptsigo)
690! IF(J==jj.and.i==ii)print*,'Debug:L,OLD T= ',
691! + L,TSL(I,J)
692 IF(.NOT.tsldone) THEN
693 tsl(i,j)=(fsl1(i,j)-fsl(i,j))*psigo/(rd*dpsig)
694 ENDIF
695 fsl1(i,j)=fsl(i,j)
696 IF(.NOT.tsldone) THEN
697 tv=tsl(i,j)
698 tsl(i,j)=tv/(h1+d608*qsl(i,j))
699 ENDIF
700 qsat=pq0/psigo *exp(a2*(tsl(i,j)-a3)/(tsl(i,j)-a4))
701!
702 rhl=qsl(i,j)/qsat
703!
704 IF(rhl>1.) qsl(i,j)=qsat
705 IF(rhl<0.01) qsl(i,j)=0.01*qsat
706 END IF
707! IF(J==jj.and.i==ii)print*,'Debug:L,T,Q,Q2,FSL=',
708! + L,TSL(I,J),QSL(I,J),Q2SL(I,J),FSL(I,J)
709
710! COMPUTE EXCHANGE COEFFICIENT ON INTERFACES
711 IF(nl1xf(i,j)<=2 .OR. nl1xf(i,j)>(llmh+1))THEN
712 akh(i,j)=0.0
713 ELSE
714 fact=(apfsigo-log(pint(i,j,ll)))/ &
715 & (log(pint(i,j,ll))-log(pint(i,j,ll-1)))
716! EXCH_H is on the bottom of model interfaces
717 IF(exch_h(i,j,ll-2)<spval .AND. &
718 & exch_h(i,j,ll-1)<spval) &
719 & akh(i,j)=exch_h(i,j,ll-1)+(exch_h(i,j,ll-1) &
720 & -exch_h(i,j,ll-2))*fact
721 END IF
722
723! END OF HYDROSTATIC TEMPERATURE DERIVATION
724 END DO
725 END DO
726!
727! VERTICAL INTERPOLATION FOR WIND FOR E and B GRIDS
728!
729 if(gridtype=='B' .or. gridtype=='E') &
730 call exch(pint(ista_2l:iend_2u,jsta_2l:jend_2u,lp1))
731 IF(gridtype=='E')THEN
732 DO j=jsta,jend
733! DO I=1,IM-MOD(J,2)
734 DO i=ista,iend-mod(j,2) !Jesse 20211014
735!
736!*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER FOR V POINT JUST BELOW
737!*** THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
738!
739 llmh = nint(lmh(i,j))
740
741!Jesse 20211014
742! IF(J == 1 .AND. I < IM)THEN !SOUTHERN BC
743! PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
744! ELSE IF(J==JM .AND. I<IM)THEN !NORTHERN BC
745! PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
746! ELSE IF(I == 1 .AND. MOD(J,2) == 0) THEN !WESTERN EVEN BC
747! PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
748! ELSE IF(I == IM .AND. MOD(J,2) == 0) THEN !EASTERN EVEN BC
749! PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
750! ELSE IF (MOD(J,2) < 1) THEN
751! PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I-1,J,LLMH+1) &
752! & +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
753! ELSE
754! PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1) &
755! & +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
756! END IF
757
758 IF(j == 1 .AND. i < iend)THEN !SOUTHERN BC
759 pdv=0.5*(pint(i,j,llmh+1)+pint(i+1,j,llmh+1))
760 ELSE IF(j==jm .AND. i<iend)THEN !NORTHERN BC
761 pdv=0.5*(pint(i,j,llmh+1)+pint(i+1,j,llmh+1))
762 ELSE IF(i == ista .AND. mod(j,2) == 0) THEN !WESTERN EVEN BC
763 pdv=0.5*(pint(i,j-1,llmh+1)+pint(i,j+1,llmh+1))
764 ELSE IF(i == iend .AND. mod(j,2) == 0) THEN !EASTERN EVEN BC
765 pdv=0.5*(pint(i,j-1,llmh+1)+pint(i,j+1,llmh+1))
766 ELSE IF (mod(j,2) < 1) THEN
767 pdv=0.25*(pint(i,j,llmh+1)+pint(i-1,j,llmh+1) &
768 & +pint(i,j+1,llmh+1)+pint(i,j-1,llmh+1))
769 ELSE
770 pdv=0.25*(pint(i,j,llmh+1)+pint(i+1,j,llmh+1) &
771 & +pint(i,j+1,llmh+1)+pint(i,j-1,llmh+1))
772 END IF
773
774 psigo=ptsigo+asigo(lp)*(pdv-ptsigo)
775 nl1x(i,j)=lp1
776 DO l=2,lm
777 IF(nl1x(i,j)==lp1.AND.pmidv(i,j,l)>psigo)THEN
778 nl1x(i,j)=l
779 ENDIF
780 ENDDO
781!
782! IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
783! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
784! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
785! WILL EXTRAPOLATE TO THAT POINT
786!
787 IF(nl1x(i,j)==lp1.AND. pdv>psigo)THEN
788 nl1x(i,j)=lm
789 ENDIF
790!
791 ENDDO
792 ENDDO
793!
794 DO 230 j=jsta,jend
795! DO 230 I=1,IM-MOD(j,2)
796 DO 230 i=ista,iend-mod(j,2) !Jesse 20211014
797
798 llmh = nint(lmh(i,j))
799
800!Jesse 20211014
801! IF(J == 1 .AND. I < IM)THEN !SOUTHERN BC
802! PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
803! ELSE IF(J==JM .AND. I<IM)THEN !NORTHERN BC
804! PDV=0.5*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1))
805! ELSE IF(I == 1 .AND. MOD(J,2) == 0) THEN !WESTERN EVEN BC
806! PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
807! ELSE IF(I == IM .AND. MOD(J,2) == 0) THEN !EASTERN EVEN BC
808! PDV=0.5*(PINT(I,J-1,LLMH+1)+PINT(I,J+1,LLMH+1))
809! ELSE IF (MOD(J,2) < 1) THEN
810! PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I-1,J,LLMH+1) &
811! & +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
812! ELSE
813! PDV=0.25*(PINT(I,J,LLMH+1)+PINT(I+1,J,LLMH+1) &
814! & +PINT(I,J+1,LLMH+1)+PINT(I,J-1,LLMH+1))
815! END IF
816
817 IF(j == 1 .AND. i < iend)THEN !SOUTHERN BC
818 pdv=0.5*(pint(i,j,llmh+1)+pint(i+1,j,llmh+1))
819 ELSE IF(j==jm .AND. i<iend)THEN !NORTHERN BC
820 pdv=0.5*(pint(i,j,llmh+1)+pint(i+1,j,llmh+1))
821 ELSE IF(i == ista .AND. mod(j,2) == 0) THEN !WESTERN EVEN BC
822 pdv=0.5*(pint(i,j-1,llmh+1)+pint(i,j+1,llmh+1))
823 ELSE IF(i == iend .AND. mod(j,2) == 0) THEN !EASTERN EVEN BC
824 pdv=0.5*(pint(i,j-1,llmh+1)+pint(i,j+1,llmh+1))
825 ELSE IF (mod(j,2) < 1) THEN
826 pdv=0.25*(pint(i,j,llmh+1)+pint(i-1,j,llmh+1) &
827 & +pint(i,j+1,llmh+1)+pint(i,j-1,llmh+1))
828 ELSE
829 pdv=0.25*(pint(i,j,llmh+1)+pint(i+1,j,llmh+1) &
830 & +pint(i,j+1,llmh+1)+pint(i,j-1,llmh+1))
831 END IF
832
833 psigo=ptsigo+asigo(lp)*(pdv-ptsigo)
834 apsigo=log(psigo)
835 ll=nl1x(i,j)
836!---------------------------------------------------------------------
837!*** VERTICAL INTERPOLATION OF WINDS FOR A-E GRID
838!---------------------------------------------------------------------
839!
840!HC IF(NL1X(I,J)<=LM)THEN
841 llmh = nint(lmh(i,j))
842 IF(nl1x(i,j)<=llmh)THEN
843!
844!---------------------------------------------------------------------
845! INTERPOLATE LINEARLY IN LOG(P)
846!*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
847!*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
848!*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
849!---------------------------------------------------------------------
850!
851
852 fact=(apsigo-log(pmidv(i,j,ll)))/ &
853 & (log(pmidv(i,j,ll))-log(pmidv(i,j,ll-1)))
854 IF(uh(i,j,ll)<spval .AND. uh(i,j,ll-1)<spval) &
855 & usl(i,j)=uh(i,j,ll)+(uh(i,j,ll)-uh(i,j,ll-1))*fact
856 IF(vh(i,j,ll)<spval .AND. vh(i,j,ll-1)<spval) &
857 & vsl(i,j)=vh(i,j,ll)+(vh(i,j,ll)-vh(i,j,ll-1))*fact
858!
859! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE
860! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
861! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
862! GOUND
863 ELSE
864 IF(uh(i,j,llmh)<spval)usl(i,j)=uh(i,j,llmh)
865 IF(vh(i,j,llmh)<spval)vsl(i,j)=vh(i,j,llmh)
866 END IF
867 230 CONTINUE
868 jjb=jsta
869 IF(mod(jsta,2)==0)jjb=jsta+1
870 jje=jend
871 IF(mod(jend,2)==0)jje=jend-1
872 DO j=jjb,jje,2 !chc
873 usl(iend,j)=usl(iend-1,j)
874 vsl(iend,j)=vsl(iend-1,j)
875 END DO
876
877 ELSE IF (gridtype=='B')THEN
878 DO j=jsta,jend_m
879 DO i=ista,iend_m
880!
881!*** LOCATE VERTICAL INDEX OF MODEL MIDLAYER FOR V POINT JUST BELOW
882!*** THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
883!
884 pdv=0.25*(pint(i,j,lp1)+pint(i+1,j,lp1) &
885 +pint(i,j+1,lp1)+pint(i+1,j+1,lp1))
886
887 psigo=ptsigo+asigo(lp)*(pdv-ptsigo)
888 nl1x(i,j)=lp1
889 DO l=2,lm
890 IF(nl1x(i,j)==lp1.AND.pmidv(i,j,l)>psigo)THEN
891 nl1x(i,j)=l
892 ENDIF
893 ENDDO
894!
895! IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
896! BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
897! WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
898! WILL EXTRAPOLATE TO THAT POINT
899!
900 IF(nl1x(i,j)==lp1.AND. pdv>psigo)THEN
901 nl1x(i,j)=lm
902 ENDIF
903!
904 ENDDO
905 ENDDO
906!
907 DO 231 j=jsta,jend_m
908 DO 231 i=ista,iend_m
909 pdv=0.25*(pint(i,j,lp1)+pint(i+1,j,lp1) &
910 +pint(i,j+1,lp1)+pint(i+1,j+1,lp1))
911 psigo=ptsigo+asigo(lp)*(pdv-ptsigo)
912 apsigo=log(psigo)
913 ll=nl1x(i,j)
914!---------------------------------------------------------------------
915!*** VERTICAL INTERPOLATION OF WINDS FOR A-E GRID
916!---------------------------------------------------------------------
917!
918!HC IF(NL1X(I,J)<=LM)THEN
919 llmh = nint(lmh(i,j))
920 IF(nl1x(i,j)<=llmh)THEN
921!
922!---------------------------------------------------------------------
923! INTERPOLATE LINEARLY IN LOG(P)
924!*** EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
925!*** INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
926!*** EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
927!---------------------------------------------------------------------
928!
929
930 fact=(apsigo-log(pmidv(i,j,ll)))/ &
931 & (log(pmidv(i,j,ll))-log(pmidv(i,j,ll-1)))
932 IF(uh(i,j,ll)<spval .AND. uh(i,j,ll-1)<spval) &
933 & usl(i,j)=uh(i,j,ll)+(uh(i,j,ll)-uh(i,j,ll-1))*fact
934 IF(vh(i,j,ll)<spval .AND. vh(i,j,ll-1)<spval) &
935 & vsl(i,j)=vh(i,j,ll)+(vh(i,j,ll)-vh(i,j,ll-1))*fact
936!
937! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE
938! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
939! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
940! GOUND
941 ELSE
942 IF(uh(i,j,llmh)<spval)usl(i,j)=uh(i,j,llmh)
943 IF(vh(i,j,llmh)<spval)vsl(i,j)=vh(i,j,llmh)
944 END IF
945 231 CONTINUE
946
947
948
949 END IF ! END OF WIND INTERPOLATION FOR NMM
950
951
952!
953!---------------------------------------------------------------------
954! LOAD GEOPOTENTIAL AND TEMPERATURE INTO STANDARD LEVEL
955! ARRAYS FOR THE NEXT PASS.
956!---------------------------------------------------------------------
957!
958!---------------------------------------------------------------------
959!*** CALCULATE 1000MB GEOPOTENTIALS CONSISTENT WITH SLP OBTAINED
960!*** FROM THE MESINGER OR NWS SHUELL SLP REDUCTION.
961!---------------------------------------------------------------------
962!
963!---------------------------------------------------------------------
964!---------------------------------------------------------------------
965! *** PART II ***
966!---------------------------------------------------------------------
967!---------------------------------------------------------------------
968!
969! INTERPOLATE/OUTPUT SELECTED FIELDS.
970!
971!---------------------------------------------------------------------
972!
973!*** OUTPUT GEOPOTENTIAL (SCALE BY GI)
974!
975 IF(iget(205)>0)THEN
976 IF(lvls(lp+1,iget(205))>0)THEN
977!$omp parallel do
978 DO j=jsta,jend
979 DO i=ista,iend
980 IF(fsl(i,j)<spval) THEN
981 grid1(i,j)=fsl(i,j)*gi
982 ELSE
983 grid1(i,j)=spval
984 ENDIF
985 ENDDO
986 ENDDO
987 if(grib=="grib2" )then
988 cfld=cfld+1
989 fld_info(cfld)%ifld=iavblfld(iget(205))
990 fld_info(cfld)%lvl=lvlsxml(lp+1,iget(205))
991 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
992 endif
993 ENDIF
994 ENDIF
995!
996!*** OUTPUT EXCHANGE COEEFICIENT
997!
998! OUTPUT FIRST INTERFACE KH Heat Diffusivity
999 IF (iget(243)>0) THEN !!Air Quality (Plee Oct2003) ^^^^^
1000 IF (lvls(lp+1,iget(243))>0) THEN
1001!$omp parallel do
1002 DO j=jsta,jend
1003 DO i=ista,iend
1004 grid1(i,j)=akh(i,j)
1005 IF(lp==(lsig+1))grid1(i,j)=0.0 !! NO SLIP ASSUMTION FOR CMAQ
1006 ENDDO
1007 ENDDO
1008 if(grib=="grib2" )then
1009 cfld=cfld+1
1010 fld_info(cfld)%ifld=iavblfld(iget(243))
1011 fld_info(cfld)%lvl=lvlsxml(lp+1,iget(243))
1012 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1013 endif
1014 if(me==0)print*,'output Heat Diffusivity'
1015 ENDIF
1016 ENDIF
1017!
1018!*** TEMPERATURE
1019!
1020 IF(iget(206)>0) THEN
1021 IF(lvls(lp,iget(206))>0) THEN
1022 DO j=jsta,jend
1023 DO i=ista,iend
1024 grid1(i,j)=tsl(i,j)
1025 ENDDO
1026 ENDDO
1027 if(grib=="grib2" )then
1028 cfld=cfld+1
1029 fld_info(cfld)%ifld=iavblfld(iget(206))
1030 fld_info(cfld)%lvl=lvlsxml(lp,iget(206))
1031 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1032 endif
1033 ENDIF
1034 ENDIF
1035!
1036!*** PRESSURE
1037!
1038 IF(iget(216)>0)THEN
1039 IF(lvls(lp,iget(216))>0)THEN
1040!$omp parallel do
1041 DO j=jsta,jend
1042 DO i=ista,iend
1043 llmh = nint(lmh(i,j))
1044 grid1(i,j)=ptsigo+asigo(lp)*(pint(i,j,llmh+1)-ptsigo)
1045 ENDDO
1046 ENDDO
1047 if(grib=="grib2" )then
1048 cfld=cfld+1
1049 fld_info(cfld)%ifld=iavblfld(iget(216))
1050 fld_info(cfld)%lvl=lvlsxml(lp,iget(216))
1051 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1052 endif
1053 ENDIF
1054 ENDIF
1055!
1056!*** SPECIFIC HUMIDITY.
1057!
1058 IF(iget(207)>0)THEN
1059 IF(lvls(lp,iget(207))>0)THEN
1060 DO j=jsta,jend
1061 DO i=ista,iend
1062 grid1(i,j)=qsl(i,j)
1063 ENDDO
1064 ENDDO
1065 CALL bound(grid1,h1m12,h99999)
1066 if(grib=="grib2" )then
1067 cfld=cfld+1
1068 fld_info(cfld)%ifld=iavblfld(iget(207))
1069 fld_info(cfld)%lvl=lvlsxml(lp,iget(207))
1070 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1071 endif
1072 ENDIF
1073 ENDIF
1074!
1075!*** OMEGA
1076!
1077 IF(iget(210)>0)THEN
1078 IF(lvls(lp,iget(210))>0)THEN
1079 DO j=jsta,jend
1080 DO i=ista,iend
1081 grid1(i,j)=osl(i,j)
1082 ENDDO
1083 ENDDO
1084 if(grib=="grib2" )then
1085 cfld=cfld+1
1086 fld_info(cfld)%ifld=iavblfld(iget(210))
1087 fld_info(cfld)%lvl=lvlsxml(lp,iget(210))
1088 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1089 endif
1090 ENDIF
1091 ENDIF
1092!
1093!*** U AND/OR V WIND
1094!
1095 IF(iget(208)>0.OR.iget(209)>0)THEN
1096 IF(lvls(lp,iget(208))>0.OR.lvls(lp,iget(209))>0) then
1097 DO j=jsta,jend
1098 DO i=ista,iend
1099 grid1(i,j)=usl(i,j)
1100 grid2(i,j)=vsl(i,j)
1101 ENDDO
1102 ENDDO
1103 if(grib=="grib2" )then
1104 cfld=cfld+1
1105 fld_info(cfld)%ifld=iavblfld(iget(208))
1106 fld_info(cfld)%lvl=lvlsxml(lp,iget(208))
1107 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1108 cfld=cfld+1
1109 fld_info(cfld)%ifld=iavblfld(iget(209))
1110 fld_info(cfld)%lvl=lvlsxml(lp,iget(209))
1111 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid2(ista:iend,jsta:jend)
1112 endif
1113 ENDIF
1114 ENDIF
1115!
1116!*** TURBULENT KINETIC ENERGY
1117!
1118 IF (iget(217)>0) THEN
1119 IF (lvls(lp,iget(217))>0) THEN
1120 DO j=jsta,jend
1121 DO i=ista,iend
1122 grid1(i,j)=q2sl(i,j)
1123 ENDDO
1124 ENDDO
1125 if(grib=="grib2" )then
1126 cfld=cfld+1
1127 fld_info(cfld)%ifld=iavblfld(iget(217))
1128 fld_info(cfld)%lvl=lvlsxml(lp,iget(217))
1129 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1130 endif
1131 ENDIF
1132 ENDIF
1133!
1134!*** CLOUD WATER
1135!
1136 IF (iget(211)>0) THEN
1137 IF (lvls(lp,iget(211))>0) THEN
1138 DO j=jsta,jend
1139 DO i=ista,iend
1140 grid1(i,j)=qw1(i,j)
1141 ENDDO
1142 ENDDO
1143 if(grib=="grib2" )then
1144 cfld=cfld+1
1145 fld_info(cfld)%ifld=iavblfld(iget(211))
1146 fld_info(cfld)%lvl=lvlsxml(lp,iget(211))
1147 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1148 endif
1149 ENDIF
1150 ENDIF
1151!
1152!*** CLOUD ICE
1153!
1154 IF (iget(212)>0) THEN
1155 IF (lvls(lp,iget(212))>0) THEN
1156 DO j=jsta,jend
1157 DO i=ista,iend
1158 grid1(i,j)=qi1(i,j)
1159 ENDDO
1160 ENDDO
1161 if(grib=="grib2" )then
1162 cfld=cfld+1
1163 fld_info(cfld)%ifld=iavblfld(iget(212))
1164 fld_info(cfld)%lvl=lvlsxml(lp,iget(212))
1165 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1166 endif
1167 ENDIF
1168 ENDIF
1169!
1170!--- RAIN
1171 IF (iget(213)>0) THEN
1172 IF (lvls(lp,iget(213))>0) THEN
1173 DO j=jsta,jend
1174 DO i=ista,iend
1175 grid1(i,j)=qr1(i,j)
1176 ENDDO
1177 ENDDO
1178 if(grib=="grib2" )then
1179 cfld=cfld+1
1180 fld_info(cfld)%ifld=iavblfld(iget(213))
1181 fld_info(cfld)%lvl=lvlsxml(lp,iget(213))
1182 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1183 endif
1184 ENDIF
1185 ENDIF
1186!
1187!--- SNOW
1188 IF (iget(214)>0) THEN
1189 IF (lvls(lp,iget(214))>0) THEN
1190 DO j=jsta,jend
1191 DO i=ista,iend
1192 grid1(i,j)=qs1(i,j)
1193 ENDDO
1194 ENDDO
1195 if(grib=="grib2" )then
1196 cfld=cfld+1
1197 fld_info(cfld)%ifld=iavblfld(iget(214))
1198 fld_info(cfld)%lvl=lvlsxml(lp,iget(214))
1199 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1200 endif
1201 ENDIF
1202 ENDIF
1203!
1204!--- GRAUPEL
1205 IF (iget(255)>0) THEN
1206 IF (lvls(lp,iget(255))>0) THEN
1207 DO j=jsta,jend
1208 DO i=ista,iend
1209 grid1(i,j)=qg1(i,j)
1210 ENDDO
1211 ENDDO
1212 if(grib=="grib2" )then
1213 cfld=cfld+1
1214 fld_info(cfld)%ifld=iavblfld(iget(255))
1215 fld_info(cfld)%lvl=lvlsxml(lp,iget(255))
1216 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1217 endif
1218 ENDIF
1219 ENDIF
1220!
1221!--- TOTAL CONDENSATE
1222 IF (iget(215)>0) THEN
1223 IF (lvls(lp,iget(215))>0) THEN
1224 DO j=jsta,jend
1225 DO i=ista,iend
1226 grid1(i,j)=c1d(i,j)
1227 ENDDO
1228 ENDDO
1229 if(grib=="grib2" )then
1230 cfld=cfld+1
1231 fld_info(cfld)%ifld=iavblfld(iget(215))
1232 fld_info(cfld)%lvl=lvlsxml(lp,iget(215))
1233 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1234 endif
1235 ENDIF
1236 ENDIF
1237!
1238! TOTAL CLOUD COVER
1239 IF (iget(222)>0) THEN
1240 IF (lvls(lp,iget(222))>0) THEN
1241 DO j=jsta,jend
1242 DO i=ista,iend
1243 grid1(i,j)=cfrsig(i,j)
1244 ENDDO
1245 ENDDO
1246 if(grib=="grib2" )then
1247 cfld=cfld+1
1248 fld_info(cfld)%ifld=iavblfld(iget(222))
1249 fld_info(cfld)%lvl=lvlsxml(lp,iget(222))
1250 datapd(1:iend-ista+1,1:jend-jsta+1,cfld)=grid1(ista:iend,jsta:jend)
1251 endif
1252 ENDIF
1253 ENDIF
1254!*** END OF MAIN VERTICAL LOOP
1255!
1256 310 CONTINUE
1257!*** ENDIF FOR IF TEST SEEING IF WE WANT ANY OTHER VARIABLES
1258!
1259 ENDIF
1260!
1261!
1262!
1263! END OF ROUTINE.
1264!
1265 RETURN
1266 END
bound
subroutine bound(fld, fmin, fmax)
Clips data in passed array.
Definition BOUND.f:37
exch
subroutine exch(a)
exch() Subroutine that exchanges one halo row.
Definition EXCH.f:27
mdl2sigma
subroutine mdl2sigma
SUBPROGRAM: MDL2P VERT INTRP OF MODEL LVLS TO PRESSURE PRGRMMR: BLACK ORG: W/NP22 DATE: 99-09-23
Definition MDL2SIGMA.f:55