subroutine tpause(mype,method)
!$$$ subprogram documentation block
! . . . .
! subprogram: tpause locate tropopause
! prgmmr: treadon org: np23 date: 2003-09-13
!
! abstract: locate tropopause using one of two methods. The default
! method uses the temperature lapse rate to identify the
! tropopause. Method 'pvoz' uses a combination of the
! potential voriticy and ozone mixing ration to locate
! the tropopause.
!
! program history log:
! 2003-09-13 treadon - initial routine
! 2003-12-23 kleist - generalized to use guess pressure
! 2004-06-15 treadon - reformat documentation
! 2004-07-26 treadon - remove call smooth121 (leads to different
! results when running code with different
! number of mpi tasks)
! 2004-07-27 treadon - add use only; add intent in/out
! 2005-11-29 derber - remove psfcg and use ges_lnps instead
! 2006-02-02 treadon - rename prsl as ges_prsl
! 2006-07-28 derber - use r1000 from constants
! - use sensible temperature rather than virtual
! 2006-07-31 kleist - change to ges_ps from ln(ps)
! 2008-04-03 safford - rm unused vars and uses
! 2013-10-19 todling - metguess now holds background
! revised how method is chosen based on guess fields
!
! input argument list:
! mype - mpi task id
! method - method used to locate tropopause
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm RS/6000 SP
!
!$$$
use kinds, only: r_kind,i_kind
use constants, only: rd_over_cp,grav,rad2deg,one,r1000,r0_01
use guess_grids, only: tropprs,geop_hgtl,&
ntguessig,ges_prsl,ges_tsen
use gridmod, only: istart,nlat,rlats,nsig,lat2,lon2
use gsi_bundlemod, only : gsi_bundlegetpointer
use gsi_metguess_mod, only : gsi_metguess_bundle
implicit none
! Declare passed variables
character(4) ,intent(in ) :: method
integer(i_kind),intent(in ) :: mype
! Declare local parameters
character(len=*),parameter::myname='tpause'
real(r_kind),parameter:: r3em7=3.0e-7_r_kind
real(r_kind),parameter:: r2em6=2.0e-6_r_kind
real(r_kind),parameter:: r0_001=0.001_r_kind
real(r_kind),parameter:: r0_7=0.7_r_kind
real(r_kind),parameter:: r20=20.0_r_kind
real(r_kind),parameter:: r40=40.0_r_kind
real(r_kind),parameter:: r1e5=1.0e5_r_kind
! Declare local variables
logical t_method,pvoz_capable
integer(i_kind) i,j,k,mm1,nt,istatus
integer(i_kind) latrad
integer(i_kind) ifound_pv,ifound_oz,itrp_pv,itrp_oz,itrop_k
real(r_kind) pm1,pp1
real(r_kind) thetam1,thetap1,pv,wgt1,ptrop
real(r_kind),dimension(nsig):: prs,tdry,hgt,pvort
real(r_kind),dimension(lat2):: slatd
real(r_kind),dimension(lat2,lon2):: trop_t,trop_pv,trop_oz
real(r_kind),dimension(lat2,lon2):: trop_pvoz
real(r_kind) psi
real(r_kind),dimension(:,:) ,pointer :: ges_ps_nt=>NULL()
real(r_kind),dimension(:,:,:),pointer :: ges_tv_nt=>NULL()
real(r_kind),dimension(:,:,:),pointer :: ges_oz_nt=>NULL()
real(r_kind),dimension(:,:,:),pointer :: ges_vor_nt=>NULL()
!================================================================================
! Set local constants
t_method = .false.
if (index(method,'pvoz') == 0) t_method = .true.
nt=ntguessig
call gsi_bundlegetpointer (gsi_metguess_bundle(nt),'ps',ges_ps_nt,istatus)
if(istatus/=0) return ! if ps not defined forget it ...
call gsi_bundlegetpointer (gsi_metguess_bundle(nt),'tv',ges_tv_nt,istatus)
pvoz_capable=istatus==0
call gsi_bundlegetpointer (gsi_metguess_bundle(nt),'oz',ges_oz_nt,istatus)
pvoz_capable=pvoz_capable.and.istatus==0
call gsi_bundlegetpointer (gsi_metguess_bundle(nt),'vor',ges_vor_nt,istatus)
pvoz_capable=pvoz_capable.and.istatus==0
if (trim(method)=='pvoz') then
if(.not.pvoz_capable) then
if (mype==0) then
write(6,*) trim(myname), ': Warning, user request pv-based method to',&
'identify tropaupose,'
write(6,*) 'but not all fields needed are available, resetting',&
'method to temperature-based one'
endif
t_method=.true.
endif
endif
! Locate tropopause based on temperature profile (WMO approach)
if (t_method) then
do j=1,lon2
do i=1,lat2
do k=1,nsig
prs(k) = r1000*ges_prsl(i,j,k,ntguessig)
tdry(k)= ges_tsen(i,j,k,ntguessig)
hgt(k) = geop_hgtl(i,j,k,ntguessig)
end do
call tpause_t(nsig,prs,tdry,hgt,ptrop)
trop_t(i,j) = ptrop*r0_01 !hPa
end do
end do
! Load tropopause pressure (hPa) into output array (passed through module guess_grids)
do j=1,lon2
do i=1,lat2
tropprs(i,j) = trop_t(i,j)
end do
end do
! Locate tropopause using combination of potential vorticity (pv) and ozone.
else
! Compute latitudes on subdomain
mm1=mype+1
do i=1,lat2
latrad=min(max(1,istart(mm1)+i-2),nlat)
slatd(i)=abs(rlats(latrad))*rad2deg
end do
! Compute pv. Use for locating tropopause poleward 30S/N
do j=1,lon2
do i = 1,lat2
psi=one/ges_ps_nt(i,j)
do k=1,nsig
prs(k) = r1000*ges_prsl(i,j,k,ntguessig)
end do
! Compute pv
do k = 2,nsig-1
pm1 = prs(k-1)
pp1 = prs(k+1)
thetam1 = ges_tv_nt(i,j,k-1)*(r1e5/pm1)**(rd_over_cp)
thetap1 = ges_tv_nt(i,j,k+1)*(r1e5/pp1)**(rd_over_cp)
pv = grav*ges_vor_nt(i,j,k)*(thetam1-thetap1)/(pm1-pp1)
pvort(k) = abs(pv)
end do
pvort(1) = pvort(2)
pvort(nsig) = pvort(nsig-1)
! Locate tropopause
ifound_pv=0; ifound_oz=0
itrp_pv=nsig; itrp_oz=nsig
! Search upward (decreasing pressure) for tropopause above sigma 0.7
do k=2,nsig
if ((prs(k)*(r0_001)*psi) < r0_7) then
! Trop at level where pv=2e-6
if (pvort(k)>r2em6 .and. ifound_pv==0) then
ifound_pv=1
itrp_pv = k
endif
! Trop at level where ozone greater than 3e-7
if (ges_oz_nt(i,j,k)>r3em7 .and. ifound_oz==0) then
ifound_oz=1
itrp_oz = k
endif
endif
end do
! Merge pv and ozone tropopause levels between 20 and 40 deg latitude
if (slatd(i) >= r40 ) then
itrop_k = itrp_pv
elseif (slatd(i) >= r20) then
wgt1 = (slatd(i)-r20)/r20
itrop_k = wgt1*itrp_pv + (one-wgt1)*itrp_oz
else
itrop_k = itrp_oz
endif
itrop_k = max(1,min(itrop_k,nsig))
trop_pv(i,j) = prs(itrp_pv)*r0_01 !hPa
trop_oz(i,j) = prs(itrp_oz)*r0_01 !hPa
trop_pvoz(i,j) = prs(itrop_k)*r0_01 !hPa
end do
end do
! Load tropopause pressure (hPa) into output array
do j=1,lon2
do i=1,lat2
tropprs(i,j) = trop_pvoz(i,j)
end do
end do
! End of tropopause location method blocks
endif
! *** NOTE ***
! The tropopause pressures are used to deflate the
! moisture sensitivity vectors for satellite radiance
! data and for IR quality control;
! here we are setting bounds on the tropopause
! pressure to make sure we are deflating at the very
! minimum above 150 mb, and nowhere below 350 mb
do j=1,lon2
do i=1,lat2
tropprs(i,j)=max(150.0_r_kind,min(350.0_r_kind,tropprs(i,j)))
end do
end do
! End of routine
return
end subroutine tpause