UPP  11.0.0
 All Data Structures Files Functions Variables Pages
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 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
339 end 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 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
578 end 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 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
702 end subroutine
703 !-----------------------------------------------------------------------
mxwind
subroutine mxwind(km, p, u, v, t, h, pmw, umw, vmw, tmw, hmw)
mxwind() computes maximum wind level fields.
Definition: GFSPOST.F:459
mptgen
subroutine mptgen(mpirank, mpisize, nd, jt1, jt2, j1, j2, jx, jm, jn)
mptgen() generates grid decomposition dimensions.
Definition: GFSPOST.F:547
p2th
subroutine p2th(km, theta, u, v, h, t, pvu, sigma, rh, omga, kth, th, lth, uth, vth, hth, tth, zth, sigmath, rhth, oth)
p2th() interpolates from pressure level to isentropic (theta) level.
Definition: GFSPOST.F:137
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:649
pvetc
subroutine pvetc(km, p, px, py, t, tx, ty, h, u, v, av, hm, s, bvf2, pvn, theta, sigma, pvu)
pvetc() computes potential vorticity, etc.
Definition: GFSPOST.F:62
physcons_post
Definition: physcons.f:1
p2pv
subroutine p2pv(km, pvu, h, t, p, u, v, kpv, pv, pvpt, pvpb, lpv, upv, vpv, hpv, tpv, ppv, spv)
p2pv() interpolates to potential vorticity level.
Definition: GFSPOST.F:203
tpause
subroutine tpause(km, p, u, v, t, h, ptp, utp, vtp, ttp, htp, shrtp)
tpause() computes tropopause level fields.
Definition: GFSPOST.F:372
machine_post
Definition: machine.f:1
rsearch1
subroutine rsearch1(km1, z1, km2, z2, l2)
rsearch1() searches for a surrounding real interval.
Definition: GFSPOST.F:297