UPP (develop)
Loading...
Searching...
No Matches
GFSPOST.F
Go to the documentation of this file.
1
40!-------------------------------------------------------------------------------
61!
62 subroutine pvetc(km,p,px,py,t,tx,ty,h,u,v,av,hm,s,bvf2,pvn,theta,sigma,pvu)
63
64 use physcons_post, only: con_cp, con_g, con_rd, con_rocp
65!
66 implicit none
67 integer,intent(in):: km
68 real,intent(in), dimension(km):: p,px,py,t,tx,ty,h,u,v,av
69 real,intent(out),dimension(km):: hm,s,bvf2,pvn,theta,sigma,pvu
70! real,parameter:: hhmin=500.,t0=2.e2,p0=1.e5
71 real,parameter:: hhmin=5.,t0=2.e2,p0=1.e5
72 integer k,kd,ku,k2(2)
73 real cprho,sx,sy,sz,uz,vz
74! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75 do k=1,km
76 hm(k) = con_cp*t(k) + con_g*h(k)
77 s(k) = con_cp*log(t(k)/t0) - con_rd*log(p(k)/p0)
78 enddo
79 do k=1,km
80 call rsearch1(km,h,2,(/h(k)-hhmin,h(k)+hhmin/),k2)
81! kd = max(k2(1),1)
82! ku = min(k2(2)+1,km)
83! kd = min(k2(1),km) ! Chuang: post counts from top down, redefine lower bound
84 kd = min(k2(1)+1,km) ! Chuang: post counts from top down,
85! ku = max(k2(2)-1,1)
86 ku = max(k2(2),1)
87 if(ku==1) kd=2 ! Chuang: make sure ku ne kd at model top
88 cprho = p(k)/(con_rocp*t(k))
89 sx = con_cp*tx(k) / t(k)-con_rd*px(k)/p(k)
90 sy = con_cp*ty(k) / t(k)-con_rd*py(k)/p(k)
91 sz = (s(ku)-s(kd)) / (h(ku)-h(kd))
92 uz = (u(ku)-u(kd)) / (h(ku)-h(kd))
93 vz = (v(ku)-v(kd)) / (h(ku)-h(kd))
94 bvf2(k) = con_g/con_cp*sz
95 pvn(k) = (av(k)*sz - vz*sx + uz*sy) / cprho
96 theta(k) = t0*exp(s(k)/con_cp)
97 sigma(k) = t(k)/con_g*bvf2(k)
98 pvu(k) = 1.e6*theta(k)*pvn(k)
99 enddo
100 end subroutine
101! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137 subroutine p2th(km,theta,u,v,h,t,pvu,sigma,rh,omga,kth,th &
138 ,lth,uth,vth,hth,tth,zth,sigmath,rhth,oth)
139 implicit none
140 integer,intent(in):: km,kth
141 real,intent(in),dimension(km):: theta,u,v,h,t,pvu,sigma,rh,omga
142 real,intent(in):: th(kth)
143 logical*1,intent(out),dimension(kth):: lth
144 real,intent(out),dimension(kth):: uth,vth,hth,tth,zth &
145 ,sigmath,rhth,oth
146 real w
147 integer loc(kth),l
148 integer k
149! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150 call rsearch1(km,theta(1),kth,th(1),loc(1))
151 do k=1,kth
152 l = loc(k)
153 lth(k) = l > 0 .and.l < km
154 if(lth(k)) then
155 w = log(th(k)/theta(l)) / log(theta(l+1)/theta(l))
156 uth(k) = u(l) + w*(u(l+1)-u(l))
157 vth(k) = v(l) + w*(v(l+1)-v(l))
158 hth(k) = h(l) + w*(h(l+1)-h(l))
159 tth(k) = t(l) + w*(t(l+1)-t(l))
160 zth(k) = pvu(l) + w*(pvu(l+1)-pvu(l))
161 sigmath(k) = sigma(l) + w*(sigma(l+1)-sigma(l))
162 rhth(k) = rh(l) + w*(rh(l+1)-rh(l))
163! pth(k) = p(l) + w*(p(l+1)-p(l))
164 oth(k) = omga(l) + w*(omga(l+1)-omga(l))
165 endif
166 enddo
167 end subroutine
168! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203 subroutine p2pv(km,pvu,h,t,p,u,v,kpv,pv,pvpt,pvpb,&
204 lpv,upv,vpv,hpv,tpv,ppv,spv)
205 use physcons_post, only: con_rog
206 implicit none
207 integer,intent(in):: km,kpv
208 real,intent(in),dimension(km):: pvu,h,t,p,u,v
209 real,intent(in):: pv(kpv),pvpt(kpv),pvpb(kpv)
210 logical*1,intent(out),dimension(kpv):: lpv
211 real,intent(out),dimension(kpv):: upv,vpv,hpv,tpv,ppv,spv
212! real,parameter:: pd=2500.
213! Increase depth criteria for identifying PV layer from 25 to 50
214! to avoid finding shallow high level PV layer in high latitudes
215 real,parameter:: pd=5000.
216 real w,spdu,spdd
217 integer k,l1,l2,lu,ld,l
218! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
219 do k=1,kpv
220 call rsearch1(km,p,1,pvpb(k),l1)
221 call rsearch1(km,p,1,pvpt(k),l2)
222! l1=l1+1
223 l = 0
224 if(pv(k) >= 0.) then
225! do lu=l2-1,l1,-1
226! do lu=l2,l1-1 ! Chuang: post counts top down
227 do lu=l2+2,l1 ! Chuang: post counts top down
228! if(pv(k)<pvu(lu+1).and.pv(k)>=pvu(lu)) then
229 if(pv(k) >= pvu(lu+1).and.pv(k) < pvu(lu)) then
230 call rsearch1(km,p,1,p(lu)+pd,ld)
231! if(all(pv(k)>=pvu(ld:lu-1))) then
232 if(all(pv(k) >= pvu(lu+1:ld))) then
233 l = lu
234 exit
235 endif
236 endif
237 enddo
238 else
239! do lu=l2-1,l1,-1
240! do lu=l2,l1-1 ! Chuang: post counts top down
241 do lu=l2+2,l1 ! Chuang: post counts top down
242! if(pv(k)>pvu(lu+1).and.pv(k)<=pvu(lu)) then
243 if(pv(k) <= pvu(lu+1).and.pv(k) > pvu(lu)) then
244 call rsearch1(km,p,1,p(lu)+pd,ld)
245! if(all(pv(k)<=pvu(ld:lu-1))) then
246 if(all(pv(k) <= pvu(lu+1:ld))) then
247 l = lu
248 exit
249 endif
250 endif
251 enddo
252 endif
253 lpv(k) = l > 0
254 if(lpv(k)) then
255 w = (pv(k)-pvu(l))/(pvu(l+1)-pvu(l))
256 upv(k) = u(l) + w*(u(l+1)-u(l))
257 vpv(k) = v(l) + w*(v(l+1)-v(l))
258 hpv(k) = h(l) + w*(h(l+1)-h(l))
259 tpv(k) = t(l) + w*(t(l+1)-t(l))
260 ppv(k) = p(l)*exp((h(l)-hpv(k))*(1-0.5*(tpv(k)/t(l)-1))/(con_rog*t(l)))
261
262 spdu = sqrt(u(l+1)*u(l+1) + v(l+1)*v(l+1))
263 spdd = sqrt(u(l)*u(l) + v(l)*v(l))
264 spv(k) = (spdu-spdd) / (h(l+1)-h(l))
265 endif
266 enddo
267 end subroutine
268!-------------------------------------------------------------------------------
269
270!-------------------------------------------------------------------------------
297 subroutine rsearch1(km1,z1,km2,z2,l2)
298 implicit none
299 integer,intent(in):: km1,km2
300 real,intent(in):: z1(km1),z2(km2)
301 integer,intent(out):: l2(km2)
302 integer k1,k2
303! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
304! find the surrounding input interval for each output point.
305 if(z1(1) <= z1(km1)) then
306! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
307! input coordinate is monotonically ascending.
308 do k2=1,km2
309 if (z1(1) >= z2(k2)) then
310 l2(k2) = 1
311 else
312 l2(k2)=km1
313 do k1=1,km1-1
314 if(z1(k1) <= z2(k2) .and. z1(k1+1) > z2(k2)) then
315 l2(k2) = k1
316 exit
317 endif
318 enddo
319 endif
320 enddo
321 else
322! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
323! input coordinate is monotonically descending.
324 do k2=1,km2
325 if (z1(1) <= z2(k2)) then
326 l2(k2) = 1
327 else
328 l2(k2)=km1
329 do k1=km1,2,-1
330 if(z2(k2) >= z1(k1) .and. z2(k2) < z1(k1-1)) then
331 l2(k2) = k1-1
332 exit
333 endif
334 enddo
335 endif
336 enddo
337 endif
338! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
339end subroutine
340! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
372 subroutine tpause(km,p,u,v,t,h,ptp,utp,vtp,ttp,htp,shrtp)
373 use physcons_post, only: con_rog
374 implicit none
375 integer,intent(in):: km
376 real,intent(in),dimension(km):: p,u,v,t,h
377 real,intent(out):: ptp,utp,vtp,ttp,htp,shrtp
378 real,parameter:: ptplim(2)=(/500.e+2,50.e+2/),gamtp=2.e-3,hd=2.e+3
379 real gamu,gamd,td,gami,wtp,spdu,spdd
380 integer klim(2),k,kd,ktp
381! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
382! find tropopause level
383 call rsearch1(km-2,p(2),2,ptplim,klim)
384 klim(1)=klim(1)+1
385 klim(2)=klim(2)+2
386 ! klim(1) > klim(2) or loops does not run ; klim(2) has a
387 ! minimum value of 3 to insure k-2 != 0 in called subprogram
388 gamd=1.e+9
389 ktp=klim(2)
390 wtp=0
391! do k=klim(1),klim(2)
392 do k=klim(1),klim(2),-1
393! gamu=(t(k-1)-t(k+1))/(h(k+1)-h(k-1))
394 gamu=(t(k+1)-t(k-1))/(h(k-1)-h(k+1))
395 if(gamu<=gamtp) then
396! call rsearch1(km-k-1,h(k+1),1,h(k)+hd,kd)
397 call rsearch1(k-2,h(2),1,h(k)+hd,kd)
398! td=t(k+kd)+(h(k)+hd-h(k+kd))/(h(k+kd+1)-h(k+kd))*(t(k+kd+1)-t(k+kd))
399 td=t(kd+2)+(h(k)+hd-h(2+kd))/(h(kd+1)-h(2+kd))*(t(kd+1)-t(2+kd))
400 gami=(t(k)-td)/hd
401 if(gami<=gamtp) then
402 ktp=k
403 wtp=(gamtp-gamu)/(max(gamd,gamtp+0.1e-3)-gamu)
404 exit
405 endif
406 endif
407 gamd=gamu
408 enddo
409! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
410! compute tropopause level fields
411 utp=u(ktp)-wtp*(u(ktp)-u(ktp-1))
412 vtp=v(ktp)-wtp*(v(ktp)-v(ktp-1))
413 ttp=t(ktp)-wtp*(t(ktp)-t(ktp-1))
414 htp=h(ktp)-wtp*(h(ktp)-h(ktp-1))
415 ptp=p(ktp)*exp((h(ktp)-htp)*(1-0.5*(ttp/t(ktp)-1))/(con_rog*t(ktp)))
416 spdu=sqrt(u(ktp)**2+v(ktp)**2)
417 spdd=sqrt(u(ktp-1)**2+v(ktp-1)**2)
418 shrtp=(spdu-spdd)/(h(ktp)-h(ktp-1))
419
420 utp=u(ktp)-wtp*(u(ktp)-u(ktp+1))
421 vtp=v(ktp)-wtp*(v(ktp)-v(ktp+1))
422 ttp=t(ktp)-wtp*(t(ktp)-t(ktp+1))
423 htp=h(ktp)-wtp*(h(ktp)-h(ktp+1))
424 ptp=p(ktp)*exp((h(ktp)-htp)*(1-0.5*(ttp/t(ktp)-1))/(con_rog*t(ktp)))
425 spdu=sqrt(u(ktp)**2+v(ktp)**2)
426 spdd=sqrt(u(ktp+1)**2+v(ktp+1)**2)
427 shrtp=(spdu-spdd)/(h(ktp)-h(ktp+1))
428 end subroutine
429! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
459 subroutine mxwind(km,p,u,v,t,h,pmw,umw,vmw,tmw,hmw)
460 use physcons_post, only: con_rog
461 implicit none
462 integer,intent(in):: km
463 real,intent(in),dimension(km):: p,u,v,t,h
464 real,intent(out):: pmw,umw,vmw,tmw,hmw
465 real,parameter:: pmwlim(2)=(/500.e+2,100.e+2/)
466 integer klim(2),k,kmw
467 real spd(km),spdmw,wmw,dhd,dhu,shrd,shru,dhmw,ub,vb,spdb
468! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
469! find maximum wind level
470 call rsearch1(km,p(1),2,pmwlim,klim)
471! klim(1)=klim(1)+1
472 klim(2)=klim(2)+1
473! spd(klim(1):klim(2))=sqrt(u(klim(1):klim(2))**2+v(klim(1):klim(2))**2)
474 spd(klim(2):klim(1))=sqrt(u(klim(2):klim(1))**2+v(klim(2):klim(1))**2)
475 spdmw=spd(klim(1))
476 kmw=klim(1)
477! do k=klim(1)+1,klim(2)
478 do k=klim(1)-1,klim(2),-1
479 if(spd(k)>spdmw) then
480 spdmw=spd(k)
481 kmw=k
482 endif
483 enddo
484! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
485! find speed and height at the maximum wind level
486 if(kmw==klim(1).or.kmw==klim(2)) then
487 hmw=h(kmw)
488 spdmw=spd(kmw)
489 wmw=0.
490 else
491! dhd=h(kmw)-h(kmw-1)
492 dhd=h(kmw)-h(kmw+1) !post counts top down
493! dhu=h(kmw+1)-h(kmw)
494 dhu=h(kmw-1)-h(kmw)
495! shrd=(spd(kmw)-spd(kmw-1))/(h(kmw)-h(kmw-1))
496 shrd=(spd(kmw)-spd(kmw+1))/(h(kmw)-h(kmw+1))
497! shru=(spd(kmw)-spd(kmw+1))/(h(kmw+1)-h(kmw))
498 shru=(spd(kmw)-spd(kmw-1))/(h(kmw-1)-h(kmw))
499 dhmw=(shrd*dhu-shru*dhd)/(2*(shrd+shru))
500 hmw=h(kmw)+dhmw
501 spdmw=spd(kmw)+dhmw**2*(shrd+shru)/(dhd+dhu)
502! if(dhmw>0) kmw=kmw+1
503 if(dhmw>0) kmw=kmw-1
504! wmw=(h(kmw)-hmw)/(h(kmw)-h(kmw-1))
505 wmw=(h(kmw)-hmw)/(h(kmw)-h(kmw+1))
506 endif
507! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
508! compute maximum wind level fields
509! ub=u(kmw)-wmw*(u(kmw)-u(kmw-1))
510 ub=u(kmw)-wmw*(u(kmw)-u(kmw+1))
511! vb=v(kmw)-wmw*(v(kmw)-v(kmw-1))
512 vb=v(kmw)-wmw*(v(kmw)-v(kmw+1))
513 spdb=max(sqrt(ub**2+vb**2),1.e-6)
514 umw=ub*spdmw/spdb
515 vmw=vb*spdmw/spdb
516! tmw=t(kmw)-wmw*(t(kmw)-t(kmw-1))
517 tmw=t(kmw)-wmw*(t(kmw)-t(kmw+1))
518 pmw=p(kmw)*exp((h(kmw)-hmw)*(1-0.5*(tmw/t(kmw)-1))/(con_rog*t(kmw)))
519 end subroutine
520
521! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
547 subroutine mptgen(mpirank,mpisize,nd,jt1,jt2,j1,j2,jx,jm,jn)
548 use machine_post,only:kint_mpi
549 implicit none
550 integer(kint_mpi),intent(in):: mpirank,mpisize,nd,jt1(nd),jt2(nd)
551 integer(kint_mpi),intent(out):: j1(nd),j2(nd),jx(nd),jm(nd),jn(nd)
552 integer msize,mrank,msn,mrn,n
553! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
554 msize=mpisize
555 mrank=mpirank
556 do n=nd,1,-1
557 if(jt2(n)>=jt1(n)) then
558 jm(n)=(jt2(n)-jt1(n))/msize+1
559 msn=max(msize/(jt2(n)-jt1(n)+1),1)
560 if(n==1) msn=1
561 jn(n)=msize/msn
562 mrn=mrank/msn
563 j1(n)=min(jt1(n)+jm(n)*mrn,jt2(n)+1)
564 j2(n)=min(jt1(n)+jm(n)*mrn+jm(n)-1,jt2(n))
565 jx(n)=j2(n)-j1(n)+1
566 msize=msn
567 mrank=mod(mrank,msn)
568 write(*,*)' mrank=',mrank,' j1=',j1(n),' j2=',j2(n),' jx=',jx(n),' jm=',jm(n)
569 else
570 jm(n)=0
571 jn(n)=1
572 j1(n)=jt1(n)
573 j2(n)=jt2(n)
574 jx(n)=0
575 endif
576 enddo
577! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
578end subroutine
579!-------------------------------------------------------------------------------
649 subroutine mptranr4(mpicomm,mpisize,im,ida,idb,&
650 jm,jma,jmb,jda,km,kma,kmb,kdb,a,b,ta,tb)
651 use machine_post,only:kint_mpi
652 implicit none
653 include 'mpif.h'
654 integer(kint_mpi),intent(in):: mpicomm,mpisize
655 integer(kint_mpi),intent(in):: im,ida,idb
656 integer(kint_mpi),intent(in):: jm,jma,jmb,jda
657 integer(kint_mpi),intent(in):: km,kma,kmb,kdb
658 real(4),dimension(ida,jda,kma),intent(in):: a
659 real(4),dimension(idb,kdb,jmb),intent(out):: b
660 real(4),dimension(im,jm,km,mpisize),intent(inout):: ta,tb
661 integer(4) jmb1(1),jmbf(mpisize),kma1(1),kmaf(mpisize)
662 integer(kint_mpi)::i,j,k,l,ierr,ja,kb
663! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
664 jmb1(1) = jmb
665 call mpi_allgather(jmb1,1,mpi_integer4,jmbf,1,mpi_integer4,mpicomm,ierr)
666 kma1(1) = kma
667 call mpi_allgather(kma1,1,mpi_integer4,kmaf,1,mpi_integer4,mpicomm,ierr)
668! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
669! internally transpose input array
670!$omp parallel do private(i,j,k,l,ja)
671 do l=1,mpisize
672 do k=1,kma
673 do j=1,jm
674 ja = j + sum(jmbf(1:l-1))
675 if(ja <= jma) then
676 do i=1,im
677 ta(i,j,k,l) = a(i,ja,k)
678 enddo
679 endif
680 enddo
681 enddo
682 enddo
683! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
684! externally transpose data
685 call mpi_alltoall(ta,im*jm*km,mpi_real4,tb,im*jm*km,mpi_real4,mpicomm,ierr)
686! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
687! internally transpose output array
688!$omp parallel do private(i,j,k,l,kb)
689 do l=1,mpisize
690 do k=1,km
691 kb = k + sum(kmaf(1:l-1))
692 if(kb <= kmb) then
693 do j=1,jmb
694 do i=1,im
695 b(i,kb,j) = tb(i,j,k,l)
696 enddo
697 enddo
698 endif
699 enddo
700 enddo
701! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
702end subroutine
703!-----------------------------------------------------------------------
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:651
mxwind
subroutine mxwind(km, p, u, v, t, h, pmw, umw, vmw, tmw, hmw)
mxwind() computes maximum wind level fields.
Definition GFSPOST.F:460
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
rsearch1
subroutine rsearch1(km1, z1, km2, z2, l2)
rsearch1() searches for a surrounding real interval.
Definition GFSPOST.F:298
tpause
subroutine tpause(km, p, u, v, t, h, ptp, utp, vtp, ttp, htp, shrtp)
tpause() computes tropopause level fields.
Definition GFSPOST.F:373
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
mptgen
subroutine mptgen(mpirank, mpisize, nd, jt1, jt2, j1, j2, jx, jm, jn)
mptgen() generates grid decomposition dimensions.
Definition GFSPOST.F:548