scanslt.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine scanslt (ztodt   ,lat     ,dtr     ,iter    ,pmap    , &
                    kdpmpf  ,kdpmph  ,lam     ,phi     ,dphi    , &
                    lbasdy  ,lbasdz  ,lbasiy  ,lbasiz  ,lbassi  , &
                    detam   ,detai   ,dlam    ,cwava   ,etamid  , &
                    etaint  ,grfu    ,grfv    ,grlps1  ,grlps2  , &
                    grt1    ,grt2    ,grq1    ,grq2    ,grfu1   , &
                    grfu2   ,grfv1   ,grfv2   ,ps      ,u3      , &
                    v3      ,t3      ,q3      ,lnpssld ,prhssld , &
                    tarrsld ,parrsld ,n3      ,n3m1    ,u3sld   , &
                    v3sld   ,etadot  ,nlon    )
!-----------------------------------------------------------------------
!
! Purpose:
! Interpolate terms for semi-lagrangian transport and SLD dynamics.
! One latitude slice only
!
! Author:  J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: scanslt.F90,v 1.14.2.1 2002/04/22 19:09:51 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid,      only: plon, plond, plev, plevp, plat, platd, beglat, endlat, beglatex, &
                         endlatex, plndlv, i1, j1
  use constituents,only: pcnst, pnats
  use comslt,      only: qfcst, gamma, hw1lat
  use rgrid,       only: nmmax
  use pspect,      only: pmmax
  use commap,      only: w, clat, t0
  use prognostics, only: ptimelevels
  use physconst,   only: cappa
  use dynconst,    only: ra

  implicit none

#include <comctl.h>
#include <comfft.h>
#include <comhyb.h>

!------------------------------Arguments--------------------------------
!
  real(r8), intent(in)   :: ztodt                 ! twice the time step unless nstep = 0
  integer , intent(in)   :: lat                   ! latitude index
  real(r8), intent(in)   :: dtr                   ! 1/dt
  integer , intent(in)   :: iter                  ! number of iterations for trajectory
  integer , intent(in)   :: pmap                  ! dimension of artificial array
!                                                 ! used to locate vertical interval
!                                                 ! in which departure point falls
  integer , intent(in)   :: kdpmpf  (pmap)        ! mapping from artificial array to
!                                                 ! model levels
  integer , intent(in)   :: kdpmph  (pmap)        ! mapping from artificial array to
!                                                 ! model interfaces
  real(r8), intent(in)   :: lam     (plond,platd) ! longitude coordinates of model grid
  real(r8), intent(in)   :: phi     (platd)       ! latitude  coordinates of model grid
  real(r8), intent(in)   :: dphi    (platd)       ! latitudinal grid increments
  real(r8), intent(in)   :: lbasdy  (4,2,platd)   ! basis functions for lat deriv est.
  real(r8), intent(in)   :: lbasdz  (4,2,plev)    ! basis functions for vert deriv est.
  real(r8), intent(in)   :: lbasiy  (4,2,platd)   ! basis functions for Lagrange interp
  real(r8), intent(in)   :: lbasiz  (4,2,plev)    ! Lagrange cubic interp wghts (vert)
  real(r8), intent(in)   :: lbassi  (4,2,plevp)   ! Lagrange cubic interp wghts (vert)
  real(r8), intent(in)   :: detam   (plev)        ! delta eta at levels
  real(r8), intent(in)   :: detai   (plevp)       ! delta eta at interfaces
  real(r8), intent(in)   :: dlam    (platd)       ! longitudinal grid increment
  real(r8), intent(in)   :: cwava   (plat)        ! weight for global water vapor int.
  real(r8), intent(in)   :: etamid  (plev)        ! eta at levels
  real(r8), intent(in)   :: etaint  (plevp)       ! eta at interfaces
  real(r8), intent(in)   :: grfu    (plond,plev,beglat:endlat) ! nonlinear term - u momentum eqn
  real(r8), intent(in)   :: grfv    (plond,plev,beglat:endlat) ! nonlinear term - v momentum eqn
