#include <misc.h>
#include <params.h>
subroutine quad(lm ,zdt ,ztdtsq ,grlps1 ,grlps2 ,& 1,7
grt1 ,grz1 ,grd1 ,grfu1 ,grfv1 ,&
grvt1 ,grrh1 ,grt2 ,grz2 ,grd2 ,&
grfu2 ,grfv2 ,grvt2 ,grrh2 )
!-----------------------------------------------------------------------
!
! Perform gaussian quadrature for 1 Fourier wavenumber (m) to obtain the
! spectral coefficients of ln(ps), temperature, vorticity, and divergence.
! Add the tendency terms requiring meridional derivatives during the
! transform.
!
!---------------------------Code history--------------------------------
!
! Original version: J. Rosinski
! Standardized: J. Rosinski, June 1992
! Reviewed: B. Boville, D. Williamson, J. Hack, August 1992
! Reviewed: B. Boville, D. Williamson, April 1996
! Modified: P. Worley, September 2002
! Modified: NEC, April 2004
!
!-----------------------------------------------------------------------
use shr_kind_mod, only: r8 => shr_kind_r8
use pmgrid
use pspect
use comspe
use rgrid
use commap
use dynconst, only: rearth
implicit none
!
! Input arguments
!
integer, intent(in) :: lm ! local Fourier wavenumber index
real(r8), intent(in) :: zdt ! timestep(dt) unless nstep = 0
real(r8), intent(in) :: ztdtsq(pnmax) ! 2*zdt*n(n+1)/(a^2)
! where n IS the 2-d wavenumber
!
! Fourier coefficient arrays which have a latitude index on them for
! multitasking. These arrays are defined in LINEMS and and in QUAD
! to compute spectral coefficients. They contain a latitude index so
! that the sums over latitude can be performed in a specified order.
!
! Suffixes 1 and 2 refer to symmetric and antisymmetric components
! respectively.
!
real(r8), intent(in) :: grlps1(2*maxm,plat/2) ! ln(ps) - symmetric
real(r8), intent(in) :: grlps2(2*maxm,plat/2) ! ln(ps) - antisymmetric
!
! symmetric components
!
real(r8), intent(in) :: grt1(2*maxm,plev,plat/2) ! temperature
real(r8), intent(in) :: grz1(2*maxm,plev,plat/2) ! vorticity
real(r8), intent(in) :: grd1(2*maxm,plev,plat/2) ! divergence
real(r8), intent(in) :: grfu1(2*maxm,plev,plat/2) ! partial u momentum tendency (fu)
real(r8), intent(in) :: grfv1(2*maxm,plev,plat/2) ! partial v momentum tendency (fv)
real(r8), intent(in) :: grvt1(2*maxm,plev,plat/2) ! heat flux
real(r8), intent(in) :: grrh1(2*maxm,plev,plat/2) ! rhs of div eqn (del^2 term)
!
! antisymmetric components
!
real(r8), intent(in) :: grt2(2*maxm,plev,plat/2) ! temperature
real(r8), intent(in) :: grz2(2*maxm,plev,plat/2) ! vorticity
real(r8), intent(in) :: grd2(2*maxm,plev,plat/2) ! divergence
real(r8), intent(in) :: grfu2(2*maxm,plev,plat/2) ! partial u momentum tend (fu)
real(r8), intent(in) :: grfv2(2*maxm,plev,plat/2) ! partial v momentum tend (fv)
real(r8), intent(in) :: grvt2(2*maxm,plev,plat/2) ! heat flux
real(r8), intent(in) :: grrh2(2*maxm,plev,plat/2) ! rhs of div eqn (del^2 term)
!
!---------------------------Local workspace-----------------------------
!
integer j ! latitude pair index
integer m ! global wavenumber index
integer n ! total wavenumber index
integer ir,ii ! spectral indices
integer lmr,lmc ! spectral indices
integer k ! level index
integer kv ! index for vectorization
real(r8) zcsj ! cos**2(lat)*radius of earth
real(r8) zrcsj ! 1./(a*cos^2(lat))
real(r8) zdtrc ! dt/(a*cos^2(lat))
real(r8) ztdtrc ! 2dt/(a*cos^2(lat))
real(r8) zw(plat/2) ! 2*w
real(r8) ztdtrw(plat/2) ! 2w*2dt/(a*cos^2(lat))
real(r8) zwalp ! zw*alp
real(r8) zwdalp ! zw*dalp
real(r8) sqzwalp ! ztdtsq*zw*alp
real(r8) tmpGR1odd(plev*6,plat/2) ! temporary space for Fourier coeffs
real(r8) tmpGR2odd(plev*6,plat/2) !
real(r8) tmpGR3odd(plev*6,plat/2) !
real(r8) tmpGR1evn(plev*6,plat/2) !
real(r8) tmpGR2evn(plev*6,plat/2) !
real(r8) tmpGR3evn(plev*6,plat/2) !
real(r8) tmpSPEodd(plev*6,2*ptrn) ! temporary space for spectral coeffs
real(r8) tmpSPEevn(plev*6,2*ptrn) !
!
!-----------------------------------------------------------------------
!
! Compute constants
!
do j=1,plat/2
zcsj = cs(j)*rearth
zrcsj = 1./zcsj
zdtrc = zdt*zrcsj
ztdtrc = 2.*zdtrc
zw(j) = w(j)*2.
ztdtrw(j) = ztdtrc*zw(j)
end do
!
! Accumulate contributions to spectral coefficients of ln(p*), the only
! single level field. Use symmetric or antisymmetric fourier cofficients
! depending on whether the total wavenumber is even or odd.
!
m = locm(lm,iam)
lmr = lnstart(lm)
lmc = 2*lmr
do n=1,2*nlen(m)
alps(lmc+n) = 0.
end do
do n=1,nlen(m),2
do j=beglatpair(m),plat/2
ir = lmc + 2*n - 1
ii = ir + 1
zwalp = zw(j)*lalp(lmr+n,j)
alps(ir) = alps(ir) + grlps1(2*lm-1,j)*zwalp
alps(ii) = alps(ii) + grlps1(2*lm ,j)*zwalp
end do
end do
do n=2,nlen(m),2
do j=beglatpair(m),plat/2
ir = lmc + 2*n - 1
ii = ir + 1
zwalp = zw(j)*lalp(lmr+n,j)
alps(ir) = alps(ir) + grlps2(2*lm-1,j)*zwalp
alps(ii) = alps(ii) + grlps2(2*lm ,j)*zwalp
end do
end do
!
! Accumulate contributions to spectral coefficients of the multilevel fields.
! Use symmetric or antisymmetric fourier coefficients depending on whether
! the total wavenumber is even or odd.
!
!
! Initialize temporary storage for spectral coefficients
!
do n=1,nlen(m)
do kv=1,plev*6
tmpSPEodd(kv,n) = 0.
tmpSPEevn(kv,n) = 0.
end do
end do
!
! Rearrange Fourier coefficients to temporal storage
!
do j = beglatpair(m),plat/2
do k=1,plev
tmpGR1odd(k ,j) = grt1 (2*lm-1,k,j) ! first term for odd n
tmpGR1odd(k+plev ,j) = grt1 (2*lm ,k,j)
tmpGR1odd(k+plev*2,j) = grd1 (2*lm-1,k,j)
tmpGR1odd(k+plev*3,j) = grd1 (2*lm ,k,j)
tmpGR1odd(k+plev*4,j) = grz1 (2*lm-1,k,j)
tmpGR1odd(k+plev*5,j) = grz1 (2*lm ,k,j)
tmpGR2odd(k ,j) = grvt2(2*lm-1,k,j) ! second term for odd n
tmpGR2odd(k+plev ,j) = grvt2(2*lm ,k,j)
tmpGR2odd(k+plev*2,j) = -grfv2(2*lm-1,k,j)
tmpGR2odd(k+plev*3,j) = -grfv2(2*lm ,k,j)
tmpGR2odd(k+plev*4,j) = grfu2(2*lm-1,k,j)
tmpGR2odd(k+plev*5,j) = grfu2(2*lm ,k,j)
tmpGR3odd(k+plev*2,j) = grrh1(2*lm-1,k,j) ! additional term for odd n
tmpGR3odd(k+plev*3,j) = grrh1(2*lm ,k,j)
tmpGR1evn(k ,j) = grt2 (2*lm-1,k,j) ! first term for even n
tmpGR1evn(k+plev ,j) = grt2 (2*lm ,k,j)
tmpGR1evn(k+plev*2,j) = grd2 (2*lm-1,k,j)
tmpGR1evn(k+plev*3,j) = grd2 (2*lm ,k,j)
tmpGR1evn(k+plev*4,j) = grz2 (2*lm-1,k,j)
tmpGR1evn(k+plev*5,j) = grz2 (2*lm ,k,j)
tmpGR2evn(k ,j) = grvt1(2*lm-1,k,j) ! first term for even n
tmpGR2evn(k+plev ,j) = grvt1(2*lm ,k,j)
tmpGR2evn(k+plev*2,j) = -grfv1(2*lm-1,k,j)
tmpGR2evn(k+plev*3,j) = -grfv1(2*lm ,k,j)
tmpGR2evn(k+plev*4,j) = grfu1(2*lm-1,k,j)
tmpGR2evn(k+plev*5,j) = grfu1(2*lm ,k,j)
tmpGR3evn(k+plev*2,j) = grrh2(2*lm-1,k,j) ! additional term for even n
tmpGR3evn(k+plev*3,j) = grrh2(2*lm ,k,j)
end do
end do
!
! Accumulate first and second terms
!
do j=beglatpair(m),plat/2
do n=1,nlen(m),2
zwdalp = ztdtrw(j)*ldalp(lmr+n,j)
zwalp = zw(j) *lalp (lmr+n,j)
do kv=1,plev*6
tmpSPEodd(kv,n) = tmpSPEodd(kv,n) + &
zwalp*tmpGR1odd(kv,j) + zwdalp*tmpGR2odd(kv,j)
end do
end do
end do
do j=beglatpair(m),plat/2
do n=2,nlen(m),2
zwdalp = ztdtrw(j)*ldalp(lmr+n,j)
zwalp = zw(j) *lalp (lmr+n,j)
do kv=1,plev*6
tmpSPEevn(kv,n) = tmpSPEevn(kv,n) + &
zwalp*tmpGR1evn(kv,j) + zwdalp*tmpGR2evn(kv,j)
end do
end do
end do
!
! Add additional term for divergence
!
do n=1,nlen(m),2
do j=beglatpair(m),plat/2
sqzwalp = ztdtsq(n+m-1)*zw(j)*lalp (lmr+n,j)
do kv=plev*2+1,plev*4
tmpSPEodd(kv,n) = tmpSPEodd(kv,n) + sqzwalp*tmpGR3odd(kv,j)
end do
end do
end do
do n=2,nlen(m),2
do j=beglatpair(m),plat/2
sqzwalp = ztdtsq(n+m-1)*zw(j)*lalp (lmr+n,j)
do kv=plev*2+1,plev*4
tmpSPEevn(kv,n) = tmpSPEevn(kv,n) + sqzwalp*tmpGR3evn(kv,j)
end do
end do
end do
!
! Save accumulated results
!
do n=1,nlen(m),2
ir = lmc+2*n-1
ii = ir+1
do k=1,plev
t (ir,k) = tmpSPEodd(k ,n)
t (ii,k) = tmpSPEodd(k+plev ,n)
d (ir,k) = tmpSPEodd(k+plev*2,n)
d (ii,k) = tmpSPEodd(k+plev*3,n)
vz(ir,k) = tmpSPEodd(k+plev*4,n)
vz(ii,k) = tmpSPEodd(k+plev*5,n)
end do
end do
do n=2,nlen(m),2
ir = lmc+2*n-1
ii = ir+1
do k=1,plev
t (ir,k) = tmpSPEevn(k ,n)
t (ii,k) = tmpSPEevn(k+plev ,n)
d (ir,k) = tmpSPEevn(k+plev*2,n)
d (ii,k) = tmpSPEevn(k+plev*3,n)
vz(ir,k) = tmpSPEevn(k+plev*4,n)
vz(ii,k) = tmpSPEevn(k+plev*5,n)
end do
end do
!
return
end subroutine quad