grmult.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

           vpdsk(i) = vpdsk(i) + vkdp(i,k)*hybd(k)
#ifdef HADVTEST
!
!jr Set etadot to zero for horizontal advection test
!
           etadot(i,k+1) = 0.
#else
           etadot(i,k+1) = etadot(i,k+1) - vpdsk(i) + hybi(k+1)*(ddpn(i)+vpdsn(i))
#endif
        end do
     end if
!
! Convert eta-dot(dp/deta) to eta-dot
!
     tmp = etamid(k+1) - etamid(k)
     do i = 1,nlon
        etadot(i,k+1) = etadot(i,k+1)*tmp/(pmid(i,k+1) - pmid(i,k))
     end do
  end do
!
! Zonal and meridional derivatives of the hydrostatic term in the
! momentum equations.
!
  do k = 1,plev
     do i = 1,nlon
        tmp1     = (1. + zvir*q3(i,k))
        tmp2     = t3(i,k)*zvir
        rtvl     = rair*( tmp1*tl(i,k) + tmp2*ql(i,k) )
        rtvm     = rair*( tmp1*tm(i,k) + tmp2*qm(i,k) )
!
        tmp      = rpmid(i,k)*pdel(i,k)
        hsl(i,k) = -rtvl*tmp
        hsm(i,k) = -rtvm*tmp
     end do
     if(k .ge. nprlev) then
        abp0 = ps0*(hyam(k)*(hybi(k+1) - hybi(k)) - hybm(k)*(hyai(k+1) - hyai(k)))
        do i = 1,nlon
           tmp      = rtv(i,k)*abp0/(pmid(i,k)*pmid(i,k))
           hsl(i,k) = hsl(i,k) - pspsl(i)*tmp
           hsm(i,k) = hsm(i,k) - pspsm(i)*tmp
        end do
     endif
  end do
!
! Calculate RHS of all prognostic eqns
!
  do k = 1,plev
     tmp  = cappa*t0(k)*hypi(plevp)/hypm(k)
     tmp1 = hypd(k)/hypi(plevp)
     do i = 1,nlon
!
! Surface Pressure eqn
!
        prhsl(i,k) =  div(i,k)*(tmp1 - pdel(i,k)/ps(i))
        psldm(i,k) = -div(i,k)*tmp1 - ed1p(i,k+1) + ed1p(i,k)
!
! Temperature eqn
!
        trhsl (i,k) = cappa*tv(i,k)/(1. + cpvir*q3(i,k))*omga(i,k)*rpmid(i,k)
        trhsl (i,k) = trhsl(i,k) - omga(i,k)/ps(i)*tmp
        tsldm (i,k) = 0.5*(ed1p(i,k+1) + ed1p(i,k))*tmp
!
! ... horizontal advection portion of Hortal Temperature correction
!
        trhsl (i,k) = trhsl(i,k) - &
                      rcoslat*( u3(i,k)*lpsstarl(i) + v3(i,k)*lpsstarm(i) )*hortalc(k)
!
! ... Ritchie damping term for Temperature eqn.
!
        tsld0a(i,k) = rcoslat*( u3l(i,k)*lpsstarl(i) + v3l(i,k)*lpsstarm(i) )
!
! U/V eqns (includes only the diagonal portion of the hydrostatic term)
!
        urhsl(i,k) = 0.5*hsl(i,k) + href(k,k)*tl(i,k) + bps(k)*dpsl(i)
        vrhsl(i,k) = 0.5*hsm(i,k) + href(k,k)*tm(i,k) + bps(k)*dpsm(i)
        usldm(i,k) = -phisl(i) + v3(i,k)*coriol*coslat -href(k,k)*tl(i,k) - bps(k)*dpsl(i)
        vsldm(i,k) = -phism(i) - u3(i,k)*coriol*coslat -href(k,k)*tm(i,k) - bps(k)*dpsm(i)
     end do
!
! Add pressure gradient terms to momentum tendencies
!
     if (k.ge.nprlev) then
        do i=1,nlon
           tmp        = rtv(i,k)*hybm(k)*rpmid(i,k)
           dpslon(i,k) = rcoslat*tmp*pspsl(i)
           dpslat(i,k) = rcoslat*tmp*pspsm(i)
           urhsl(i,k) = urhsl(i,k) - dpslon(i,k)*coslat
           vrhsl(i,k) = vrhsl(i,k) - dpslat(i,k)*coslat
        end do
     else
        do i = 1,nlon
           dpslon(i,k) = 0.
           dpslat(i,k) = 0.
        end do
     end if
  end do
!
! Interior levels of the hydrostatic term
!
  do k=1,plev-1
     do l=k+1,plev
        do i=1,nlon
           urhsl(i,k) = urhsl(i,k) + href(l,k)*tl(i,l) + hsl(i,l)
           vrhsl(i,k) = vrhsl(i,k) + href(l,k)*tm(i,l) + hsm(i,l)
!
           usldm(i,k) = usldm(i,k) - href(l,k)*tl(i,l)
           vsldm(i,k) = vsldm(i,k) - href(l,k)*tm(i,l)
        end do
     end do
  end do

  if(is_first_step()) then
     do k = 1,plevp
        do i = 1,nlon
           etadotm1(i,k) = etadot(i,k)
        end do
     end do
  endif
