#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