#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       ,grlps1  ,grlps2  ,grt1    ,grq1    , &                grz1    ,grd1    ,grfu1   ,grfv1   ,grt2    , &                grq2    ,grz2    ,grd2    ,grfu2   ,grfv2   )!-----------------------------------------------------------------------!! Purpose:! Perform gaussian quadrature to obtain the spectral coefficients of ! ln(ps), temperature, vorticity, divergence, and specific humidity! (to calculate derivatives only). Add the tendency terms requiring! meridional derivatives during the transform.!! 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:  J. Rosinski!!-----------------------------------------------------------------------!! $Id: quad.F90,v 1.4 2001/04/13 22:40:42 rosinski Exp $! $Author: rosinski $!!-----------------------------------------------------------------------  use precision  use pmgrid  use pspect  use comspe  use rgrid  use dynconst, only: rearth  use commap  implicit none!------------------------------Arguments--------------------------------!! 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.!  integer , intent(in)   :: n                           ! spectral space "diagonal" index  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)   :: grq1  (2*pmmax,plev,plat/2) ! moisture  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)!! antisymmetric components!  real(r8), intent(in)   :: grt2  (2*pmmax,plev,plat/2) ! temperature  real(r8), intent(in)   :: grq2  (2*pmmax,plev,plat/2) ! moisture  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)!!---------------------------Local workspace-----------------------------!  real(r8) alp2  (2*pmax,plat/2)  ! expand polynomials to complex  real(r8) dalp2 (2*pmax,plat/2)  ! expand polynomials to complex  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 ialp                    ! index into polynomials  real(r8) zcsj                   ! cos**2(lat)*radius of earth  real(r8) zrcsj                  ! 1./(a*cos^2(lat))  real(r8) zw    (plat/2)         ! 2*w  real(r8) ztdtrw(plat/2)         ! 2w*2dt/(a*cos^2(lat))!!-----------------------------------------------------------------------!! Compute constants!  do j = 1,plat/2     zcsj      = cs(j)*rearth     zrcsj     = 1./zcsj     zw(j)     = w(j)*2.     ztdtrw(j) = zrcsj*zw(j)  end do  ialp = nalp(n)  isp  = nco2(n) - 2  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.           q (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)                 q (isp+m,k) = q (isp+m,k) + grq1 (m,k,j)*alp2 (m,j)                 d (isp+m,k) = d (isp+m,k) + grd1 (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)                 q (isp+m,k) = q (isp+m,k) + grq2 (m,k,j)*alp2 (m,j)                 d (isp+m,k) = d (isp+m,k) + grd2 (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     end do  else                       ! vectorize over levels     do m = 1,2*nm(n)        do k = 1,plev           t (isp+m,k) = 0.           q (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 = beglatpair((m+1)/2),plat/2              do k = 1,plev                 t (isp+m,k) = t (isp+m,k) + grt1 (m,k,j)*alp2 (m,j)                 q (isp+m,k) = q (isp+m,k) + grq1 (m,k,j)*alp2 (m,j)                 d (isp+m,k) = d (isp+m,k) + grd1 (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=beglatpair((m+1)/2),plat/2              do k = 1,plev                 t (isp+m,k) = t (isp+m,k) + grt2 (m,k,j)*alp2 (m,j)                 q (isp+m,k) = q (isp+m,k) + grq2 (m,k,j)*alp2 (m,j)                 d (isp+m,k) = d (isp+m,k) + grd2 (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     end do  end if!  returnend subroutine quad#elsesubroutine quad(m       ,grlps1  ,grlps2  ,grt1    ,grq1    , &                grz1    ,grd1    ,grfu1   ,grfv1   ,grt2    , &                grq2    ,grz2    ,grd2    ,grfu2   ,grfv2   )!-----------------------------------------------------------------------!! Purpose:! Perform gaussian quadrature for 1 Fourier wavenumber (m) to obtain the! spectral coefficients of ln(ps), temperature, vorticity, divergence! and specific humidity.  Add the tendency terms requiring meridional! derivatives during the transform.!! Author:  J. Rosinski!!-----------------------------------------------------------------------  use precision  use pmgrid  use pspect  use comspe  use rgrid  use commap  use dynconst, only: rearth  implicit none!------------------------------Arguments--------------------------------!! 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.!  integer , intent(in)   :: m                           ! Fourier wavenumber  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  (plev,2*pmmax,plat/2) ! temperature  real(r8), intent(in)   :: grq1  (plev,2*pmmax,plat/2) ! moisture  real(r8), intent(in)   :: grz1  (plev,2*pmmax,plat/2) ! vorticity  real(r8), intent(in)   :: grd1  (plev,2*pmmax,plat/2) ! divergence  real(r8), intent(in)   :: grfu1 (plev,2*pmmax,plat/2) ! partial u momentum tendency (fu)  real(r8), intent(in)   :: grfv1 (plev,2*pmmax,plat/2) ! partial v momentum tendency (fv)!! antisymmetric components!  real(r8), intent(in)   :: grt2  (plev,2*pmmax,plat/2) ! temperature  real(r8), intent(in)   :: grq2  (plev,2*pmmax,plat/2) ! moisture  real(r8), intent(in)   :: grz2  (plev,2*pmmax,plat/2) ! vorticity  real(r8), intent(in)   :: grd2  (plev,2*pmmax,plat/2) ! divergence  real(r8), intent(in)   :: grfu2 (plev,2*pmmax,plat/2) ! partial u momentum tend (fu)  real(r8), intent(in)   :: grfv2 (plev,2*pmmax,plat/2) ! partial v momentum tend (fv)!!---------------------------Local workspace-----------------------------!  integer j                 ! latitude pair index  integer n                 ! index  integer ir,ii             ! indices  integer mr,mc             ! indices  integer k                 ! level index  real(r8) zcsj             ! cos**2(lat)*radius of earth  real(r8) zrcsj            ! 1./(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  real(r8) zwdalp!!-----------------------------------------------------------------------!! Compute constants!  do j = 1,plat/2     zcsj      = cs(j)*rearth     zrcsj     = 1./zcsj     zw(j)     = w(j)*2.     ztdtrw(j) = zrcsj*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.!  mr = nstart(m)  mc = 2*mr  do n = 1,2*nlen(m)     alps(mc+n) = 0.  end do  do j=beglatpair(m),plat/2     do n=1,nlen(m),2        ir = mc + 2*n - 1        ii = ir + 1        zwalp    = zw(j)*alp(mr+n,j)        alps(ir) = alps(ir) + grlps1(2*m-1,j)*zwalp        alps(ii) = alps(ii) + grlps1(2*m  ,j)*zwalp     end do     do n = 2,nlen(m),2        ir = mc + 2*n - 1        ii = ir + 1        zwalp    = zw(j)*alp(mr+n,j)        alps(ir) = alps(ir) + grlps2(2*m-1,j)*zwalp        alps(ii) = alps(ii) + grlps2(2*m  ,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.!  do k = 1,plev     do n = 1,2*nlen(m)        t (mc+n,k) = 0.        q (mc+n,k) = 0.        d (mc+n,k) = 0.        vz(mc+n,k) = 0.     end do     do j=beglatpair(m),plat/2        do n=1,nlen(m),2           zwdalp = ztdtrw(j)*dalp(mr+n,j)           zwalp  = zw(j)    *alp (mr+n,j)           ir = mc + 2*n - 1           ii = ir + 1           t (ir,k) = t (ir,k) + zwalp*grt1 (k,2*m-1,j)           t (ii,k) = t (ii,k) + zwalp*grt1 (k,2*m  ,j)           q (ir,k) = q (ir,k) + zwalp*grq1 (k,2*m-1,j)           q (ii,k) = q (ii,k) + zwalp*grq1 (k,2*m  ,j)           d (ir,k) = d (ir,k) + grd1 (k,2*m-1,j)*zwalp - grfv2(k,2*m-1,j)*zwdalp           d (ii,k) = d (ii,k) + grd1 (k,2*m  ,j)*zwalp - grfv2(k,2*m  ,j)*zwdalp           vz(ir,k) = vz(ir,k) + grz1 (k,2*m-1,j)*zwalp + grfu2(k,2*m-1,j)*zwdalp           vz(ii,k) = vz(ii,k) + grz1 (k,2*m  ,j)*zwalp + grfu2(k,2*m  ,j)*zwdalp        end do     end do     do j=beglatpair(m),plat/2        do n = 2,nlen(m),2           zwdalp = ztdtrw(j)*dalp(mr+n,j)           zwalp  = zw(j)    *alp (mr+n,j)           ir = mc + 2*n - 1           ii = ir + 1           t (ir,k) = t (ir,k) + zwalp*grt2(k,2*m-1,j)           t (ii,k) = t (ii,k) + zwalp*grt2(k,2*m  ,j)           q (ir,k) = q (ir,k) + zwalp*grq2(k,2*m-1,j)           q (ii,k) = q (ii,k) + zwalp*grq2(k,2*m  ,j)           d (ir,k) = d (ir,k) + grd2 (k,2*m-1,j)*zwalp - grfv1(k,2*m-1,j)*zwdalp           d (ii,k) = d (ii,k) + grd2 (k,2*m  ,j)*zwalp - grfv1(k,2*m  ,j)*zwdalp           vz(ir,k) = vz(ir,k) + grz2 (k,2*m-1,j)*zwalp + grfu1(k,2*m-1,j)*zwdalp           vz(ii,k) = vz(ii,k) + grz2 (k,2*m  ,j)*zwalp + grfu1(k,2*m  ,j)*zwdalp        end do     end do  end do!  returnend subroutine quad#endif