#if ( defined PVP )
  real(r8), intent(out)  :: grlps1(2*pmmax     ,plat/2) ! ------------------------------
  real(r8), intent(out)  :: grlps2(2*pmmax     ,plat/2) ! |
  real(r8), intent(out)  :: grt1  (2*pmmax,plev,plat/2) ! |
  real(r8), intent(out)  :: grt2  (2*pmmax,plev,plat/2) ! |
  real(r8), intent(out)  :: grq1  (2*pmmax,plev,plat/2) ! |- see quad for definitions
  real(r8), intent(out)  :: grq2  (2*pmmax,plev,plat/2) ! |
  real(r8), intent(out)  :: grfu1 (2*pmmax,plev,plat/2) ! |
  real(r8), intent(out)  :: grfu2 (2*pmmax,plev,plat/2) ! |
  real(r8), intent(out)  :: grfv1 (2*pmmax,plev,plat/2) ! |
  real(r8), intent(out)  :: grfv2 (2*pmmax,plev,plat/2) ! ------------------------------
#else
  real(r8), intent(out)  :: grlps1(     2*pmmax,plat/2) ! ------------------------------
  real(r8), intent(out)  :: grlps2(     2*pmmax,plat/2) ! |
  real(r8), intent(out)  :: grt1  (plev,2*pmmax,plat/2) ! |
  real(r8), intent(out)  :: grt2  (plev,2*pmmax,plat/2) ! |
  real(r8), intent(out)  :: grq1  (plev,2*pmmax,plat/2) ! |
  real(r8), intent(out)  :: grq2  (plev,2*pmmax,plat/2) ! |- see quad for definitions
  real(r8), intent(out)  :: grfu1 (plev,2*pmmax,plat/2) ! |
  real(r8), intent(out)  :: grfu2 (plev,2*pmmax,plat/2) ! |
  real(r8), intent(out)  :: grfv1 (plev,2*pmmax,plat/2) ! |
  real(r8), intent(out)  :: grfv2 (plev,2*pmmax,plat/2) ! ------------------------------
#endif
  real(r8), intent(in)   :: ps      (plond,beglat:endlat,ptimelevels)
  real(r8), intent(in)   :: u3      (plond,plev,beglatex:endlatex,ptimelevels) ! u-wind com
  real(r8), intent(in)   :: v3      (plond,plev,beglatex:endlatex,ptimelevels) ! v-wind comp
  real(r8), intent(in)   :: t3      (plond,plev,beglatex:endlatex,ptimelevels) ! temperature
  real(r8), intent(in)   :: q3      (plond,plev,pcnst+pnats,beglatex:endlatex,ptimelevels)
!                                                                    ! q and const
  real(r8), intent(in)   :: lnpssld (plond,plev,beglatex:endlatex)   ! RHS Ps term for SLD
  real(r8), intent(in)   :: prhssld (plond,plev,beglatex:endlatex)   ! RHS Ps term for SLD
  real(r8), intent(in)   :: tarrsld (plond,plev,beglatex:endlatex)   ! T  at arr. pt.(SLD)
  real(r8), intent(inout):: parrsld (plond,plev,beglatex:endlatex)   ! Ps at arr. pt.(SLD)
  integer , intent(in)   :: n3                                       ! time index
  integer , intent(in)   :: n3m1                                     ! time index
  real(r8), intent(in)   :: u3sld   (plond,plev ,beglatex:endlatex)  ! u3 inpt for SLD int
  real(r8), intent(in)   :: v3sld   (plond,plev ,beglatex:endlatex)  ! v3 inpt for SLD int
  real(r8), intent(in)   :: etadot  (plond,plevp,beglatex:endlatex,ptimelevels)! Vertical motion

  integer , intent(in)   :: nlon                                     ! # of longitudes
!
!---------------------------Local workspace-----------------------------
!
  integer i                    ! index
  integer k                    ! index
  integer l                    ! index
  integer m                    ! constituent index
  integer inc                  ! increment for fft991
  integer ntr                  ! number of fft's to perform
  integer isign                ! flag indicates fft transform directn
  integer irow                 ! N/S latitude pair index
  integer jcen                 ! lat index (extended grid)
