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