sphdep.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine sphdep(jcen    ,jgc     ,dt      ,ra      ,iterdp  , &
                  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.1 2002/05/13 17:57:36 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid
#if (!defined CRAY)
  use srchutil, only: ismax,ismin
#endif
  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    )
!
! 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
!
! 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  , &
                  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.1 2002/05/13 17:57:36 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  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.