#include <misc.h>#include <params.h>subroutine spetru(ps      ,phis    ,u3      ,v3      ,t3      , &                  vort    ,div     ,dpsl    ,dpsm    ,phis_hires)!----------------------------------------------------------------------- ! ! Purpose: ! ! Method: ! 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: ! Original version:  J. Rosinski! Standardized:      J. Rosinski, June 1992! Reviewed:          B. Boville, J. Hack, August 1992!!-----------------------------------------------------------------------   use precision   use pmgrid,   only: plon, plond, plev, plat   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>!! Input/Output 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 temperature   logical, intent(in) :: phis_hires          ! true => PHIS came from hi res topo file!! Output arguments!   real(r8), intent(out) :: vort(plond,plev,plat)    ! Spectrally truncated vorticity   real(r8), intent(out) :: div(plond,plev,plat)     ! Spectrally truncated divergence   real(r8), intent(out) :: dpsl(plond,plat)         ! Spectrally trunc d(ln(ps))/d(longitude)   real(r8), intent(out) :: dpsm(plond,plat)         ! Spectrally trunc d(ln(ps))/d(latitude)!!---------------------------Local workspace-----------------------------!   real(r8) phi(2,psp/2)       ! used in spectral truncation of phis   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) zcor               ! correction for absolute vorticity   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   integer m                   ! longitudinal wavenumber index (non-PVP)!                                   along-diagonal index (PVP)   integer n                   ! latitudinal wavenumber index (non-PVP)!                                   diagonal index (PVP)   integer nspec#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!   vz(:,:) = 0.   d(:,:) = 0.   t(:,:) = 0.   alps(:) = 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: note contiguity assumptions! 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(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 160 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         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 150 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)               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)               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) + &