#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 quad(n       ,zdt     ,ztdtsq  ,grlps1  ,grlps2  ,&                grt1    ,grz1    ,grd1    ,grfu1   ,grfv1   ,&                grvt1   ,grrh1   ,grt2    ,grz2    ,grd2    ,&                grfu2   ,grfv2   ,grvt2   ,grrh2   )!----------------------------------------------------------------------- ! ! Purpose: ! Perform gaussian quadrature to obtain the spectral coefficients of ! ln(ps), temperature, vorticity, and divergence. Add the tendency terms! requiring meridional derivatives during the transform.! ! Method: ! Computational note: This routine is multitasked over each diagonal in! the spectral (i.e., [m,n]) data structure since each of these vector! elements are computationally independent of each other. Care should be ! taken in interpreting the code since the diagonals are indexed using the ! variable "n" (which is often used to denote 2-dimensional wavenumber).! ! Author: ! 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!!-----------------------------------------------------------------------   use precision   use pmgrid   use pspect   use comspe   use rgrid   use commap   use dynconst, only: rearth!-----------------------------------------------------------------------   implicit none!-----------------------------------------------------------------------!! Input arguments!   integer, intent(in) :: n                          ! spectral space "diagonal" index   real(r8), intent(in) :: zdt                           ! timestep(dt) unless nstep = 0   real(r8), intent(in) :: ztdtsq(2*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 used 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*pmmax,plat/2)        ! ln(ps) - symmetric   real(r8), intent(in) :: grlps2(2*pmmax,plat/2)        ! ln(ps) - antisymmetric!! symmetric components!   real(r8), intent(in) :: grt1(2*pmmax,plev,plat/2)     ! temperature   real(r8), intent(in) :: grz1(2*pmmax,plev,plat/2)     ! vorticity   real(r8), intent(in) :: grd1(2*pmmax,plev,plat/2)     ! divergence   real(r8), intent(in) :: grfu1(2*pmmax,plev,plat/2)    ! partial u momentum tendency (fu)   real(r8), intent(in) :: grfv1(2*pmmax,plev,plat/2)    ! partial v momentum tendency (fv)   real(r8), intent(in) :: grvt1(2*pmmax,plev,plat/2)    ! heat flux   real(r8), intent(in) :: grrh1(2*pmmax,plev,plat/2)    ! rhs of div eqn (del^2 term)!! antisymmetric components!   real(r8), intent(in) :: grt2(2*pmmax,plev,plat/2)     ! temperature   real(r8), intent(in) :: grz2(2*pmmax,plev,plat/2)     ! vorticity   real(r8), intent(in) :: grd2(2*pmmax,plev,plat/2)     ! divergence   real(r8), intent(in) :: grfu2(2*pmmax,plev,plat/2)    ! partial u momentum tend (fu)   real(r8), intent(in) :: grfv2(2*pmmax,plev,plat/2)    ! partial v momentum tend (fv)   real(r8), intent(in) :: grvt2(2*pmmax,plev,plat/2)    ! heat flux   real(r8), intent(in) :: grrh2(2*pmmax,plev,plat/2)    ! rhs of div eqn (del^2 term)!!---------------------------Local workspace-----------------------------!   real(r8) alp2(2*pmax,plat/2)           ! expand polynomials to complex   real(r8) dalp2(2*pmax,plat/2)          ! expand polynomials to complex   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))   integer j                          ! latitude pair index   integer m                          ! diagonal element(index)of sp.arr.   integer isp                        ! index into levls of sp. arrs.   integer k                          ! level index   integer ne                         ! index into ztdtsq   integer ialp                       ! index into polynomials!!-----------------------------------------------------------------------!! 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   ialp = nalp(n)   isp = nco2(n) - 2   ne = 2*(n-1)   if (nm(n).gt.plat/2) then   ! vectorize over diagonals      do j=1,plat/2!! Expand polynomials to complex form to allow largest possible vector length!!DIR$ IVDEP         do m=1,nmreduced(n,j)            alp2 (2*m-1,j) = alp(ialp+m,j)*zw(j)            alp2 (2*m  ,j) = alp(ialp+m,j)*zw(j)            dalp2(2*m-1,j) = dalp(ialp+m,j)*ztdtrw(j)            dalp2(2*m  ,j) = dalp(ialp+m,j)*ztdtrw(j)         end do      end do   else                         ! vectorize over latitude      do m=1,nm(n)!! Expand polynomials to complex form to allow largest possible vector length!!DIR$ IVDEP         do j=1,plat/2            alp2 (2*m-1,j) = alp(ialp+m,j)*zw(j)            alp2 (2*m  ,j) = alp(ialp+m,j)*zw(j)            dalp2(2*m-1,j) = dalp(ialp+m,j)*ztdtrw(j)            dalp2(2*m  ,j) = dalp(ialp+m,j)*ztdtrw(j)         end do      end do   end if!! 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.! Vectorize over diagonals or latitude pair index depending on vector length.! Comparison is 2*nm(n) vs. plat/8 instead of plat/2 since latitude sum is! into a scalar variable.!   do m=1,2*nm(n)      alps(isp+m) = 0.   end do   if (2*nm(n).gt.plat/8) then     ! vectorize over diagonals      if (mod(n,2).ne.0) then         do j=1,plat/2            do m=1,2*nmreduced(n,j)               alps(isp+m) = alps(isp+m) + grlps1(m,j)*alp2(m,j)            end do         end do      else         do j=1,plat/2            do m=1,2*nmreduced(n,j)               alps(isp+m) = alps(isp+m) + grlps2(m,j)*alp2(m,j)            end do         end do      end if   else                            ! vectorize over latitude      if (mod(n,2).ne.0) then         do m=1,2*nm(n)            do j=beglatpair((m+1)/2),plat/2               alps(isp+m) = alps(isp+m) + grlps1(m,j)*alp2(m,j)            end do         end do      else         do m=1,2*nm(n)            do j=beglatpair((m+1)/2),plat/2               alps(isp+m) = alps(isp+m) + grlps2(m,j)*alp2(m,j)            end do         end do      end if   end if!! Accumulate contributions to spectral coefficients of the multilevel fields! (temperature, divergence, vorticity).! Use symmetric or antisymmetric fourier coefficients depending on whether! the total wavenumber is even or odd.!   if (2*nm(n).gt.plev) then     ! vectorize over diagonals      do k=1,plev         do m=1,2*nm(n)            t(isp+m,k) = 0.            d(isp+m,k) = 0.            vz(isp+m,k) = 0.         end do         if (mod(n,2).ne.0) then ! n is odd            do j=1,plat/2               do m=1,2*nmreduced(n,j)                  t(isp+m,k) = t(isp+m,k) + grt1(m,k,j)* alp2(m,j) + &                     grvt2(m,k,j)*dalp2(m,j)                  d(isp+m,k) = d(isp+m,k) + (grd1(m,k,j) + &                     ztdtsq(ne+m)*grrh1(m,k,j))*alp2(m,j) - &                     grfv2(m,k,j)*dalp2(m,j)                  vz(isp+m,k) = vz(isp+m,k) + grz1(m,k,j)* alp2(m,j) + &                     grfu2(m,k,j)*dalp2(m,j)               end do            end do         else                    ! n is even            do j=1,plat/2               do m=1,2*nmreduced(n,j)                  t(isp+m,k) = t(isp+m,k) + grt2(m,k,j)* alp2(m,j) + &                     grvt1(m,k,j)*dalp2(m,j)                  d(isp+m,k) = d(isp+m,k) + (grd2(m,k,j) + &                     ztdtsq(ne+m)*grrh2(m,k,j))*alp2(m,j) - &                     grfv1(m,k,j)*dalp2(m,j)                  vz(isp+m,k) = vz(isp+m,k) + grz2(m,k,j)* alp2(m,j) + &                     grfu1(m,k,j)*dalp2(m,j)               end do            end do         end if