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