!                              ! of forecast
  real(r8) fdp  (plon,plev,2)  ! interpolant
  real(r8) pmid (plond,plev)   ! pressure at model levels
  real(r8) pint (plond,plevp)  ! pressure at interfaces
  real(r8) pdel (plond,plev)   ! pressure difference between
  real(r8) lamdp(plon,plev)    ! x-coord of dep pt
  real(r8) phidp(plon,plev)    ! y-coord of dep pt
  real(r8) sigdp(plon,plev)    ! z-coord of dep pt

  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
  integer kkdp  (plon,plev)    ! index of z-coordinate of dep pt (alt)

  real(r8) xl   (plon,plev,4)  ! weight for x-interpolants (left)
  real(r8) xr   (plon,plev,4)  ! weight for x-interpolants (right)
  real(r8) wgt1x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)
  real(r8) wgt2x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)
  real(r8) wgt3x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)
  real(r8) wgt4x(plon,plev,4)  ! weight for x-interpolants (Lag Cubic)
  real(r8) hl   (plon,plev,4)  ! weight for x-interpolants (Hermite)
  real(r8) hr   (plon,plev,4)  ! weight for x-interpolants (Hermite)
  real(r8) dhl  (plon,plev,4)  ! weight for x-interpolants (Hermite)
  real(r8) dhr  (plon,plev,4)  ! weight for x-interpolants (Hermite)

  real(r8) ys   (plon,plev)    ! weight for y-interpolants (south)
  real(r8) yn   (plon,plev)    ! weight for y-interpolants (north)
  real(r8) wgt1y(plon,plev)    ! weight for y-interpolants (Lag Cubic)
  real(r8) wgt2y(plon,plev)    ! weight for y-interpolants (Lag Cubic)
  real(r8) wgt3y(plon,plev)    ! weight for y-interpolants (Lag Cubic)
  real(r8) wgt4y(plon,plev)    ! weight for y-interpolants (Lag Cubic)
  real(r8) hs   (plon,plev)    ! weight for y-interpolants (Hermite)
  real(r8) hn   (plon,plev)    ! weight for y-interpolants (Hermite)
  real(r8) dhs  (plon,plev)    ! weight for y-interpolants (Hermite)
  real(r8) dhn  (plon,plev)    ! weight for y-interpolants (Hermite)
  real(r8) rdphi(plon,plev)    ! reciprocal of y-interval

  real(r8) wgt1z(plon,plev)    ! weight for z-interpolants (Lag Cubic)
  real(r8) wgt2z(plon,plev)    ! weight for z-interpolants (Lag Cubic)
  real(r8) wgt3z(plon,plev)    ! weight for z-interpolants (Lag Cubic)
  real(r8) wgt4z(plon,plev)    ! weight for z-interpolants (Lag Cubic)
  real(r8) hb   (plon,plev)    ! weight for z-interpolants (Hermite)
  real(r8) ht   (plon,plev)    ! weight for z-interpolants (Hermite)
  real(r8) dhb  (plon,plev)    ! weight for z-interpolants (Hermite)
  real(r8) dht  (plon,plev)    ! weight for z-interpolants (Hermite)
  real(r8) rdz  (plon,plev)    ! reciprocal of z-interval
  real(r8) zt   (plon,plev)    ! top vertical interpolation weight 
  real(r8) zb   (plon,plev)    ! bot vertical interpolation weight 

  real(r8) lampr(plon,plev)    ! trajectory increment (x-direction)
  real(r8) phipr(plon,plev)    ! trajectory increment (y-direction)
  real(r8) upr  (plon,plev)    ! interpolated u field (local geodesic)
  real(r8) vpr  (plon,plev)    ! interpolated v field (local geodesic)

#if ( ! defined USEFFTLIB )
  real(r8) work((plon+1)*5*plev)   ! workspace array for fft991
#else
  real(r8) work((plon+1)*pcray)    ! workspace array for fft991
#endif
  real(r8) xnlin(plndlv*4 + plond) ! non-linear terms (equivalence 
!                                  ! region for following arrays to 
!                                  ! optimize fft performance)
  real(r8) grfulat(plond,plev)     ! non-linear terms for u-momentum 
  real(r8) grfvlat(plond,plev)     ! non-linear terms for u-momentum 
  real(r8) grtlat (plond,plev)     ! RHS of T-eqn
  real(r8) grqlat (plond,plev)     ! q
  real(r8) grpslat(plond)          ! RHS of Ps-eqn
!
! The following equivalences are for optimal fft performance
!
#if ( defined SUNOS )
!
! Stupid Sun compiler hoop-jumping again
!
  save grfulat, grfvlat, grtlat, grqlat, grpslat, xnlin
#endif
  equivalence   (grfulat,xnlin(1))
  equivalence   (grfvlat,xnlin(1+1*plndlv))
  equivalence   (grtlat ,xnlin(1+2*plndlv))
  equivalence   (grqlat ,xnlin(1+3*plndlv))
  equivalence   (grpslat,xnlin(1+4*plndlv))
