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!   returnend subroutine grcalcssubroutine 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!   returnend subroutine grcalca#endif