grcalc.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

      grvh1s(:) = 0.
!
! Loop over n for t,q,d,and end of u and v
!
      do m=1,nmmax(irow)
         mr = nstart(m)
         mc = 2*mr
         do n=1,nlen(m),2
            ir = mc + 2*n - 1
            ii = ir + 1
            grts (2*m-1,k) = grts (2*m-1,k) + t(ir,k)*alp(mr+n,irow)
            grts (2*m  ,k) = grts (2*m  ,k) + t(ii,k)*alp(mr+n,irow)
!
            tmp = alp(mr+n,irow)*hdiftq(n+m-1,k)
            grths(2*m-1,k) = grths(2*m-1,k) - t(ir,k)*tmp
            grths(2*m  ,k) = grths(2*m  ,k) - t(ii,k)*tmp
!
            grds(2*m-1,k) = grds(2*m-1,k) + d(ir,k)*alp(mr+n,irow)
            grds(2*m  ,k) = grds(2*m  ,k) + d(ii,k)*alp(mr+n,irow)
!
            grzs(2*m-1,k) = grzs(2*m-1,k) + vz(ir,k)*alp(mr+n,irow)
            grzs(2*m  ,k) = grzs(2*m  ,k) + vz(ii,k)*alp(mr+n,irow)
!
            gru1s (2*m-1) = gru1s (2*m-1) + d(ir,k)*alpn(mr+n)
            gru1s (2*m  ) = gru1s (2*m  ) + d(ii,k)*alpn(mr+n)
!
            tmp = alpn(mr+n)*hdifzd(n+m-1,k)
            gruh1s(2*m-1) = gruh1s(2*m-1) - d(ir,k)*tmp
            gruh1s(2*m  ) = gruh1s(2*m  ) - d(ii,k)*tmp
!
            grv1s (2*m-1) = grv1s (2*m-1) + vz(ir,k)*alpn(mr+n)
            grv1s (2*m  ) = grv1s (2*m  ) + vz(ii,k)*alpn(mr+n)
!
            grvh1s(2*m-1) = grvh1s(2*m-1) - vz(ir,k)*tmp
            grvh1s(2*m  ) = grvh1s(2*m  ) - vz(ii,k)*tmp
         end do
      end do

      do m=1,nmmax(irow)
         mr = nstart(m)
         mc = 2*mr
         do n=2,nlen(m),2
            ir = mc + 2*n - 1
            ii = ir + 1
!
            grus (2*m-1,k) = grus (2*m-1,k) + vz(ir,k)*dalpn(mr+n)
            grus (2*m  ,k) = grus (2*m  ,k) + vz(ii,k)*dalpn(mr+n)
!
            tmp = dalpn(mr+n)*hdifzd(n+m-1,k)
            gruhs(2*m-1,k) = gruhs(2*m-1,k) - vz(ir,k)*tmp
            gruhs(2*m  ,k) = gruhs(2*m  ,k) - vz(ii,k)*tmp
!
            grvs (2*m-1,k) = grvs (2*m-1,k) - d(ir,k)*dalpn(mr+n)
            grvs (2*m  ,k) = grvs (2*m  ,k) - d(ii,k)*dalpn(mr+n)
!
            grvhs(2*m-1,k) = grvhs(2*m-1,k) + d(ir,k)*tmp
            grvhs(2*m  ,k) = grvhs(2*m  ,k) + d(ii,k)*tmp
         end do
      end do
!
! Combine the two parts of u(m) and v(m)
!
      do m=1,nmmax(irow)
         grus (2*m-1,k) = grus (2*m-1,k) + gru1s (2*m  )
         gruhs(2*m-1,k) = gruhs(2*m-1,k) + gruh1s(2*m  )
         grus (2*m  ,k) = grus (2*m  ,k) - gru1s (2*m-1)
         gruhs(2*m  ,k) = gruhs(2*m  ,k) - gruh1s(2*m-1)
         grvs (2*m-1,k) = grvs (2*m-1,k) + grv1s (2*m  )
         grvhs(2*m-1,k) = grvhs(2*m-1,k) + grvh1s(2*m  )
         grvs (2*m  ,k) = grvs (2*m  ,k) - grv1s (2*m-1)
         grvhs(2*m  ,k) = grvhs(2*m  ,k) - grvh1s(2*m-1)
      end do
!
! Remove Coriolis contribution to absolute vorticity from u(m)
! Correction for u:zeta=vz-ez=(zeta+f)-f
!
      grus(1,k) = grus(1,k) - zurcor
   end do
