#include <misc.h>
#include <params.h>
subroutine sphdep(jcen ,jgc ,dt ,ra ,iterdp , & 1,20
locgeo ,ub ,uxl ,uxr ,lam , &
phib ,lbasiy ,lammp ,phimp ,lamdp , &
phidp ,idp ,jdp ,vb ,vxl , &
vxr ,nlon ,nlonex )
!-----------------------------------------------------------------------
!
! Purpose:
! Compute departure points for semi-Lagrangian transport on surface of
! sphere using midpoint quadrature. Computations are done in:
!
! 1) "local geodesic" coordinates for "locgeo" = .true.
! 2) "global spherical" coordinates for "locgeo" = .false.
!
! Method:
!
! Author: J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: sphdep.F90,v 1.1.44.5.4.1 2004/05/20 18:36:06 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
use abortutils, only: endrun
implicit none
#include <parslt.h>
!------------------------------Arguments--------------------------------
integer , intent(in) :: nlon ! longitude dimension
integer , intent(in) :: nlonex(platd) ! extended longitude dimension
integer , intent(in) :: jcen ! index of lat slice (extnd)
integer , intent(in) :: jgc ! index of lat slice (model)
real(r8), intent(in) :: dt ! time step (seconds)
real(r8), intent(in) :: ra ! 1./(radius of earth)
integer , intent(in) :: iterdp ! number of iterations
logical , intent(in) :: locgeo ! computation type flag
real(r8), intent(in) :: ub (plond,plev,beglatex:endlatex) ! x-deriv
real(r8), intent(in) :: vb (plond,plev,beglatex:endlatex) ! x-deriv
real(r8), intent(in) :: uxl(plond,plev,beglatex:endlatex) ! left x-deriv (u)
real(r8), intent(in) :: uxr(plond,plev,beglatex:endlatex) ! right x-deriv
real(r8), intent(in) :: vxl(plond,plev,beglatex:endlatex) ! left x-deriv (v)
real(r8), intent(in) :: vxr(plond,plev,beglatex:endlatex) ! right x-deriv
real(r8), intent(in) :: lam(plond,platd) ! long. coord. of model grid
real(r8), intent(in) :: phib(platd) ! lat. coord. of model grid
real(r8), intent(in) :: lbasiy(4,2,platd) ! lat interpolation weights
real(r8), intent(inout) :: lammp(plon,plev) ! long coord of midpoint
real(r8), intent(inout) :: phimp(plon,plev) ! lat coord of midpoint
real(r8), intent(out) :: lamdp(plon,plev) ! long coord of dep. point
real(r8), intent(out) :: phidp(plon,plev) ! lat coord of dep. point
integer , intent(out) :: idp(plon,plev,4) ! long index of dep. point
integer , intent(out) :: jdp(plon,plev) ! lat index of dep. point
!
! jcen Index in extended grid corresponding to latitude being
! forecast.
! jgc Index in model grid corresponding to latitude being
! forecast.
! dt Time interval that parameterizes the parcel trajectory.
! ra Reciprocal of radius of earth.
! iterdp Number of iterations used for departure point calculation.
! locgeo Logical flag to indicate computation in "local geodesic" or
! "global spherical" space.
! ub Longitudinal velocity components in spherical coordinates.
! uxl x-derivatives of u at the left (west) edge of given interval
! vxl x-derivatives of v at the left (west) edge of given interval
! uxr x-derivatives of u at the right (east) edge of given interval
! vxr x-derivatives of v at the right (east) edge of given interval
! lam Longitude values for the extended grid.
! phib Latitude values for the extended grid.
! lbasiy Weights for Lagrange cubic interpolation on the unequally
! spaced latitude grid.
! lammp Longitude coordinates of the trajectory mid-points of the
! parcels that correspond to the global grid points contained
! in the latitude slice being forecast. On entry lammp
! is an initial guess.
! phimp Latitude coordinates of the trajectory mid-points of the
! parcels that correspond to the global grid points contained
! in the latitude slice being forecast. On entry phimp
! is an initial guess.
! lamdp Longitude coordinates of the departure points that correspond
! to the global grid points contained in the latitude slice
! being forecast. lamdp is constrained so that
! 0.0 .le. lamdp(i) .lt. 2*pi .
! phidp Latitude coordinates of the departure points that correspond
! to the global grid points contained in the latitude slice
! being forecast. If phidp is computed outside the latitudinal
! domain of the extended grid, then an abort will be called by
! subroutine "trjgl".
! idp Longitude index of departure points. This index points into
! the extended arrays, e.g.,
! lam (idp(i,k)) .le. lamdp(i,k) .lt. lam (idp(i,k)+1).
! jdp Latitude index of departure points. This index points into
! the extended arrays, e.g.,
! phib(jdp(i,k)) .le. phidp(i,k) .lt. phib(jdp(i,k)+1).
!-----------------------------------------------------------------------
!------------------------ local variables ------------------------------
integer iter ! index
integer i, j, k ! indices
integer imax, imin, kmin, kmax ! indices
real(r8) finc ! time step factor
real(r8) dttmp ! time step (seconds)
real(r8) dlam(platd) ! increment of grid in x-direction
real(r8) phicen ! latitude coord of current lat slice
real(r8) cphic ! cos(phicen)
real(r8) sphic ! sin(phicen)
real(r8) upr (plon,plev) ! u in local geodesic coords
real(r8) vpr (plon,plev) ! v in local geodesic coords
real(r8) lampr(plon,plev) ! relative long coord of dep pt
real(r8) phipr(plon,plev) ! relative lat coord of dep pt
real(r8) uvmp (plon,plev,2) ! u/v (spherical) interpltd to dep pt
real(r8) fint (plon,plev,ppdy,2) ! u/v x-interpolants
real(r8) phidpmax
real(r8) phidpmin
real(r8) phimpmax
real(r8) phimpmin
!-----------------------------------------------------------------------
!
do j=1,platd
dlam(j) = lam(nxpt+2,j) - lam(nxpt+1,j)
end do
phicen = phib(jcen)
cphic = cos( phicen )
sphic = sin( phicen )
!
! Convert latitude coordinates of trajectory midpoints from spherical
! to local geodesic basis.
!
if( locgeo ) call s2gphi(lam(i1,jcen) ,cphic ,sphic ,lammp ,phimp, &
phipr ,nlon )
!
! Loop over departure point iterates.
!
do 30 iter = 1,iterdp
!
! Compute midpoint indicies.
!
call bandij(dlam ,phib ,lammp ,phimp ,idp , &
jdp ,nlon )
!
! Hermite cubic interpolation to the x-coordinate of each
! departure point at each y-coordinate required to compute the
! y-interpolants.
!
call herxin(1 ,1 ,ub ,uxl ,uxr , &
lam ,lammp ,idp ,jdp ,fint(1,1,1,1), &
nlon ,nlonex )
call herxin(1 ,1 ,vb ,vxl ,vxr , &
lam ,lammp ,idp ,jdp ,fint(1,1,1,2), &
nlon ,nlonex )
call lagyin(2 ,fint ,lbasiy ,phimp ,jdp , &
jcen ,uvmp ,nlon )
#if ( defined SCAM )
!
! turn off u/v for scam
!
uvmp(:,:,1:2) = 0
#else
!
! Put u/v on unit sphere
!
do k = 1,plev
do i = 1,nlon
uvmp(i,k,1) = uvmp(i,k,1)*ra
uvmp(i,k,2) = uvmp(i,k,2)*ra
end do
end do
#endif
!
! For local geodesic:
!
! a) Convert velocity coordinates at trajectory midpoints from
! spherical coordinates to local geodesic coordinates,
! b) Estimate midpoint parcel trajectory,
! c) Convert back to spherical coordinates
!
! Else, for global spherical
!
! Estimate midpoint trajectory with no conversions
!
if ( locgeo ) then
call s2gvel(uvmp(1,1,1),uvmp(1,1,2) ,lam(i1,jcen) ,cphic ,sphic , &
lammp ,phimp ,upr ,vpr ,nlon )
call trajmp(dt ,upr ,vpr ,phipr ,lampr , &
nlon )
dttmp = 0.5*dt
call g2spos(dttmp ,lam(i1,jcen) ,phib ,phicen ,cphic , &
sphic ,upr ,vpr ,lampr ,phipr , &
lammp ,phimp ,nlon )
else
call trjmps(dt ,uvmp(1,1,1) ,uvmp(1,1,2), phimp ,lampr , &
phipr ,nlon )
finc = 1.
call trjgl (finc ,phicen ,lam(i1,jcen) ,lampr ,phipr , &
lammp ,phimp ,nlon )
end if
!
! Test that the latitudinal extent of trajectory is NOT over the poles
! Distributed memory case: check that the latitudinal extent of the
! trajectory is not more than "jintmx" gridpoints away.
!
phimpmax = -1.e36
phimpmin = 1.e36
do k=1,plev
do i=1,nlon
if (phimp(i,k)>phimpmax) then
phimpmax = phimp(i,k)
imax = i
kmax = k
end if
if (phimp(i,k)<phimpmin) then
phimpmin = phimp(i,k)
imin = i
kmin = k
end if
end do
end do
#if ( defined SPMD )
if ( phimp(imax,kmax) >= phib(endlatex-nxpt) ) then
#else
if ( phimp(imax,kmax) >= phib(j1+plat) ) then
#endif
write(6,*)'SPHDEP: ****** MODEL IS BLOWING UP *********'
write(6,9000) imax,kmax,jgc
call endrun
#if ( defined SPMD )
else if( phimp(imin,kmin) < phib(beglatex+nxpt) ) then
#else
else if( phimp(imin,kmin) < phib(j1-1) ) then
#endif
write(6,*)'SPHDEP: ****** MODEL IS BLOWING UP *********'
write(6,9000) imin,kmin,jgc
call endrun
end if
30 continue ! End of iter=1,iterdp loop
!
! Compute departure points in geodesic coordinates, and convert back
! to spherical coordinates.
!
! Else, compute departure points directly in spherical coordinates
!
if (locgeo) then
do k = 1,plev
do i = 1,nlon
lampr(i,k) = 2.*lampr(i,k)
phipr(i,k) = 2.*phipr(i,k)
end do
end do
dttmp = dt
call g2spos(dttmp ,lam(i1,jcen) ,phib ,phicen ,cphic , &
sphic ,upr ,vpr ,lampr ,phipr , &
lamdp ,phidp ,nlon )
else
finc = 2.
call trjgl (finc ,phicen ,lam(i1,jcen) ,lampr ,phipr , &
lamdp ,phidp ,nlon )
end if
!
! Test that the latitudinal extent of trajectory is NOT over the poles
! Distributed memory case: check that the latitudinal extent of the
! trajectory is not more than "jintmx" gridpoints away.
!
phidpmax = -1.e36
phidpmin = 1.e36
do k=1,plev
do i=1,nlon
if (phidp(i,k)>phidpmax) then
phidpmax = phidp(i,k)
imax = i
kmax = k
end if
if (phidp(i,k)<phidpmin) then
phidpmin = phidp(i,k)
imin = i
kmin = k
end if
end do
end do
#if ( defined SPMD )
if ( phidp(imax,kmax) >= phib(endlatex-nxpt) ) then
#else
if ( phidp(imax,kmax) >= phib(j1+plat) ) then
#endif
write(6,*)'SPHDEP: ****** MODEL IS BLOWING UP *********'
write(6,9000) imax,kmax,jgc
call endrun
#if ( defined SPMD )
else if( phidp(imin,kmin) < phib(beglatex+nxpt) ) then
#else
else if( phidp(imin,kmin) < phib(j1-1) ) then
#endif
write(6,*)'SPHDEP: ****** MODEL IS BLOWING UP *********'
write(6,9000) imin,kmin,jgc
call endrun
end if
!
! Compute departure point indicies.
!
call bandij(dlam ,phib ,lamdp ,phidp ,idp , &
jdp ,nlon )
9000 format(//'Parcel associated with longitude ',i5,', level ',i5, &
' and latitude ',i5,' is outside the model domain.')
return
end subroutine sphdep
!============================================================================================
subroutine g2spos(dttmp ,lam ,phib ,phi ,cosphi , & 2,2
sinphi ,upr ,vpr ,lamgc ,phigc , &
lamsc ,phisc ,nlon )
!-----------------------------------------------------------------------
!
! Purpose:
! Transform position coordinates for a set of points, each of which is
! associated with a grid point in a global latitude slice, from local
! geodesic to spherical coordinates.
!
! Method:
!
! Author: J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: sphdep.F90,v 1.1.44.5.4.1 2004/05/20 18:36:06 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
implicit none
!------------------------------Arguments--------------------------------
real(r8), intent(in) :: dttmp ! time step
real(r8), intent(in) :: lam(plond1) ! model longitude coordinates
real(r8), intent(in) :: phib(platd) ! extended grid latitude coordinates
real(r8), intent(in) :: phi ! current latitude coordinate (radians)
real(r8), intent(in) :: cosphi ! cos of current latitude
real(r8), intent(in) :: sinphi ! sin of current latitude
real(r8), intent(in) :: upr (plon,plev) ! u-wind in geodesic coord
real(r8), intent(in) :: vpr (plon,plev) ! v-wind in geodesic coord
real(r8), intent(in) :: lamgc(plon,plev) ! geodesic long coord. of dep. point
real(r8), intent(in) :: phigc(plon,plev) ! geodesic lat coord. of dep. point
integer , intent(in) :: nlon ! longitude dimension
real(r8), intent(out):: lamsc(plon,plev) ! spherical long coord. of dep. point
real(r8), intent(out):: phisc(plon,plev) ! spherical lat coord. of dep. point
!
!
! dttmp Time step over which midpoint/endpoint trajectory is
! calculated (seconds).
! lam Longitude coordinates of the global grid points in spherical
! system. The grid points in the global array are the reference
! points for the local geodesic systems.
! phib Latitude values for the extended grid.
! phi Latitude coordinate (in the global grid) of the current
! latitude slice.
! cosphi cos( phi )
! sinphi sin( phi )
! upr zonal velocity at departure point in local geodesic coord
! vpr Meridional velocity at departure point in local geodesic coord
! lamgc Longitude coordinate of points in geodesic coordinates.
! phigc Latitude coordinate of points in geodesic coordinates.
! lamsc Longitude coordinate of points in spherical coordinates.
! phisc Latitude coordinate of points in spherical coordinates.
!-----------------------------------------------------------------------
!---------------------------Local variables-----------------------------
integer i,ii,k ! indices
integer nval(plev) ! number of values returned from whenfgt
integer indx(plon,plev) ! index holder
real(r8) pi ! 4.*atan(1.)
real(r8) twopi ! 2.*pi
real(r8) pi2 ! pi/2
real(r8) sgnphi ! holds sign of phi
real(r8) sphigc ! sin(phigc)
real(r8) cphigc ! cos(phigc)
real(r8) clamgc ! cos(lamgc)
real(r8) slam2 ! sin(lamgc)**2
real(r8) phipi2 ! tmp variable
real(r8) slamgc(plon,plev) ! sin(lamgc)
real(r8) dlam(plon,plev) ! zonal extent of trajectory
real(r8) coeff ! tmp variable
real(r8) distmx ! max distance
real(r8) dist(plon,plev) ! approx. distance traveled along traj.
real(r8) fac ! 1. - 10*eps, eps from mach. precision
integer s_nval
!-----------------------------------------------------------------------
!
fac = 1. - 10.*epsilon (fac)
pi = 4.*atan(1.)
twopi = pi*2.
pi2 = pi/2.
coeff = (1.1*dttmp)**2
distmx = (sign(pi2,phi) - phi)**2/coeff
sgnphi = sign( 1._r8, phi )
do k=1,plev
do i=1,nlon
sphigc = sin( phigc(i,k) )
cphigc = cos( phigc(i,k) )
slamgc(i,k) = sin( lamgc(i,k) )
clamgc = cos( lamgc(i,k) )
phisc(i,k) = asin((sphigc*cosphi + cphigc*sinphi*clamgc)*fac)
if ( abs(phisc(i,k)) .ge. phib(j1+plat)*fac ) then
phisc(i,k) = sign( phib(j1+plat),phisc(i,k) )*fac
end if
dlam(i,k) = asin((slamgc(i,k)*cphigc/cos(phisc(i,k)))*fac)
!
! Compute estimated trajectory distance based upon winds alone
!
dist(i,k) = upr(i,k)**2 + vpr(i,k)**2
end do
!
! Determine which trajectories may have crossed over pole
!
s_nval = 0
do i=1,nlon
if (dist(i,k) > distmx) then
s_nval = s_nval + 1
indx(s_nval,k) = i
end if
end do
nval(k) = s_nval
end do
!
! Check that proper branch of arcsine is used for calculation of
! dlam for those trajectories which may have crossed over pole.
!
do k=1,plev
!cdir nodep
do ii=1,nval(k)
i = indx(ii,k)
slam2 = slamgc(i,k)**2
phipi2 = asin((sqrt((slam2 - 1.)/(slam2 - 1./cosphi**2)))*fac)
if (sgnphi*phigc(i,k) > phipi2) then
dlam(i,k) = sign(pi,lamgc(i,k)) - dlam(i,k)
end if
end do
do i=1,nlon
lamsc(i,k) = lam(i) + dlam(i,k)
!
! Restrict longitude to be in the range [0, twopi).
!
if( lamsc(i,k) >= twopi ) lamsc(i,k) = lamsc(i,k) - twopi
if( lamsc(i,k) < 0.0 ) lamsc(i,k) = lamsc(i,k) + twopi
end do
end do
return
end subroutine g2spos
!============================================================================================
subroutine s2gphi(lam ,cosphi ,sinphi ,lamsc ,phisc , & 1,2
phigc ,nlon )
!-----------------------------------------------------------------------
!
! Purpose:
! Calculate transformed local geodesic latitude coordinates for a set
! of points, each of which is associated with a grid point in a global
! latitude slice. Transformation is spherical to local geodesic.
! (Williamson and Rasch, 1991)
!
! Method:
!
! Author: J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: sphdep.F90,v 1.1.44.5.4.1 2004/05/20 18:36:06 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
implicit none
!------------------------------Arguments--------------------------------
real(r8), intent(in) :: lam(plond1) ! long coordinates of model grid
real(r8), intent(in) :: cosphi ! cos(latitude)
real(r8), intent(in) :: sinphi ! sin(latitude)
real(r8), intent(in) :: lamsc(plon,plev) ! spher. long coords of dep points
real(r8), intent(in) :: phisc(plon,plev) ! spher. lat coords of dep points
integer , intent(in) :: nlon ! longitude dimension
real(r8), intent(out) :: phigc(plon,plev) ! loc geod. lat coords of dep points
!
! lam longitude coordinates of the global grid points in spherical
! system. The grid points in the global array are the reference
! points for the local geodesic systems.
! cosphi cosine of the latitude of the global latitude slice.
! sinphi sine of the latitude of the global latitude slice.
! lamsc longitude coordinate of dep. points in spherical coordinates.
! phisc latitude coordinate of dep. points in spherical coordinates.
! phigc latitude coordinate of dep. points in local geodesic coords.
!-----------------------------------------------------------------------
!---------------------------Local variables-----------------------------
integer i,k ! longitude, level indices
real(r8) sphisc ! |
real(r8) cphisc ! | -- temporary variables
real(r8) clamsc ! |
!-----------------------------------------------------------------------
!
do k = 1,plev
do i = 1,nlon
sphisc = sin( phisc(i,k) )
cphisc = cos( phisc(i,k) )
clamsc = cos( lam(i) - lamsc(i,k) )
phigc(i,k) = asin( sphisc*cosphi - cphisc*sinphi*clamsc )
end do
end do
return
end subroutine s2gphi
!============================================================================================
subroutine s2gvel(udp ,vdp ,lam ,cosphi ,sinphi , & 1,2
lamdp ,phidp ,upr ,vpr ,nlon )
!-----------------------------------------------------------------------
!
! Purpose:
! Transform velocity components at departure points associated with a
! single latitude slice from spherical coordinates to local geodesic
! coordinates. (Williamson and Rasch, 1991)
!
! Method:
!
! Author: J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: sphdep.F90,v 1.1.44.5.4.1 2004/05/20 18:36:06 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
implicit none
!------------------------------Arguments--------------------------------
integer , intent(in) :: nlon ! longitude dimension
real(r8), intent(in) :: udp(plon,plev) ! u in spherical coords.
real(r8), intent(in) :: vdp(plon,plev) ! v in spherical coords.
real(r8), intent(in) :: lam(plond1) ! x-coordinates of model grid
real(r8), intent(in) :: cosphi ! cos(latitude)
real(r8), intent(in) :: sinphi ! sin(latitude)
real(r8), intent(in) :: lamdp(plon,plev) ! spherical longitude coord of dep pt.
real(r8), intent(in) :: phidp(plon,plev) ! spherical latitude coord of dep pt.
real(r8), intent(out) :: upr(plon,plev) ! u in local geodesic coords.
real(r8), intent(out) :: vpr(plon,plev) ! v in local geodesic coords.
!
! udp u-component of departure point velocity in spherical coords.
! vdp v-component of departure point velocity in spherical coords.
! lam Longitude of arrival point position (model grid point) in spherical coordinates.
! cosphi Cos of latitude of arrival point positions (model grid pt).
! sinphi Sin of latitude of arrival point positions (model grid pt).
! lamdp Longitude of departure point position in spherical coordinates.
! phidp Latitude of departure point position in spherical coordinates.
! upr u-component of departure point velocity in geodesic coords.
! vpr v-component of departure point velocity in geodesic coords.
!-----------------------------------------------------------------------
!---------------------------Local variables-----------------------------
integer i,k ! longitude, level indices
real(r8) cdlam ! |
real(r8) clamp ! |
real(r8) cphid ! |
real(r8) cphip ! |
real(r8) dlam ! | -- temporary variables
real(r8) sdlam ! |
real(r8) slamp ! |
real(r8) sphid ! |
real(r8) sphip ! |
!-----------------------------------------------------------------------
!
do k = 1,plev
do i = 1,nlon
dlam = lam(i) - lamdp(i,k)
sdlam = sin( dlam )
cdlam = cos( dlam )
sphid = sin( phidp(i,k) )
cphid = cos( phidp(i,k) )
sphip = sphid*cosphi - cphid*sinphi*cdlam
cphip = cos( asin( sphip ) )
slamp = -sdlam*cphid/cphip
clamp = cos( asin( slamp ) )
vpr(i,k) = (vdp(i,k)*(cphid*cosphi + sphid*sinphi*cdlam) - &
udp(i,k)*sinphi*sdlam)/cphip
upr(i,k) = (udp(i,k)*cdlam + vdp(i,k)*sphid*sdlam + &
vpr(i,k)*slamp*sphip)/clamp
end do
end do
return
end subroutine s2gvel
!============================================================================================
subroutine trajmp(dt ,upr ,vpr ,phipr ,lampr , & 1,2
nlon )
!-----------------------------------------------------------------------
!
! Purpose:
! Estimate mid-point of parcel trajectory (geodesic coordinates) based
! upon horizontal wind field.
!
! Method:
!
! Author: J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: sphdep.F90,v 1.1.44.5.4.1 2004/05/20 18:36:06 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
implicit none
!------------------------------Arguments--------------------------------
integer , intent(in) :: nlon ! longitude dimension
real(r8), intent(in) :: dt ! time step (seconds)
real(r8), intent(in) :: upr(plon,plev) ! u-component of wind in local geodesic
real(r8), intent(in) :: vpr(plon,plev) ! v-component of wind in local geodesic
real(r8), intent(inout) :: phipr(plon,plev) ! latitude coord of trajectory mid-point
real(r8), intent(out) :: lampr(plon,plev) ! longitude coord of traj. mid-point
!
! dt Time interval that corresponds to the parcel trajectory.
! upr u-coordinate of velocity corresponding to the most recent
! estimate of the trajectory mid-point (in geodesic system).
! vpr v-coordinate of velocity corresponding to the most recent
! estimate of the trajectory mid-point (in geodesic system).
! phipr Phi value at trajectory mid-point (geodesic coordinates).
! On entry this is the most recent estimate.
! lampr Lambda value at trajectory mid-point (geodesic coordinates).
!-----------------------------------------------------------------------
!---------------------------Local variables-----------------------------
integer i,k ! index
!-----------------------------------------------------------------------
!
do k=1,plev
do i = 1,nlon
lampr(i,k) = -.5*dt* upr(i,k) / cos( phipr(i,k) )
phipr(i,k) = -.5*dt* vpr(i,k)
end do
end do
return
end subroutine trajmp
!============================================================================================
subroutine trjgl(finc ,phicen ,lam ,lampr ,phipr , & 2,2
lamp ,phip ,nlon )
!-----------------------------------------------------------------------
!
! Purpose:
! Map relative trajectory mid/departure point coordinates to global
! latitude/longitude coordinates and test limits
!
! Method:
!
! Author: J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: sphdep.F90,v 1.1.44.5.4.1 2004/05/20 18:36:06 mvr Exp $
! $Author: mvr $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
implicit none
!------------------------------Arguments--------------------------------
integer , intent(in) :: nlon ! longitude dimension
real(r8), intent(in) :: finc ! number of time increments
real(r8), intent(in) :: phicen ! current latitude value in extnded grid
real(r8), intent(in) :: lam(plond1) ! longitude values for the extended grid
real(r8), intent(in) :: lampr(plon,plev) ! relative x-coordinate of departure pt.
real(r8), intent(in) :: phipr(plon,plev) ! relative y-coordinate of departure pt.
real(r8), intent(out) :: lamp (plon,plev) ! long coords of traj midpoints
real(r8), intent(out) :: phip (plon,plev) ! lat coords of traj midpoints
!
! finc Time step factor (1. for midpoint, 2. for dep. point)
! phicen Latitude value for current latitude being forecast.
! lam Longitude values for the extended grid.
! lampr Longitude coordinates (relative to the arrival point) of the
! trajectory mid-points of the parcels that correspond to the
! global grid points contained in the latitude slice being forecast.
! phipr Latitude coordinates (relative to the arrival point) of the
! trajectory mid-points of the parcels that correspond to the
! global grid points contained in the latitude slice being forecast.
! lamp Longitude coordinates of the trajectory mid-points of the
! parcels that correspond to the global grid points contained
! in the latitude slice being forecast.
! phip Latitude coordinates of the trajectory mid-points of the
! parcels that correspond to the global grid points contained
! in the latitude slice being forecast.
!-----------------------------------------------------------------------
!--------------------------Local variables------------------------------
integer i ! longitude index
integer k ! level index
real(r8) pi ! 3.14.......
real(r8) twopi ! 2*pi
!-----------------------------------------------------------------------
!
pi = 4.*atan(1.)
twopi = pi*2.
do k = 1,plev
do i = 1,nlon
lamp(i,k) = lam(i) + finc*lampr(i,k)
phip(i,k) = phicen + finc*phipr(i,k)
if(lamp(i,k) >= twopi) lamp(i,k) = lamp(i,k) - twopi
if(lamp(i,k) < 0.0) lamp(i,k) = lamp(i,k) + twopi
end do
end do
return
end subroutine trjgl