spetru.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#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)