UPP  V11.0.0
 All Data Structures Files Functions Pages
MDL2THANDPV.f
Go to the documentation of this file.
1 
22  SUBROUTINE mdl2thandpv(kth,kpv,th,pv)
23 
24 !
25 !
26  use vrbls3d, only: pmid, t, uh, q, vh, zmid, omga, pint
27  use vrbls2d, only: f
28  use masks, only: gdlat, gdlon, dx, dy
29  use physcons_post, only: con_eps, con_epsm1
30  use params_mod, only: dtr, small, erad, d608, rhmin
31  use ctlblk_mod, only: spval, lm, jsta_2l, jend_2u, grib, cfld, datapd, fld_info,&
32  im, jm, jsta, jend, jsta_m, jend_m, modelname, global,gdsdegr,me,&
33  ista, iend, ista_m, iend_m, ista_2l, iend_2u
34  use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml
35  use gridspec_mod, only: gridtype,dyval
36  use upp_physics, only: fpvsnew
37  use upp_math, only: dvdxdudy, ddvdx, ddudy, uuavg, h2u
38 !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39  implicit none
40 !
41 ! DECLARE VARIABLES.
42 !
43  integer,intent(in) :: kth, kpv
44  real, intent(in) :: th(kth), pv(kpv)
45  real, dimension(ista:iend,jsta:jend) :: grid1, grid2
46  real, dimension(kpv) :: pvpt, pvpb
47 
48  LOGICAL ioomg,ioall
49  LOGICAL lth(kth), lpv(kpv)
50 
51  REAL, allocatable:: dum1d1(:), dum1d2(:), dum1d3(:),dum1d4(:) &
52  , dum1d5(:), dum1d6(:), dum1d7(:),dum1d8(:) &
53  , dum1d9(:), dum1d10(:),dum1d11(:) &
54  , dum1d12(:),dum1d13(:),dum1d14(:)
55 !
56  real, dimension(ISTA:IEND,JSTA:JEND,KTH) :: uth, vth, hmth, tth, pvth, &
57  sigmath, rhth, oth
58  real, dimension(ISTA:IEND,JSTA:JEND,KPV) :: upv, vpv, hpv, tpv, ppv, spv
59  real, dimension(IM,2) :: glatpoles, coslpoles, pvpoles
60  real, dimension(IM,2,LM) :: upoles, tpoles, ppoles
61  real, dimension(IM,JSTA:JEND) :: cosltemp, pvtemp
62 !
63  real, allocatable :: wrk1(:,:), wrk2(:,:), wrk3(:,:), wrk4(:,:), cosl(:,:), dum2d(:,:)
64  real, allocatable :: tuv(:,:,:),pmiduv(:,:,:)
65 !
66  integer, dimension(im) :: iw, ie
67  integer i,j,l,k,lp,imb2,ip1,im1,ii,jj,jmt2,ihw,ihe
68  real dvdx,dudy,uavg,tphi, es, qstl, eradi, tem
69  real, allocatable :: dvdxl(:,:,:), dudyl(:,:,:), uavgl(:,:,:)
70 !
71 !
72 !******************************************************************************
73 !
74 ! START MDL2TH.
75 !
76  if(me==0) write(0,*) 'MDL2THANDPV starts'
77 !
78 ! SET TOTAL NUMBER OF POINTS ON OUTPUT GRID.
79 !
80 !---------------------------------------------------------------
81 !
82 ! *** PART I ***
83 !
84 ! VERTICAL INTERPOLATION OF EVERYTHING ELSE. EXECUTE ONLY
85 ! IF THERE'S SOMETHING WE WANT.
86 !
87  IF((iget(332) > 0).OR.(iget(333) > 0).OR. &
88  (iget(334) > 0).OR.(iget(335) > 0).OR. &
89  (iget(336) > 0).OR.(iget(337) > 0).OR. &
90  (iget(338) > 0).OR.(iget(339) > 0).OR. &
91  (iget(340) > 0).OR.(iget(341) > 0).OR. &
92  (iget(351) > 0).OR.(iget(352) > 0).OR. &
93  (iget(353) > 0).OR.(iget(378) > 0) ) THEN
94 !
95 !---------------------------------------------------------------------
96 !***
97 !*** BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND, DO ALL INTERPOLATION ABOVE GROUND NOW.
98 !***
99 !
100  do k=1,kpv
101  pvpt(k) = 5000.0 ! top limit for PV search
102  pvpb(k) = 0.8 ! Bottome limit for PV search in sigma
103  enddo
104 
105  do k=1,kth
106 !$omp parallel do private(i,j)
107  do j=jsta,jend
108  do i=ista,iend
109  uth(i,j,k) = spval
110  vth(i,j,k) = spval
111  hmth(i,j,k) = spval
112  tth(i,j,k) = spval
113  pvth(i,j,k) = spval
114  sigmath(i,j,k) = spval
115  rhth(i,j,k) = spval
116  oth(i,j,k) = spval
117  enddo
118  enddo
119  enddo
120  do k=1,kpv
121 !$omp parallel do private(i,j)
122  do j=jsta,jend
123  do i=ista,iend
124  upv(i,j,k) = spval
125  vpv(i,j,k) = spval
126  hpv(i,j,k) = spval
127  tpv(i,j,k) = spval
128  ppv(i,j,k) = spval
129  spv(i,j,k) = spval
130  enddo
131  enddo
132  enddo
133  ALLOCATE(dum1d1(lm), dum1d2(lm), dum1d3(lm),dum1d4(lm))
134  ALLOCATE(dum1d5(lm), dum1d6(lm)) !TV and Vorticity
135  ALLOCATE(dum1d7(lm), dum1d8(lm), dum1d9(lm),dum1d10(lm))
136  ALLOCATE(dum1d11(lm),dum1d12(lm),dum1d13(lm))
137  ALLOCATE(dum1d14(lm))
138 !
139  DO l=1,lm
140  CALL exch(pmid(ista_2l:iend_2u,jsta_2l:jend_2u,l))
141  CALL exch(t(ista_2l:iend_2u,jsta_2l:jend_2u,l))
142  CALL exch(uh(ista_2l:iend_2u,jsta_2l:jend_2u,l))
143  CALL exch(vh(ista_2l:iend_2u,jsta_2l:jend_2u,l))
144  END DO
145  CALL exch(gdlat(ista_2l,jsta_2l))
146  CALL exch(gdlon(ista_2l,jsta_2l))
147 
148 ! print *,' JSTA_2L=',JSTA_2L,' JSTA=',JSTA_2L,' JEND_2U=', &
149 ! &JEND_2U,' JEND=',JEND,' IM=',IM
150 ! print *,' GDLATa=',gdlat(1,:)
151 ! print *,' GDLATb=',gdlat(im,:)
152 !
153  allocate (wrk1(ista:iend,jsta:jend), wrk2(ista:iend,jsta:jend), &
154  & wrk3(ista:iend,jsta:jend), cosl(ista_2l:iend_2u,jsta_2l:jend_2u))
155  allocate (dum2d(ista_2l:iend_2u,jsta_2l:jend_2u))
156  allocate (wrk4(ista:iend,jsta:jend))
157 
158  imb2 = im /2
159  eradi = 1.0 / erad
160 
161 !! IF(MODELNAME == 'GFS' .or. global) THEN
162  IF(gridtype == 'A')THEN
163 !$omp parallel do private(i)
164  do i=1,im
165  ie(i) = i + 1
166  iw(i) = i - 1
167  enddo
168 ! iw(1) = im
169 ! ie(im) = 1
170 !
171 !$omp parallel do private(i,j,ip1,im1)
172  DO j=jsta,jend
173  do i=ista,iend
174  ip1 = ie(i)
175  im1 = iw(i)
176  cosl(i,j) = cos(gdlat(i,j)*dtr)
177  IF(cosl(i,j) >= small) then
178  wrk1(i,j) = eradi / cosl(i,j)
179  else
180  wrk1(i,j) = 0.
181  end if
182  if(i == im .or. i == 1)then
183  wrk2(i,j) = 1.0 / ((360.+gdlon(ip1,j)-gdlon(im1,j))*dtr)!1/dlam
184  else
185  wrk2(i,j) = 1.0 / ((gdlon(ip1,j)-gdlon(im1,j))*dtr) !1/dlam
186  end if
187  wrk4(i,j) = wrk1(i,j) * wrk2(i,j) ! 1/dx
188  enddo
189  enddo
190  CALL exch(cosl)
191 
192  call fullpole(cosl,coslpoles)
193  call fullpole(gdlat(ista_2l:iend_2u,jsta_2l:jend_2u),glatpoles)
194 
195 !$omp parallel do private(i,j,ii,tem)
196  DO j=jsta,jend
197  if (j == 1) then
198  do i=ista,iend
199  ii = i + imb2
200  if (ii > im) ii = ii - im
201  ! wrk3(i,j) = 1.0 / ((180.-GDLAT(i,J+1)-GDLAT(II,J))*DTR) !1/dphi
202  wrk3(i,j) = 1.0 / ((180.-gdlat(i,j+1)-glatpoles(ii,1))*dtr) !1/dphi
203  enddo
204  elseif (j == jm) then
205  do i=ista,iend
206  ii = i + imb2
207  if (ii > im) ii = ii - im
208  ! wrk3(i,j) = 1.0 / ((180.+GDLAT(i,J-1)+GDLAT(II,J))*DTR) !1/dphi
209  wrk3(i,j) = 1.0 / ((180.+gdlat(i,j-1)+glatpoles(ii,2))*dtr) !1/dphi
210  enddo
211  else
212  !print *,' j=',j,' GDLATJm1=',gdlat(:,j-1)
213  !print *,' j=',j,' GDLATJp1=',gdlat(:,j+1)
214  do i=ista,iend
215  tem = gdlat(i,j-1) - gdlat(i,j+1)
216  if (abs(tem) > small) then
217  wrk3(i,j) = 1.0 / (tem*dtr) !1/dphi
218  else
219  wrk3(i,j) = 0.0
220  endif
221  enddo
222  endif
223  !if (j == 181) print*,' wrk3=',wrk3(126,j),' gdlat=',&
224  ! GDLAT(126,J-1), gdlat(126,j+1)
225  enddo
226  else !!global?
227 !$omp parallel do private(i,j)
228  DO j=jsta_m,jend_m
229  DO i=ista_m,iend_m
230  wrk2(i,j) = 0.5 / dx(i,j)
231  wrk3(i,j) = 0.5 / dy(i,j)
232  END DO
233  END DO
234  endif
235 
236 ! need to put T and P on V points for computing dp/dx for e grid
237  IF(gridtype == 'E')THEN
238  allocate(tuv(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
239  allocate(pmiduv(ista_2l:iend_2u,jsta_2l:jend_2u,lm))
240  do l=1,lm
241  call h2u(t(ista_2l:iend_2u,jsta_2l:jend_2u,l),tuv(ista_2l:iend_2u,jsta_2l:jend_2u,l))
242  call h2u(pmid(ista_2l:iend_2u,jsta_2l:jend_2u,l),pmiduv(ista_2l:iend_2u,jsta_2l:jend_2u,l))
243  end do
244  end if
245 
246 !add A-grid regional models
247  IF(gridtype == 'A')THEN
248  IF(modelname == 'GFS' .or. global) THEN
249 
250  DO l=1,lm
251  CALL fullpole(pmid(ista_2l:iend_2u,jsta_2l:jend_2u,l),ppoles(:,:,l))
252  CALL fullpole( t(ista_2l:iend_2u,jsta_2l:jend_2u,l),tpoles(:,:,l))
253  CALL fullpole( uh(ista_2l:iend_2u,jsta_2l:jend_2u,l),upoles(:,:,l))
254  ENDDO
255 !!$omp parallel do private(i,j,ip1,im1,ii,jj,l,es,dum1d1,dum1d2,dum1d3,dum1d4,dum1d5,dum1d6,dum1d14,tem)
256  DO j=jsta,jend
257  DO i=ista,iend
258  ip1 = ie(i)
259  im1 = iw(i)
260  ii = i + imb2
261  if (ii > im) ii = ii - im
262 
263  IF(j == 1) then ! Near North pole
264  IF(cosl(i,j) >= small) THEN !not a pole point
265  tem = wrk3(i,j) * eradi
266 !$omp parallel do private(l,es)
267  DO l=1,lm
268  dum1d5(l) = t(i,j,l)*(1.+d608*q(i,j,l)) !Tv
269  es = min(fpvsnew(t(i,j,l)),pmid(i,j,l))
270  dum1d14(l) = q(i,j,l) * (pmid(i,j,l)+con_epsm1*es)/(con_eps*es) ! RH
271  dum1d1(l) = (pmid(ip1,j,l)- pmid(im1,j,l)) * wrk4(i,j) !dp/dx
272  dum1d3(l) = (t(ip1,j,l) - t(im1,j,l)) * wrk4(i,j) !dt/dx
273  ! DUM1D2(L) = (PMID(II,J,L) - PMID(I,J+1,L)) * tem !dp/dy
274  dum1d2(l) = (ppoles(ii,1,l) - pmid(i,j+1,l)) * tem !dp/dy
275  ! DUM1D4(L) = (T(II,J,L) - T(I,J+1,L)) * tem !dt/dy
276  dum1d4(l) = (tpoles(ii,1,l) - t(i,j+1,l)) * tem !dt/dy
277  dum1d6(l) = ((vh(ip1,j,l)-vh(im1,j,l))*wrk2(i,j) &
278  !& ! + (UH(II,J,L)*COSL(II,J) &
279  & + (upoles(ii,1,l)*coslpoles(ii,1) &
280  & + uh(i,j+1,l)*cosl(i,j+1))*wrk3(i,j))*wrk1(i,j) &
281  & + f(i,j)
282  END DO
283  ELSE !pole point, compute at j=2
284  jj = 2
285  tem = wrk3(i,jj) * eradi
286 !$omp parallel do private(l,es)
287  DO l=1,lm
288  dum1d5(l) = t(i,j,l)*(1.+d608*q(i,j,l)) !Tv
289  es = min(fpvsnew(t(i,j,l)),pmid(i,j,l))
290  dum1d14(l) = q(i,j,l) * (pmid(i,j,l)+con_epsm1*es)/(con_eps*es) ! RH
291  dum1d1(l) = (pmid(ip1,jj,l)- pmid(im1,jj,l)) * wrk4(i,jj) !dp/dx
292  dum1d3(l) = (t(ip1,jj,l) - t(im1,jj,l)) * wrk4(i,jj) !dt/dx
293  dum1d2(l) = (pmid(i,j,l)-pmid(i,jj+1,l)) * tem !dp/dy
294  dum1d4(l) = (t(i,j,l) - t(i,jj+1,l)) * tem !dt/dy
295  dum1d6(l) = ((vh(ip1,jj,l)-vh(im1,jj,l))*wrk2(i,jj) &
296  & + (uh(i,j,l)*cosl(i,j) &
297  + uh(i,jj+1,l)*cosl(i,jj+1))*wrk3(i,jj))*wrk1(i,jj) &
298  & + f(i,jj)
299  END DO
300  END IF
301  ELSE IF(j == jm) THEN ! Near South Pole
302  IF(cosl(i,j) >= small) THEN !not a pole point
303  tem = wrk3(i,j) * eradi
304 !$omp parallel do private(l,es)
305  DO l=1,lm
306  dum1d5(l) = t(i,j,l)*(1.+d608*q(i,j,l)) !Tv
307  es = min(fpvsnew(t(i,j,l)),pmid(i,j,l))
308  dum1d14(l) = q(i,j,l) * (pmid(i,j,l)+con_epsm1*es)/(con_eps*es) ! RH
309  dum1d1(l) = (pmid(ip1,j,l)- pmid(im1,j,l)) * wrk4(i,j) !dp/dx
310  dum1d3(l) = (t(ip1,j,l) - t(im1,j,l)) * wrk4(i,j) !dt/dx
311  ! DUM1D2(L) = (PMID(I,J-1,L)-PMID(II,J,L)) * tem !dp/dy
312  dum1d2(l) = (pmid(i,j-1,l)-ppoles(ii,2,l)) * tem !dp/dy
313  ! DUM1D4(L) = (T(I,J-1,L)-T(II,J,L)) * tem !dt/dy
314  dum1d4(l) = (t(i,j-1,l)-tpoles(ii,2,l)) * tem !dt/dy
315  dum1d6(l) = ((vh(ip1,j,l)-vh(im1,j,l))* wrk2(i,j) &
316  & + (uh(i,j-1,l)*cosl(i,j-1) &
317  !& ! + UH(II,J,L)*COSL(II,J))*wrk3(i,j))*wrk1(i,j) &
318  & + upoles(ii,2,l)*coslpoles(ii,2))*wrk3(i,j))*wrk1(i,j) &
319  & + f(i,j)
320  END DO
321  ELSE !pole point, compute at j=jm-1
322  jj = jm-1
323  tem = wrk3(i,jj) * eradi
324 !$omp parallel do private(l,es)
325  DO l=1,lm
326  dum1d5(l) = t(i,j,l)*(1.+d608*q(i,j,l)) !Tv
327  es = min(fpvsnew(t(i,j,l)),pmid(i,j,l))
328  dum1d14(l) = q(i,j,l) * (pmid(i,j,l)+con_epsm1*es)/(con_eps*es) ! RH
329  dum1d1(l) = (pmid(ip1,jj,l)- pmid(im1,jj,l)) * wrk4(i,jj) !dp/dx
330  dum1d3(l) = (t(ip1,jj,l) - t(im1,jj,l)) * wrk4(i,jj) !dt/dx
331  dum1d2(l) = (pmid(i,jj-1,l)-pmid(i,j,l)) * tem !dp/dy
332  dum1d4(l) = (t(i,jj-1,l)-t(i,j,l)) * tem !dt/dy
333  dum1d6(l) = ((vh(ip1,jj,l)-vh(im1,jj,l))*wrk2(i,jj) &
334  & + (uh(i,jj-1,l)*cosl(i,jj-1) &
335  & + uh(i,j,l)*cosl(i,j))*wrk3(i,jj))*wrk1(i,jj) &
336  & + f(i,jj)
337  END DO
338  END IF
339  ELSE
340  tem = wrk3(i,j) * eradi
341 !$omp parallel do private(l,es)
342  DO l=1,lm
343  dum1d5(l) = t(i,j,l)*(1.+d608*q(i,j,l)) !Tv
344  es = min(fpvsnew(t(i,j,l)),pmid(i,j,l))
345  dum1d14(l) = q(i,j,l) * (pmid(i,j,l)+con_epsm1*es)/(con_eps*es) ! RH
346  dum1d1(l) = (pmid(ip1,j,l)- pmid(im1,j,l)) * wrk4(i,j) !dp/dx
347  dum1d3(l) = (t(ip1,j,l) - t(im1,j,l)) * wrk4(i,j) !dt/dx
348 ! if (j >= 181) print *,' i=',i,' tem=',tem,' pmid=',pmid(i,j-1,l)&
349 ! ,pmid(i,j-1,l),' l=',l,' j=',j
350  dum1d2(l) = (pmid(i,j-1,l)-pmid(i,j+1,l)) * tem !dp/dy
351  dum1d4(l) = (t(i,j-1,l)-t(i,j+1,l)) * tem !dt/dy
352  dum1d6(l) = ((vh(ip1,j,l)-vh(im1,j,l))* wrk2(i,j) &
353  & - (uh(i,j-1,l)*cosl(i,j-1) &
354  & - uh(i,j+1,l)*cosl(i,j+1))*wrk3(i,j))*wrk1(i,j) &
355  & + f(i,j)
356  END DO
357  END IF
358 
359 
360  IF(i==im/2 .AND. j==jm/2)then
361  print*,'SAMPLE PVETC INPUT ', &
362  'p,dpdx,dpdy,tv,dtdx,dtdy,h,u,v,vort= '
363  DO l=1,lm
364  print*,pmid(i,j,l),dum1d1(l),dum1d2(l),dum1d5(l) &
365  ,dum1d3(l),dum1d4(l),zmid(i,j,l),uh(i,j,l),vh(i,j,l) &
366  ,dum1d6(l),l
367  end do
368  end if
369 
370  CALL pvetc(lm,pmid(i,j,1:lm),dum1d1,dum1d2 &
371  ,dum1d5,dum1d3,dum1d4,zmid(i,j,1:lm),uh(i,j,1:lm) &
372  ,vh(i,j,1:lm),dum1d6 &
373  ,dum1d7,dum1d8,dum1d9,dum1d10,dum1d11,dum1d12,dum1d13)!output
374 
375  IF(i==im/2 .AND. j==jm/2)then
376  print*,'SAMPLE PVETC OUTPUT ' &
377  ,'hm,s,bvf2,pvn,theta,sigma,pvu= '
378  DO l=1,lm
379  print*,dum1d7(l),dum1d8(l),dum1d9(l),dum1d10(l),dum1d11(l) &
380  ,dum1d12(l),dum1d13(l),l
381  end do
382  end if
383 
384  IF((iget(332) > 0).OR.(iget(333) > 0).OR. &
385  (iget(334) > 0).OR.(iget(335) > 0).OR. &
386  (iget(351) > 0).OR.(iget(352) > 0).OR. &
387  (iget(353) > 0).OR.(iget(378) > 0))THEN
388 ! interpolate to isentropic levels
389  CALL p2th(lm,dum1d11,uh(i,j,1:lm),vh(i,j,1:lm) &
390  ,dum1d7,t(i,j,1:lm),dum1d13,dum1d12,dum1d14 &
391  ,omga(i,j,1:lm),kth,th &
392  ,lth,uth(i,j,1:kth),vth(i,j,1:kth) &
393 !output
394  ,hmth(i,j,1:kth) &
395  ,tth(i,j,1:kth),pvth(i,j,1:kth) &
396  ,sigmath(i,j,1:kth),rhth(i,j,1:kth) &
397  ,oth(i,j,1:kth))!output
398  END IF
399 ! interpolate to PV levels
400  IF((iget(336) > 0).OR.(iget(337) > 0).OR. &
401  (iget(338) > 0).OR.(iget(339) > 0).OR. &
402  (iget(340) > 0).OR.(iget(341) > 0)) THEN
403  CALL p2pv(lm,dum1d13,zmid(i,j,1:lm),t(i,j,1:lm),pmid(i,j,1:lm) &
404  ,uh(i,j,1:lm),vh(i,j,1:lm),kpv,pv,pvpt,pvpb*pint(i,j,lm+1)&
405  ,lpv,upv(i,j,1:kpv),vpv(i,j,1:kpv),hpv(i,j,1:kpv) &
406 !output
407  ,tpv(i,j,1:kpv),ppv(i,j,1:kpv),spv(i,j,1:kpv) ) !output
408  END IF
409 
410  END DO
411  END DO
412 
413 !-----------------------------------------------------------------
414 !add A-grid regional models
415  ELSE !regional models start here (GFS or global ends here)
416  DO j=jsta_m,jend_m
417  jmt2=jm/2+1
418  tphi=(j-jmt2)*(dyval/gdsdegr)*dtr
419  DO i=ista_m,iend_m
420  ip1 = i + 1
421  im1 = i - 1
422  tem = wrk3(i,j) * eradi
423 !$omp parallel do private(l,es)
424  DO l=1,lm
425  dum1d5(l) = t(i,j,l)*(1.+d608*q(i,j,l)) !Tv
426  es = min(fpvsnew(t(i,j,l)),pmid(i,j,l))
427  dum1d14(l) = q(i,j,l) * (pmid(i,j,l)+con_epsm1*es)/(con_eps*es) ! RH
428  dum1d1(l) = (pmid(ip1,j,l)- pmid(im1,j,l)) * wrk4(i,j) !dp/dx
429  dum1d3(l) = (t(ip1,j,l) - t(im1,j,l)) * wrk4(i,j) !dt/dx
430  dum1d2(l) = (pmid(i,j+1,l)-pmid(i,j-1,l)) * tem !dp/dy
431  dum1d4(l) = (t(i,j+1,l)-t(i,j-1,l)) * tem !dt/dy
432  dum1d6(l) = ((vh(ip1,j,l)-vh(im1,j,l))* wrk2(i,j) &
433  & - (uh(i,j+1,l)*cosl(i,j+1) &
434  & - uh(i,j-1,l)*cosl(i,j-1))*wrk3(i,j))*wrk1(i,j) &
435  & + f(i,j)
436  END DO
437 
438  IF(i==im/2 .AND. j==jm/2)then
439  print*,'SAMPLE PVETC INPUT for regional ', &
440  'p,dpdx,dpdy,tv,dtdx,dtdy,h,u,v,vort ', &
441  'JSTA_m,Jend_m, L= '
442  DO l=1,lm
443  print*,pmid(i,j,l),dum1d1(l),dum1d2(l),dum1d5(l) &
444  ,dum1d3(l),dum1d4(l),zmid(i,j,l),uh(i,j,l),vh(i,j,l) &
445  ,dum1d6(l),jsta_m,jend_m,l
446  end do
447  end if
448 
449  CALL pvetc(lm,pmid(i,j,1:lm),dum1d1,dum1d2 &
450  ,dum1d5,dum1d3,dum1d4,zmid(i,j,1:lm),uh(i,j,1:lm) &
451  ,vh(i,j,1:lm),dum1d6 &
452  ,dum1d7,dum1d8,dum1d9,dum1d10,dum1d11,dum1d12,dum1d13)!output
453 
454  IF(i==im/2 .AND. j==jm/2)then
455  print*,'SAMPLE PVETC OUTPUT ' &
456  ,'hm,s,bvf2,pvn,theta,sigma,pvu,pvort= '
457  DO l=1,lm
458  print*,dum1d7(l),dum1d8(l),dum1d9(l),dum1d10(l),dum1d11(l) &
459  ,dum1d12(l),dum1d13(l),dum1d6(l),l
460  end do
461  end if
462 
463  IF((iget(332) > 0).OR.(iget(333) > 0).OR. &
464  (iget(334) > 0).OR.(iget(335) > 0).OR. &
465  (iget(351) > 0).OR.(iget(352) > 0).OR. &
466  (iget(353) > 0).OR.(iget(378) > 0))THEN
467 ! interpolate to isentropic levels
468  CALL p2th(lm,dum1d11,uh(i,j,1:lm),vh(i,j,1:lm) &
469  ,dum1d7,t(i,j,1:lm),dum1d13,dum1d12,dum1d14 &
470  ,omga(i,j,1:lm),kth,th &
471  ,lth,uth(i,j,1:kth),vth(i,j,1:kth) &
472 !output
473  ,hmth(i,j,1:kth) &
474  ,tth(i,j,1:kth),pvth(i,j,1:kth) &
475  ,sigmath(i,j,1:kth),rhth(i,j,1:kth) &
476  ,oth(i,j,1:kth))!output
477  END IF
478 ! interpolate to PV levels
479  IF((iget(336) > 0).OR.(iget(337) > 0).OR. &
480  (iget(338) > 0).OR.(iget(339) > 0).OR. &
481  (iget(340) > 0).OR.(iget(341) > 0)) THEN
482  CALL p2pv(lm,dum1d13,zmid(i,j,1:lm),t(i,j,1:lm),pmid(i,j,1:lm) &
483  ,uh(i,j,1:lm),vh(i,j,1:lm),kpv,pv,pvpt,pvpb*pint(i,j,lm+1)&
484  ,lpv,upv(i,j,1:kpv),vpv(i,j,1:kpv),hpv(i,j,1:kpv) &
485 !output
486  ,tpv(i,j,1:kpv),ppv(i,j,1:kpv),spv(i,j,1:kpv) ) !output
487  END IF
488 
489  ENDDO
490  ENDDO
491 
492  ENDIF !regional models and A-grid end here
493 !-----------------------------------------------------------------
494  ELSE IF (gridtype == 'B')THEN
495  allocate(dvdxl(ista_m:iend_m,jsta_m:jend_m,lm))
496  allocate(dudyl(ista_m:iend_m,jsta_m:jend_m,lm))
497  allocate(uavgl(ista_m:iend_m,jsta_m:jend_m,lm))
498  DO l=1,lm
499  CALL exch(vh(ista_2l:iend_2u,jsta_2l:jend_2u,l))
500  CALL exch(uh(ista_2l:iend_2u,jsta_2l:jend_2u,l))
501  CALL dvdxdudy(uh(:,:,l),vh(:,:,l))
502  DO j=jsta_m,jend_m
503  DO i=ista_m,iend_m
504  dvdxl(i,j,l) = ddvdx(i,j)
505  dudyl(i,j,l) = ddudy(i,j)
506  uavgl(i,j,l) = uuavg(i,j)
507  END DO
508  END DO
509  END DO
510  DO j=jsta_m,jend_m
511  jmt2=jm/2+1
512  tphi=(j-jmt2)*(dyval/gdsdegr)*dtr
513  DO i=ista_m,iend_m
514  ip1 = i + 1
515  im1 = i - 1
516  DO l=1,lm
517  dum1d5(l) = t(i,j,l)*(1.+d608*q(i,j,l)) !TV
518  es = min(fpvsnew(t(i,j,l)),pmid(i,j,l))
519  qstl = con_eps*es/(pmid(i,j,l)+con_epsm1*es)
520  dum1d14(l) = q(i,j,l)/qstl !RH
521  dum1d1(l) = (pmid(ip1,j,l)- pmid(im1,j,l)) * wrk2(i,j) !dp/dx
522  dum1d3(l) = (t(ip1,j,l) - t(im1,j,l)) * wrk2(i,j) !dt/dx
523  dum1d2(l) = (pmid(i,j+1,l)-pmid(i,j-1,l)) * wrk3(i,j) !dp/dy
524  dum1d4(l) = (t(i,j+1,l)-t(i,j-1,l)) * wrk3(i,j) !dt/dy
525  dvdx = dvdxl(i,j,l)
526  dudy = dudyl(i,j,l)
527  uavg = uavgl(i,j,l)
528 ! is there a (f+tan(phi)/erad)*u term?
529  dum1d6(l) = dvdx - dudy + f(i,j) + uavg*tan(tphi)/erad !vort
530 
531  END DO
532 
533 ! IF(I==IM/2 .AND. J==JM/2)then
534 ! PRINT*,'SAMPLE PVETC INPUT ' &
535 ! ,'p,dpdx,dpdy,tv,dtdx,dtdy,h,u,v,vort= '
536 ! DO L=1,LM
537 ! print*,pmid(i,j,l),dum1d1(l),dum1d2(l),dum1d5(l) &
538 ! ,dum1d3(l),dum1d4(l),zmid(i,j,l),uh(i,j,l),vh(i,j,l) &
539 ! ,dum1d6(l)
540 ! end do
541 ! end if
542 
543  CALL pvetc(lm,pmid(i,j,1:lm),dum1d1,dum1d2 &
544  ,dum1d5,dum1d3,dum1d4,zmid(i,j,1:lm),uh(i,j,1:lm) &
545  ,vh(i,j,1:lm),dum1d6 &
546  ,dum1d7,dum1d8,dum1d9,dum1d10,dum1d11,dum1d12,dum1d13)!output
547 
548 ! IF(I==IM/2 .AND. J==JM/2)then
549 ! PRINT*,'SAMPLE PVETC OUTPUT ' &
550 ! ,'hm,s,bvf2,pvn,theta,sigma,pvu= '
551 ! DO L=1,LM
552 ! print*,dum1d7(l),dum1d8(l),dum1d9(l),dum1d10(l),dum1d11(l) &
553 ! ,dum1d12(l),dum1d13(l)
554 ! end do
555 ! end if
556  IF((iget(332)>0).OR.(iget(333)>0).OR. &
557  (iget(334)>0).OR.(iget(335)>0).OR. &
558  (iget(351)>0).OR.(iget(352)>0).OR. &
559  (iget(353)>0).OR.(iget(378)>0))THEN
560 ! interpolate to isentropic levels
561  CALL p2th(lm,dum1d11,uh(i,j,1:lm),vh(i,j,1:lm) &
562  ,dum1d7,t(i,j,1:lm),dum1d13,dum1d12,dum1d14 &
563  ,omga(i,j,1:lm),kth,th &
564  ,lth,uth(i,j,1:kth),vth(i,j,1:kth) &
565 !output
566  ,hmth(i,j,1:kth) &
567  ,tth(i,j,1:kth),pvth(i,j,1:kth) &
568  ,sigmath(i,j,1:kth),rhth(i,j,1:kth) &
569  ,oth(i,j,1:kth))!output
570  END IF
571 ! interpolate to PV levels
572  IF((iget(336)>0).OR.(iget(337)>0).OR. &
573  (iget(338)>0).OR.(iget(339)>0).OR. &
574  (iget(340)>0).OR.(iget(341)>0))THEN
575  CALL p2pv(lm,dum1d13,zmid(i,j,1:lm),t(i,j,1:lm),pmid(i,j,1:lm) &
576  ,uh(i,j,1:lm),vh(i,j,1:lm),kpv,pv,pvpt,pvpb*pint(i,j,lm+1) &
577  ,lpv,upv(i,j,1:kpv),vpv(i,j,1:kpv),hpv(i,j,1:kpv) &
578 !output
579  ,tpv(i,j,1:kpv),ppv(i,j,1:kpv),spv(i,j,1:kpv) ) !output
580  END IF
581  END DO
582  END DO
583  deallocate(dvdxl,dudyl,uavgl)
584  ELSE IF (gridtype == 'E')THEN
585  DO j=jsta_m,jend_m
586  jmt2 = jm/2+1
587  tphi = (j-jmt2)*(dyval/gdsdegr)*dtr
588  ihw= - mod(j,2)
589  ihe = ihw + 1
590  DO i=ista_m,iend_m
591  ip1 = i + 1
592  im1 = i - 1
593  DO l=1,lm
594  dum1d5(l)=t(i,j,l)*(1.+d608*q(i,j,l)) !TV
595  es=fpvsnew(t(i,j,l))
596  es=min(es,pmid(i,j,l))
597  qstl=con_eps*es/(pmid(i,j,l)+con_epsm1*es)
598  dum1d14(l)=q(i,j,l)/qstl !RH
599  dum1d1(l) = (pmiduv(i+ihe,j,l)- pmiduv(i+ihw,j,l))*wrk2(i,j) !dp/dx
600  dum1d3(l) = (tuv(i+ihe,j,l) - tuv(i+ihw,j,l)) * wrk2(i,j) !dt/dx
601  dum1d2(l) = (pmiduv(i,j+1,l)- pmiduv(i,j-1,l))*wrk3(i,j)!dp/dy
602  dum1d4(l)= (tuv(i,j+1,l)- tuv(i,j-1,l))*wrk3(i,j)!dt/dy
603  dvdx=(vh(i+ihe,j,l)-vh(i+ihw,j,l))* wrk2(i,j)
604  dudy=(uh(i,j+1,l)-uh(i,j-1,l))* wrk3(i,j)
605  uavg=0.25*(uh(i+ihe,j,l)+uh(i+ihw,j,l)+uh(i,j-1,l)+uh(i,j+1,l))
606 ! is there a (f+tan(phi)/erad)*u term?
607  dum1d6(l)=dvdx-dudy+f(i,j)+uavg*tan(tphi)/erad !vort
608  END DO
609 
610  IF(i==im/2 .AND. j==jm/2)then
611  print*,'SAMPLE PVETC INPUT ' &
612  ,'p,dpdx,dpdy,tv,dtdx,dtdy,h,u,v,vort= '
613  DO l=1,lm
614  print*,pmid(i,j,l),dum1d1(l),dum1d2(l),dum1d5(l) &
615  ,dum1d3(l),dum1d4(l),zmid(i,j,l),uh(i,j,l),vh(i,j,l) &
616  ,dum1d6(l)
617  end do
618  end if
619 
620  CALL pvetc(lm,pmid(i,j,1:lm),dum1d1,dum1d2 &
621  ,dum1d5,dum1d3,dum1d4,zmid(i,j,1:lm),uh(i,j,1:lm) &
622  ,vh(i,j,1:lm),dum1d6 &
623  ,dum1d7,dum1d8,dum1d9,dum1d10,dum1d11,dum1d12,dum1d13)!output
624 
625  IF(i==im/2 .AND. j==jm/2)then
626  print*,'SAMPLE PVETC OUTPUT ' &
627  ,'hm,s,bvf2,pvn,theta,sigma,pvu= '
628  DO l=1,lm
629  print*,dum1d7(l),dum1d8(l),dum1d9(l),dum1d10(l),dum1d11(l) &
630  ,dum1d12(l),dum1d13(l)
631  end do
632  end if
633  IF((iget(332) > 0).OR.(iget(333) > 0).OR. &
634  (iget(334) > 0).OR.(iget(335) > 0).OR. &
635  (iget(351) > 0).OR.(iget(352) > 0).OR. &
636  (iget(353) > 0).OR.(iget(378) > 0))THEN
637 ! interpolat e to isentropic levels
638  CALL p2th(lm,dum1d11,uh(i,j,1:lm),vh(i,j,1:lm) &
639  ,dum1d7,t(i,j,1:lm),dum1d13,dum1d12,dum1d14 &
640  ,omga(i,j,1:lm),kth,th &
641  ,lth,uth(i,j,1:kth),vth(i,j,1:kth) &
642 !output
643  ,hmth(i,j,1:kth) &
644  ,tth(i,j,1:kth),pvth(i,j,1:kth) &
645  ,sigmath(i,j,1:kth),rhth(i,j,1:kth) &
646  ,oth(i,j,1:kth))!output
647  END IF
648 ! interpolate to PV levels
649  IF((iget(336) > 0) .OR. (iget(337) > 0).OR. &
650  (iget(338) > 0) .OR. (iget(339) > 0).OR. &
651  (iget(340) > 0) .OR. (iget(341) > 0)) THEN
652  CALL p2pv(lm,dum1d13,zmid(i,j,1:lm),t(i,j,1:lm),pmid(i,j,1:lm) &
653  ,uh(i,j,1:lm),vh(i,j,1:lm),kpv,pv,pvpt,pvpb*pint(i,j,lm+1) &
654  ,lpv,upv(i,j,1:kpv),vpv(i,j,1:kpv),hpv(i,j,1:kpv) &
655 !output
656  ,tpv(i,j,1:kpv),ppv(i,j,1:kpv),spv(i,j,1:kpv) ) !output
657  END IF
658  END DO
659  END DO
660 
661  END IF ! for different grids
662 
663 
664 ! print *,'in mdlthpv,af P2PV,tpv=',maxval(TPV(I,J,1:KPV)), &
665 ! minval(TPV(I,J,1:KPV)),'kpv=',kpv,'zmid=',maxval(ZMID(1:im,jsta:jend,1:LM)), &
666 ! minval(ZMID(1:im,jsta:jend,1:LM)),'uth=',maxval(uth(1:im,jsta:jend,1:kth)), &
667 ! minval(uth(1:im,jsta:jend,1:kth)),'vth=',maxval(vth(1:im,jsta:jend,1:kth)), &
668 ! minval(vth(1:im,jsta:jend,1:kth)),'uth(1,1,1)=',uth(1,1,1:kth)
669 !---------------------------------------------------------------------
670 !---------------------------------------------------------------------
671 ! *** PART II *** Write out fields
672 !---------------------------------------------------------------------
673 !
674 !
675  DO lp=1,kth
676 !*** ISENTROPIC U AND/OR V WIND
677 !
678 
679  IF(iget(332) > 0 .OR. iget(333) > 0)THEN
680  IF(lvls(lp,iget(332)) > 0 .OR. lvls(lp,iget(333)) > 0)THEN
681 !$omp parallel do private(i,j)
682  DO j=jsta,jend
683  DO i=ista,iend
684  grid1(i,j) = uth(i,j,lp)
685  grid2(i,j) = vth(i,j,lp)
686  ENDDO
687  ENDDO
688  if(grib=='grib2')then
689  cfld = cfld + 1
690  fld_info(cfld)%ifld = iavblfld(iget(332))
691  fld_info(cfld)%lvl = lvlsxml(lp,iget(332))
692 !$omp parallel do private(i,j,ii,jj)
693  do j=1,jend-jsta+1
694  jj = jsta+j-1
695  do i=1,iend-ista+1
696  ii = ista+i-1
697  datapd(i,j,cfld) = grid1(ii,jj)
698  enddo
699  enddo
700  cfld = cfld + 1
701  fld_info(cfld)%ifld = iavblfld(iget(333))
702  fld_info(cfld)%lvl = lvlsxml(lp,iget(333))
703 !$omp parallel do private(i,j,ii,jj)
704  do j=1,jend-jsta+1
705  jj = jsta+j-1
706  do i=1,iend-ista+1
707  ii = ista+i-1
708  datapd(i,j,cfld) = grid2(ii,jj)
709  enddo
710  enddo
711  endif
712  ENDIF
713  ENDIF
714 
715 !*** OUTPUT ISENTROPIC T
716 !
717  IF(iget(334) > 0)THEN
718  IF(lvls(lp,iget(334)) > 0)THEN
719 !!$omp parallel do
720 ! GFS use lon avg as one scaler value for pole point
721 ! JJ=1
722 ! IF(JJ>=jsta .and. JJ<=jend .and. cosl(1,JJ)<SMALL)then
723 ! WORK=0.
724 ! DO I=1,IM
725 ! WORK=TTH(I,JJ,LP)+WORK
726 ! END DO
727 ! DO I=1,IM
728 ! TTH(I,JJ,LP)=WORK/IM
729 ! END DO
730 ! END IF
731 ! JJ=JM
732 ! IF(JJ>=jsta .and. JJ<=jend .and. cosl(1,JJ)<SMALL)then
733 ! WORK=0.
734 ! DO I=1,IM
735 ! WORK=TTH(I,JJ,LP)+WORK
736 ! END DO
737 ! DO I=1,IM
738 ! TTH(I,JJ,LP)=WORK/IM
739 ! END DO
740 ! END IF
741 !$omp parallel do private(i,j)
742  DO j=jsta,jend
743  DO i=ista,iend
744  grid1(i,j) = tth(i,j,lp)
745  ENDDO
746  ENDDO
747  if(grib=='grib2')then
748  cfld = cfld + 1
749  fld_info(cfld)%ifld=iavblfld(iget(334))
750  fld_info(cfld)%lvl=lvlsxml(lp,iget(334))
751 !$omp parallel do private(i,j,ii,jj)
752  do j=1,jend-jsta+1
753  jj = jsta+j-1
754  do i=1,iend-ista+1
755  ii = ista+i-1
756  datapd(i,j,cfld) = grid1(ii,jj)
757  enddo
758  enddo
759  endif
760  ENDIF
761  ENDIF
762 !
763 !*** ISENTROPIC PV
764 !
765  IF(iget(335) > 0) THEN
766  IF(lvls(lp,iget(335)) > 0)THEN
767  ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1:IM,JSTA:JEND) &
768  ! ,SPVAL,PVTH(1:IM,JSTA:JEND,LP))
769  ! IF(1>=jsta .and. 1<=jend)print*,'PVTH at N POLE= ' &
770  ! ,pvth(1,1,lp),pvth(im/2,1,lp) &
771  ! ,pvth(10,10,lp),pvth(im/2,10,lp),SPVAL,grib,LP
772  dum2d(ista:iend,jsta:jend)=pvth(ista:iend,jsta:jend,lp)
773  CALL exch(dum2d)
774  CALL fullpole(dum2d,pvpoles)
775  cosltemp=spval
776  IF(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
777  IF(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
778  pvtemp=spval
779  IF(jsta== 1) pvtemp(1:im, 1)=pvpoles(1:im,1)
780  IF(jend==jm) pvtemp(1:im,jm)=pvpoles(1:im,2)
781 
782  call poleavg(im,jm,jsta,jend,small,cosltemp(1:im,jsta:jend) &
783  ,spval,pvtemp(1:im,jsta:jend))
784 
785  IF(jsta== 1) pvth(ista:iend, 1,lp)=pvtemp(ista:iend, 1)
786  IF(jend==jm) pvth(ista:iend,jm,lp)=pvtemp(ista:iend,jm)
787 
788 !$omp parallel do private(i,j)
789  DO j=jsta,jend
790  DO i=ista,iend
791  IF(pvth(i,j,lp) /= spval)THEN
792  grid1(i,j) = pvth(i,j,lp)*1.0e-6
793  ELSE
794  grid1(i,j) = pvth(i,j,lp)
795  END IF
796  ENDDO
797  ENDDO
798  if(grib=='grib2')then
799  cfld=cfld+1
800  fld_info(cfld)%ifld=iavblfld(iget(335))
801  fld_info(cfld)%lvl=lvlsxml(lp,iget(335))
802 !$omp parallel do private(i,j,ii,jj)
803  do j=1,jend-jsta+1
804  jj = jsta+j-1
805  do i=1,iend-ista+1
806  ii = ista+i-1
807  datapd(i,j,cfld) = grid1(ii,jj)
808  enddo
809  enddo
810  endif
811  ENDIF
812  ENDIF
813 !
814 !*** ISENTROPIC Montgomery function
815 !
816  IF(iget(353) > 0) THEN
817  IF(lvls(lp,iget(353)) > 0)THEN
818 !$omp parallel do private(i,j)
819  DO j=jsta,jend
820  DO i=ista,iend
821  grid1(i,j) = hmth(i,j,lp)
822  ENDDO
823  ENDDO
824  if(grib=='grib2')then
825  cfld=cfld+1
826  fld_info(cfld)%ifld=iavblfld(iget(353))
827  fld_info(cfld)%lvl=lvlsxml(lp,iget(353))
828 !$omp parallel do private(i,j,ii,jj)
829  do j=1,jend-jsta+1
830  jj = jsta+j-1
831  do i=1,iend-ista+1
832  ii = ista+i-1
833  datapd(i,j,cfld) = grid1(ii,jj)
834  enddo
835  enddo
836  endif
837  ENDIF
838  ENDIF
839 !
840 !*** ISENTROPIC static stability
841 !
842  IF(iget(351) > 0) THEN
843  IF(lvls(lp,iget(351)) > 0)THEN
844 !$omp parallel do private(i,j)
845  DO j=jsta,jend
846  DO i=ista,iend
847  grid1(i,j) = sigmath(i,j,lp)
848  ENDDO
849  ENDDO
850  if(grib=='grib2') then
851  cfld=cfld+1
852  fld_info(cfld)%ifld=iavblfld(iget(351))
853  fld_info(cfld)%lvl=lvlsxml(lp,iget(351))
854 !$omp parallel do private(i,j,ii,jj)
855  do j=1,jend-jsta+1
856  jj = jsta+j-1
857  do i=1,iend-ista+1
858  ii = ista+i-1
859  datapd(i,j,cfld) = grid1(ii,jj)
860  enddo
861  enddo
862  endif
863  ENDIF
864  ENDIF
865 !
866 !*** ISENTROPIC RH
867 !
868  IF(iget(352) > 0) THEN
869  IF(lvls(lp,iget(352)) > 0)THEN
870 !$omp parallel do private(i,j)
871  DO j=jsta,jend
872  DO i=ista,iend
873  IF(rhth(i,j,lp) /= spval) THEN
874  grid1(i,j) = 100.0 * min(1.,max(rhmin,rhth(i,j,lp)))
875  ELSE
876  grid1(i,j) = spval
877  END IF
878  ENDDO
879  ENDDO
880  if(grib=='grib2') then
881  cfld=cfld+1
882  fld_info(cfld)%ifld=iavblfld(iget(352))
883  fld_info(cfld)%lvl=lvlsxml(lp,iget(352))
884 !$omp parallel do private(i,j,ii,jj)
885  do j=1,jend-jsta+1
886  jj = jsta+j-1
887  do i=1,iend-ista+1
888  ii = ista+i-1
889  datapd(i,j,cfld) = grid1(ii,jj)
890  enddo
891  enddo
892  endif
893  ENDIF
894  ENDIF
895 !
896 !*** ISENTROPIC OMEGA
897 !
898  IF(iget(378) > 0) THEN
899  IF(lvls(lp,iget(378)) > 0)THEN
900 !$omp parallel do private(i,j)
901  DO j=jsta,jend
902  DO i=ista,iend
903  grid1(i,j) = oth(i,j,lp)
904  ENDDO
905  ENDDO
906  if(grib=='grib2') then
907  cfld=cfld+1
908  fld_info(cfld)%ifld=iavblfld(iget(378))
909  fld_info(cfld)%lvl=lvlsxml(lp,iget(378))
910 !$omp parallel do private(i,j,ii,jj)
911  do j=1,jend-jsta+1
912  jj = jsta+j-1
913  do i=1,iend-ista+1
914  ii = ista+i-1
915  datapd(i,j,cfld) = grid1(ii,jj)
916  enddo
917  enddo
918  endif
919  ENDIF
920  ENDIF
921  END DO ! end loop for isentropic levels
922 ! Lopp through PV levels now
923  DO lp=1,kpv
924 !*** U AND/OR V WIND on PV surface
925 !
926  IF(iget(336) > 0.OR.iget(337) > 0)THEN
927  IF(lvls(lp,iget(336)) > 0.OR.lvls(lp,iget(337)) > 0)THEN
928 ! GFS use lon avg as one scaler value for pole point
929  ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1:IM,JSTA:JEND) &
930  ! ,SPVAL,VPV(1:IM,JSTA:JEND,LP))
931  dum2d(ista:iend,jsta:jend)=vpv(ista:iend,jsta:jend,lp)
932  CALL exch(dum2d)
933  CALL fullpole(dum2d,pvpoles)
934  cosltemp=spval
935  IF(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
936  IF(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
937  pvtemp=spval
938  IF(jsta== 1) pvtemp(1:im, 1)=pvpoles(1:im,1)
939  IF(jend==jm) pvtemp(1:im,jm)=pvpoles(1:im,2)
940 
941  call poleavg(im,jm,jsta,jend,small,cosltemp(1:im,jsta:jend) &
942  ,spval,pvtemp(1:im,jsta:jend))
943 
944  IF(jsta== 1) vpv(ista:iend, 1,lp)=pvtemp(ista:iend, 1)
945  IF(jend==jm) vpv(ista:iend,jm,lp)=pvtemp(ista:iend,jm)
946 
947 !$omp parallel do private(i,j)
948  DO j=jsta,jend
949  DO i=ista,iend
950  grid1(i,j) = upv(i,j,lp)
951  grid2(i,j) = vpv(i,j,lp)
952  ENDDO
953  ENDDO
954  if(grib=='grib2') then
955  cfld=cfld+1
956  fld_info(cfld)%ifld=iavblfld(iget(336))
957  fld_info(cfld)%lvl=lvlsxml(lp,iget(336))
958 !$omp parallel do private(i,j,ii,jj)
959  do j=1,jend-jsta+1
960  jj = jsta+j-1
961  do i=1,iend-ista+1
962  ii = ista+i-1
963  datapd(i,j,cfld) = grid1(ii,jj)
964  enddo
965  enddo
966  cfld=cfld+1
967  fld_info(cfld)%ifld=iavblfld(iget(337))
968  fld_info(cfld)%lvl=lvlsxml(lp,iget(337))
969 !$omp parallel do private(i,j,ii,jj)
970  do j=1,jend-jsta+1
971  jj = jsta+j-1
972  do i=1,iend-ista+1
973  ii = ista+i-1
974  datapd(i,j,cfld) = grid2(ii,jj)
975  enddo
976  enddo
977  endif
978  ENDIF
979  ENDIF
980 
981 !*** T on constant PV
982 !
983 
984  IF(iget(338) > 0)THEN
985  IF(lvls(lp,iget(338)) > 0)THEN
986 ! GFS use lon avg as one scaler value for pole point
987  ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1:IM,JSTA:JEND) &
988  ! ,SPVAL,TPV(1:IM,JSTA:JEND,LP))
989  dum2d(ista:iend,jsta:jend)=tpv(ista:iend,jsta:jend,lp)
990  CALL exch(dum2d)
991  CALL fullpole(dum2d,pvpoles)
992  cosltemp=spval
993  IF(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
994  IF(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
995  pvtemp=spval
996  IF(jsta== 1) pvtemp(1:im, 1)=pvpoles(1:im,1)
997  IF(jend==jm) pvtemp(1:im,jm)=pvpoles(1:im,2)
998 
999  call poleavg(im,jm,jsta,jend,small,cosltemp(1:im,jsta:jend) &
1000  ,spval,pvtemp(1:im,jsta:jend))
1001 
1002  IF(jsta== 1) tpv(ista:iend, 1,lp)=pvtemp(ista:iend, 1)
1003  IF(jend==jm) tpv(ista:iend,jm,lp)=pvtemp(ista:iend,jm)
1004 
1005 !$omp parallel do private(i,j)
1006  DO j=jsta,jend
1007  DO i=ista,iend
1008  grid1(i,j) = tpv(i,j,lp)
1009  ENDDO
1010  ENDDO
1011  if(grib=='grib2') then
1012  cfld=cfld+1
1013  fld_info(cfld)%ifld=iavblfld(iget(338))
1014  fld_info(cfld)%lvl=lvlsxml(lp,iget(338))
1015 !$omp parallel do private(i,j,ii,jj)
1016  do j=1,jend-jsta+1
1017  jj = jsta+j-1
1018  do i=1,iend-ista+1
1019  ii = ista+i-1
1020  datapd(i,j,cfld) = grid1(ii,jj)
1021  enddo
1022  enddo
1023  endif
1024  ENDIF
1025  ENDIF
1026 !
1027 !*** Height on constant PV
1028 !
1029  IF(iget(339) > 0) THEN
1030  IF(lvls(lp,iget(339)) > 0)THEN
1031 ! GFS use lon avg as one scaler value for pole point
1032  ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1:IM,JSTA:JEND) &
1033  ! ,SPVAL,HPV(1:IM,JSTA:JEND,LP))
1034  dum2d(ista:iend,jsta:jend)=hpv(ista:iend,jsta:jend,lp)
1035  CALL exch(dum2d)
1036  CALL fullpole(dum2d,pvpoles)
1037  cosltemp=spval
1038  IF(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
1039  IF(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
1040  pvtemp=spval
1041  IF(jsta== 1) pvtemp(1:im, 1)=pvpoles(1:im,1)
1042  IF(jend==jm) pvtemp(1:im,jm)=pvpoles(1:im,2)
1043 
1044  call poleavg(im,jm,jsta,jend,small,cosltemp(1:im,jsta:jend) &
1045  ,spval,pvtemp(1:im,jsta:jend))
1046 
1047  IF(jsta== 1) hpv(ista:iend, 1,lp)=pvtemp(ista:iend, 1)
1048  IF(jend==jm) hpv(ista:iend,jm,lp)=pvtemp(ista:iend,jm)
1049 
1050 !$omp parallel do private(i,j)
1051  DO j=jsta,jend
1052  DO i=ista,iend
1053  grid1(i,j) = hpv(i,j,lp)
1054  ENDDO
1055  ENDDO
1056  if(grib=='grib2') then
1057  cfld=cfld+1
1058  fld_info(cfld)%ifld=iavblfld(iget(339))
1059  fld_info(cfld)%lvl=lvlsxml(lp,iget(339))
1060 !$omp parallel do private(i,j,ii,jj)
1061  do j=1,jend-jsta+1
1062  jj = jsta+j-1
1063  do i=1,iend-ista+1
1064  ii = ista+i-1
1065  datapd(i,j,cfld) = grid1(ii,jj)
1066  enddo
1067  enddo
1068  endif
1069  ENDIF
1070  ENDIF
1071 !
1072 !*** Pressure on constant PV
1073 !
1074  IF(iget(340) > 0) THEN
1075  IF(lvls(lp,iget(340)) > 0)THEN
1076 ! GFS use lon avg as one scaler value for pole point
1077  ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1:IM,JSTA:JEND) &
1078  ! ,SPVAL,PPV(1:IM,JSTA:JEND,LP))
1079  dum2d(ista:iend,jsta:jend)=ppv(ista:iend,jsta:jend,lp)
1080  CALL exch(dum2d)
1081  CALL fullpole(dum2d,pvpoles)
1082  cosltemp=spval
1083  IF(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
1084  IF(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
1085  pvtemp=spval
1086  IF(jsta== 1) pvtemp(1:im, 1)=pvpoles(1:im,1)
1087  IF(jend==jm) pvtemp(1:im,jm)=pvpoles(1:im,2)
1088 
1089  call poleavg(im,jm,jsta,jend,small,cosltemp(1:im,jsta:jend) &
1090  ,spval,pvtemp(1:im,jsta:jend))
1091 
1092  IF(jsta== 1) ppv(ista:iend, 1,lp)=pvtemp(ista:iend, 1)
1093  IF(jend==jm) ppv(ista:iend,jm,lp)=pvtemp(ista:iend,jm)
1094 
1095 !$omp parallel do private(i,j)
1096  DO j=jsta,jend
1097  DO i=ista,iend
1098  grid1(i,j) = ppv(i,j,lp)
1099  ENDDO
1100  ENDDO
1101  if(grib=='grib2') then
1102  cfld=cfld+1
1103  fld_info(cfld)%ifld=iavblfld(iget(340))
1104  fld_info(cfld)%lvl=lvlsxml(lp,iget(340))
1105 !$omp parallel do private(i,j,ii,jj)
1106  do j=1,jend-jsta+1
1107  jj = jsta+j-1
1108  do i=1,iend-ista+1
1109  ii = ista+i-1
1110  datapd(i,j,cfld) = grid1(ii,jj)
1111  enddo
1112  enddo
1113  endif
1114  ENDIF
1115  ENDIF
1116 !
1117 !*** Wind Shear on constant PV
1118 !
1119  IF(iget(341) > 0) THEN
1120  IF(lvls(lp,iget(341)) > 0)THEN
1121 ! GFS use lon avg as one scaler value for pole point
1122  ! call poleavg(IM,JM,JSTA,JEND,SMALL,COSL(1:IM,JSTA:JEND) &
1123  ! ,SPVAL,SPV(1:IM,JSTA:JEND,LP))
1124  dum2d(ista:iend,jsta:jend)=spv(ista:iend,jsta:jend,lp)
1125  CALL exch(dum2d)
1126  CALL fullpole(dum2d,pvpoles)
1127  cosltemp=spval
1128  IF(jsta== 1) cosltemp(1:im, 1)=coslpoles(1:im,1)
1129  IF(jend==jm) cosltemp(1:im,jm)=coslpoles(1:im,2)
1130  pvtemp=spval
1131  IF(jsta== 1) pvtemp(1:im, 1)=pvpoles(1:im,1)
1132  IF(jend==jm) pvtemp(1:im,jm)=pvpoles(1:im,2)
1133 
1134  call poleavg(im,jm,jsta,jend,small,cosltemp(1:im,jsta:jend) &
1135  ,spval,pvtemp(1:im,jsta:jend))
1136 
1137  IF(jsta== 1) spv(ista:iend, 1,lp)=pvtemp(ista:iend, 1)
1138  IF(jend==jm) spv(ista:iend,jm,lp)=pvtemp(ista:iend,jm)
1139 
1140 !$omp parallel do private(i,j)
1141  DO j=jsta,jend
1142  DO i=ista,iend
1143  grid1(i,j) = spv(i,j,lp)
1144  ENDDO
1145  ENDDO
1146  if(grib=='grib2') then
1147  cfld=cfld+1
1148  fld_info(cfld)%ifld=iavblfld(iget(341))
1149  fld_info(cfld)%lvl=lvlsxml(lp,iget(341))
1150 !$omp parallel do private(i,j,ii,jj)
1151  do j=1,jend-jsta+1
1152  jj = jsta+j-1
1153  do i=1,iend-ista+1
1154  ii = ista+i-1
1155  datapd(i,j,cfld) = grid1(ii,jj)
1156  enddo
1157  enddo
1158  endif
1159  ENDIF
1160  ENDIF
1161 
1162  END DO ! end loop for constant PV levels
1163 
1164  DEALLOCATE(dum1d1,dum1d2,dum1d3,dum1d4,dum1d5,dum1d6,dum1d7, &
1165  dum1d8,dum1d9,dum1d10,dum1d11,dum1d12,dum1d13, &
1166  dum1d14,wrk1, wrk2, wrk3, wrk4, cosl, dum2d)
1167 
1168  END IF ! end of selection for isentropic and constant PV fields
1169  if(me==0) write(0,*) 'MDL2THANDPV ends'
1170 !
1171 !
1172 ! END OF ROUTINE.
1173 !
1174  RETURN
1175  END
p2th
subroutine p2th(km, theta, u, v, h, t, pvu, sigma, rh, omga, kth, th, lth, uth, vth, hth, tth, zth, sigmath, rhth, oth)
p2th() interpolates to isentropic level.
Definition: GFSPOST.F:109
masks
Definition: MASKS_mod.f:1
physcons_post
Definition: physcons.f:1
params_mod
Definition: params.F:1
upp_math
dvdxdudy() computes dudy, dvdx, uwnd
Definition: UPP_MATH.f:17
p2pv
subroutine p2pv(km, pvu, h, t, p, u, v, kpv, pv, pvpt, pvpb, lpv, upv, vpv, hpv, tpv, ppv, spv)
p2pv() interpolates to potential vorticity level.
Definition: GFSPOST.F:175
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition: UPP_PHYSICS.f:345
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition: UPP_PHYSICS.f:27