sphdep.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

!  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
!-----------------------------------------------------------------------
!
  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
!
     nval(k) = 0
     do i=1,nlon
        if (dist(i,k) > distmx) then
           nval(k) = nval(k) + 1
           indx(nval(k),k) = i
        end if
     end do
  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$ IVDEP
     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   , &
                  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.1 2002/05/13 17:57:36 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  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  , &
                  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.1 2002/05/13 17:57:36 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  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   , &
                  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.1 2002/05/13 17:57:36 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  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   , &
                 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.1 2002/05/13 17:57:36 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

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