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