!
!-----------------------------------------------------------------------
!
! Computation for 1-level variables (ln(p*) and derivatives).
!
   do m=1,nmmax(irow)
      mr = nstart(m)
      mc = 2*mr
      do n=1,nlen(m),2
         ir = mc + 2*n - 1
         ii = ir + 1
!
         tmpr = alps(ir)*alp(mr+n,irow)
         tmpi = alps(ii)*alp(mr+n,irow)
         
         grpss(2*m-1) = grpss(2*m-1) + tmpr
         grpss(2*m  ) = grpss(2*m  ) + tmpi
!
         grdpss(2*m-1) = grdpss(2*m-1) + tmpr*hdfst4(m+n-1)*ztodt
         grdpss(2*m  ) = grdpss(2*m  ) + tmpi*hdfst4(m+n-1)*ztodt
      end do
   end do

   do m=1,nmmax(irow)
      mr = nstart(m)
      mc = 2*mr
      do n=2,nlen(m),2
         ir = mc + 2*n - 1
         ii = ir + 1
!
         grpms(2*m-1) = grpms(2*m-1) + alps(ir)*dalp(mr+n,irow)*ra
         grpms(2*m  ) = grpms(2*m  ) + alps(ii)*dalp(mr+n,irow)*ra
      end do
!
! Multiply by m/a to get d(ln(p*))/dlamda
! and by 1/a to get (1-mu**2)d(ln(p*))/dmu
!
      grpls(2*m-1) = -grpss(2*m  )*ra*xm(m)
      grpls(2*m  ) =  grpss(2*m-1)*ra*xm(m)
   end do
!
   return
end subroutine grcalcs


subroutine grcalca (irow    ,ztodt   ,grta    ,grtha   ,grda    ,&
                    grza    ,grua    ,gruha   ,grva    ,grvha   ,&
                    grpsa   ,grdpsa  ,grpma   ,grpla   )

!-----------------------------------------------------------------------
!
! Complete inverse Legendre transforms from spectral to Fourier space at 
! the the given latitude. Only positive latitudes are considered and 
! symmetric and antisymmetric (about equator) components are computed. 
! The sum and difference of these components give the actual fourier 
! coefficients for the latitude circle in the northern and southern 
! hemispheres respectively.
!
! The naming convention is as follows:
!  - The fourier coefficient arrays all begin with "gr";
!  - "t, q, d, z, ps" refer to temperature, specific humidity, 
!     divergence, vorticity, and surface pressure;
!  - "h" refers to the horizontal diffusive tendency for the field.
!  - "s" suffix to an array => symmetric component;
!  - "a" suffix to an array => antisymmetric component.
! Thus "grts" contains the symmetric Fourier coeffs of temperature and
! "grtha" contains the antisymmetric Fourier coeffs of the temperature
! tendency due to horizontal diffusion.
! Three additional surface pressure related quantities are returned:
!  1. "grdpss" and "grdpsa" contain the surface pressure factor
!      (proportional to del^4 ps) used for the partial correction of 
!      the horizontal diffusion to pressure surfaces.
!  2. "grpms" and "grpma" contain the longitudinal component of the 
!      surface pressure gradient.
!  3. "grpls" and "grpla" contain the latitudinal component of the 
!      surface pressure gradient.
!
!---------------------------Code history--------------------------------
!
! Original version:  CCM1
! Standardized:      J. Rosinski, June 1992
! Reviewed:          B. Boville, D. Williamson, J. Hack, August 1992
! Reviewed:          B. Boville, D. Williamson, April 1996
!
!-----------------------------------------------------------------------
!
! $Id: grcalc.F90,v 1.5 2001/09/16 22:13:25 rosinski Exp $
! $Author: rosinski $
!
   use precision
   use pmgrid
   use pspect
   use comspe
   use rgrid
   use commap
   use dynconst, only: ra

   implicit none

#include <comhd.h>
!
! Input arguments
!
   integer, intent(in) :: irow         ! latitude pair index
   real(r8), intent(in) :: ztodt       ! twice the timestep unless nstep = 0
