#include <misc.h>#include <params.h>! Note that this routine has 2 complete blocks of code for PVP vs. non-PVP.! Make sure to make appropriate coding changes where necessary.#if ( defined PVP )subroutine grcalcs (irow    ,ztodt   ,grts    ,grths   ,grds    ,&                    grzs    ,grus    ,gruhs   ,grvs    ,grvhs   ,&                    grpss   ,grdpss  ,grpms   ,grpls   )!-----------------------------------------------------------------------!! 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, ez   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: symmetric fourier coefficients!   real(r8), intent(out) :: grts(plond,plev)    ! sum(n) of t(n,m)*P(n,m)   real(r8), intent(out) :: grths(plond,plev)   ! sum(n) of K(2i)*t(n,m)*P(n,m)   real(r8), intent(out) :: grds(plond,plev)    ! sum(n) of d(n,m)*P(n,m)   real(r8), intent(out) :: grzs(plond,plev)    ! sum(n) of z(n,m)*P(n,m)   real(r8), intent(out) :: grus(plond,plev)    ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))   real(r8), intent(out) :: gruhs(plond,plev)   ! sum(n) of K(2i)*z(n,m)*H(n,m)*a/(n(n+1))    real(r8), intent(out) :: grvs(plond,plev)    ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))   real(r8), intent(out) :: grvhs(plond,plev)   ! sum(n) of K(2i)*d(n,m)*H(n,m)*a/(n(n+1))   real(r8), intent(out) :: grpss(plond)        ! sum(n) of lnps(n,m)*P(n,m)   real(r8), intent(out) :: grdpss(plond)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)   real(r8), intent(out) :: grpms(plond)        ! sum(n) of lnps(n,m)*H(n,m)   real(r8), intent(out) :: grpls(plond)        ! sum(n) of lnps(n,m)*P(n,m)*m/a!!---------------------------Local workspace-----------------------------!   real(r8) gru1s(plond,plev)   ! sum(n) of d(n,m)*P(n,m)*m*a/(n(n+1))   real(r8) gruh1s(plond,plev)  ! sum(n) of K(2i)*d(n,m)*P(n,m)*m*a/(n(n+1))   real(r8) grv1s(plond,plev)   ! sum(n) of z(n,m)*P(n,m)*m*a/(n(n+1))   real(r8) grvh1s(plond,plev)  ! sum(n) of K(2i)*z(n,m)*P(n,m)*m*a/(n(n+1))   real(r8) zdfac(2*pnmax,plev) ! horiz. diffusion factor (vort,div) (complex)   real(r8) tqfac(2*pnmax,plev) ! horiz. diffusion factor (t,q) (complex)   real(r8) alp2(2*pspt)        ! Legendre functions (complex)   real(r8) dalp2(2*pspt)       ! derivative of Legendre functions (complex)   real(r8) alpn2(2*pspt)       ! (a*m/(n(n+1)))*Legendre functions (complex)   real(r8) dalpn2(2*pspt)      ! (a/(n(n+1)))*derivative of Legendre functions (complex)   real(r8) dlpnz               ! (a/(n(n+1)))*H(0,1) for conversion bet abs and rel vort   real(r8) zurcor              ! conversion term relating abs. & rel. vort.   integer k                ! level index   integer m                ! diagonal element(index) of spectral array   integer n                ! meridional wavenumber index   integer ne               ! index into spectral arrays   integer mn               ! index into spectral arrays   integer mnc              ! index into spectral arrays   integer mnev             ! index into spectral arrays!!-----------------------------------------------------------------------!! Compute alpn and dalpn! Expand polynomials and derivatives to complex form to allow largest ! possible vector length and multiply by appropriate factors!   do n=1,pmax      ne = n - 1!dir$ ivdep      do m=1,nmreduced(n,irow)         mnc = 2*(m+nalp(n))         mn = m + nalp(n)         alp2(mnc-1) = alp(mn,irow)         alp2(mnc  ) = alp(mn,irow)         dalp2(mnc-1) = dalp(mn,irow)*ra         dalp2(mnc  ) = dalp(mn,irow)*ra         alpn2(mnc-1) = alp(mn,irow)*(rsq(m+ne)*ra)*xm(m)         alpn2(mnc  ) = alp(mn,irow)*(rsq(m+ne)*ra)*xm(m)         dalpn2(mnc-1) = dalp(mn,irow)*(rsq(m+ne)*ra)         dalpn2(mnc  ) = dalp(mn,irow)*(rsq(m+ne)*ra)      end do   end do   dlpnz = dalpn2(2*nalp(2)+1)   zurcor = ez*dlpnz!! Initialize sums!   grzs(:,:)  = 0.   grds(:,:)  = 0.   gruhs(:,:) = 0.   grvhs(:,:) = 0.   grths(:,:) = 0.   grpss(:)   = 0.   grus(:,:)  = 0.   grvs(:,:)  = 0.   grts(:,:)  = 0.   grpls(:)   = 0.   grpms(:)   = 0.   grdpss(:)   = 0.   do k=1,plev!! Diffusion factors: expand for longest possible vectors!!dir$ ivdep      do n=1,pnmax         zdfac(n*2-1,k) = -hdifzd(n,k)         zdfac(n*2  ,k) = -hdifzd(n,k)         tqfac(n*2-1,k) = -hdiftq(n,k)         tqfac(n*2  ,k) = -hdiftq(n,k)      end do      gru1s(:,k) = 0.      gruh1s(:,k) = 0.      grv1s(:,k) = 0.      grvh1s(:,k) = 0.!! Evaluate symmetric components involving P and antisymmetric involving ! H. Loop over n for t(m), q(m), d(m),and the two parts of u(m) and v(m).! The inner (vector) loop accumulates sums over n along the diagonals! of the spectral truncation to obtain the maximum length vectors.!! "ncutoff" is used to switch to vectorization in the vertical when the ! length of the diagonal is less than the number of levels.!      do n=1,ncutoff,2         ne = 2*(n-1)         do m=1,2*nmreduced(n,irow)            mnev = m + nco2(n) - 2            mn = m + 2*nalp(n)            grts (m,k) = grts (m,k) + t(mnev,k)*alp2(mn)            grths(m,k) = grths(m,k) + t(mnev,k)*alp2(mn)*tqfac(m+ne,k)            grds(m,k) = grds(m,k) + d(mnev,k)*alp2(mn)            grzs(m,k) = grzs(m,k) + vz(mnev,k)*alp2(mn)            gru1s (m,k) = gru1s (m,k) + d(mnev,k)*alpn2(mn)            gruh1s(m,k) = gruh1s(m,k) + d(mnev,k)*alpn2(mn)*zdfac(m+ne,k)            grv1s (m,k) = grv1s (m,k) + vz(mnev,k)*alpn2(mn)            grvh1s(m,k) = grvh1s(m,k) + vz(mnev,k)*alpn2(mn)*zdfac(m+ne,k)         end do      end do!! Evaluate antisymmetric components involving P and symmetric involving ! H. Loop over n for t(m), q(m), d(m),and the two parts of u(m) and v(m).! The inner (vector) loop accumulates sums over n along the diagonals! of the spectral truncation to obtain the maximum length vectors.!! "ncutoff" is used to switch to vectorization in the vertical when the ! length of the diagonal is less than the number of levels.!      do n=2,ncutoff,2         ne = 2*(n-1)         do m=1,2*nmreduced(n,irow)            mnev = m + nco2(n) - 2            mn = m + 2*nalp(n)            grus (m,k) = grus (m,k) + vz(mnev,k)*dalpn2(mn)            gruhs(m,k) = gruhs(m,k) + vz(mnev,k)*dalpn2(mn)*zdfac(m+ne,k)            grvs (m,k) = grvs (m,k) - d(mnev,k)*dalpn2(mn)            grvhs(m,k) = grvhs(m,k) - d(mnev,k)*dalpn2(mn)*zdfac(m+ne,k)         end do      end do   end do!! For short diagonals, repeat above loops with vectorization in vertical,! instead of along diagonals, to keep vector lengths from getting too short.!   if (ncutoff.lt.pmax) then      do n=ncutoff+1,pmax,2   ! ncutoff guaranteed even         ne = 2*(n-1)         do m=1,2*nmreduced(n,irow)            mnev = m + nco2(n) - 2            mn = m + 2*nalp(n)            do k=1,plev               grts (m,k) = grts (m,k) + t(mnev,k)*alp2(mn)               grths(m,k) = grths(m,k) + t(mnev,k)*alp2(mn)*tqfac(m+ne,k)               grds(m,k) = grds(m,k) + d(mnev,k)*alp2(mn)               grzs(m,k) = grzs(m,k) + vz(mnev,k)*alp2(mn)               gru1s (m,k) = gru1s (m,k) + d(mnev,k)*alpn2(mn)               gruh1s(m,k) = gruh1s(m,k) + d(mnev,k)*alpn2(mn)*zdfac(m+ne,k)               grv1s (m,k) = grv1s (m,k) + vz(mnev,k)*alpn2(mn)               grvh1s(m,k) = grvh1s(m,k) + vz(mnev,k)*alpn2(mn)*zdfac(m+ne,k)            end do         end do      end do      do n=ncutoff+2,pmax,2         ne = 2*(n-1)         do m=1,2*nmreduced(n,irow)            mnev = m + nco2(n) - 2            mn = m + 2*nalp(n)            do k=1,plev               grus (m,k) = grus (m,k) + vz(mnev,k)*dalpn2(mn)               gruhs(m,k) = gruhs(m,k) + vz(mnev,k)*dalpn2(mn)*zdfac(m+ne,k)               grvs (m,k) = grvs (m,k) - d(mnev,k)*dalpn2(mn)               grvhs(m,k) = grvhs(m,k) - d(mnev,k)*dalpn2(mn)*zdfac(m+ne,k)            end do         end do      end do   end if                      ! ncutoff.lt.pmax   do k=1,plev!! Combine the two parts of u(m) and v(m)!!dir$ ivdep      do m=1,nmmax(irow)         grus (2*m-1,k) = grus (2*m-1,k) + gru1s (2*m  ,k)         gruhs(2*m-1,k) = gruhs(2*m-1,k) + gruh1s(2*m  ,k)         grus (2*m  ,k) = grus (2*m  ,k) - gru1s (2*m-1,k)         gruhs(2*m  ,k) = gruhs(2*m  ,k) - gruh1s(2*m-1,k)         grvs (2*m-1,k) = grvs (2*m-1,k) + grv1s (2*m  ,k)         grvhs(2*m-1,k) = grvhs(2*m-1,k) + grvh1s(2*m  ,k)         grvs (2*m  ,k) = grvs (2*m  ,k) - grv1s (2*m-1,k)         grvhs(2*m  ,k) = grvhs(2*m  ,k) - grvh1s(2*m-1,k)      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 single level variables.!! Evaluate symmetric components involving P and antisymmetric involving ! H.  Loop over n for lnps(m) and derivatives.! The inner loop accumulates over n along diagonal of the truncation.!   do n=1,pmax,2      ne = n - 1      do m=1,2*nmreduced(n,irow)         mnev = m + nco2(n) - 2         mn = m + 2*nalp(n)         grpss (m) = grpss (m) + alps(mnev)*alp2(mn)         grdpss(m) = grdpss(m) + alps(mnev)*alp2(mn)*hdfst4(ne+(m+1)/2)*ztodt      end do   end do!! Evaluate antisymmetric components involving P and symmetric involving ! H.  Loop over n for lnps(m) and derivatives.! The inner loop accumulates over n along diagonal of the truncation.!   do n=2,pmax,2      ne = n - 1      do m=1,2*nmreduced(n,irow)         mnev = m + nco2(n) - 2         mn = m + 2*nalp(n)         grpms(m) = grpms(m) + alps(mnev)*dalp2(mn)      end do   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!   do m=1,nmmax(irow)      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, ez   implicit none#include <comhd.h>!! Input arguments!