#include <misc.h>
#include <params.h>
subroutine grcalcs (irow ,ztodt ,grts ,grths ,grds ,& 2,8
grzs ,grus ,gruhs ,grvs ,grvhs ,&
grpss ,grdpss ,grpms ,grpls ,tmpSPEcoef)
!-----------------------------------------------------------------------
!
! 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
! Modified: P. Worley, October 2002
!
!-----------------------------------------------------------------------
!
! $Id: grcalc.F90,v 1.5.6.9 2004/04/28 23:43:54 eaton Exp $
! $Author: eaton $
!
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
use pspect
use comspe
use rgrid
use commap
use dynconst, only: ra, ez
use comhd
implicit none
!
! Input arguments
!
integer, intent(in) :: irow ! latitude pair index
real(r8), intent(in) :: ztodt ! twice the timestep unless nstep = 0
real(r8), intent(in) :: tmpSPEcoef(plev*24,pnmax,maxm) ! rearranged variables array
!
! 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) :: 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) :: grzs(2*maxm,plev) ! sum(n) of z(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
!
!---------------------------Local workspace-----------------------------
!
real(r8) dalpn(pspt) ! (a/(n(n+1)))*derivative of Legendre functions (complex)
real(r8) zurcor ! conversion term relating abs. & rel. vort.
real(r8) tmpGRcoef(plev*24,maxm) ! temporal storage for Fourier coeffs
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 lmr,lmc ! spectral indices
integer lmwave0 ! local index for wavenumber 0
integer lmrwave0 ! local offset for wavenumber 0
integer kv ! level x variable index
!
!-----------------------------------------------------------------------
!
! Compute alpn and dalpn
!
!DIR$ NOSTREAM
lmwave0 = -1
lmrwave0 = 0
dalpn(2) = 0.0
mlength = numm(iam)
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
if (m .eq. 1) then
lmwave0 = lm
lmrwave0 = lmr
endif
do n=1,nlen(m)
dalpn(lmr+n) = ldalp(lmr+n,irow)*rsq(m+n-1)*ra
end do
end do
zurcor = ez*dalpn(lmrwave0 + 2)
!
! Initialize sums
!
grpss (:) = 0.
grpls (:) = 0.
grpms (:) = 0.
grdpss(:) = 0.
!cdir collapse
tmpGRcoef (:,:) = 0.
!
! Loop over n for t,q,d,and end of u and v
!
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
do n=2,nlen(m),2
do kv=1,plev*8
tmpGRcoef(kv,lm) = tmpGRcoef(kv,lm) + tmpSPEcoef(kv,n,lm)*dalpn(lmr+n)
end do
end do
end do
!
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
do n=1,nlen(m),2
do kv=plev*8+1,plev*24
tmpGRcoef(kv,lm) = tmpGRcoef(kv,lm) + tmpSPEcoef(kv,n,lm)*lalp(lmr+n,irow)
end do
end do
end do
!
! Combine the two parts of u(m) and v(m)
!
do lm=1,mlength
do kv=1,plev*8
tmpGRcoef(kv,lm) = tmpGRcoef(kv,lm) + tmpGRcoef(kv+plev*16,lm)
end do
end do
!
! Save accumulated results to gr* arrays
!
do lm=1,mlength
!DIR$ PREFERVECTOR
do k=1,plev
grus (2*lm-1,k) = tmpGRcoef(k ,lm)
grus (2*lm ,k) = tmpGRcoef(k+plev ,lm)
grvs (2*lm-1,k) = tmpGRcoef(k+plev*2 ,lm)
grvs (2*lm ,k) = tmpGRcoef(k+plev*3 ,lm)
gruhs(2*lm-1,k) = tmpGRcoef(k+plev*4 ,lm)
gruhs(2*lm ,k) = tmpGRcoef(k+plev*5 ,lm)
grvhs(2*lm-1,k) = tmpGRcoef(k+plev*6 ,lm)
grvhs(2*lm ,k) = tmpGRcoef(k+plev*7 ,lm)
grts (2*lm-1,k) = tmpGRcoef(k+plev*8 ,lm)
grts (2*lm ,k) = tmpGRcoef(k+plev*9 ,lm)
grths(2*lm-1,k) = tmpGRcoef(k+plev*10,lm)
grths(2*lm ,k) = tmpGRcoef(k+plev*11,lm)
grds (2*lm-1,k) = tmpGRcoef(k+plev*12,lm)
grds (2*lm ,k) = tmpGRcoef(k+plev*13,lm)
grzs (2*lm-1,k) = tmpGRcoef(k+plev*14,lm)
grzs (2*lm ,k) = tmpGRcoef(k+plev*15,lm)
end do
end do
!
! Remove Coriolis contribution to absolute vorticity from u(m)
! Correction for u:zeta=vz-ez=(zeta+f)-f
!
if (lmwave0 .ne. -1) then
do k=1,plev
! grus(1,k) = grus(1,k) - zurcor
grus(2*lmwave0-1,k) = grus(2*lmwave0-1,k) - zurcor
end do
endif
!
!-----------------------------------------------------------------------
!
! Computation for 1-level variables (ln(p*) and derivatives).
!
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
lmc = 2*lmr
!cdir shortloop
do n=1,nlen(m),2
ir = lmc + 2*n - 1
ii = ir + 1
!
grpss (2*lm-1) = grpss (2*lm-1) + alps(ir)*lalp(lmr+n,irow)
grpss (2*lm ) = grpss (2*lm ) + alps(ii)*lalp(lmr+n,irow)
!
grdpss(2*lm-1) = grdpss(2*lm-1) + alps(ir)*lalp(lmr+n,irow)*hdfst4(m+n-1)*ztodt
grdpss(2*lm ) = grdpss(2*lm ) + alps(ii)*lalp(lmr+n,irow)*hdfst4(m+n-1)*ztodt
end do
end do
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
lmc = 2*lmr
!cdir shortloop
do n=2,nlen(m),2
ir = lmc + 2*n - 1
ii = ir + 1
!
grpms(2*lm-1) = grpms(2*lm-1) + alps(ir)*ldalp(lmr+n,irow)*ra
grpms(2*lm ) = grpms(2*lm ) + alps(ii)*ldalp(lmr+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 ,grtha ,grda ,& 2,8
grza ,grua ,gruha ,grva ,grvha ,&
grpsa ,grdpsa ,grpma ,grpla ,tmpSPEcoef)
!-----------------------------------------------------------------------
!
! 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
! Modified: P. Worley, October 2002
!
!-----------------------------------------------------------------------
!
! $Id: grcalc.F90,v 1.5.6.9 2004/04/28 23:43:54 eaton Exp $
! $Author: eaton $
!
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
real(r8), intent(in) :: tmpSPEcoef(plev*24,pnmax,maxm) ! array for rearranged variables
!
!
! Output arguments: antisymmetric fourier coefficients
!
real(r8), intent(out) :: grta(2*maxm,plev) ! sum(n) of t(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) :: grza(2*maxm,plev) ! sum(n) of z(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
!
!---------------------------Local workspace-----------------------------
!
real(r8) dalpn(pspt) ! (a/(n(n+1)))*derivative of Legendre functions (complex)
real(r8) tmpGRcoef(plev*24,maxm) ! temporal storage for Fourier coefficients
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 lmr,lmc ! spectral indices
integer kv ! level x variable index
!
!-----------------------------------------------------------------------
!
! Compute alpn and dalpn
!
!DIR$ NOSTREAM
mlength = numm(iam)
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
do n=1,nlen(m)
dalpn(lmr+n) = ldalp(lmr+n,irow)*rsq(m+n-1)*ra
end do
end do
!
! Initialize sums
!
grpsa (:) = 0.
grpla (:) = 0.
grpma (:) = 0.
grdpsa(:) = 0.
!cdir collapse
tmpGRcoef(:,:) = 0.
!
! Loop over n for t,q,d,and end of u and v
!
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
do n=1,nlen(m),2
do kv=1,plev*8
tmpGRcoef(kv,lm) = tmpGRcoef(kv,lm) + tmpSPEcoef(kv,n,lm)*dalpn(lmr+n)
end do
end do
end do
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
do n=2,nlen(m),2
do kv=plev*8+1,plev*24
tmpGRcoef(kv,lm) = tmpGRcoef(kv,lm) + tmpSPEcoef(kv,n,lm)*lalp(lmr+n,irow)
end do
end do
end do
!
! Combine the two parts of u(m) and v(m)
!
do lm=1,mlength
do kv=1,plev*8
tmpGRcoef(kv,lm) = tmpGRcoef(kv,lm) + tmpGRcoef(kv+plev*16,lm)
end do
end do
!
! Save accumulated results to gr* arrays
!
do lm=1,mlength
!DIR$ PREFERVECTOR
do k=1,plev
grua (2*lm-1,k) = tmpGRcoef(k ,lm)
grua (2*lm ,k) = tmpGRcoef(k+plev ,lm)
grva (2*lm-1,k) = tmpGRcoef(k+plev*2 ,lm)
grva (2*lm ,k) = tmpGRcoef(k+plev*3 ,lm)
gruha(2*lm-1,k) = tmpGRcoef(k+plev*4 ,lm)
gruha(2*lm ,k) = tmpGRcoef(k+plev*5 ,lm)
grvha(2*lm-1,k) = tmpGRcoef(k+plev*6 ,lm)
grvha(2*lm ,k) = tmpGRcoef(k+plev*7 ,lm)
grta (2*lm-1,k) = tmpGRcoef(k+plev*8 ,lm)
grta (2*lm ,k) = tmpGRcoef(k+plev*9 ,lm)
grtha(2*lm-1,k) = tmpGRcoef(k+plev*10,lm)
grtha(2*lm ,k) = tmpGRcoef(k+plev*11,lm)
grda (2*lm-1,k) = tmpGRcoef(k+plev*12,lm)
grda (2*lm ,k) = tmpGRcoef(k+plev*13,lm)
grza (2*lm-1,k) = tmpGRcoef(k+plev*14,lm)
grza (2*lm ,k) = tmpGRcoef(k+plev*15,lm)
end do
end do
!
!-----------------------------------------------------------------------
!
! Computation for 1-level variables (ln(p*) and derivatives).
!
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
lmc = 2*lmr
!cdir shortloop
do n=1,nlen(m),2
ir = lmc + 2*n - 1
ii = ir + 1
grpma(2*lm-1) = grpma(2*lm-1) + alps(ir)*ldalp(lmr+n,irow)*ra
grpma(2*lm ) = grpma(2*lm ) + alps(ii)*ldalp(lmr+n,irow)*ra
end do
end do
!cdir novector
do lm=1,mlength
m = locm(lm,iam)
lmr = lnstart(lm)
lmc = 2*lmr
!cdir shortloop
do n=2,nlen(m),2
ir = lmc + 2*n - 1
ii = ir + 1
!
grpsa (2*lm-1) = grpsa (2*lm-1) + alps(ir)*lalp(lmr+n,irow)
grpsa (2*lm ) = grpsa (2*lm ) + alps(ii)*lalp(lmr+n,irow)
!
grdpsa(2*lm-1) = grdpsa(2*lm-1) + alps(ir)*lalp(lmr+n,irow)*hdfst4(m+n-1)*ztodt
grdpsa(2*lm ) = grdpsa(2*lm ) + alps(ii)*lalp(lmr+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
subroutine prepGRcalc(tmpSPEcoef) 1,8
!-----------------------------------------------------------------------
!
! Rearrange multi-level spectral coefficients for vectorization.
! The results are saved to "tmpSPEcoef" and will be used in
! "grcalcs" and "grcalca".
!
!-----------------------------------------------------------------------
!
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, only: hdiftq, hdifzd
!
implicit none
!
!
!---------------------------Output argument-----------------------------
!
real(r8), intent(out) :: tmpSPEcoef(plev*24,pnmax,maxm) ! array for rearranged variables
!
!---------------------------Local workspace-----------------------------
!
real(r8) raxm
!
integer lm, m, n, k
integer lmr, lmc
integer ir ,ii
!
!-----------------------------------------------------------------------
!
do lm=1,numm(iam)
m = locm(lm,iam)
lmr = lnstart(lm)
lmc = 2*lmr
raxm = ra*xm(m)
do n=1,nlen(m)
ir = lmc + 2*n - 1
ii = ir + 1
do k=1,plev
tmpSPEcoef(k ,n,lm) = vz(ir,k)
tmpSPEcoef(k+plev ,n,lm) = vz(ii,k)
tmpSPEcoef(k+plev*2 ,n,lm) = -d(ir,k)
tmpSPEcoef(k+plev*3 ,n,lm) = -d(ii,k)
tmpSPEcoef(k+plev*4 ,n,lm) = -vz(ir,k)*hdifzd(n+m-1,k)
tmpSPEcoef(k+plev*5 ,n,lm) = -vz(ii,k)*hdifzd(n+m-1,k)
tmpSPEcoef(k+plev*6 ,n,lm) = d(ir,k)*hdifzd(n+m-1,k)
tmpSPEcoef(k+plev*7 ,n,lm) = d(ii,k)*hdifzd(n+m-1,k)
tmpSPEcoef(k+plev*8 ,n,lm) = t(ir,k)
tmpSPEcoef(k+plev*9 ,n,lm) = t(ii,k)
tmpSPEcoef(k+plev*10,n,lm) = -t(ir,k)*hdiftq(n+m-1,k)
tmpSPEcoef(k+plev*11,n,lm) = -t(ii,k)*hdiftq(n+m-1,k)
tmpSPEcoef(k+plev*12,n,lm) = d(ir,k)
tmpSPEcoef(k+plev*13,n,lm) = d(ii,k)
tmpSPEcoef(k+plev*14,n,lm) = vz(ir,k)
tmpSPEcoef(k+plev*15,n,lm) = vz(ii,k)
tmpSPEcoef(k+plev*16,n,lm) = d (ii,k)*rsq(m+n-1)*raxm
tmpSPEcoef(k+plev*17,n,lm) = -d (ir,k)*rsq(m+n-1)*raxm
tmpSPEcoef(k+plev*18,n,lm) = vz(ii,k)*rsq(m+n-1)*raxm
tmpSPEcoef(k+plev*19,n,lm) = -vz(ir,k)*rsq(m+n-1)*raxm
tmpSPEcoef(k+plev*20,n,lm) = -d (ii,k)*hdifzd(n+m-1,k)*rsq(m+n-1)*raxm
tmpSPEcoef(k+plev*21,n,lm) = d (ir,k)*hdifzd(n+m-1,k)*rsq(m+n-1)*raxm
tmpSPEcoef(k+plev*22,n,lm) = -vz(ii,k)*hdifzd(n+m-1,k)*rsq(m+n-1)*raxm
tmpSPEcoef(k+plev*23,n,lm) = vz(ir,k)*hdifzd(n+m-1,k)*rsq(m+n-1)*raxm
end do
end do
end do
!
return
end subroutine prepGRcalc