#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    ,grqs    ,grths   , &                    grds    ,grus    ,gruhs   ,grvs    ,grvhs   , &                    grpss   ,grdpss  ,grpms   ,grpls   ,grtms   , &                    grtls   ,grqms   ,grqls   )!-----------------------------------------------------------------------!! Purpose:! 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.!! Original version:  CCM1!!-----------------------------------------------------------------------!! $Id: grcalc.F90,v 1.5 2001/04/13 22:40:37 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: symmetric fourier coefficients!   real(r8), intent(out) :: grts(plond,plev)  ! sum(n) of t(n,m)*P(n,m)   real(r8), intent(out) :: grqs(plond,plev)  ! sum(n) of q(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) :: 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   real(r8), intent(out) :: grtms (plond,plev)   real(r8), intent(out) :: grtls (plond,plev)   real(r8), intent(out) :: grqms (plond,plev)   real(r8), intent(out) :: grqls (plond,plev)!!---------------------------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)!                                ! betwn absolute and relative vorticity   integer k                     ! level index   integer m                     ! diagonal element(index) of spec. 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!! Initialize sums!   grts(:,:) = 0.   grqs(:,:) = 0.   grths(:,:) = 0.   grds(:,:)  = 0.   grus(:,:)  = 0.   gruhs(:,:) = 0.   grvs(:,:)  = 0.   grvhs(:,:) = 0.   grpss(:)   = 0.   grdpss(:)   = 0.   grpms(:)   = 0.   grpls(:)   = 0.   grtms(:,:)   = 0.   grtls(:,:)   = 0.   grqms(:,:)   = 0.   grqls(:,:)   = 0.!!-----------------------------------------------------------------------!! Computation for multilevel variables!   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!! Initialize local sums!      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)            grqs  (m,k) = grqs  (m,k) + q (mnev,k)*alp2  (mn)            grds  (m,k) = grds  (m,k) + d (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)            grtms (m,k) = grtms (m,k) + t (mnev,k)*dalp2 (mn)            grqms (m,k) = grqms (m,k) + q (mnev,k)*dalp2 (mn)            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                    ! k=1,plev!! 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)               grqs  (m,k) = grqs  (m,k) + q (mnev,k)*alp2  (mn)               grds  (m,k) = grds  (m,k) + d (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               grtms (m,k) = grtms (m,k) + t (mnev,k)*dalp2 (mn)               grqms (m,k) = grqms (m,k) + q (mnev,k)*dalp2 (mn)               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) and compute derivatives!!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)!! Derivatives!         grtls(2*m-1,k) = -grts(2*m  ,k)*ra*xm(m)         grtls(2*m  ,k) =  grts(2*m-1,k)*ra*xm(m)         grqls(2*m-1,k) = -grqs(2*m  ,k)*ra*xm(m)         grqls(2*m  ,k) =  grqs(2*m-1,k)*ra*xm(m)      end do   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    ,grqa    ,grtha   , &                    grda    ,grua    ,gruha   ,grva    ,grvha   , &                    grpsa   ,grdpsa  ,grpma   ,grpla   ,grtma   , &                    grtla   ,grqma   ,grqla   )!-----------------------------------------------------------------------!! Purpose:! 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.!! Original version:  CCM1!!-----------------------------------------------------------------------!! $Id: grcalc.F90,v 1.5 2001/04/13 22:40:37 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: anti-symmetric fourier coefficients!   real(r8), intent(out) :: grta(plond,plev)  ! sum(n) of t(n,m)*P(n,m)   real(r8), intent(out) :: grqa(plond,plev)  ! sum(n) of q(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) :: grua(plond,plev)  ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))