UPP v11.0.0
Loading...
Searching...
No Matches
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:111
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:177
upp_math
dvdxdudy() computes dudy, dvdx, uwnd
Definition UPP_MATH.f:17
upp_physics
calcape() computes CAPE/CINS and other storm related variables.
Definition UPP_PHYSICS.f:27
upp_physics::fpvsnew
elemental real function, public fpvsnew(t)
Definition UPP_PHYSICS.f:346