end do   else                       ! vectorize over levels      do m=1,2*nm(n)         do k=1,plev            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=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) + &                     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=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) + &                     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      end do   end if!   returnend subroutine quad#elsesubroutine quad(m       ,zdt     ,ztdtsq  ,grlps1  ,grlps2  ,&                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!!-----------------------------------------------------------------------   use precision   use pmgrid   use pspect   use comspe   use rgrid   use commap   use dynconst, only: rearth   implicit none!! Input arguments!   integer, intent(in) :: m                          ! Fourier wavenumber   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 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(plev,2*pmmax,plat/2)     ! temperature   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)   real(r8), intent(in) :: grvt1(plev,2*pmmax,plat/2)    ! heat flux   real(r8), intent(in) :: grrh1(plev,2*pmmax,plat/2)    ! rhs of div eqn (del^2 term)!! antisymmetric components!   real(r8), intent(in) :: grt2(plev,2*pmmax,plat/2)     ! temperature   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)   real(r8), intent(in) :: grvt2(plev,2*pmmax,plat/2)    ! heat flux   real(r8), intent(in) :: grrh2(plev,2*pmmax,plat/2)    ! rhs of div eqn (del^2 term)!!---------------------------Local workspace-----------------------------!   integer j                          ! latitude pair index   integer n                          ! total wavenumber index   integer ir,ii                      ! spectral indices   integer mr,mc                      ! spectral indices   integer k                          ! level index   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!!-----------------------------------------------------------------------!! 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.!   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.         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) + zwdalp*grvt2(k,2*m-1,j)            t(ii,k) = t(ii,k) + zwalp*grt1 (k,2*m  ,j) + zwdalp*grvt2(k,2*m  ,j)            d(ir,k) = d(ir,k) + (grd1(k,2*m-1,j) + &               ztdtsq(n+m-1)*grrh1(k,2*m-1,j))*zwalp - &               grfv2(k,2*m-1,j)*zwdalp            d(ii,k) = d(ii,k) + (grd1(k,2*m  ,j) + &               ztdtsq(n+m-1)*grrh1(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) + &               zwdalp*grvt1(k,2*m-1,j)            t(ii,k) = t(ii,k) + zwalp*grt2(k,2*m  ,j) + &               zwdalp*grvt1(k,2*m  ,j)            d(ir,k) = d(ir,k) + (grd2(k,2*m-1,j) + &               ztdtsq(n+m-1)*grrh2(k,2*m-1,j))*zwalp - &               grfv1(k,2*m-1,j)*zwdalp            d(ii,k) = d(ii,k) + (grd2(k,2*m  ,j) + &               ztdtsq(n+m-1)*grrh2(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