#include <misc.h>#include <params.h>subroutine spetru (ps      ,phis    ,u3      ,v3      ,t3      , &                   q3      ,div     ,dpsl    ,dpsm    ,tl      , &                   tm      ,ql      ,qm      ,phi     ,phisl   , &                   phism   ,phis_hires)!-----------------------------------------------------------------------!! Purpose:! Spectrally truncate input fields which have already been transformed into ! fourier space.  Some arrays are dimensioned (2,...), where (1,...) is the! real part of the complex fourier coefficient, and (2,...) is the imaginary.! Any array dimensioned (plond,...) *cannot* be dimensioned (2,plond/2,...) ! because plond *may* be (and in fact currently is) odd. In these cases ! reference to real and imaginary parts is by (2*m-1,...) and (2*m  ,...) ! respectively.!! Author:  J. Rosinski!!-----------------------------------------------------------------------  use precision  use pmgrid  use constituents, only: pcnst, pnats  use pspect  use comspe  use rgrid,        only: nlon, nmmax  use commap,       only: w, xm, rsq, cs  use dynconst,     only: ez, ra, rearth  implicit none#include <comctl.h>#include <comfft.h>!------------------------------Arguments--------------------------------!  real(r8), intent(inout):: ps   (plond,plat)      ! Fourier -> spec. coeffs. for ln(ps)  real(r8), intent(inout):: phis (plond,plat)      ! Fourier -> spec. coeffs. for sfc geo.  real(r8), intent(inout):: u3   (plond,plev,plat) ! Fourier -> spec. coeffs. for u-wind  real(r8), intent(inout):: v3   (plond,plev,plat) ! Fourier -> spec. coeffs. for v-wind  real(r8), intent(inout):: t3   (plond,plev,plat) ! Fourier -> spec. coeffs. for temp  real(r8), intent(inout):: q3   (plond,plev*(pcnst+pnats),plat)!                                                  ! Fourier -> spec. coeffs. for q  real(r8), intent(out)  :: div  (plond,plev,plat) ! Spectrally truncated divergence  real(r8), intent(out)  :: dpsl (plond,plat)      ! Spectrally trunc d(ln(ps))/d(long)  real(r8), intent(out)  :: dpsm (plond,plat)      ! Spectrally trunc d(ln(ps))/d(lat )  real(r8), intent(out)  :: tl   (plond,plev,plat) ! Spectrally trunc d(T)/d(longitude)  real(r8), intent(out)  :: tm   (plond,plev,plat) ! Spectrally trunc d(T)/d(latitude)  real(r8), intent(out)  :: ql   (plond,plev,plat) ! Spectrally trunc d(q)/d(longitude)  real(r8), intent(out)  :: qm   (plond,plev,plat) ! Spectrally trunc d(q)/d(latitude)  real(r8), intent(out)  :: phi  (2,psp/2)         ! used in spectral truncation of phis  real(r8), intent(out)  :: phisl(plond,plat)      ! Spectrally trunc d(phis)/d(longitude)  real(r8), intent(out)  :: phism(plond,plat)      ! Spectrally trunc d(phis)/d(latitude)  logical , intent(in)   :: phis_hires             ! true => PHIS came from hi res topog!!---------------------------Local workspace-----------------------------!  real(r8) alpn (pspt)          ! alp*rsq*xm*ra  real(r8) dalpn(pspt)          ! dalp*rsq*ra  real(r8) tmp1                 ! vector temporary  real(r8) tmp2                 ! vector temporary  real(r8) tmpr                 ! vector temporary (real)  real(r8) tmpi                 ! vector temporary (imaginary)  real(r8) phialpr,phialpi      ! phi*alp (real and imaginary)  real(r8) phdalpr,phdalpi      ! phi*dalp (real and imaginary)  real(r8) zwalp                ! zw*alp  real(r8) zwdalp               ! zw*dalp  real(r8) psdalpr,psdalpi      ! alps (real and imaginary)*dalp  real(r8) psalpr,psalpi        ! alps (real and imaginary)*alp  real(r8) zrcsj                ! ra/(cos**2 latitude)  real(r8) zw                   ! w**2  real(r8) filtlim              ! filter function  real(r8) ft                   ! filter multiplier for spectral coefficients#if ( ! defined USEFFTLIB )  real(r8) work((plon+1)*plev)  ! Workspace for fft#else  real(r8) work((plon+1)*pcray) ! Workspace for fft#endif  real(r8) zsqcs  integer ir,ii                 ! indices complex coeffs. of spec. arrs.  integer i,k                   ! longitude, level indices  integer irow                  ! latitude pair index  integer latm,latp             ! symmetric latitude indices  integer lat                   ! index  integer m                     ! longitudinal wavenumber index (non-PVP)!                               ! along-diagonal index (PVP)  integer n                     ! latitudinal wavenumber index (non-PVP)!                               ! diagonal index (PVP)  integer nspec                 ! index#if ( defined PVP )                integer ne                    ! index into rsq  integer ic                    ! complex coeff. index  integer ialp                  ! index into legendre polynomials#else  integer mr,mc                 ! spectral indices#endif!!-----------------------------------------------------------------------!! Zero spectral arrays!  alps(:)   = 0.  vz  (:,:) = 0.  d   (:,:) = 0.  t   (:,:) = 0.  q   (:,:) = 0.  phi (:,:) = 0.!! Compute the quantities which are transformed to spectral space:!   1. u = u*sqrt(1-mu*mu),   u * cos(phi)!   2. v = v*sqrt(1-mu*mu),   v * cos(phi)!   3. t = t                  T!   4. ps= ln(ps). !  do lat=1,plat     irow = lat     if (lat.gt.plat/2) irow = plat - lat + 1     zsqcs = sqrt(cs(irow))     do k=1,plev        do i=1,nlon(lat)           u3(i,k,lat) = u3(i,k,lat)*zsqcs           v3(i,k,lat) = v3(i,k,lat)*zsqcs        end do     end do     do i=1,nlon(lat)        ps(i,lat) = log(ps(i,lat))     end do!! Transform grid -> fourier! 1st transform: U,V,T and Q! 2nd transform: LN(PS).  3rd transform: surface geopotential!     call fft991 (u3(1,1,lat),work     ,trig(1,irow),ifax(1,irow),1       , &                     plond   ,nlon(lat),plev        ,-1         )     call fft991 (v3(1,1,lat),work     ,trig(1,irow),ifax(1,irow),1       , &                     plond, nlon(lat)  ,plev        ,-1         )     call fft991 (t3(1,1,lat),work     ,trig(1,irow),ifax(1,irow),1       , &                     plond, nlon(lat)  ,plev        ,-1         )     call fft991 (q3(1,1,lat),work     ,trig(1,irow),ifax(1,irow),1       , &                     plond, nlon(lat)  ,plev        ,-1         )     call fft991 (ps(1,lat),work       ,trig(1,irow),ifax(1,irow),1       , &                     plond, nlon(lat)  ,1           ,-1         )     call fft991 (phis(1,lat),work     ,trig(1,irow),ifax(1,irow),1       , &                     plond, nlon(lat)  ,1           ,-1         )  end do                    ! lat=1,plat!! Loop over latitude pairs!  do irow=1,plat/2     latp = irow     latm = plat - irow + 1     zrcsj = ra/cs(irow)     zw = w(irow)*2.     do i=1,2*nmmax(irow)!! Compute symmetric and antisymmetric components: ps first, then phis!        tmp1 = 0.5*(ps(i,latm) - ps(i,latp))        tmp2 = 0.5*(ps(i,latm) + ps(i,latp))        ps(i,latm) = tmp1        ps(i,latp) = tmp2        tmp1 = 0.5*(phis(i,latm) - phis(i,latp))        tmp2 = 0.5*(phis(i,latm) + phis(i,latp))        phis(i,latm) = tmp1        phis(i,latp) = tmp2     end do!! Multi-level fields: U, V, T!     do k=1,plev        do i=1,2*nmmax(irow)           tmp1 = 0.5*(u3(i,k,latm) - u3(i,k,latp))           tmp2 = 0.5*(u3(i,k,latm) + u3(i,k,latp))           u3(i,k,latm) = tmp1           u3(i,k,latp) = tmp2           tmp1 = 0.5*(v3(i,k,latm) - v3(i,k,latp))           tmp2 = 0.5*(v3(i,k,latm) + v3(i,k,latp))           v3(i,k,latm) = tmp1           v3(i,k,latp) = tmp2           tmp1 = 0.5*(t3(i,k,latm) - t3(i,k,latp))           tmp2 = 0.5*(t3(i,k,latm) + t3(i,k,latp))           t3(i,k,latm) = tmp1           t3(i,k,latp) = tmp2           tmp1 = 0.5*(q3(i,k,latm) - q3(i,k,latp))           tmp2 = 0.5*(q3(i,k,latm) + q3(i,k,latp))           q3(i,k,latm) = tmp1           q3(i,k,latp) = tmp2        end do     end do!     ! Compute vzmn,dmn and ln(p*)mn and also phi*mn,tmn and qmn!#if ( defined PVP )     do n=1,pmax,2        ic = ncoefi(n) - 1        ialp = nalp(n)        do m=1,nmreduced(n,irow)           zwalp = zw*alp(ialp+m,irow)           phi(1,ic+m) = phi(1,ic+m) + zwalp*phis(2*m-1,latp)           phi(2,ic+m) = phi(2,ic+m) + zwalp*phis(2*m  ,latp)           ir = 2*(ic+m) - 1           ii = ir + 1           alps(ir) = alps(ir) + zwalp*ps(2*m-1,latp)           alps(ii) = alps(ii) + zwalp*ps(2*m  ,latp)        end do     end do!     do n=2,pmax,2        ic = ncoefi(n) - 1        ialp = nalp(n)        do m=1,nmreduced(n,irow)           zwalp = zw*alp(ialp+m,irow)           phi(1,ic+m) = phi(1,ic+m) + zwalp*phis(2*m-1,latm)           phi(2,ic+m) = phi(2,ic+m) + zwalp*phis(2*m  ,latm)           ir = 2*(ic+m) - 1           ii = ir + 1           alps(ir) = alps(ir) + zwalp*ps(2*m-1,latm)           alps(ii) = alps(ii) + zwalp*ps(2*m  ,latm)        end do     end do#else     do m=1,nmmax(irow)        mr = nstart(m)        mc = 2*mr        do n=1,nlen(m),2           zwalp = zw*alp(mr+n,irow)           phi(1,mr+n) = phi(1,mr+n) + zwalp*phis(2*m-1,latp)           phi(2,mr+n) = phi(2,mr+n) + zwalp*phis(2*m  ,latp)           ir = mc + 2*n - 1           ii = ir + 1           alps(ir) = alps(ir) + zwalp*ps(2*m-1,latp)           alps(ii) = alps(ii) + zwalp*ps(2*m  ,latp)        end do        do n=2,nlen(m),2           zwalp = zw*alp(mr+n,irow)           phi(1,mr+n) = phi(1,mr+n) + zwalp*phis(2*m-1,latm)           phi(2,mr+n) = phi(2,mr+n) + zwalp*phis(2*m  ,latm)           ir = mc + 2*n - 1           ii = ir + 1           alps(ir) = alps(ir) + zwalp*ps(2*m-1,latm)           alps(ii) = alps(ii) + zwalp*ps(2*m  ,latm)        end do     end do#endif     do k=1,plev#if ( defined PVP )        do n=1,pmax,2           ic = ncoefi(n) - 1           ialp = nalp(n)           do m=1,nmreduced(n,irow)              zwdalp   = zw*dalp(ialp+m,irow)              zwalp    = zw*alp (ialp+m,irow)              ir       = 2*(ic+m) - 1              ii       = ir + 1              d(ir,k)  = d(ir,k) - (zwdalp*v3(2*m-1,k,latm) + &                         xm(m)*zwalp*u3(2*m  ,k,latp))*zrcsj              d(ii,k)  = d(ii,k) - (zwdalp*v3(2*m  ,k,latm) - &                         xm(m)*zwalp*u3(2*m-1,k,latp))*zrcsj              t(ir,k)  = t(ir,k) + zwalp*t3(2*m-1,k,latp)              t(ii,k)  = t(ii,k) + zwalp*t3(2*m  ,k,latp)              q(ir,k)  = q(ir,k) + zwalp*q3(2*m-1,k,latp)              q(ii,k)  = q(ii,k) + zwalp*q3(2*m  ,k,latp)              vz(ir,k) = vz(ir,k) + (zwdalp*u3(2*m-1,k,latm) - &                         xm(m)*zwalp*v3(2*m  ,k,latp))*zrcsj              vz(ii,k) = vz(ii,k) + (zwdalp*u3(2*m  ,k,latm) + &                         xm(m)*zwalp*v3(2*m-1,k,latp))*zrcsj           end do        end do!        do n=2,pmax,2           ic = ncoefi(n) - 1           ialp = nalp(n)           do m=1,nmreduced(n,irow)              zwdalp   = zw*dalp(ialp+m,irow)              zwalp    = zw*alp (ialp+m,irow)              ir       = 2*(ic+m) - 1              ii       = ir + 1              d(ir,k)  = d(ir,k) - (zwdalp*v3(2*m-1,k,latp) + &                         xm(m)*zwalp*u3(2*m  ,k,latm))*zrcsj              d(ii,k)  = d(ii,k) - (zwdalp*v3(2*m  ,k,latp) - &                         xm(m)*zwalp*u3(2*m-1,k,latm))*zrcsj              t(ir,k)  = t(ir,k) + zwalp*t3(2*m-1,k,latm)              t(ii,k)  = t(ii,k) + zwalp*t3(2*m  ,k,latm)              q(ir,k)  = q(ir,k) + zwalp*q3(2*m-1,k,latm)              q(ii,k)  = q(ii,k) + zwalp*q3(2*m  ,k,latm)              vz(ir,k) = vz(ir,k) + (zwdalp*u3(2*m-1,k,latp) - &                         xm(m)*zwalp*v3(2*m  ,k,latm))*zrcsj              vz(ii,k) = vz(ii,k) + (zwdalp*u3(2*m  ,k,latp) + &                         xm(m)*zwalp*v3(2*m-1,k,latm))*zrcsj           end do        end do#else        do m=1,nmmax(irow)           mr = nstart(m)           mc = 2*mr           do n=1,nlen(m),2              zwdalp   = zw*dalp(mr+n,irow)              zwalp    = zw*alp (mr+n,irow)              ir       = mc + 2*n - 1              ii       = ir + 1              d(ir,k)  = d(ir,k) - (zwdalp*v3(2*m-1,k,latm) + &                         xm(m)*zwalp*u3(2*m  ,k,latp))*zrcsj              d(ii,k)  = d(ii,k) - (zwdalp*v3(2*m  ,k,latm) - &                         xm(m)*zwalp*u3(2*m-1,k,latp))*zrcsj              t(ir,k)  = t(ir,k) + zwalp*t3(2*m-1,k,latp)              t(ii,k)  = t(ii,k) + zwalp*t3(2*m  ,k,latp)              q(ir,k)  = q(ir,k) + zwalp*q3(2*m-1,k,latp)              q(ii,k)  = q(ii,k) + zwalp*q3(2*m  ,k,latp)              vz(ir,k) = vz(ir,k) + (zwdalp*u3(2*m-1,k,latm) - &                         xm(m)*zwalp*v3(2*m  ,k,latp))*zrcsj              vz(ii,k) = vz(ii,k) + (zwdalp*u3(2*m  ,k,latm) + &                         xm(m)*zwalp*v3(2*m-1,k,latp))*zrcsj           end do        end do        do m=1,nmmax(irow)           mr = nstart(m)           mc = 2*mr           do n=2,nlen(m),2              zwdalp   = zw*dalp(mr+n,irow)              zwalp    = zw*alp (mr+n,irow)              ir       = mc + 2*n - 1              ii       = ir + 1              d(ir,k)  = d(ir,k) - (zwdalp*v3(2*m-1,k,latp) + &                         xm(m)*zwalp*u3(2*m  ,k,latm))*zrcsj              d(ii,k)  = d(ii,k) - (zwdalp*v3(2*m  ,k,latp) - &                         xm(m)*zwalp*u3(2*m-1,k,latm))*zrcsj              t(ir,k)  = t(ir,k) + zwalp*t3(2*m-1,k,latm)              t(ii,k)  = t(ii,k) + zwalp*t3(2*m  ,k,latm)