subroutine pvetc(km,p,px,py,t,tx,ty,h,u,v,av,hm,s,bvf2,pvn,theta,sigma,pvu)
!$$$ Subprogram documentation block
!
! Subprogram: pvetc Compute potential vorticity, etc
! Prgmmr: Iredell Org: np23 Date: 1999-10-18
!
! Abstract: This subprogram computes
! the Montgomery streamfunction
! hm=cp*t+g*z
! the specific entropy
! s=cp*log(t/t0)-r*log(p/p0)
! the Brunt-Vaisala frequency squared
! bvf2=g/cp*ds/dz
! the potential vorticity defined as
! pvn=(av*ds/dz-dv/dz*ds/dx+du/dz*ds/dy)/rho/cp
! the potential temperature
! theta=t0*exp(s/cp)
! the static stability
! sigma=t/g*bvf2
! and the potential vorticity in PV units
! pvu=10**-6*theta*pvn
!
! Program history log:
! 1999-10-18 Mark Iredell
!
! Usage: call pvetc(km,p,px,py,t,tx,ty,h,u,v,av,s,bvf2,pvn,theta,sigma,pvu)
! Input argument list:
! km integer number of levels
! p real (km) pressure (Pa)
! px real (km) pressure x-gradient (Pa/m)
! py real (km) pressure y-gradient (Pa/m)
! t real (km) (virtual) temperature (K)
! tx real (km) (virtual) temperature x-gradient (K/m)
! ty real (km) (virtual) temperature y-gradient (K/m)
! h real (km) height (m)
! u real (km) x-component wind (m/s)
! v real (km) y-component wind (m/s)
! av real (km) absolute vorticity (1/s)
! Output argument list:
! hm real (km) Montgomery streamfunction (m**2/s**2)
! s real (km) specific entropy (J/K/kg)
! bvf2 real (km) Brunt-Vaisala frequency squared (1/s**2)
! pvn real (km) potential vorticity (m**2/kg/s)
! theta real (km) (virtual) potential temperature (K)
! sigma real (km) static stability (K/m)
! pvu real (km) potential vorticity (10**-6*K*m**2/kg/s)
!
! Modules used:
! physcons Physical constants
!
! Attributes:
! Language: Fortran 90
!
!$$$
use physcons_post, only: con_cp, con_g, con_rd, con_rocp
!
implicit none
integer,intent(in):: km
real,intent(in), dimension(km):: p,px,py,t,tx,ty,h,u,v,av
real,intent(out),dimension(km):: hm,s,bvf2,pvn,theta,sigma,pvu
! real,parameter:: hhmin=500.,t0=2.e2,p0=1.e5
real,parameter:: hhmin=5.,t0=2.e2,p0=1.e5
integer k,kd,ku,k2(2)
real cprho,sx,sy,sz,uz,vz
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
do k=1,km
hm(k) = con_cp*t(k) + con_g*h(k)
s(k) = con_cp*log(t(k)/t0) - con_rd*log(p(k)/p0)
enddo
do k=1,km
call rsearch1(km,h,2,(/h(k)-hhmin,h(k)+hhmin/),k2)
! kd = max(k2(1),1)
! ku = min(k2(2)+1,km)
! kd = min(k2(1),km) ! Chuang: post counts from top down, redefine lower bound
kd = min(k2(1)+1,km) ! Chuang: post counts from top down,
! ku = max(k2(2)-1,1)
ku = max(k2(2),1)
if(ku==1) kd=2 ! Chuang: make sure ku ne kd at model top
cprho = p(k)/(con_rocp*t(k))
sx = con_cp*tx(k) / t(k)-con_rd*px(k)/p(k)
sy = con_cp*ty(k) / t(k)-con_rd*py(k)/p(k)
sz = (s(ku)-s(kd)) / (h(ku)-h(kd))
uz = (u(ku)-u(kd)) / (h(ku)-h(kd))
vz = (v(ku)-v(kd)) / (h(ku)-h(kd))
bvf2(k) = con_g/con_cp*sz
pvn(k) = (av(k)*sz - vz*sx + uz*sy) / cprho
theta(k) = t0*exp(s(k)/con_cp)
sigma(k) = t(k)/con_g*bvf2(k)
pvu(k) = 1.e6*theta(k)*pvn(k)
enddo
end subroutine
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine p2th(km,theta,u,v,h,t,pvu,sigma,rh,omga,kth,th &
,lth,uth,vth,hth,tth,zth,sigmath,rhth,oth)
!$$$ Subprogram documentation block
!
! Subprogram: p2th Interpolate to isentropic level
! Prgmmr: Iredell Org: np23 Date: 1999-10-18
!
! Abstract: This subprogram interpolates fields to given isentropic levels.
! The interpolation is linear in entropy.
! Outside the domain the bitmap is set to false.
!
! Program history log:
! 1999-10-18 Mark Iredell
!
! Usage: call p2th(km,theta,u,v,h,t,puv,kth,th,uth,vth,tth)
! Input argument list:
! km integer number of levels
! theta real (km) potential temperature (K)
! u real (km) x-component wind (m/s)
! v real (km) y-component wind (m/s)
! h real (km) height (m)
! t real (km) temperature (K)
! pvu real (km) potential vorticity in PV units (10**-6*K*m**2/kg/s)
! kth integer number of isentropic levels
! th real (kth) isentropic levels (K)
! Output argument list:
! lpv logical*1 (kth) bitmap
! uth real (kth) x-component wind (m/s)
! vth real (kth) y-component wind (m/s)
! hth real (kth) height (m)
! tth real (kth) temperature (K)
! zth real (kth) potential vorticity in PV units (10**-6*K*m**2/kg/s)
!
! Subprograms called:
! rsearch1 search for a surrounding real interval
!
! Attributes:
! Language: Fortran 90
!
!$$$
implicit none
integer,intent(in):: km,kth
real,intent(in),dimension(km):: theta,u,v,h,t,pvu,sigma,rh,omga
real,intent(in):: th(kth)
logical*1,intent(out),dimension(kth):: lth
real,intent(out),dimension(kth):: uth,vth,hth,tth,zth &
,sigmath,rhth,oth
real w
integer loc(kth),l
integer k
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
call rsearch1(km,theta(1),kth,th(1),loc(1))
do k=1,kth
l = loc(k)
lth(k) = l > 0 .and.l < km
if(lth(k)) then
w = log(th(k)/theta(l)) / log(theta(l+1)/theta(l))
uth(k) = u(l) + w*(u(l+1)-u(l))
vth(k) = v(l) + w*(v(l+1)-v(l))
hth(k) = h(l) + w*(h(l+1)-h(l))
tth(k) = t(l) + w*(t(l+1)-t(l))
zth(k) = pvu(l) + w*(pvu(l+1)-pvu(l))
sigmath(k) = sigma(l) + w*(sigma(l+1)-sigma(l))
rhth(k) = rh(l) + w*(rh(l+1)-rh(l))
! pth(k) = p(l) + w*(p(l+1)-p(l))
oth(k) = omga(l) + w*(omga(l+1)-omga(l))
endif
enddo
end subroutine
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine p2pv(km,pvu,h,t,p,u,v,kpv,pv,pvpt,pvpb,&
lpv,upv,vpv,hpv,tpv,ppv,spv)
!$$$ Subprogram documentation block
!
! Subprogram: p2pv Interpolate to potential vorticity level
! Prgmmr: Iredell Org: np23 Date: 1999-10-18
!
! Abstract: This subprogram interpolates fields to given potential vorticity
! levels within given pressure limits.
! The output level is the first encountered from the top pressure limit.
! If the given potential vorticity level is not found, the outputs are zero
! and the bitmap is false. The interpolation is linear in potential vorticity.
!
! Program history log:
! 1999-10-18 Mark Iredell
!
! Usage: call p2pv(km,pvu,h,t,p,u,v,kpv,pv,pvpt,pvpb,&
! lpv,upv,vpv,hpv,tpv,ppv,spv)
! Input argument list:
! km integer number of levels
! pvu real (km) potential vorticity in PV units (10**-6*K*m**2/kg/s)
! h real (km) height (m)
! t real (km) temperature (K)
! p real (km) pressure (Pa)
! u real (km) x-component wind (m/s)
! v real (km) y-component wind (m/s)
! kpv integer number of potential vorticity levels
! pv real (kpv) potential vorticity levels (10**-6*K*m**2/kg/s)
! pvpt real (kpv) top pressures for PV search (Pa)
! pvpb real (kpv) bottom pressures for PV search (Pa)
! Output argument list:
! lpv logical*1 (kpv) bitmap
! upv real (kpv) x-component wind (m/s)
! vpv real (kpv) y-component wind (m/s)
! hpv real (kpv) temperature (K)
! tpv real (kpv) temperature (K)
! ppv real (kpv) pressure (Pa)
! spv real (kpv) wind speed shear (1/s)
!
! Subprograms called:
! rsearch1 search for a surrounding real interval
!
! Attributes:
! Language: Fortran 90
!
!$$$
use physcons_post, only: con_rog
implicit none
integer,intent(in):: km,kpv
real,intent(in),dimension(km):: pvu,h,t,p,u,v
real,intent(in):: pv(kpv),pvpt(kpv),pvpb(kpv)
logical*1,intent(out),dimension(kpv):: lpv
real,intent(out),dimension(kpv):: upv,vpv,hpv,tpv,ppv,spv
real,parameter:: pd=2500.
real w,spdu,spdd
integer k,l1,l2,lu,ld,l
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
do k=1,kpv
call rsearch1(km,p,1,pvpb(k),l1)
call rsearch1(km,p,1,pvpt(k),l2)
! l1=l1+1
l = 0
if(pv(k) >= 0.) then
! do lu=l2-1,l1,-1
! do lu=l2,l1-1 ! Chuang: post counts top down
do lu=l2+2,l1 ! Chuang: post counts top down
! if(pv(k).lt.pvu(lu+1).and.pv(k).ge.pvu(lu)) then
if(pv(k) >= pvu(lu+1).and.pv(k) < pvu(lu)) then
call rsearch1(km,p,1,p(lu)+pd,ld)
! if(all(pv(k).ge.pvu(ld:lu-1))) then
if(all(pv(k) >= pvu(lu+1:ld))) then
l = lu
exit
endif
endif
enddo
else
! do lu=l2-1,l1,-1
! do lu=l2,l1-1 ! Chuang: post counts top down
do lu=l2+2,l1 ! Chuang: post counts top down
! if(pv(k).gt.pvu(lu+1).and.pv(k).le.pvu(lu)) then
if(pv(k) <= pvu(lu+1).and.pv(k) > pvu(lu)) then
call rsearch1(km,p,1,p(lu)+pd,ld)
! if(all(pv(k).le.pvu(ld:lu-1))) then
if(all(pv(k) <= pvu(lu+1:ld))) then
l = lu
exit
endif
endif
enddo
endif
lpv(k) = l > 0
if(lpv(k)) then
w = (pv(k)-pvu(l))/(pvu(l+1)-pvu(l))
upv(k) = u(l) + w*(u(l+1)-u(l))
vpv(k) = v(l) + w*(v(l+1)-v(l))
hpv(k) = h(l) + w*(h(l+1)-h(l))
tpv(k) = t(l) + w*(t(l+1)-t(l))
ppv(k) = p(l)*exp((h(l)-hpv(k))*(1-0.5*(tpv(k)/t(l)-1))/(con_rog*t(l)))
spdu = sqrt(u(l+1)*u(l+1) + v(l+1)*v(l+1))
spdd = sqrt(u(l)*u(l) + v(l)*v(l))
spv(k) = (spdu-spdd) / (h(l+1)-h(l))
endif
enddo
end subroutine
!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
subroutine rsearch1(km1,z1,km2,z2,l2)
!$$$ subprogram documentation block
!
! subprogram: rsearch1 search for a surrounding real interval
! prgmmr: iredell org: w/nmc23 date: 98-05-01
!
! abstract: this subprogram searches a monotonic sequences of real numbers
! for intervals that surround a given search set of real numbers.
! the sequences may be monotonic in either direction; the real numbers
! may be single or double precision.
!
! program history log:
! 1999-01-05 mark iredell
!
! usage: call rsearch1(km1,z1,km2,z2,l2)
! input argument list:
! km1 integer number of points in the sequence
! z1 real (km1) sequence values to search
! (z1 must be monotonic in either direction)
! km2 integer number of points to search for
! z2 real (km2) set of values to search for
! (z2 need not be monotonic)
!
! output argument list:
! l2 integer (km2) interval locations from 0 to km1
! (z2 will be between z1(l2) and z1(l2+1))
!
! subprograms called:
! sbsrch essl binary search
! dbsrch essl binary search
!
! remarks:
! returned values of 0 or km1 indicate that the given search value
! is outside the range of the sequence.
!
! if a search value is identical to one of the sequence values
! then the location returned points to the identical value.
! if the sequence is not strictly monotonic and a search value is
! identical to more than one of the sequence values, then the
! location returned may point to any of the identical values.
!
! if l2(k)=0, then z2(k) is less than the start point z1(1)
! for ascending sequences (or greater than for descending sequences).
! if l2(k)=km1, then z2(k) is greater than or equal to the end point
! z1(km1) for ascending sequences (or less than or equal to for
! descending sequences). otherwise z2(k) is between the values
! z1(l2(k)) and z1(l2(k+1)) and may equal the former.
!
! attributes:
! language: fortran
!
!$$$
implicit none
integer,intent(in):: km1,km2
real,intent(in):: z1(km1),z2(km2)
integer,intent(out):: l2(km2)
integer k1,k2
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! find the surrounding input interval for each output point.
if(z1(1) <= z1(km1)) then
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! input coordinate is monotonically ascending.
do k2=1,km2
if (z1(1) >= z2(k2)) then
l2(k2) = 1
else
l2(k2)=km1
do k1=1,km1-1
if(z1(k1) <= z2(k2) .and. z1(k1+1) > z2(k2)) then
l2(k2) = k1
exit
endif
enddo
endif
enddo
else
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! input coordinate is monotonically descending.
do k2=1,km2
if (z1(1) <= z2(k2)) then
l2(k2) = 1
else
l2(k2)=km1
do k1=km1,2,-1
if(z2(k2) >= z1(k1) .and. z2(k2) < z1(k1-1)) then
l2(k2) = k1-1
exit
endif
enddo
endif
enddo
endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end subroutine
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine tpause(km,p,u,v,t,h,ptp,utp,vtp,ttp,htp,shrtp)
!$$$ Subprogram documentation block
!
! Subprogram: tpause Compute tropopause level fields
! Prgmmr: Iredell Org: np23 Date: 1999-10-18
!
! Abstract: This subprogram finds the tropopause level and computes fields
! at the tropopause level. The tropopause is defined as the lowest level
! above 500 mb which has a temperature lapse rate of less than 2 K/km.
! The lapse rate must average less than 2 K/km over a 2 km depth.
! If no such level is found below 50 mb, the tropopause is set to 50 mb.
! The tropopause fields are interpolated linearly in lapse rate.
! The tropopause pressure is found hydrostatically.
! The tropopause wind shear is computed as the partial derivative
! of wind speed with respect to height at the tropopause level.
!
! Program history log:
! 1999-10-18 Mark Iredell
!
! Usage: call tpause(km,p,u,v,t,h,ptp,utp,vtp,ttp,htp,shrtp)
! Input argument list:
! km integer number of levels
! p real (km) pressure (Pa)
! u real (km) x-component wind (m/s)
! v real (km) y-component wind (m/s)
! t real (km) temperature (K)
! h real (km) height (m)
! Output argument list:
! ptp real tropopause pressure (Pa)
! utp real tropopause x-component wind (m/s)
! vtp real tropopause y-component wind (m/s)
! ttp real tropopause temperature (K)
! htp real tropopause height (m)
! shrtp real tropopause wind shear (1/s)
!
! Files included:
! physcons.h Physical constants
!
! Subprograms called:
! rsearch1 search for a surrounding real interval
!
! Attributes:
! Language: Fortran 90
!
!$$$
use physcons_post, only: con_rog
implicit none
integer,intent(in):: km
real,intent(in),dimension(km):: p,u,v,t,h
real,intent(out):: ptp,utp,vtp,ttp,htp,shrtp
real,parameter:: ptplim(2)=(/500.e+2,50.e+2/),gamtp=2.e-3,hd=2.e+3
real gamu,gamd,td,gami,wtp,spdu,spdd
integer klim(2),k,kd,ktp
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! find tropopause level
call rsearch1(km-2,p(2),2,ptplim(1),klim(1))
klim(1)=klim(1)+1
klim(2)=klim(2)+2
! klim(1) > klim(2) or loops does not run ; klim(2) has a
! minimum value of 3 to insure k-2 != 0 in called subprogram
gamd=1.e+9
ktp=klim(2)
wtp=0
! do k=klim(1),klim(2)
do k=klim(1),klim(2),-1
! gamu=(t(k-1)-t(k+1))/(h(k+1)-h(k-1))
gamu=(t(k+1)-t(k-1))/(h(k-1)-h(k+1))
if(gamu.le.gamtp) then
! call rsearch1(km-k-1,h(k+1),1,h(k)+hd,kd)
call rsearch1(k-2,h(2),1,h(k)+hd,kd)
! td=t(k+kd)+(h(k)+hd-h(k+kd))/(h(k+kd+1)-h(k+kd))*(t(k+kd+1)-t(k+kd))
td=t(kd+2)+(h(k)+hd-h(2+kd))/(h(kd+1)-h(2+kd))*(t(kd+1)-t(2+kd))
gami=(t(k)-td)/hd
if(gami.le.gamtp) then
ktp=k
wtp=(gamtp-gamu)/(max(gamd,gamtp+0.1e-3)-gamu)
exit
endif
endif
gamd=gamu
enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! compute tropopause level fields
utp=u(ktp)-wtp*(u(ktp)-u(ktp-1))
vtp=v(ktp)-wtp*(v(ktp)-v(ktp-1))
ttp=t(ktp)-wtp*(t(ktp)-t(ktp-1))
htp=h(ktp)-wtp*(h(ktp)-h(ktp-1))
ptp=p(ktp)*exp((h(ktp)-htp)*(1-0.5*(ttp/t(ktp)-1))/(con_rog*t(ktp)))
spdu=sqrt(u(ktp)**2+v(ktp)**2)
spdd=sqrt(u(ktp-1)**2+v(ktp-1)**2)
shrtp=(spdu-spdd)/(h(ktp)-h(ktp-1))
utp=u(ktp)-wtp*(u(ktp)-u(ktp+1))
vtp=v(ktp)-wtp*(v(ktp)-v(ktp+1))
ttp=t(ktp)-wtp*(t(ktp)-t(ktp+1))
htp=h(ktp)-wtp*(h(ktp)-h(ktp+1))
ptp=p(ktp)*exp((h(ktp)-htp)*(1-0.5*(ttp/t(ktp)-1))/(con_rog*t(ktp)))
spdu=sqrt(u(ktp)**2+v(ktp)**2)
spdd=sqrt(u(ktp+1)**2+v(ktp+1)**2)
shrtp=(spdu-spdd)/(h(ktp)-h(ktp+1))
end subroutine
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
subroutine mxwind(km,p,u,v,t,h,pmw,umw,vmw,tmw,hmw)
!$$$ Subprogram documentation block
!
! Subprogram: mxwind Compute maximum wind level fields
! Prgmmr: Iredell Org: np23 Date: 1999-10-18
!
! Abstract: This subprogram finds the maximum wind level and computes fields
! at the maximum wind level. The maximum wind level is searched for
! between 500 mb and 100 mb. The height and wind speed at the maximum wind
! speed level is calculated by assuming the wind speed varies quadratically
! in height in the neighborhood of the maximum wind level. The other fields
! are interpolated linearly in height to the maximum wind level.
! The maximum wind level pressure is found hydrostatically.
!
! Program history log:
! 1999-10-18 Mark Iredell
! 2005-02-02 Mark Iredell changed upper limit to 100 mb
!
! Usage: call mxwind(km,p,u,v,t,h,pmw,umw,vmw,tmw,hmw)
! Input argument list:
! km integer number of levels
! p real (km) pressure (Pa)
! u real (km) x-component wind (m/s)
! v real (km) y-component wind (m/s)
! t real (km) temperature (K)
! h real (km) height (m)
! Output argument list:
! pmw real maximum wind level pressure (Pa)
! umw real maximum wind level x-component wind (m/s)
! vmw real maximum wind level y-component wind (m/s)
! tmw real maximum wind level temperature (K)
! hmw real maximum wind level height (m)
!
! Files included:
! physcons.h Physical constants
!
! Subprograms called:
! rsearch1 search for a surrounding real interval
!
! Attributes:
! Language: Fortran 90
!
!$$$
use physcons_post, only: con_rog
implicit none
integer,intent(in):: km
real,intent(in),dimension(km):: p,u,v,t,h
real,intent(out):: pmw,umw,vmw,tmw,hmw
real,parameter:: pmwlim(2)=(/500.e+2,100.e+2/)
integer klim(2),k,kmw
real spd(km),spdmw,wmw,dhd,dhu,shrd,shru,dhmw,ub,vb,spdb
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! find maximum wind level
call rsearch1(km,p(1),2,pmwlim(1),klim(1))
! klim(1)=klim(1)+1
klim(2)=klim(2)+1
! spd(klim(1):klim(2))=sqrt(u(klim(1):klim(2))**2+v(klim(1):klim(2))**2)
spd(klim(2):klim(1))=sqrt(u(klim(2):klim(1))**2+v(klim(2):klim(1))**2)
spdmw=spd(klim(1))
kmw=klim(1)
! do k=klim(1)+1,klim(2)
do k=klim(1)-1,klim(2),-1
if(spd(k).gt.spdmw) then
spdmw=spd(k)
kmw=k
endif
enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! find speed and height at the maximum wind level
if(kmw.eq.klim(1).or.kmw.eq.klim(2)) then
hmw=h(kmw)
spdmw=spd(kmw)
wmw=0.
else
! dhd=h(kmw)-h(kmw-1)
dhd=h(kmw)-h(kmw+1) !post counts top down
! dhu=h(kmw+1)-h(kmw)
dhu=h(kmw-1)-h(kmw)
! shrd=(spd(kmw)-spd(kmw-1))/(h(kmw)-h(kmw-1))
shrd=(spd(kmw)-spd(kmw+1))/(h(kmw)-h(kmw+1))
! shru=(spd(kmw)-spd(kmw+1))/(h(kmw+1)-h(kmw))
shru=(spd(kmw)-spd(kmw-1))/(h(kmw-1)-h(kmw))
dhmw=(shrd*dhu-shru*dhd)/(2*(shrd+shru))
hmw=h(kmw)+dhmw
spdmw=spd(kmw)+dhmw**2*(shrd+shru)/(dhd+dhu)
! if(dhmw.gt.0) kmw=kmw+1
if(dhmw.gt.0) kmw=kmw-1
! wmw=(h(kmw)-hmw)/(h(kmw)-h(kmw-1))
wmw=(h(kmw)-hmw)/(h(kmw)-h(kmw+1))
endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! compute maximum wind level fields
! ub=u(kmw)-wmw*(u(kmw)-u(kmw-1))
ub=u(kmw)-wmw*(u(kmw)-u(kmw+1))
! vb=v(kmw)-wmw*(v(kmw)-v(kmw-1))
vb=v(kmw)-wmw*(v(kmw)-v(kmw+1))
spdb=max(sqrt(ub**2+vb**2),1.e-6)
umw=ub*spdmw/spdb
vmw=vb*spdmw/spdb
! tmw=t(kmw)-wmw*(t(kmw)-t(kmw-1))
tmw=t(kmw)-wmw*(t(kmw)-t(kmw+1))
pmw=p(kmw)*exp((h(kmw)-hmw)*(1-0.5*(tmw/t(kmw)-1))/(con_rog*t(kmw)))
end subroutine
subroutine mptgen(mpirank,mpisize,nd,jt1,jt2,j1,j2,jx,jm,jn)
!$$$ Subprogram documentation block
!
! Subprogram: mptgen Generate grid decomposition dimensions
! Prgmmr: Iredell Org: W/NP23 Date: 1999-02-12
!
! Abstract: This subprogram decomposes total dimensions of a problem
! into smaller domains to be managed on a distributed memory system.
! The last dimension given is decomposed first. If more decompositions
! are possible, the next to last dimension is decomposed next, and so on.
! The transpositions between decompositions should be done by mptran*.
!
! Program history log:
! 1999-02-12 Mark Iredell
!
! Usage: call mptgen(mpirank,mpisize,nd,jt1,jt2,j1,j2,jx,jm,jn)
! Input argument list:
! mpirank integer(kint_mpi) rank of the process (from mpi_comm_rank)
! mpisize integer(kint_mpi) size of the process (from mpi_comm_size)
! nd integer(kint_mpi) number of dimensions to decompose
! jt1 integer(kint_mpi) (nd) lower bounds of total dimensions
! jt2 integer(kint_mpi) (nd) upper bounds of total dimensions
! Output argument list:
! j1 integer(kint_mpi) (nd) lower bounds of local decompositions
! j2 integer(kint_mpi) (nd) upper bounds of local decompositions
! jx integer(kint_mpi) (nd) local size of decompositions
! jm integer(kint_mpi) (nd) maximum size of decompositions
! jn integer(kint_mpi) (nd) number of decompositions
!
! Attributes:
! Language: Fortran 90
!
!$$$
use machine_post,only:kint_mpi
implicit none
integer(kint_mpi),intent(in):: mpirank,mpisize,nd,jt1(nd),jt2(nd)
integer(kint_mpi),intent(out):: j1(nd),j2(nd),jx(nd),jm(nd),jn(nd)
integer msize,mrank,msn,mrn,n
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
msize=mpisize
mrank=mpirank
do n=nd,1,-1
if(jt2(n).ge.jt1(n)) then
jm(n)=(jt2(n)-jt1(n))/msize+1
msn=max(msize/(jt2(n)-jt1(n)+1),1)
if(n.eq.1) msn=1
jn(n)=msize/msn
mrn=mrank/msn
j1(n)=min(jt1(n)+jm(n)*mrn,jt2(n)+1)
j2(n)=min(jt1(n)+jm(n)*mrn+jm(n)-1,jt2(n))
jx(n)=j2(n)-j1(n)+1
msize=msn
mrank=mod(mrank,msn)
write(0,*)' mrank=',mrank,' j1=',j1(n),' j2=',j2(n),' jx=',jx(n),' jm=',jm(n)
else
jm(n)=0
jn(n)=1
j1(n)=jt1(n)
j2(n)=jt2(n)
jx(n)=0
endif
enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end subroutine
!-------------------------------------------------------------------------------
subroutine mptranr4(mpicomm,mpisize,im,ida,idb,&
jm,jma,jmb,jda,km,kma,kmb,kdb,a,b,ta,tb)
!$$$ Subprogram documentation block
!
! Subprogram: mptranr4 Transpose grid decompositions
! Prgmmr: Iredell Org: W/NP23 Date: 1999-02-12
!
! Abstract: This subprogram transposes an array of data from one
! grid decomposition to another by using message passing.
! The grid decompositions should be generated by mptgen.
!
! Program history log:
! 1999-02-12 Mark Iredell
!
! Usage: call mptranr4(mpicomm,mpisize,im,ida,idb,&
! jm,jma,jmb,jda,km,kma,kmb,kdb,a,b)
! Input argument list:
! mpicomm integer(kint_mpi) mpi communicator
! mpisize integer(kint_mpi) size of the process (from mpi_comm_size)
! im integer(kint_mpi) undecomposed range
! ida integer(kint_mpi) undecomposed input dimension
! idb integer(kint_mpi) undecomposed output dimension
! jm integer(kint_mpi) output grid decomposition size
! jma integer(kint_mpi) input grid undecomposed range
! jmb integer(kint_mpi) output grid decomposed range
! jda integer(kint_mpi) input grid undecomposed dimension
! km integer(kint_mpi) input grid decomposition size
! kma integer(kint_mpi) input grid decomposed range
! kmb integer(kint_mpi) output grid undecomposed range
! kdb integer(kint_mpi) output grid undecomposed dimension
! a real(4) (ida,jda,kma) input array
! Output argument list:
! b real(4) (idb,kdb,jmb) output array
! ta,tb real(4) (im,jm,km,mpisize) work arrays
!
! Subprograms called:
! mpi_alltoall mpi exchange of data between every process pair
!
! Remarks:
! While this routine serves a wide variety of scalable transpose functions
! for multidimensional grids,
! (a) it does not work with nonrectanguloid grids;
! (b) it does not do any load balancing;
! (c) it does not do any communication hiding.
!
! This subprogram must be used rather than mpi_alltoall
! in any of the following cases:
! (a) The undecomposed range is less than the respective dimension
! (either im.lt.ida or im.lt.idb)
! (b) The decomposition size is greater than one
! (either km.gt.1 or jm.gt.1)
! (c) The decomposed range is ever zero
! (either kma.eq.0 or jmb.eq.0 for any process)
! (d) The output grid range is not the full extent
! (either kmb.lt.mpisize or kmb.lt.kda or jma.lt.mpisize or jma.lt.jda)
! If none of these conditions apply, mpi_alltoall could be used directly
! rather than this subprogram and would be more efficient.
!
! Example 1. Transpose a 1000 x 10000 matrix.
!
!!! include 'mpif.h' ! use mpi
!!! parameter(jt=1000,kt=10000) ! set problem size
!!! real,allocatable:: a(:,:),b(:,:) ! declare arrays
!!! call mpi_init(ierr) ! initialize mpi
!!! call mpi_comm_rank(MPI_COMM_WORLD,mpirank,ierr) ! get mpi rank
!!! call mpi_comm_size(MPI_COMM_WORLD,mpisize,ierr) ! get mpi size
!!! call mptgen(mpirank,mpisize,1,1,jt,j1,j2,jx,jm,jn) ! decompose output
!!! call mptgen(mpirank,mpisize,1,1,kt,k1,k2,kx,km,kn) ! decompose input
!!! allocate(a(jt,k1:k2),b(kt,j1:j2)) ! allocate arrays
!!! a=reshape((/((j+k,j=1,jt),k=k1,k2)/),(/jt,k2-k1+1/)) ! initialize input
!!! call mptranr4(MPI_COMM_WORLD,mpisize,1,1,1, ! transpose arrays
!!! & jm,jt,j2-j1+1,jt,km,k2-k1+1,kt,kt,a,b)
!!! print '(2i8,f16.1)',((k,j,b(k,j),k=2000,kt,2000), ! print some values
!!! & j=((j1-1)/200+1)*200,j2,200)
!!! call mpi_finalize(ierr) ! finalize mpi
!!! end
!
! This transpose took 0.6 seconds on 4 2-way winterhawk nodes.
! A 20000x10000 transpose took 3.4 seconds on 16 2-way winterhawk nodes.
! Thus a transpose may take about 1 second for every 16 Mb per node.
!
! Attributes:
! language: fortran 90
!
!$$$
use machine_post,only:kint_mpi
implicit none
include 'mpif.h'
integer(kint_mpi),intent(in):: mpicomm,mpisize
integer(kint_mpi),intent(in):: im,ida,idb
integer(kint_mpi),intent(in):: jm,jma,jmb,jda
integer(kint_mpi),intent(in):: km,kma,kmb,kdb
real(4),dimension(ida,jda,kma),intent(in):: a
real(4),dimension(idb,kdb,jmb),intent(out):: b
real(4),dimension(im,jm,km,mpisize),intent(inout):: ta,tb
integer(4) jmb1(1),jmbf(mpisize),kma1(1),kmaf(mpisize)
integer(kint_mpi)::i,j,k,l,ierr,ja,kb
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
jmb1(1) = jmb
call mpi_allgather(jmb1,1,MPI_INTEGER4,jmbf,1,MPI_INTEGER4,mpicomm,ierr)
kma1(1) = kma
call mpi_allgather(kma1,1,MPI_INTEGER4,kmaf,1,MPI_INTEGER4,mpicomm,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! internally transpose input array
!$omp parallel do private(i,j,k,l,ja)
do l=1,mpisize
do k=1,kma
do j=1,jm
ja = j + sum(jmbf(1:l-1))
if(ja <= jma) then
do i=1,im
ta(i,j,k,l) = a(i,ja,k)
enddo
endif
enddo
enddo
enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! externally transpose data
call mpi_alltoall(ta,im*jm*km,MPI_REAL4,tb,im*jm*km,MPI_REAL4,mpicomm,ierr)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! internally transpose output array
!$omp parallel do private(i,j,k,l,kb)
do l=1,mpisize
do k=1,km
kb = k + sum(kmaf(1:l-1))
if(kb <= kmb) then
do j=1,jmb
do i=1,im
b(i,kb,j) = tb(i,j,k,l)
enddo
enddo
endif
enddo
enddo
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
end subroutine
!-----------------------------------------------------------------------