#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