!
! Output arguments: antisymmetric fourier coefficients
!
   real(r8), intent(out) :: grta(plond,plev)    ! sum(n) of t(n,m)*P(n,m)
   real(r8), intent(out) :: grtha(plond,plev)   ! sum(n) of K(2i)*t(n,m)*P(n,m)
   real(r8), intent(out) :: grda(plond,plev)    ! sum(n) of d(n,m)*P(n,m)
   real(r8), intent(out) :: grza(plond,plev)    ! sum(n) of z(n,m)*P(n,m)
   real(r8), intent(out) :: grua(plond,plev)    ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(out) :: gruha(plond,plev)   ! sum(n) of K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(out) :: grva(plond,plev)    ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(out) :: grvha(plond,plev)   ! sum(n) of K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(out) :: grpsa(plond)        ! sum(n) of lnps(n,m)*P(n,m)
   real(r8), intent(out) :: grdpsa(plond)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
   real(r8), intent(out) :: grpma(plond)        ! sum(n) of lnps(n,m)*H(n,m)
   real(r8), intent(out) :: grpla(plond)        ! sum(n) of lnps(n,m)*P(n,m)*m/a
!
!---------------------------Local workspace-----------------------------
!
   real(r8) gru1a(plond)        ! sum(n) of d(n,m)*P(n,m)*m*a/(n(n+1))
   real(r8) gruh1a(plond)       ! sum(n) of K(2i)*d(n,m)*P(n,m)*m*a/(n(n+1))
   real(r8) grv1a(plond)        ! sum(n) of z(n,m)*P(n,m)*m*a/(n(n+1))
   real(r8) grvh1a(plond)       ! sum(n) of K(2i)*z(n,m)*P(n,m)*m*a/(n(n+1))
   real(r8) alpn(pspt)          ! (a*m/(n(n+1)))*Legendre functions (complex)
   real(r8) dalpn(pspt)         ! (a/(n(n+1)))*derivative of Legendre functions (complex)

   integer k                ! level index
   integer m                ! Fourier wavenumber index of spectral array
   integer n                ! meridional wavenumber index
   integer ir,ii            ! spectral indices
   integer mr,mc            ! spectral indices

   real(r8) tmp,tmpr,tmpi,raxm  ! temporary workspace
!
!-----------------------------------------------------------------------
!
! Compute alpn and dalpn
!
   do m=1,nmmax(irow)
      mr = nstart(m)
      raxm = ra*xm(m)
      do n=1,nlen(m)
         alpn(mr+n) = alp(mr+n,irow)*rsq(m+n-1)*raxm
         dalpn(mr+n) = dalp(mr+n,irow)*rsq(m+n-1)*ra
      end do
   end do
!
! Initialize sums
!
   grza(:,:)  = 0.
   grda(:,:)  = 0.
   gruha(:,:) = 0.
   grvha(:,:) = 0.
   grtha(:,:) = 0.
   grpsa(:)   = 0.
   grua(:,:)  = 0.
   grva(:,:)  = 0.
   grta(:,:)  = 0.
   grpla(:)   = 0.
   grpma(:)   = 0.
   grdpsa(:)   = 0.

   do k=1,plev
      gru1a(:) = 0.
      gruh1a(:) = 0.
      grv1a(:) = 0.
      grvh1a(:) = 0.
!
! Loop over n for t,q,d,and end of u and v
!
      do m=1,nmmax(irow)
         mr = nstart(m)
         mc = 2*mr
         do n=1,nlen(m),2
            ir = mc + 2*n - 1
            ii = ir + 1
!
            grua (2*m-1,k) = grua (2*m-1,k) + vz(ir,k)*dalpn(mr+n)
            grua (2*m  ,k) = grua (2*m  ,k) + vz(ii,k)*dalpn(mr+n)
!
            tmp = dalpn(mr+n)*hdifzd(n+m-1,k)
            gruha(2*m-1,k) = gruha(2*m-1,k) - vz(ir,k)*tmp
            gruha(2*m  ,k) = gruha(2*m  ,k) - vz(ii,k)*tmp
!
            grva (2*m-1,k) = grva (2*m-1,k) - d(ir,k)*dalpn(mr+n)
            grva (2*m  ,k) = grva (2*m  ,k) - d(ii,k)*dalpn(mr+n)
!
            grvha(2*m-1,k) = grvha(2*m-1,k) + d(ir,k)*tmp
            grvha(2*m  ,k) = grvha(2*m  ,k) + d(ii,k)*tmp
         end do
      end do

      do m=1,nmmax(irow)
         mr = nstart(m)
         mc = 2*mr
         do n=2,nlen(m),2
            ir = mc + 2*n - 1
            ii = ir + 1
            grta (2*m-1,k) = grta (2*m-1,k) + t(ir,k)*alp(mr+n,irow)
            grta (2*m  ,k) = grta (2*m  ,k) + t(ii,k)*alp(mr+n,irow)
