UPP v11.0.0
Loading...
Searching...
No Matches
GFSPOST.F
Go to the documentation of this file.
1
40 subroutine pvetc(km,p,px,py,t,tx,ty,h,u,v,av,hm,s,bvf2,pvn,theta,sigma,pvu)
41
42 use physcons_post, only: con_cp, con_g, con_rd, con_rocp
43!
44 implicit none
45 integer,intent(in):: km
46 real,intent(in), dimension(km):: p,px,py,t,tx,ty,h,u,v,av
47 real,intent(out),dimension(km):: hm,s,bvf2,pvn,theta,sigma,pvu
48! real,parameter:: hhmin=500.,t0=2.e2,p0=1.e5
49 real,parameter:: hhmin=5.,t0=2.e2,p0=1.e5
50 integer k,kd,ku,k2(2)
51 real cprho,sx,sy,sz,uz,vz
52! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53 do k=1,km
54 hm(k) = con_cp*t(k) + con_g*h(k)
55 s(k) = con_cp*log(t(k)/t0) - con_rd*log(p(k)/p0)
56 enddo
57 do k=1,km
58 call rsearch1(km,h,2,(/h(k)-hhmin,h(k)+hhmin/),k2)
59! kd = max(k2(1),1)
60! ku = min(k2(2)+1,km)
61! kd = min(k2(1),km) ! Chuang: post counts from top down, redefine lower bound
62 kd = min(k2(1)+1,km) ! Chuang: post counts from top down,
63! ku = max(k2(2)-1,1)
64 ku = max(k2(2),1)
65 if(ku==1) kd=2 ! Chuang: make sure ku ne kd at model top
66 cprho = p(k)/(con_rocp*t(k))
67 sx = con_cp*tx(k) / t(k)-con_rd*px(k)/p(k)
68 sy = con_cp*ty(k) / t(k)-con_rd*py(k)/p(k)
69 sz = (s(ku)-s(kd)) / (h(ku)-h(kd))
70 uz = (u(ku)-u(kd)) / (h(ku)-h(kd))
71 vz = (v(ku)-v(kd)) / (h(ku)-h(kd))
72 bvf2(k) = con_g/con_cp*sz
73 pvn(k) = (av(k)*sz - vz*sx + uz*sy) / cprho
74 theta(k) = t0*exp(s(k)/con_cp)
75 sigma(k) = t(k)/con_g*bvf2(k)
76 pvu(k) = 1.e6*theta(k)*pvn(k)
77 enddo
78 end subroutine
79! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 subroutine p2th(km,theta,u,v,h,t,pvu,sigma,rh,omga,kth,th &
110 ,lth,uth,vth,hth,tth,zth,sigmath,rhth,oth)
111 implicit none
112 integer,intent(in):: km,kth
113 real,intent(in),dimension(km):: theta,u,v,h,t,pvu,sigma,rh,omga
114 real,intent(in):: th(kth)
115 logical*1,intent(out),dimension(kth):: lth
116 real,intent(out),dimension(kth):: uth,vth,hth,tth,zth &
117 ,sigmath,rhth,oth
118 real w
119 integer loc(kth),l
120 integer k
121! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122 call rsearch1(km,theta(1),kth,th(1),loc(1))
123 do k=1,kth
124 l = loc(k)
125 lth(k) = l > 0 .and.l < km
126 if(lth(k)) then
127 w = log(th(k)/theta(l)) / log(theta(l+1)/theta(l))
128 uth(k) = u(l) + w*(u(l+1)-u(l))
129 vth(k) = v(l) + w*(v(l+1)-v(l))
130 hth(k) = h(l) + w*(h(l+1)-h(l))
131 tth(k) = t(l) + w*(t(l+1)-t(l))
132 zth(k) = pvu(l) + w*(pvu(l+1)-pvu(l))
133 sigmath(k) = sigma(l) + w*(sigma(l+1)-sigma(l))
134 rhth(k) = rh(l) + w*(rh(l+1)-rh(l))
135! pth(k) = p(l) + w*(p(l+1)-p(l))
136 oth(k) = omga(l) + w*(omga(l+1)-omga(l))
137 endif
138 enddo
139 end subroutine
140! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175 subroutine p2pv(km,pvu,h,t,p,u,v,kpv,pv,pvpt,pvpb,&
176 lpv,upv,vpv,hpv,tpv,ppv,spv)
177 use physcons_post, only: con_rog
178 implicit none
179 integer,intent(in):: km,kpv
180 real,intent(in),dimension(km):: pvu,h,t,p,u,v
181 real,intent(in):: pv(kpv),pvpt(kpv),pvpb(kpv)
182 logical*1,intent(out),dimension(kpv):: lpv
183 real,intent(out),dimension(kpv):: upv,vpv,hpv,tpv,ppv,spv
184! real,parameter:: pd=2500.
185! Increase depth criteria for identifying PV layer from 25 to 50
186! to avoid finding shallow high level PV layer in high latitudes
187 real,parameter:: pd=5000.
188 real w,spdu,spdd
189 integer k,l1,l2,lu,ld,l
190! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
191 do k=1,kpv
192 call rsearch1(km,p,1,pvpb(k),l1)
193 call rsearch1(km,p,1,pvpt(k),l2)
194! l1=l1+1
195 l = 0
196 if(pv(k) >= 0.) then
197! do lu=l2-1,l1,-1
198! do lu=l2,l1-1 ! Chuang: post counts top down
199 do lu=l2+2,l1 ! Chuang: post counts top down
200! if(pv(k)<pvu(lu+1).and.pv(k)>=pvu(lu)) then
201 if(pv(k) >= pvu(lu+1).and.pv(k) < pvu(lu)) then
202 call rsearch1(km,p,1,p(lu)+pd,ld)
203! if(all(pv(k)>=pvu(ld:lu-1))) then
204 if(all(pv(k) >= pvu(lu+1:ld))) then
205 l = lu
206 exit
207 endif
208 endif
209 enddo
210 else
211! do lu=l2-1,l1,-1
212! do lu=l2,l1-1 ! Chuang: post counts top down
213 do lu=l2+2,l1 ! Chuang: post counts top down
214! if(pv(k)>pvu(lu+1).and.pv(k)<=pvu(lu)) then
215 if(pv(k) <= pvu(lu+1).and.pv(k) > pvu(lu)) then
216 call rsearch1(km,p,1,p(lu)+pd,ld)
217! if(all(pv(k)<=pvu(ld:lu-1))) then
218 if(all(pv(k) <= pvu(lu+1:ld))) then
219 l = lu
220 exit
221 endif
222 endif
223 enddo
224 endif
225 lpv(k) = l > 0
226 if(lpv(k)) then
227 w = (pv(k)-pvu(l))/(pvu(l+1)-pvu(l))
228 upv(k) = u(l) + w*(u(l+1)-u(l))
229 vpv(k) = v(l) + w*(v(l+1)-v(l))
230 hpv(k) = h(l) + w*(h(l+1)-h(l))
231 tpv(k) = t(l) + w*(t(l+1)-t(l))
232 ppv(k) = p(l)*exp((h(l)-hpv(k))*(1-0.5*(tpv(k)/t(l)-1))/(con_rog*t(l)))
233
234 spdu = sqrt(u(l+1)*u(l+1) + v(l+1)*v(l+1))
235 spdd = sqrt(u(l)*u(l) + v(l)*v(l))
236 spv(k) = (spdu-spdd) / (h(l+1)-h(l))
237 endif
238 enddo
239 end subroutine
240!-------------------------------------------------------------------------------
241
242!-------------------------------------------------------------------------------
269 subroutine rsearch1(km1,z1,km2,z2,l2)
270 implicit none
271 integer,intent(in):: km1,km2
272 real,intent(in):: z1(km1),z2(km2)
273 integer,intent(out):: l2(km2)
274 integer k1,k2
275! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276! find the surrounding input interval for each output point.
277 if(z1(1) <= z1(km1)) then
278! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
279! input coordinate is monotonically ascending.
280 do k2=1,km2
281 if (z1(1) >= z2(k2)) then
282 l2(k2) = 1
283 else
284 l2(k2)=km1
285 do k1=1,km1-1
286 if(z1(k1) <= z2(k2) .and. z1(k1+1) > z2(k2)) then
287 l2(k2) = k1
288 exit
289 endif
290 enddo
291 endif
292 enddo
293 else
294! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
295! input coordinate is monotonically descending.
296 do k2=1,km2
297 if (z1(1) <= z2(k2)) then
298 l2(k2) = 1
299 else
300 l2(k2)=km1
301 do k1=km1,2,-1
302 if(z2(k2) >= z1(k1) .and. z2(k2) < z1(k1-1)) then
303 l2(k2) = k1-1
304 exit
305 endif
306 enddo
307 endif
308 enddo
309 endif
310! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
311end subroutine
312! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
344 subroutine tpause(km,p,u,v,t,h,ptp,utp,vtp,ttp,htp,shrtp)
345 use physcons_post, only: con_rog
346 implicit none
347 integer,intent(in):: km
348 real,intent(in),dimension(km):: p,u,v,t,h
349 real,intent(out):: ptp,utp,vtp,ttp,htp,shrtp
350 real,parameter:: ptplim(2)=(/500.e+2,50.e+2/),gamtp=2.e-3,hd=2.e+3
351 real gamu,gamd,td,gami,wtp,spdu,spdd
352 integer klim(2),k,kd,ktp
353! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
354! find tropopause level
355 call rsearch1(km-2,p(2),2,ptplim,klim)
356 klim(1)=klim(1)+1
357 klim(2)=klim(2)+2
358 ! klim(1) > klim(2) or loops does not run ; klim(2) has a
359 ! minimum value of 3 to insure k-2 != 0 in called subprogram
360 gamd=1.e+9
361 ktp=klim(2)
362 wtp=0
363! do k=klim(1),klim(2)
364 do k=klim(1),klim(2),-1
365! gamu=(t(k-1)-t(k+1))/(h(k+1)-h(k-1))
366 gamu=(t(k+1)-t(k-1))/(h(k-1)-h(k+1))
367 if(gamu<=gamtp) then
368! call rsearch1(km-k-1,h(k+1),1,h(k)+hd,kd)
369 call rsearch1(k-2,h(2),1,h(k)+hd,kd)
370! td=t(k+kd)+(h(k)+hd-h(k+kd))/(h(k+kd+1)-h(k+kd))*(t(k+kd+1)-t(k+kd))
371 td=t(kd+2)+(h(k)+hd-h(2+kd))/(h(kd+1)-h(2+kd))*(t(kd+1)-t(2+kd))
372 gami=(t(k)-td)/hd
373 if(gami<=gamtp) then
374 ktp=k
375 wtp=(gamtp-gamu)/(max(gamd,gamtp+0.1e-3)-gamu)
376 exit
377 endif
378 endif
379 gamd=gamu
380 enddo
381! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
382! compute tropopause level fields
383 utp=u(ktp)-wtp*(u(ktp)-u(ktp-1))
384 vtp=v(ktp)-wtp*(v(ktp)-v(ktp-1))
385 ttp=t(ktp)-wtp*(t(ktp)-t(ktp-1))
386 htp=h(ktp)-wtp*(h(ktp)-h(ktp-1))
387 ptp=p(ktp)*exp((h(ktp)-htp)*(1-0.5*(ttp/t(ktp)-1))/(con_rog*t(ktp)))
388 spdu=sqrt(u(ktp)**2+v(ktp)**2)
389 spdd=sqrt(u(ktp-1)**2+v(ktp-1)**2)
390 shrtp=(spdu-spdd)/(h(ktp)-h(ktp-1))
391
392 utp=u(ktp)-wtp*(u(ktp)-u(ktp+1))
393 vtp=v(ktp)-wtp*(v(ktp)-v(ktp+1))
394 ttp=t(ktp)-wtp*(t(ktp)-t(ktp+1))
395 htp=h(ktp)-wtp*(h(ktp)-h(ktp+1))
396 ptp=p(ktp)*exp((h(ktp)-htp)*(1-0.5*(ttp/t(ktp)-1))/(con_rog*t(ktp)))
397 spdu=sqrt(u(ktp)**2+v(ktp)**2)
398 spdd=sqrt(u(ktp+1)**2+v(ktp+1)**2)
399 shrtp=(spdu-spdd)/(h(ktp)-h(ktp+1))
400 end subroutine
401! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
431 subroutine mxwind(km,p,u,v,t,h,pmw,umw,vmw,tmw,hmw)
432 use physcons_post, only: con_rog
433 implicit none
434 integer,intent(in):: km
435 real,intent(in),dimension(km):: p,u,v,t,h
436 real,intent(out):: pmw,umw,vmw,tmw,hmw
437 real,parameter:: pmwlim(2)=(/500.e+2,100.e+2/)
438 integer klim(2),k,kmw
439 real spd(km),spdmw,wmw,dhd,dhu,shrd,shru,dhmw,ub,vb,spdb
440! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
441! find maximum wind level
442 call rsearch1(km,p(1),2,pmwlim,klim)
443! klim(1)=klim(1)+1
444 klim(2)=klim(2)+1
445! spd(klim(1):klim(2))=sqrt(u(klim(1):klim(2))**2+v(klim(1):klim(2))**2)
446 spd(klim(2):klim(1))=sqrt(u(klim(2):klim(1))**2+v(klim(2):klim(1))**2)
447 spdmw=spd(klim(1))
448 kmw=klim(1)
449! do k=klim(1)+1,klim(2)
450 do k=klim(1)-1,klim(2),-1
451 if(spd(k)>spdmw) then
452 spdmw=spd(k)
453 kmw=k
454 endif
455 enddo
456! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
457! find speed and height at the maximum wind level
458 if(kmw==klim(1).or.kmw==klim(2)) then
459 hmw=h(kmw)
460 spdmw=spd(kmw)
461 wmw=0.
462 else
463! dhd=h(kmw)-h(kmw-1)
464 dhd=h(kmw)-h(kmw+1) !post counts top down
465! dhu=h(kmw+1)-h(kmw)
466 dhu=h(kmw-1)-h(kmw)
467! shrd=(spd(kmw)-spd(kmw-1))/(h(kmw)-h(kmw-1))
468 shrd=(spd(kmw)-spd(kmw+1))/(h(kmw)-h(kmw+1))
469! shru=(spd(kmw)-spd(kmw+1))/(h(kmw+1)-h(kmw))
470 shru=(spd(kmw)-spd(kmw-1))/(h(kmw-1)-h(kmw))
471 dhmw=(shrd*dhu-shru*dhd)/(2*(shrd+shru))
472 hmw=h(kmw)+dhmw
473 spdmw=spd(kmw)+dhmw**2*(shrd+shru)/(dhd+dhu)
474! if(dhmw>0) kmw=kmw+1
475 if(dhmw>0) kmw=kmw-1
476! wmw=(h(kmw)-hmw)/(h(kmw)-h(kmw-1))
477 wmw=(h(kmw)-hmw)/(h(kmw)-h(kmw+1))
478 endif
479! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
480! compute maximum wind level fields
481! ub=u(kmw)-wmw*(u(kmw)-u(kmw-1))
482 ub=u(kmw)-wmw*(u(kmw)-u(kmw+1))
483! vb=v(kmw)-wmw*(v(kmw)-v(kmw-1))
484 vb=v(kmw)-wmw*(v(kmw)-v(kmw+1))
485 spdb=max(sqrt(ub**2+vb**2),1.e-6)
486 umw=ub*spdmw/spdb
487 vmw=vb*spdmw/spdb
488! tmw=t(kmw)-wmw*(t(kmw)-t(kmw-1))
489 tmw=t(kmw)-wmw*(t(kmw)-t(kmw+1))
490 pmw=p(kmw)*exp((h(kmw)-hmw)*(1-0.5*(tmw/t(kmw)-1))/(con_rog*t(kmw)))
491 end subroutine
492
493! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
519 subroutine mptgen(mpirank,mpisize,nd,jt1,jt2,j1,j2,jx,jm,jn)
520 use machine_post,only:kint_mpi
521 implicit none
522 integer(kint_mpi),intent(in):: mpirank,mpisize,nd,jt1(nd),jt2(nd)
523 integer(kint_mpi),intent(out):: j1(nd),j2(nd),jx(nd),jm(nd),jn(nd)
524 integer msize,mrank,msn,mrn,n
525! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
526 msize=mpisize
527 mrank=mpirank
528 do n=nd,1,-1
529 if(jt2(n)>=jt1(n)) then
530 jm(n)=(jt2(n)-jt1(n))/msize+1
531 msn=max(msize/(jt2(n)-jt1(n)+1),1)
532 if(n==1) msn=1
533 jn(n)=msize/msn
534 mrn=mrank/msn
535 j1(n)=min(jt1(n)+jm(n)*mrn,jt2(n)+1)
536 j2(n)=min(jt1(n)+jm(n)*mrn+jm(n)-1,jt2(n))
537 jx(n)=j2(n)-j1(n)+1
538 msize=msn
539 mrank=mod(mrank,msn)
540 write(0,*)' mrank=',mrank,' j1=',j1(n),' j2=',j2(n),' jx=',jx(n),' jm=',jm(n)
541 else
542 jm(n)=0
543 jn(n)=1
544 j1(n)=jt1(n)
545 j2(n)=jt2(n)
546 jx(n)=0
547 endif
548 enddo
549! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
550end subroutine
551!-------------------------------------------------------------------------------
621 subroutine mptranr4(mpicomm,mpisize,im,ida,idb,&
622 jm,jma,jmb,jda,km,kma,kmb,kdb,a,b,ta,tb)
623 use machine_post,only:kint_mpi
624 implicit none
625 include 'mpif.h'
626 integer(kint_mpi),intent(in):: mpicomm,mpisize
627 integer(kint_mpi),intent(in):: im,ida,idb
628 integer(kint_mpi),intent(in):: jm,jma,jmb,jda
629 integer(kint_mpi),intent(in):: km,kma,kmb,kdb
630 real(4),dimension(ida,jda,kma),intent(in):: a
631 real(4),dimension(idb,kdb,jmb),intent(out):: b
632 real(4),dimension(im,jm,km,mpisize),intent(inout):: ta,tb
633 integer(4) jmb1(1),jmbf(mpisize),kma1(1),kmaf(mpisize)
634 integer(kint_mpi)::i,j,k,l,ierr,ja,kb
635! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
636 jmb1(1) = jmb
637 call mpi_allgather(jmb1,1,mpi_integer4,jmbf,1,mpi_integer4,mpicomm,ierr)
638 kma1(1) = kma
639 call mpi_allgather(kma1,1,mpi_integer4,kmaf,1,mpi_integer4,mpicomm,ierr)
640! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
641! internally transpose input array
642!$omp parallel do private(i,j,k,l,ja)
643 do l=1,mpisize
644 do k=1,kma
645 do j=1,jm
646 ja = j + sum(jmbf(1:l-1))
647 if(ja <= jma) then
648 do i=1,im
649 ta(i,j,k,l) = a(i,ja,k)
650 enddo
651 endif
652 enddo
653 enddo
654 enddo
655! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
656! externally transpose data
657 call mpi_alltoall(ta,im*jm*km,mpi_real4,tb,im*jm*km,mpi_real4,mpicomm,ierr)
658! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
659! internally transpose output array
660!$omp parallel do private(i,j,k,l,kb)
661 do l=1,mpisize
662 do k=1,km
663 kb = k + sum(kmaf(1:l-1))
664 if(kb <= kmb) then
665 do j=1,jmb
666 do i=1,im
667 b(i,kb,j) = tb(i,j,k,l)
668 enddo
669 enddo
670 endif
671 enddo
672 enddo
673! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
674end subroutine
675!-----------------------------------------------------------------------
mptranr4
subroutine mptranr4(mpicomm, mpisize, im, ida, idb, jm, jma, jmb, jda, km, kma, kmb, kdb, a, b, ta, tb)
mptranr4() transposes grid decompositions.
Definition GFSPOST.F:623
mxwind
subroutine mxwind(km, p, u, v, t, h, pmw, umw, vmw, tmw, hmw)
mxwind() computes maximum wind level fields.
Definition GFSPOST.F:432
rsearch1
subroutine rsearch1(km1, z1, km2, z2, l2)
rsearch1() searches for a surrounding real interval.
Definition GFSPOST.F:270
tpause
subroutine tpause(km, p, u, v, t, h, ptp, utp, vtp, ttp, htp, shrtp)
tpause() computes tropopause level fields.
Definition GFSPOST.F:345
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
mptgen
subroutine mptgen(mpirank, mpisize, nd, jt1, jt2, j1, j2, jx, jm, jn)
mptgen() generates grid decomposition dimensions.
Definition GFSPOST.F:520