!
! The modified etadotm1 will be used later for trajectory calculation in SCANSLT
!
  do k = 1,plevp
     do i = 1,nlon
        etadotm1(i,k) = 2.*etadot(i,k) - etadotm1(i,k)
     end do
  end do
!
! Compute vertical advection portion of Hortal Temperature correction
!
  npr = nprlev - 1
  if(npr .lt. 1) npr = 1
  do k = npr,plev-1
     tmpk   = 0.5*hdel(k-1)/detai(k)
     tmpkp1 = 0.5*hdel(k  )/detai(k)
     do i = 1,nlon
        trhsl(i,k) = trhsl(i,k) - (etadot(i,k  )*tmpk + etadot(i,k+1)*tmpkp1)*lpsstar(i)
     end do
  end do
!
! ... bottom level
!
  tmpk = 0.5*hdel(plev-1)/detai(plev)
  do i = 1,nlon
     trhsl(i,plev) = trhsl(i,plev) - etadot(i,plev)*tmpk*lpsstar(i)
  end do
!
! Compute (by extrapolation) RHS terms for time n + 1/2
!
  if(is_first_step()) then
     do k = 1,plev
        do i = 1,nlon
           urhs  (i,k) = urhsl (i,k)
           vrhs  (i,k) = vrhsl (i,k)
           trhs  (i,k) = trhsl (i,k)
           prhs  (i,k) = prhsl (i,k)
        end do
     end do
  endif
!
  do k = 1,plev
     do i = 1,nlon
        tmp         =       urhsl (i,k)
        urhsl (i,k) = facm1*urhsl (i,k) + facm2*urhs  (i,k)
        urhs  (i,k) = tmp
        tmp         =       vrhsl (i,k)
        vrhsl (i,k) = facm1*vrhsl (i,k) + facm2*vrhs  (i,k)
        vrhs  (i,k) = tmp
        tmp         =       trhsl (i,k)
        trhsl (i,k) = facm1*trhsl (i,k) + facm2*trhs  (i,k)
        trhs  (i,k) = tmp
        tmp         =       prhsl (i,k)
        prhsl (i,k) = facm1*prhsl (i,k) + facm2*prhs  (i,k)
        prhs  (i,k) = tmp
     end do
  end do
!
! Prepare SLD arrays for interpolation by SCANSLT
!
! Append appropriate coefficients (including decentering epsilon values -- See
! Notes 13,14,15)
!
  dt   = 0.5*ztodt
  tmp1 = 0.5*ztodt/coslat
  do k = 1,plev
     tmp = cappa*t0(k)*hypi(plevp)/hypm(k)
     do i = 1,nlon
        urhsl (i,k) = tmp1*urhsl (i,k)
        vrhsl (i,k) = tmp1*vrhsl (i,k)
        trhsl (i,k) = dt  *trhsl (i,k)
        tsld0a(i,k) = dt  *tsld0a(i,k)
        prhsl (i,k) = dt  *prhsl (i,k)
        t2    (i,k) = ztodt*t2   (i,k)
        fu    (i,k) = ztodt*fu   (i,k)
        fv    (i,k) = ztodt*fv   (i,k)
!
! Combine terms.
! (Time split the midpoint results between arrival and departure
!  points)
!
        tmp2          = lpsstar(i)*hortalc(k)
        u3m1    (i,k) = (u3  (i,k)*coslat + onemeps*usldm(i,k)*dt)/coslat + &
                                            onemeps*urhsl(i,k) + &
                                            onemeps*fu   (i,k)*0.5
        v3m1    (i,k) = (v3  (i,k)*coslat + onemeps*vsldm(i,k)*dt)/coslat + &
                                            onemeps*vrhsl(i,k) + &
                                            onemeps*fv   (i,k)*0.5
#ifdef HADVTEST
        t3m1    (i,k) = t3  (i,k)
#else
        t3m1    (i,k) = t3  (i,k) +         onemeps*tsldm(i,k)*dt + &
                                            onemeps*trhsl(i,k) - tmp2 + &
                                            onemeps*t2(i,k)*0.5
#endif
        prhssld (i,k) = (psldm(i,k)*dt + prhsl(i,k))*onemeps
        tarrsld (i,k) = (trhsl(i,k) + &
                         hybm(k)*tmp*tsld0a(i,k))*onepeps + tmp2 + onepeps*t2(i,k)*0.5
        parrsld (i,k) = (prhsl(i,k) - hybd(k)*tsld0a(i,k))*onepeps
        lnpssld (i,k) = logps(i) - lpsstar(i) - tsld0a(i,k)*onemeps
        fu      (i,k) = urhsl(i,k)*coslat*onepeps + fu(i,k)*coslat*onepeps*0.5
        fv      (i,k) = vrhsl(i,k)*coslat*onepeps + fv(i,k)*coslat*onepeps*0.5
     end do
  end do
#ifdef VADVTEST
  do k=2,plev
     do i=1,nlon
        etadot(i,k) = -0.5/(12.*86400.)
     end do
  end do
#endif
!
  return
end subroutine grmult