#include <misc.h>
#include <params.h>


subroutine sltb1(pmap    ,jcen    ,jgc     ,dt      ,ra      , & 1,11
                 iterdp  ,uxl     ,uxr     ,vxl     ,vxr     , &
                 wb      ,fxl     ,fxr     ,lam     ,phib    , &
                 dphib   ,sig     ,sigh    ,dsig    ,dsigh   , &
                 lbasdy  ,lbasdz  ,lbassd  ,lbasiy  ,kdpmpf  , &
                 kdpmph  ,lammp   ,phimp   ,sigmp   ,fbout   , &
                 u3      ,v3      ,qminus  ,n3m1    , &
                 nlon    ,nlonex  )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Drive the slt algorithm on a given latitude slice in the extended
! data arrays using information from the entire latitudinal extent
! of the arrays.
! 
! Method: 
! Compute departure points and corresponding indices.
! Poleward of latitude phigs (radians), perform the computation in
! local geodesic coordinates.
! Equatorward of latitude phigs, perform the computation in global
! spherical coordinates
! 
! Author: J. Olson
! 
!-----------------------------------------------------------------------
!
! $Id: 
!
!-----------------------------------------------------------------------

  use shr_kind_mod, only: r8 => shr_kind_r8
  use pmgrid
  use prognostics, only: ptimelevels
  use constituents,only: pcnst

  implicit none

#include <parslt.h>

!------------------------------Parameters-------------------------------
  real(r8), parameter :: phigs = 1.221730 ! cut-off latitude: about 70 degrees
!-----------------------------------------------------------------------

!------------------------------Arguments--------------------------------
  integer , intent(in) :: nlon                ! longitude dimension
  integer , intent(in) :: nlonex(platd)       ! extended longitude dimension
  integer , intent(in) :: pmap                ! artificial vert grid dim.
  integer , intent(in) :: jcen                ! index of lat slice(extend)
  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              ! iteration count
  real(r8), intent(in) :: uxl(plond,plev,beglatex:endlatex) ! left  x-deriv of ub
  real(r8), intent(in) :: uxr(plond,plev,beglatex:endlatex) ! right x-deriv of ub
  real(r8), intent(in) :: vxl(plond,plev,beglatex:endlatex) ! left  x-deriv of vb
  real(r8), intent(in) :: vxr(plond,plev,beglatex:endlatex) ! right x-deriv of vb
  real(r8), intent(in) :: wb(plon,plevp)                    ! eta-dot
  real(r8), intent(in) :: fxl(plond,plev,  pcnst,beglatex:endlatex) ! left  fb x-deriv
  real(r8), intent(in) :: fxr(plond,plev,  pcnst,beglatex:endlatex) ! right fb 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) :: dphib(platd)        ! increment between lats.
  real(r8), intent(in) :: sig  (plev)         ! vertical full levels
  real(r8), intent(in) :: sigh (plevp)        ! vertical half levels
  real(r8), intent(in) :: dsig (plev)         ! inc. between full levs
  real(r8), intent(in) :: dsigh(plevp)        ! inc. between half levs
  real(r8), intent(in) :: lbasdy(4,2,platd)   ! lat deriv weights
  real(r8), intent(in) :: lbasdz(4,2,plev)    ! vert full level deriv wts
  real(r8), intent(in) :: lbassd(4,2,plevp)   ! vert half level deriv wts
  real(r8), intent(in) :: lbasiy(4,2,platd)   ! lat interp wts(lagrng)
  integer , intent(in) :: kdpmpf(pmap)        ! artificial vert grid index
  integer , intent(in) :: kdpmph(pmap)        ! artificial vert grid index
  real(r8), intent(inout) ::  u3(plond, plev, beglatex:endlatex,ptimelevels)     ! u wind vel
  real(r8), intent(inout) ::  v3(plond, plev, beglatex:endlatex,ptimelevels)     ! v wind vel
  real(r8), intent(inout) ::  qminus(plond, plev, pcnst, beglatex:endlatex) !moist
  integer , intent(inout) ::  n3m1               ! time indicies
  real(r8), intent(inout) ::  lammp(plon,plev)       ! long coord of mid-point
  real(r8), intent(inout) ::  phimp(plon,plev)       ! lat  coord of mid-point
  real(r8), intent(inout) ::  sigmp(plon,plev)       ! vert coord of mid-point
  real(r8), intent(out) :: fbout(plond,plev,pcnst)   ! advected constituents
