grmult.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine grmult (ztodt   ,rcoslat ,coslat  ,coriol  ,div     , &
                   q3      ,t3      ,u3      ,v3      ,u3sld   , &
                   v3sld   ,ed1     ,ql      ,qm      ,tl      , &
                   tm      ,phis    ,phisl   ,phism   ,dpsl    , &
                   dpsm    ,omga    ,pmid    ,pdel    ,ps      , &
                   logps   ,rpmid   ,etamid  ,fu      ,fv      , &
                   t2      ,t3m1    ,u3m1    ,v3m1    ,etadot  , &
                   dpslon  ,dpslat  ,urhs    ,vrhs    ,trhs    , &
                   prhs    ,etadotm1,lnpssld ,prhssld ,tarrsld , &
                   parrsld ,nlon    )
!-----------------------------------------------------------------------
!
! Purpose:
! Compute non-linear dynamics terms in grid point space (in preparation
! for SLD interpolation)
!
!---------------------------Code history--------------------------------
!
! Author:  J. Olson
!
!-----------------------------------------------------------------------
!
! $Id: grmult.F90,v 1.6.2.1 2002/04/22 19:09:50 erik Exp $
! $Author: erik $
!
!-----------------------------------------------------------------------

  use precision
  use pmgrid
  use pspect
  use comslt
  use commap
  use physconst, only: rair, zvir, cappa, cpvir
  use time_manager, only: is_first_step

  implicit none

#include <comhyb.h>
!------------------------------Arguments--------------------------------
!
  real(r8), intent(in)   :: ztodt                 ! delta-t
  real(r8), intent(in)   :: rcoslat               ! 1/(cos lat)
  real(r8), intent(in)   :: coslat                ! cos(lat)
  real(r8), intent(in)   :: coriol                ! Coriolis parameter
  real(r8), intent(in)   :: div     (plond,plev)  ! divergence
  real(r8), intent(in)   :: q3      (plond,plev)  ! specific humidity
  real(r8), intent(in)   :: t3      (plond,plev)  ! temperature
  real(r8), intent(in)   :: u3      (plond,plev)  ! zonal wind
  real(r8), intent(in)   :: v3      (plond,plev)  ! meridional wind
  real(r8), intent(out)  :: u3sld   (plond,plev)  ! u-wind (time n) (used for advection)
  real(r8), intent(out)  :: v3sld   (plond,plev)  ! v-wind (time n) (used for advection)
  real(r8), intent(in)   :: ed1     (plond,plev)  ! (1/ps)etadot(dp/deta) (time n-1)
  real(r8), intent(in)   :: ql      (plond,plev)  ! longitudinal derivative of q
  real(r8), intent(in)   :: qm      (plond,plev)  ! latitudinal  derivative of q
  real(r8), intent(in)   :: tl      (plond,plev)  ! zonal derivative of T
  real(r8), intent(in)   :: tm      (plond,plev)  ! meridional derivative of T
  real(r8), intent(in)   :: phis    (plond)       ! Phi at surface
  real(r8), intent(in)   :: phisl   (plond)       ! longitudinal derivative of phis
  real(r8), intent(in)   :: phism   (plond)       ! latitudinal  devivative of phis
  real(r8), intent(in)   :: dpsl    (plond)       ! longitudinal component of grad ln(ps)
  real(r8), intent(in)   :: dpsm    (plond)       ! latitudinal  component of grad ln(ps)
  real(r8), intent(in)   :: omga    (plond,plev)  ! vertical pressure velocity
  real(r8), intent(in)   :: pmid    (plond,plev)  ! pressure at full levels
  real(r8), intent(in)   :: pdel    (plond,plev)  ! layer thicknesses (pressure)
  real(r8), intent(in)   :: ps      (plond)       ! Surface pressure (n)
  real(r8), intent(in)   :: logps   (plond)       ! log(ps)
  real(r8), intent(in)   :: rpmid   (plond,plev)  ! 1./pmid
  real(r8), intent(in)   :: etamid  (plev)        ! midpoint values of eta (a+b)
  real(r8), intent(inout):: fu      (plond,plev)  ! nonlinear term - u momentum eqn
  real(r8), intent(inout):: fv      (plond,plev)  ! nonlinear term - v momentum eqn
  real(r8), intent(inout):: t2      (plond,plev)  ! nonlinear term - temperature
  real(r8), intent(out)  :: t3m1    (plond,plev)  ! T at time n-1
  real(r8), intent(inout):: u3m1    (plond,plev)  ! U at previous time step
  real(r8), intent(inout):: v3m1    (plond,plev)  ! V at previous time step
  real(r8), intent(out)  :: etadot  (plond,plevp) ! vertical velocity in eta coordinates
  real(r8), intent(out)  :: dpslon  (plond,plev)  ! Pressure gradient term
  real(r8), intent(out)  :: dpslat  (plond,plev)  ! Pressure gradient term
  real(r8), intent(inout):: urhs    (plond,plev)  ! RHS of U  eqn valid for mid-point
  real(r8), intent(inout):: vrhs    (plond,plev)  ! RHS of V  eqn valid for mid-point
  real(r8), intent(inout):: trhs    (plond,plev)  ! RHS of T  eqn valid for mid-point
  real(r8), intent(inout):: prhs    (plond,plev)  ! RHS of Ps eqn valid for mid-point
  real(r8), intent(inout):: etadotm1(plond,plevp) ! etadot for time n-1
  real(r8), intent(out)  :: lnpssld (plond,plev)  ! RHS Ps term for SLD
  real(r8), intent(out)  :: prhssld (plond,plev)  ! RHS Ps term for SLD
  real(r8), intent(out)  :: tarrsld (plond,plev)  ! T  at arr. pt. (SLD)
  real(r8), intent(out)  :: parrsld (plond,plev)  ! Ps at arr. pt. (SLD)
  integer , intent(in)   :: nlon                  ! number of longitudes for this latitude