!
  real(r8) pd   (plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta) 
  real(r8) pdsum(plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta) 
  real(r8) pd1  (plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta) 
  real(r8) pdsm1(plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta) 
  real(r8) pa   (plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta) 
  real(r8) pasum(plond)            ! RHS term for Ps and (1/ps)etadot(dp/deta) 
  real(r8) coslat                  ! cos(latitude)
  real(r8) tmp1                    ! temp space
!
  logical limdrh                   ! horizontal derivative limiter flag
  logical limdrv                   ! vertical   derivative limiter flag
  logical lhrzint                  ! horizontal interp flag
  logical lvrtint                  ! vertical   interp flag
  logical lhrzwgt                  ! flag to compute horizontal weights
  logical lvrtwgt                  ! flag to compute vertical   weights
!
!-----------------------------------------------------------------------
!
  if(lat.le.plat/2) then
     irow = lat
  else
     irow = plat + 1 - lat
  end if
  jcen = j1 - 1 + lat
  coslat = cos(clat(lat))
!
! Initial guess for trajectory midpoints in spherical coords.
! Use arrival points as initial guess for trajectory midpoints.
!
  do k=1,plev
     do i=1,nlon
        phidp(i,k) = clat(lat)
        sigdp(i,k) = etamid(k)
     end do
  end do
!        
! Offset bottom level departure point first guess by epsilon
!
  do i = 1,nlon
     sigdp(i,plev) = sigdp(i,plev)*(1. - 10.*epsilon(sigdp))
  end do
!
! Loop through latitudes producing departure point calculation
!
  call slttraj(pmap    ,jcen    ,lat     ,ztodt   ,ra      , &
               iter    ,lam     ,phi     ,dphi    ,etamid  , &
               etaint  ,detam   ,detai   ,lbasiy  ,lbasiz  , &
               lbassi  ,kdpmpf  ,kdpmph  ,idp     ,jdp     , &
               kdp     ,kkdp    ,xl      ,xr      ,wgt1x   , &
               wgt2x   ,wgt3x   ,wgt4x   ,hl      ,hr      , &
               dhl     ,dhr     ,ys      ,yn      ,wgt1y   , &
               wgt2y   ,wgt3y   ,wgt4y   ,hs      ,hn      , &
               dhs     ,dhn     ,rdphi   ,wgt1z   ,wgt2z   , &
               wgt3z   ,wgt4z   ,hb      ,ht      ,dhb     , &
               dht     ,rdz     ,lampr   ,phipr   ,upr     , &
               vpr     ,lamdp   ,phidp   ,sigdp   ,u3      , &
               v3      ,u3sld   ,v3sld   ,etadot  ,n3      , &
               n3m1    ,dlam    ,nlon    )
!
! Calculate mass of moisture in field being advected by slt.
!
  call plevs0(nlon    ,plond   ,plev    ,ps(1,lat,n3),pint    ,pmid    ,pdel)
  call qmassa(cwava(lat)  ,w(irow) ,q3(i1,1,1,jcen,n3),pdel    ,hw1lat(1,lat), &
              nlon)
!
! Compute constituent forecast
!
  lhrzwgt = .true.
  lvrtwgt = .true.
  lhrzint = .true.
  lvrtint = .true.
  limdrh  = .true.
  limdrv  = .true.
  call bandij (dlam    ,phi     ,lamdp   ,phidp   ,idp     , &
               jdp     ,nlon    )
  call kdpfnd (plev    ,pmap    ,etamid  ,sigdp   ,kdpmpf  , &
               kdp     ,nlon    )
  call sltwgts(limdrh  ,limdrv  ,lhrzwgt ,lvrtwgt ,plev    , &
               idp     ,jdp     ,kdp     ,lam     ,phi     , &
               etamid  ,dphi    ,detam   ,lamdp   ,phidp   , &
               sigdp   ,lbasiy  ,lbasiz  ,kkdp    ,xl      , &
               xr      ,wgt1x   ,wgt2x   ,wgt3x   ,wgt4x   , &
               hl      ,hr      ,dhl     ,dhr     ,ys      , &
               yn      ,wgt1y   ,wgt2y   ,wgt3y   ,wgt4y   , &
               hs      ,hn      ,dhs     ,dhn     ,rdphi   , &
               wgt1z   ,wgt2z   ,wgt3z   ,wgt4z   ,hb      , &
               ht      ,dhb     ,dht     ,rdz     ,zt      , &