!
!  pmap    Dimension of kdpmpX arrays
!  jcen    Latitude index in extended grid corresponding to lat slice
!          being forecasted.
!  jgc     Latitude index in model    grid corresponding to lat slice 
!          being forecasted.
!  dt      Time interval that parameterizes the parcel trajectory.
!  ra      Reciprocal of radius of earth.
!  iterdp  Number of iterations used for departure point calculation.
!  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
!  wb      z-velocity component (eta-dot).
!  fxl     x-derivatives at the left  edge of each interval containing 
!          the departure point.
!  fxr     x-derivatives at the right edge of each interval containing 
!          the departure point.
!  lam     Longitude values for the extended grid.
!  phib    Latitude  values for the extended grid.
!  dphib   Interval between latitudes in the extended grid.
!  sig     Hybrid eta values at the "full-index" levels.
!  sigh    Half-index eta-levels including sigh(i,1) = eta(1/2) = 0.0
!          and sigh(i,plev+1) = eta(plev+1/2) = 1.  Note that in general
!          sigh(i,k) .lt. sig(i,k)  where sig(i,k) is the hybrid value
!          at the k_th full-index level.
!  dsig    Interval lengths in full-index hybrid level grid.
!  dsigh   Interval lengths in half-index hybrid level grid.
!  lbasdy  Weights for Lagrange cubic derivative estimates on the
!          unequally spaced latitude grid.
!  lbasdz  Weights for Lagrange cubic derivative estimates on the
!          unequally spaced vertical grid (full levels).
!  lbassd  Weights for Lagrange cubic derivative estimates on the
!          unequally spaced vertical grid (half levels).
!  lbasiy  Weights for Lagrange cubic interpolation on the unequally
!          spaced latitude grid.
!  kdpmpf  indices of artificial grid mapped into the full level grid
!  kdpmph  indices of artificial grid mapped into the half level 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 forecasted.  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 forecasted.  On entry phimp
!          is an initial guess.
!  sigmp   Hybrid value at the trajectory midpoint for each gridpoint
!          in a vertical slice from the global grid.  On entry sigmp is
!          an initial guess.
!  fbout   Extended array only one latitude of which, however, is filled
!          with forecasted (transported) values.  This routine must be
!          called multiple times to fill the entire array.  This is
!          done to facilitate multi-tasking.
!-----------------------------------------------------------------------

!---------------------------Local variables-----------------------------
  integer m                            ! constituent index
  integer idp(plon,plev,4)             ! zonal      dep point index
  integer jdp(plon,plev)               ! meridional dep point index
  integer kdp(plon,plev)               ! vertical   dep point index
  real(r8) fhr(plon,plev,pcnst)        ! horizontal interpolants
  real(r8) lamdp(plon,plev)            ! zonal      departure pt. coord.
  real(r8) phidp(plon,plev)            ! meridional departure pt. coord.
  real(r8) sigdp(plon,plev)            ! vertical   departure pt. coord.
  real(r8) fhst(plon,plev,pcnst)       ! derivative at top of interval
  real(r8) fhsb(plon,plev,pcnst)       ! derivative at bot of interval
  real(r8) wst(plon,plevp)             ! w derivative at top of interval
  real(r8) wsb(plon,plevp)             ! w derivative at bot of interval
  real(r8) fint(plon,plev,ppdy,pcnst)  ! work space
  real(r8) fyb(plon,plev,pcnst)        ! work space
  real(r8) fyt(plon,plev,pcnst)        ! work space
  real(r8) fbout1(plond,plev,pcnst)    ! work space to please lf95 compiler
  logical locgeo                       ! flag indicating coordinate sys
#if ( defined SCAM )
   integer k,i
#endif!-----------------------------------------------------------------------
#if ( !defined SCAM )
!
! Horizontal interpolation
!
  locgeo = abs(phib(jcen))>=phigs
!
  call sphdep(jcen    ,jgc     ,dt      ,ra      ,iterdp  ,                     &
              locgeo  ,u3(1,1,beglatex,n3m1)       ,uxl     ,uxr     ,lam     ,   &
              phib    ,lbasiy  ,lammp   ,phimp   ,lamdp   ,                     &
              phidp   ,idp     ,jdp     ,v3(1,1,beglatex,n3m1),                   &
              vxl     ,vxr     ,nlon    ,nlonex  )
