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