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