UPP  001
 All Data Structures Files Functions Pages
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 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
311 end 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 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
550 end 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 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
674 end subroutine
675 !-----------------------------------------------------------------------
mxwind
subroutine mxwind(km, p, u, v, t, h, pmw, umw, vmw, tmw, hmw)
mxwind() computes maximum wind level fields.
Definition: GFSPOST.F:431
mptgen
subroutine mptgen(mpirank, mpisize, nd, jt1, jt2, j1, j2, jx, jm, jn)
mptgen() generates grid decomposition dimensions.
Definition: GFSPOST.F:519
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:109
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:621
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:175
tpause
subroutine tpause(km, p, u, v, t, h, ptp, utp, vtp, ttp, htp, shrtp)
tpause() computes tropopause level fields.
Definition: GFSPOST.F:344
rsearch1
subroutine rsearch1(km1, z1, km2, z2, l2)
rsearch1() searches for a surrounding real interval.
Definition: GFSPOST.F:269