!
! Interpolate scalar fields to the departure points.
!
  call hrintp(pcnst   ,pcnst   ,qminus(1,1,1,beglatex), fxl     ,fxr     , &
              lam     ,phib    ,dphib   ,lbasdy  ,lamdp   ,                    &
              phidp   ,idp     ,jdp     ,jcen    ,plimdr  ,                    &
              fint    ,fyb     ,fyt     ,fhr     ,nlon    ,                    &   
              nlonex  )
#else
!
! fill in fhr in leiu of horizontal interpolation
!
  do m = 1,pcnst
     do k = 1,plev
        do i = 1,nlon
           fhr(i,k,m) = qminus(i1+i-1,k,m,jcen)
        end do
     end do
  end do
#endif
!
! Vertical interpolation.
! Compute vertical derivatives of vertical wind
!
  call cubzdr(nlon    ,plevp   ,wb      ,lbassd  ,wst     , &
              wsb     )
!
! Compute departure points and corresponding indices.
!
  call vrtdep(pmap    ,dt      ,iterdp  ,wb      ,wst     , &
              wsb     ,sig     ,sigh    ,dsigh   ,kdpmpf  , &
              kdpmph  ,sigmp   ,sigdp   ,kdp     ,nlon    )
!
! Vertical derivatives of scalar fields.
! Loop over constituents.
!
!CSD$ PARALLEL DO 
  do m = 1,pcnst
     call cubzdr(nlon    ,plev    ,fhr(1,1,m), lbasdz  ,fhst(1,1,m), &
                 fhsb(1,1,m) )
  end do
!CSD$ END PARALLEL DO
  if( plimdr )then
     call limdz(fhr     ,dsig    ,fhst    ,fhsb    ,nlon    )
  end if
!
! Vertical interpolation of scalar fields.
!
  call herzin(plev    ,pcnst   ,fhr     ,fhst    ,fhsb    , &
              sig     ,dsig    ,sigdp   ,kdp     ,fbout1  , &
              nlon    )

  fbout(i1:nlon+i1-1,:,:) = fbout1(1:nlon,:,:)

  return
end subroutine sltb1

!============================================================================================


subroutine vrtdep(pmap    ,dt      ,iterdp  ,wb      ,wst     , & 1,7
                  wsb     ,sig     ,sigh    ,dsigh   ,kdpmpf  , & 
                  kdpmph  ,sigmp   ,sigdp   ,kdp     ,nlon    )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Compute vertical departure point and departure point index.
! 
! Method: 
! 
! Author: J. Olson
! 
!-----------------------------------------------------------------------

  use shr_kind_mod, only: r8 => shr_kind_r8
  use pmgrid
  implicit none

!------------------------------Arguments--------------------------------
  integer , intent(in) :: nlon                ! longitude dimension
  integer , intent(in) :: pmap                ! dimension of artificial vert grid
  real(r8), intent(in) :: dt                  ! time step (seconds)
  integer , intent(in) :: iterdp              ! number of iterations
  real(r8), intent(in) :: wb (plon,plevp)     ! vertical velocity
  real(r8), intent(in) :: wst(plon,plevp)     ! z-derivative of wb at top of interval
  real(r8), intent(in) :: wsb(plon,plevp)     ! z-derivative of wb at bot of interval
  real(r8), intent(in) :: sig  (plev )        ! sigma values of model full levels
  real(r8), intent(in) :: sigh (plevp)        ! sigma values of model half levels
  real(r8), intent(in) :: dsigh(plevp)        ! increment between half levels
  integer , intent(in) :: kdpmpf(pmap)        ! artificial grid indices
  integer , intent(in) :: kdpmph(pmap)        ! artificial grid indices
  real(r8), intent(inout) :: sigmp(plon,plev) ! vert coords of traj mid-points
  real(r8), intent(out) :: sigdp(plon,plev)   ! vert coords of traj departure points
  integer , intent(out) :: kdp(plon,plev)     ! vertical departure point indices