!
            tmp = alp(mr+n,irow)*hdiftq(n+m-1,k)
            grtha(2*m-1,k) = grtha(2*m-1,k) - t(ir,k)*tmp
            grtha(2*m  ,k) = grtha(2*m  ,k) - t(ii,k)*tmp
!
            grda(2*m-1,k) = grda(2*m-1,k) + d(ir,k)*alp(mr+n,irow)
            grda(2*m  ,k) = grda(2*m  ,k) + d(ii,k)*alp(mr+n,irow)
!
            grza(2*m-1,k) = grza(2*m-1,k) + vz(ir,k)*alp(mr+n,irow)
            grza(2*m  ,k) = grza(2*m  ,k) + vz(ii,k)*alp(mr+n,irow)
!
            gru1a (2*m-1) = gru1a (2*m-1) + d(ir,k)*alpn(mr+n)
            gru1a (2*m  ) = gru1a (2*m  ) + d(ii,k)*alpn(mr+n)
!
            tmp = alpn(mr+n)*hdifzd(n+m-1,k)
            gruh1a(2*m-1) = gruh1a(2*m-1) - d(ir,k)*tmp
            gruh1a(2*m  ) = gruh1a(2*m  ) - d(ii,k)*tmp
!
            grv1a (2*m-1) = grv1a (2*m-1) + vz(ir,k)*alpn(mr+n)
            grv1a (2*m  ) = grv1a (2*m  ) + vz(ii,k)*alpn(mr+n)
!
            grvh1a(2*m-1) = grvh1a(2*m-1) - vz(ir,k)*tmp
            grvh1a(2*m  ) = grvh1a(2*m  ) - vz(ii,k)*tmp
         end do
      end do
!
! Combine the two parts of u(m) and v(m)
!
      do m=1,nmmax(irow)
         grua (2*m-1,k) = grua (2*m-1,k) + gru1a (2*m  )
         gruha(2*m-1,k) = gruha(2*m-1,k) + gruh1a(2*m  )
         grua (2*m  ,k) = grua (2*m  ,k) - gru1a (2*m-1)
         gruha(2*m  ,k) = gruha(2*m  ,k) - gruh1a(2*m-1)
         grva (2*m-1,k) = grva (2*m-1,k) + grv1a (2*m  )
         grvha(2*m-1,k) = grvha(2*m-1,k) + grvh1a(2*m  )
         grva (2*m  ,k) = grva (2*m  ,k) - grv1a (2*m-1)
         grvha(2*m  ,k) = grvha(2*m  ,k) - grvh1a(2*m-1)
      end do
   end do
!
!-----------------------------------------------------------------------
!
! Computation for 1-level variables (ln(p*) and derivatives).
!
   do m=1,nmmax(irow)
      mr = nstart(m)
      mc = 2*mr
      do n=1,nlen(m),2
         ir = mc + 2*n - 1
         ii = ir + 1

         grpma(2*m-1) = grpma(2*m-1) + alps(ir)*dalp(mr+n,irow)*ra
         grpma(2*m  ) = grpma(2*m  ) + alps(ii)*dalp(mr+n,irow)*ra
      end do
   end do

   do m=1,nmmax(irow)
      mr = nstart(m)
      mc = 2*mr
      do n=2,nlen(m),2
         ir = mc + 2*n - 1
         ii = ir + 1
!
         tmpr = alps(ir)*alp(mr+n,irow)
         tmpi = alps(ii)*alp(mr+n,irow)
         grpsa(2*m-1) = grpsa(2*m-1) + tmpr
         grpsa(2*m  ) = grpsa(2*m  ) + tmpi
!
         grdpsa(2*m-1) = grdpsa(2*m-1) + tmpr*hdfst4(m+n-1)*ztodt
         grdpsa(2*m  ) = grdpsa(2*m  ) + tmpi*hdfst4(m+n-1)*ztodt
      end do
!
! Multiply by m/a to get d(ln(p*))/dlamda
! and by 1/a to get (1-mu**2)d(ln(p*))/dmu
!
      grpla(2*m-1) = -grpsa(2*m  )*ra*xm(m)
      grpla(2*m  ) =  grpsa(2*m-1)*ra*xm(m)
   end do
!
   return
end subroutine grcalca

#endif