spetru.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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