!
!  pmap    Dimension of kdpmap arrays
!  dt      Time interval that parameterizes the parcel trajectory.
!  iterdp  Number of iterations used for departure point calculation.
!  wb      Vertical velocity component (sigma dot).
!  wst     z-derivs at the top edge of each interval contained in wb
!  wsb     z-derivs at the bot edge of each interval contained in wb
!  sig     Sigma values at the full-index levels.
!  sigh    Half-index sigma levels including sigh(1) = sigma(1/2) = 0.0
!          sigh(plev+1) = sigma(plev+1/2) = 1.0 .  Note that in general
!          sigh(k) .lt. sig(k)  where sig(k) is the sigma value at the
!          k_th full-index level.
!  dsigh   Increment in half-index sigma levels.
!  kdpmpf  Array of indices of the model full levels which are mapped
!          into an artificial evenly spaced vertical grid.  Used to aid
!          in search for vertical position of departure point 
!  kdpmph  Array of indices of the model half levels which are mapped
!          into an artificial evenly spaced vertical grid.  Used to aid
!          in search for vertical position of departure point 
!  sigmp   Sigma value at the trajectory midpoint for each gridpoint
!          in a vertical slice from the global grid.  On entry sigmp is
!          an initial guess.
!  sigdp   Sigma value at the trajectory endpoint for each gridpoint
!          in a vertical slice from the global grid.
!  kdp     Vertical index for each gridpoint.  This index points into a
!          vertical slice array whose vertical grid is given by sig.
!          E.g.,   sig(kdp(i,k)) .le. sigdp(i,k) .lt. sig(kdp(i,k)+1).
!-----------------------------------------------------------------------

!---------------------------Local variables-----------------------------
  integer i                 ! |
  integer iter              ! |-- indices
  integer k                 ! |
  real(r8) wmp(plond,plev)  ! vert vel. at midpoint
!-----------------------------------------------------------------------
!
! Loop over departure point iterates.
!
  do iter = 1,iterdp
!
! Compute midpoint indices in half-index sigma-level arrays (use kdp
! as temporary storage).
!
     call kdpfnd(plevp   ,pmap    ,sigh    ,sigmp   ,kdpmph  , &
                 kdp     ,nlon    )
!
! Interpolate sigma dot field to trajectory midpoints using Hermite
! cubic interpolant.
!
     call herzin(plevp   ,1       ,wb      ,wst     ,wsb     , &
                 sigh    ,dsigh   ,sigmp   ,kdp     ,wmp     , &
                 nlon    )
!
! Update estimate of trajectory midpoint.
!
     do k = 1,plev
        do i = 1,nlon
           sigmp(i,k) = sig(k) - .5*dt*wmp(i,k)
        end do
     end do
!
! Restrict vertical midpoints to be between the top and bottom half-
! index sigma levels.
!
     call vdplim(plevp   ,sigh    ,sigmp   ,nlon)
  end do
!
! Compute trajectory endpoints.
!
  do k = 1,plev
!DIR$ NOUNROLL
     do i = 1,nlon
        sigdp(i,k) = sig(k) - dt*wmp(i,k)
     end do
  end do
!
! Restrict vertical departure points to be between the top and bottom
! full-index sigma levels.
!
  call vdplim(plev    ,sig     ,sigdp   ,nlon)
!
! Vertical indices for trajectory endpoints that point into full-index
! sigma level arrays.
!
  call kdpfnd(plev    ,pmap    ,sig     ,sigdp   ,kdpmpf  , &
              kdp     ,nlon    )
!
  return
end subroutine vrtdep


!============================================================================================


subroutine vdplim(pkdim   ,sig     ,sigdp   ,nlon    ) 2,2

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Restrict vertical departure points to be between the top and bottom
! sigma levels of the "full-" or "half-" level grid
! 
! Method: 
! 
! Author: J. Olson
! 
!-----------------------------------------------------------------------
!
! $Id: sltb1.F90,v 1.6.4.4.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

!-----------------------------------------------------------------------
  integer , intent(in)    :: nlon               ! longitude dimension
  integer , intent(in)    :: pkdim              ! vertical dimension
  real(r8), intent(in)    :: sig(pkdim)         ! vertical coordinate of model grid
  real(r8), intent(inout) :: sigdp(plon,plev)   ! vertical coords. of departure points.
! pkdim   Vertical dimension of "sig"
! sig     Sigma values at the "full" or "half" model levels
! sigdp   Sigma value at the trajectory endpoint or midpoint for each
!         gridpoint in a vertical slice from the global grid.  This
!         routine restricts those departure points to within the
!         model's vertical grid.
!-----------------------------------------------------------------------

!---------------------------Local variables-----------------------------
  integer i,k                 ! index
!-----------------------------------------------------------------------
!
  do k=1,plev
     do i = 1,nlon
        if (sigdp(i,k) < sig(1)) then
           sigdp(i,k) = sig(1)
        end if
        if (sigdp(i,k) >= sig(pkdim)) then
           sigdp(i,k) = sig(pkdim)*(1. - 10.*epsilon(sigdp))
        end if
     end do
  end do

  return
end subroutine vdplim