#include <misc.h>
#include <params.h>
subroutine grcalcs (irow ,ztodt ,grts ,grqs ,grths , & 2,8
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.8.4 2004/03/03 19:54:10 pworley Exp $
! $Author: pworley $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
use pspect
use comspe
use rgrid
use commap
use dynconst, only: ra
use comhd
implicit none
!
! 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(2*maxm,plev) ! sum(n) of t(n,m)*P(n,m)
real(r8), intent(out) :: grqs(2*maxm,plev) ! sum(n) of q(n,m)*P(n,m)
real(r8), intent(out) :: grths(2*maxm,plev) ! sum(n) of K(2i)*t(n,m)*P(n,m)
real(r8), intent(out) :: grds(2*maxm,plev) ! sum(n) of d(n,m)*P(n,m)
real(r8), intent(out) :: grus(2*maxm,plev) ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: gruhs(2*maxm,plev) ! sum(n) of K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: grvs(2*maxm,plev) ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: grvhs(2*maxm,plev) ! sum(n) of K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: grpss(2*maxm) ! sum(n) of lnps(n,m)*P(n,m)
real(r8), intent(out) :: grdpss(2*maxm) ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
real(r8), intent(out) :: grpms(2*maxm) ! sum(n) of lnps(n,m)*H(n,m)
real(r8), intent(out) :: grpls(2*maxm) ! sum(n) of lnps(n,m)*P(n,m)*m/a
real(r8), intent(out) :: grtms (2*maxm,plev)
real(r8), intent(out) :: grtls (2*maxm,plev)
real(r8), intent(out) :: grqms (2*maxm,plev)
real(r8), intent(out) :: grqls (2*maxm,plev)
!
!---------------------------Local workspace-----------------------------
!
real(r8) gru1s (2*maxm) ! sum(n) of d(n,m)*P(n,m)*m*a/(n(n+1))
real(r8) gruh1s(2*maxm) ! sum(n) of K(2i)*d(n,m)*P(n,m)*m*a/(n(n+1))
real(r8) grv1s (2*maxm) ! sum(n) of z(n,m)*P(n,m)*m*a/(n(n+1))
real(r8) grvh1s(2*maxm) ! 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 lm, m ! local and global Fourier wavenumber indices of spectral array
integer mlength ! number of local wavenumbers
integer n ! meridional wavenumber index
integer ir,ii ! spectral indices
integer mr,mc ! spectral indices
real(r8) tmp,raxm ! temporary workspace
!
!-----------------------------------------------------------------------
!
! Compute alpn and dalpn
!
mlength = numm(iam)
do lm=1,mlength
m = locm(lm,iam)
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
!
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
!
! Initialize local sums
!
gru1s(:) = 0.
gruh1s(:) = 0.
grv1s(:) = 0.
grvh1s(:) = 0.
!
! Loop over n for t,q,d,and end of u and v
!
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
grts (2*lm-1,k) = grts (2*lm-1,k) + t(ir,k)*alp(mr+n,irow)
grts (2*lm ,k) = grts (2*lm ,k) + t(ii,k)*alp(mr+n,irow)
!
grqs (2*lm-1,k) = grqs (2*lm-1,k) + q(ir,k)*alp(mr+n,irow)
grqs (2*lm ,k) = grqs (2*lm ,k) + q(ii,k)*alp(mr+n,irow)
!
tmp = alp(mr+n,irow)*hdiftq(n+m-1,k)
grths(2*lm-1,k) = grths(2*lm-1,k) - t(ir,k)*tmp
grths(2*lm ,k) = grths(2*lm ,k) - t(ii,k)*tmp
!
grds(2*lm-1,k) = grds(2*lm-1,k) + d(ir,k)*alp(mr+n,irow)
grds(2*lm ,k) = grds(2*lm ,k) + d(ii,k)*alp(mr+n,irow)
!
gru1s (2*lm-1) = gru1s (2*lm-1) + d(ir,k)*alpn(mr+n)
gru1s (2*lm ) = gru1s (2*lm ) + d(ii,k)*alpn(mr+n)
!
tmp = alpn(mr+n)*hdifzd(n+m-1,k)
gruh1s(2*lm-1) = gruh1s(2*lm-1) - d(ir,k)*tmp
gruh1s(2*lm ) = gruh1s(2*lm ) - d(ii,k)*tmp
!
grv1s (2*lm-1) = grv1s (2*lm-1) + vz(ir,k)*alpn(mr+n)
grv1s (2*lm ) = grv1s (2*lm ) + vz(ii,k)*alpn(mr+n)
!
grvh1s(2*lm-1) = grvh1s(2*lm-1) - vz(ir,k)*tmp
grvh1s(2*lm ) = grvh1s(2*lm ) - vz(ii,k)*tmp
end do
end do
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
grtms(2*lm-1,k) = grtms(2*lm-1,k) + t(ir,k)*dalp(mr+n,irow)*ra
grtms(2*lm ,k) = grtms(2*lm ,k) + t(ii,k)*dalp(mr+n,irow)*ra
!
grqms(2*lm-1,k) = grqms(2*lm-1,k) + q(ir,k)*dalp(mr+n,irow)*ra
grqms(2*lm ,k) = grqms(2*lm ,k) + q(ii,k)*dalp(mr+n,irow)*ra
!
grus (2*lm-1,k) = grus (2*lm-1,k) + vz(ir,k)*dalpn(mr+n)
grus (2*lm ,k) = grus (2*lm ,k) + vz(ii,k)*dalpn(mr+n)
!
tmp = dalpn(mr+n)*hdifzd(n+m-1,k)
gruhs(2*lm-1,k) = gruhs(2*lm-1,k) - vz(ir,k)*tmp
gruhs(2*lm ,k) = gruhs(2*lm ,k) - vz(ii,k)*tmp
!
grvs (2*lm-1,k) = grvs (2*lm-1,k) - d(ir,k)*dalpn(mr+n)
grvs (2*lm ,k) = grvs (2*lm ,k) - d(ii,k)*dalpn(mr+n)
!
grvhs(2*lm-1,k) = grvhs(2*lm-1,k) + d(ir,k)*tmp
grvhs(2*lm ,k) = grvhs(2*lm ,k) + d(ii,k)*tmp
end do
end do
!
! Combine the two parts of u(m) and v(m)
!
do lm=1,mlength
m = locm(lm,iam)
grus (2*lm-1,k) = grus (2*lm-1,k) + gru1s (2*lm )
gruhs(2*lm-1,k) = gruhs(2*lm-1,k) + gruh1s(2*lm )
grus (2*lm ,k) = grus (2*lm ,k) - gru1s (2*lm-1)
gruhs(2*lm ,k) = gruhs(2*lm ,k) - gruh1s(2*lm-1)
grvs (2*lm-1,k) = grvs (2*lm-1,k) + grv1s (2*lm )
grvhs(2*lm-1,k) = grvhs(2*lm-1,k) + grvh1s(2*lm )
grvs (2*lm ,k) = grvs (2*lm ,k) - grv1s (2*lm-1)
grvhs(2*lm ,k) = grvhs(2*lm ,k) - grvh1s(2*lm-1)
!
! Derivatives
!
grtls(2*lm-1,k) = -grts(2*lm ,k)*ra*xm(m)
grtls(2*lm ,k) = grts(2*lm-1,k)*ra*xm(m)
grqls(2*lm-1,k) = -grqs(2*lm ,k)*ra*xm(m)
grqls(2*lm ,k) = grqs(2*lm-1,k)*ra*xm(m)
end do
end do
!
!-----------------------------------------------------------------------
!
! Computation for 1-level variables (ln(p*) and derivatives).
!
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
grpss(2*lm-1) = grpss(2*lm-1) + alps(ir)*alp(mr+n,irow)
grpss(2*lm ) = grpss(2*lm ) + alps(ii)*alp(mr+n,irow)
!
grdpss(2*lm-1) = grdpss(2*lm-1) + alps(ir)*alp(mr+n,irow)*hdfst4(m+n-1)*ztodt
grdpss(2*lm ) = grdpss(2*lm ) + alps(ii)*alp(mr+n,irow)*hdfst4(m+n-1)*ztodt
!
end do
end do
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
grpms(2*lm-1) = grpms(2*lm-1) + alps(ir)*dalp(mr+n,irow)*ra
grpms(2*lm ) = grpms(2*lm ) + 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*lm-1) = -grpss(2*lm )*ra*xm(m)
grpls(2*lm ) = grpss(2*lm-1)*ra*xm(m)
end do
!
return
end subroutine grcalcs
subroutine grcalca (irow ,ztodt ,grta ,grqa ,grtha , & 2,8
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.8.4 2004/03/03 19:54:10 pworley Exp $
! $Author: pworley $
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
use pspect
use comspe
use rgrid
use commap
use dynconst, only: ra
use comhd
implicit none
!
! 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(2*maxm,plev) ! sum(n) of t(n,m)*P(n,m)
real(r8), intent(out) :: grqa(2*maxm,plev) ! sum(n) of q(n,m)*P(n,m)
real(r8), intent(out) :: grtha(2*maxm,plev) ! sum(n) of K(2i)*t(n,m)*P(n,m)
real(r8), intent(out) :: grda(2*maxm,plev) ! sum(n) of d(n,m)*P(n,m)
real(r8), intent(out) :: grua(2*maxm,plev) ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: gruha(2*maxm,plev) ! sum(n) of K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: grva(2*maxm,plev) ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: grvha(2*maxm,plev) ! sum(n) of K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
real(r8), intent(out) :: grpsa(2*maxm) ! sum(n) of lnps(n,m)*P(n,m)
real(r8), intent(out) :: grdpsa(2*maxm) ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt
! ! *lnps(n,m)*P(n,m)
real(r8), intent(out) :: grpma(2*maxm) ! sum(n) of lnps(n,m)*H(n,m)
real(r8), intent(out) :: grpla(2*maxm) ! sum(n) of lnps(n,m)*P(n,m)*m/a
real(r8), intent(out) :: grtma (2*maxm,plev)
real(r8), intent(out) :: grtla (2*maxm,plev)
real(r8), intent(out) :: grqma (2*maxm,plev)
real(r8), intent(out) :: grqla (2*maxm,plev)
!
!---------------------------Local workspace-----------------------------
!
real(r8) gru1a (2*maxm) ! sum(n) of d(n,m)*P(n,m)*m*a/(n(n+1))
real(r8) gruh1a(2*maxm) ! sum(n) of K(2i)*d(n,m)*P(n,m)*m*a/(n(n+1))
real(r8) grv1a (2*maxm) ! sum(n) of z(n,m)*P(n,m)*m*a/(n(n+1))
real(r8) grvh1a(2*maxm) ! 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 lm, m ! local and global Fourier wavenumber indices of spectral array
integer mlength ! number of local wavenumbers
integer n ! meridional wavenumber index
integer ir,ii ! spectral indices
integer mr,mc ! spectral indices
real(r8) tmp,raxm ! temporary workspace
!
!-----------------------------------------------------------------------
!
! Compute alpn and dalpn
!
mlength = numm(iam)
do lm=1,mlength
m = locm(lm,iam)
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
!
grta(:,:) = 0.
grqa(:,:) = 0.
grtha(:,:) = 0.
grda(:,:) = 0.
grua(:,:) = 0.
gruha(:,:) = 0.
grva(:,:) = 0.
grvha(:,:) = 0.
grpsa(:) = 0.
grdpsa(:) = 0.
grpma(:) = 0.
grpla(:) = 0.
grtma(:,:) = 0.
grtla(:,:) = 0.
grqma(:,:) = 0.
grqla(:,:) = 0.
!
!-----------------------------------------------------------------------
!
! Computation for multilevel variables
!
do k=1,plev
!
! Initialize local sums
!
gru1a(:) = 0.
gruh1a(:) = 0.
grv1a(:) = 0.
grvh1a(:) = 0.
!
! Loop over n for t,q,d,and end of u and v
!
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
grtma(2*lm-1,k) = grtma(2*lm-1,k) + t(ir,k)*dalp(mr+n,irow)*ra
grtma(2*lm ,k) = grtma(2*lm ,k) + t(ii,k)*dalp(mr+n,irow)*ra
!
grqma(2*lm-1,k) = grqma(2*lm-1,k) + q(ir,k)*dalp(mr+n,irow)*ra
grqma(2*lm ,k) = grqma(2*lm ,k) + q(ii,k)*dalp(mr+n,irow)*ra
!
grua (2*lm-1,k) = grua (2*lm-1,k) + vz(ir,k)*dalpn(mr+n)
grua (2*lm ,k) = grua (2*lm ,k) + vz(ii,k)*dalpn(mr+n)
!
tmp = dalpn(mr+n)*hdifzd(n+m-1,k)
gruha(2*lm-1,k) = gruha(2*lm-1,k) - vz(ir,k)*tmp
gruha(2*lm ,k) = gruha(2*lm ,k) - vz(ii,k)*tmp
!
grva (2*lm-1,k) = grva (2*lm-1,k) - d(ir,k)*dalpn(mr+n)
grva (2*lm ,k) = grva (2*lm ,k) - d(ii,k)*dalpn(mr+n)
!
grvha(2*lm-1,k) = grvha(2*lm-1,k) + d(ir,k)*tmp
grvha(2*lm ,k) = grvha(2*lm ,k) + d(ii,k)*tmp
end do
end do
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
grta (2*lm-1,k) = grta (2*lm-1,k) + t(ir,k)*alp(mr+n,irow)
grta (2*lm ,k) = grta (2*lm ,k) + t(ii,k)*alp(mr+n,irow)
!
grqa (2*lm-1,k) = grqa (2*lm-1,k) + q(ir,k)*alp(mr+n,irow)
grqa (2*lm ,k) = grqa (2*lm ,k) + q(ii,k)*alp(mr+n,irow)
!
tmp = alp(mr+n,irow)*hdiftq(n+m-1,k)
grtha(2*lm-1,k) = grtha(2*lm-1,k) - t(ir,k)*tmp
grtha(2*lm ,k) = grtha(2*lm ,k) - t(ii,k)*tmp
!
grda(2*lm-1,k) = grda(2*lm-1,k) + d(ir,k)*alp(mr+n,irow)
grda(2*lm ,k) = grda(2*lm ,k) + d(ii,k)*alp(mr+n,irow)
!
gru1a (2*lm-1) = gru1a (2*lm-1) + d(ir,k)*alpn(mr+n)
gru1a (2*lm ) = gru1a (2*lm ) + d(ii,k)*alpn(mr+n)
!
tmp = alpn(mr+n)*hdifzd(n+m-1,k)
gruh1a(2*lm-1) = gruh1a(2*lm-1) - d(ir,k)*tmp
gruh1a(2*lm ) = gruh1a(2*lm ) - d(ii,k)*tmp
!
grv1a (2*lm-1) = grv1a (2*lm-1) + vz(ir,k)*alpn(mr+n)
grv1a (2*lm ) = grv1a (2*lm ) + vz(ii,k)*alpn(mr+n)
!
grvh1a(2*lm-1) = grvh1a(2*lm-1) - vz(ir,k)*tmp
grvh1a(2*lm ) = grvh1a(2*lm ) - vz(ii,k)*tmp
end do
end do
!
! Combine the two parts of u(m) and v(m)
!
do lm=1,mlength
m = locm(lm,iam)
grua (2*lm-1,k) = grua (2*lm-1,k) + gru1a (2*lm )
gruha(2*lm-1,k) = gruha(2*lm-1,k) + gruh1a(2*lm )
grua (2*lm ,k) = grua (2*lm ,k) - gru1a (2*lm-1)
gruha(2*lm ,k) = gruha(2*lm ,k) - gruh1a(2*lm-1)
grva (2*lm-1,k) = grva (2*lm-1,k) + grv1a (2*lm )
grvha(2*lm-1,k) = grvha(2*lm-1,k) + grvh1a(2*lm )
grva (2*lm ,k) = grva (2*lm ,k) - grv1a (2*lm-1)
grvha(2*lm ,k) = grvha(2*lm ,k) - grvh1a(2*lm-1)
!
! Derivatives
!
grtla(2*lm-1,k) = -grta(2*lm ,k)*ra*xm(m)
grtla(2*lm ,k) = grta(2*lm-1,k)*ra*xm(m)
grqla(2*lm-1,k) = -grqa(2*lm ,k)*ra*xm(m)
grqla(2*lm ,k) = grqa(2*lm-1,k)*ra*xm(m)
end do
end do
!
!-----------------------------------------------------------------------
!
! Computation for 1-level variables (ln(p*) and derivatives).
!
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=1,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
grpma(2*lm-1) = grpma(2*lm-1) + alps(ir)*dalp(mr+n,irow)*ra
grpma(2*lm ) = grpma(2*lm ) + alps(ii)*dalp(mr+n,irow)*ra
end do
end do
do lm=1,mlength
m = locm(lm,iam)
mr = nstart(m)
mc = 2*mr
do n=2,nlen(m),2
ir = mc + 2*n - 1
ii = ir + 1
!
grpsa(2*lm-1) = grpsa(2*lm-1) + alps(ir)*alp(mr+n,irow)
grpsa(2*lm ) = grpsa(2*lm ) + alps(ii)*alp(mr+n,irow)
!
grdpsa(2*lm-1) = grdpsa(2*lm-1) + alps(ir)*alp(mr+n,irow)*hdfst4(m+n-1)*ztodt
grdpsa(2*lm ) = grdpsa(2*lm ) + alps(ii)*alp(mr+n,irow)*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*lm-1) = -grpsa(2*lm )*ra*xm(m)
grpla(2*lm ) = grpsa(2*lm-1)*ra*xm(m)
end do
!
return
end subroutine grcalca