!
!---------------------------Local workspace-----------------------------
!
  real(r8) tmp1                  ! temporary workspace
  real(r8) tmp2                  ! temporary workspace
  real(r8) tmp                   ! temporary workspace
  real(r8) tmpk                  ! workspace
  real(r8) tmpkp1                ! workspace
  real(r8) u3l     (plond,plev)  ! u-wind used locally only
  real(r8) v3l     (plond,plev)  ! v-wind used locally only
  real(r8) tv      (plond,plev)  ! virtual temperature
  real(r8) ddpk    (plond)       ! partial sum of div*delta p
  real(r8) ddpn    (plond)       ! complete sum of div*delta p
  real(r8) vkdp    (plond,plev)  ! V dot grad(ln(ps))
  real(r8) vpdsk   (plond)       ! partial sum  V dot grad(ln(ps)) delta b
  real(r8) vpdsn   (plond)       ! complete sum V dot grad(ln(ps)) delta b
  real(r8) lpsstar (plond)       ! Reference ln(Ps) (used to define a new 
!                                ! perturbation Ps)
  real(r8) lpsstarl(plond)       ! long. grad of reference ln(Ps)
  real(r8) lpsstarm(plond)       ! lat.  grad of reference ln(Ps)
  real(r8) rtv     (plond,plev)  ! rair*(tv+t0)
  real(r8) dt                    ! time step
  real(r8) hsl     (plond,plev)  ! zonal      deriv of hydrostatic term
  real(r8) hsm     (plond,plev)  ! meridional deriv of hydrostatic term
  real(r8) pspsl   (plond)       ! Ps*d(lnPs)/d(long.)
  real(r8) pspsm   (plond)       ! Ps*d(lnPs)/d(lat. )
  real(r8) ed1p    (plond,plevp) ! (1/ps)etadot(dp/deta) (time n-1)
  real(r8) rtvl                  ! zonal      derivative of R*Tv
  real(r8) rtvm                  ! meridional derivative of R*Tv
  real(r8) abp0                  ! constant for grad(H(n)) matrix
!
! Arrays which hold results from the RHS of the prognostic equations.
! Results will be interpolated to trajectory departure points in the routine
! SCANSLT.
!
  real(r8) tsld0a  (plond,plev)  ! RHS of T  eqn valid for mid-point
  real(r8) usldm   (plond,plev)  ! RHS of U  eqn valid for departure pt (n)
  real(r8) vsldm   (plond,plev)  ! RHS of V  eqn valid for departure pt (n)
  real(r8) tsldm   (plond,plev)  ! RHS of T  eqn valid for departure pt (n)
  real(r8) psldm   (plond,plev)  ! RHS of Ps eqn valid for departure pt (n)
  real(r8) urhsl   (plond,plev)  ! RHS of U  eqn valid for mid-point (n+1/2)
  real(r8) vrhsl   (plond,plev)  ! RHS of V  eqn valid for mid-point (n+1/2)
  real(r8) trhsl   (plond,plev)  ! RHS of T  eqn valid for mid-point (n+1/2)
  real(r8) prhsl   (plond,plev)  ! RHS of Ps eqn valid for mid-point (n+1/2)

  real(r8) onemeps               ! 1 - epssld (SLD decentering coefficient)
  real(r8) onepeps               ! 1 + epssld (SLD decentering coefficient)
  real(r8) detai(plevp)          ! interval between interfaces
  real(r8) facm1                 ! interpolation factor for time n-1
  real(r8) facm2                 ! interpolation factor for time n-2

  integer i,l,k                  ! longitude, level indices
  integer npr                    ! index
!
!-----------------------------------------------------------------------
!
  onemeps = 1. - epssld
  onepeps = 1. + epssld
  facm1   =  3./2.
  facm2   = -1./2.
!
  do k = 1,plev
     detai  (k) = (hyai(k+1) + hybi(k+1)) - (hyai(k) + hybi(k))
  end do
!
! Compute U/V for time n + 1/2
! The first formula will be used later for trajectory calculation
! The second will be used locally in the Hortal Temperature correction
!
  do k = 1,plev
     do i = 1,nlon
        u3sld(i,k) = 2.   *u3(i,k) -       u3m1(i,k)
        v3sld(i,k) = 2.   *v3(i,k) -       v3m1(i,k)
        u3l  (i,k) = facm1*u3(i,k) + facm2*u3m1(i,k)
        v3l  (i,k) = facm1*v3(i,k) + facm2*v3m1(i,k)
     end do
  end do
!
! Zero auxiliary fields
!
  tmp = 1./(rair*t0(plev))
  do i=1,nlon
     ddpk    (i) = 0.0
     ddpn    (i) = 0.0
     vpdsk   (i) = 0.0
     vpdsn   (i) = 0.0
     lpsstar (i) = -phis (i)*tmp
     lpsstarl(i) = -phisl(i)*tmp
     lpsstarm(i) = -phism(i)*tmp
     pspsl   (i) = ps(i)*dpsl(i)
     pspsm   (i) = ps(i)*dpsm(i)
     etadot(i,1) = 0.0
     ed1p  (i,1) = 0.0
     etadot(i,plevp) = 0.0
  end do
!
! Virtual temperature
!
  call virtem(nlon, plond, plev, t3      ,q3      ,zvir    ,tv)
!
! calculate some auxiliary quantities
!
  do k=1,plev
     do i=1,nlon
        ed1p(i,k+1) = ed1(i,k)
        rtv(i,k) = rair*tv(i,k)
!
! sum(plev)(div(k)*dp(k))
!
        ddpn(i) = ddpn(i) + div(i,k)*pdel(i,k)
     end do
  end do
!
! sum(plev)(v(k)*grad(lnps)*db(k))
!
  do k=nprlev,plev
     do i=1,nlon
        vkdp(i,k)= rcoslat*(u3(i,k)*pspsl(i) + v3(i,k)*pspsm(i))
        vpdsn(i) = vpdsn(i) + vkdp(i,k)*hybd(k)
     end do
  end do
!
! Compute etadot (top and bottom = 0.)
!
  do k = 1,plev-1
!
! Compute etadot(dp/deta)(k+1/2) and sum(k)(div(j)*dp(j))
!
     do i=1,nlon
        ddpk(i) = ddpk(i) + div(i,k)*pdel(i,k)
#ifdef HADVTEST
!
!jr Set etadot to zero for horizontal advection test
!
        etadot(i,k+1) = 0.
#else
        etadot(i,k+1) = -ddpk(i)
#endif
     end do
!
! sum(k)(v(j)*grad(ps)*db(j))
!
     if (k.ge.nprlev) then
        do i=1,nlon