spegrd.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>

subroutine spegrd (ztodt   ,lat     ,cwava   ,qfcst   ,          &
                   etamid  ,ps      ,u3      ,v3      ,t3      , &
                   qminus  ,vort    ,div     ,hw2al   ,hw2bl   , &
                   hw3al   ,hw3bl   ,hwxal   ,hwxbl   ,q3m1    , &
                   grdps   ,grzs    ,grds    ,gruhs   ,grvhs   , &
                   grths   ,grpss   ,grus    ,grvs    ,grts    , &
                   grpls   ,grpms   ,grdpa   ,grza    ,grda    , &
                   gruha   ,grvha   ,grtha   ,grpsa   ,grua    , &
                   grva    ,grta    ,grpla   ,grpma   ,dps     , &
                   dpsl    ,dpsm    ,nlon)

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Transfrom variables from spherical harmonic coefficients 
! to grid point values during second gaussian latitude scan (scan2)
! 
! Method: 
! 
! Author: 
! Original version:  J. Rosinski
! Standardized:      J. Rosinski, June 1992
! Reviewed:          B. Boville, J. Hack, August 1992
! Reviewed:          B. Boville, April 1996
! 
!-----------------------------------------------------------------------
!
! $Id: spegrd.F90,v 1.12.4.1 2002/04/22 19:09:46 erik Exp $
! $Author: erik $
!
   use precision
   use pmgrid
   use pspect
   use commap
   use history, only: outfld
   use physconst, only: rga
   use constituents, only: pcnst, pnats
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comfft.h>
!---------------------------------------------------------------------
#include <comhd.h>
!---------------------------------------------------------------------
#include <comhyb.h>
!---------------------------------------------------------------------
#include <comlun.h>
!---------------------------------------------------------------------
#include <comqfl.h>
!---------------------------------------------------------------------
!
! Arguments
!
   real(r8), intent(in) :: ztodt              ! twice the timestep unles nstep=0
   real(r8), intent(in) :: cwava              ! normalization factor (1/g*plon)
   real(r8), intent(in) :: qfcst(plond,plev,pcnst)
   real(r8), intent(in) :: qminus(plond,plev,pcnst)
   real(r8), intent(in) :: etamid(plev)       ! vertical coords at midpoints

   real(r8), intent(inout) :: ps(plond)  
   real(r8), intent(inout) :: u3(plond,plev)  
   real(r8), intent(inout) :: v3(plond,plev)  
   real(r8), intent(inout) :: t3(plond,plev)  
   real(r8), intent(inout) :: vort(plond,plev)
   real(r8), intent(inout) :: div(plond,plev)
   real(r8), intent(inout) :: q3m1(plond,plev,pcnst+pnats)

   real(r8), intent(out) :: hw2al(pcnst)  ! -
   real(r8), intent(out) :: hw2bl(pcnst)  !  | lat contributions to components
   real(r8), intent(out) :: hw3al(pcnst)  !  | of slt global mass integrals 
   real(r8), intent(out) :: hw3bl(pcnst)  ! -
   real(r8), intent(out) :: hwxal(pcnst,4)
   real(r8), intent(out) :: hwxbl(pcnst,4)

   real(r8), intent(in) :: grdps(plond)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grzs(plond,plev)   ! sum(n) of z(n,m)*P(n,m) 
   real(r8), intent(in) :: grds(plond,plev)   ! sum(n) of d(n,m)*P(n,m)
   real(r8), intent(in) :: gruhs(plond,plev)  ! sum(n) of K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grvhs(plond,plev)  ! sum(n) of K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grths(plond,plev)  ! sum(n) of K(2i)*t(n,m)*P(n,m)
   real(r8), intent(in) :: grpss(plond)       ! sum(n) of lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grus(plond,plev)   ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grvs(plond,plev)   ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grts(plond,plev)   ! sum(n) of t(n,m)*P(n,m)
   real(r8), intent(in) :: grpls(plond)       ! sum(n) of lnps(n,m)*P(n,m)*m/a
   real(r8), intent(in) :: grpms(plond)       ! sum(n) of lnps(n,m)*H(n,m)
!
! Antisymmetric fourier coefficient arrays for all variables transformed
! from spherical harmonics (see grcalc)
!
   real(r8), intent(in) :: grdpa(plond)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grza(plond,plev)   ! sum(n) of z(n,m)*P(n,m)
   real(r8), intent(in) :: grda(plond,plev)   ! sum(n) of d(n,m)*P(n,m)
   real(r8), intent(in) :: gruha(plond,plev)  ! sum(n)K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grvha(plond,plev)  ! sum(n)K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grtha(plond,plev)  ! sum(n) of K(2i)*t(n,m)*P(n,m)
   real(r8), intent(in) :: grpsa(plond)       ! sum(n) of lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grua(plond,plev)   ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grva(plond,plev)   ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grta(plond,plev)   ! sum(n) of t(n,m)*P(n,m)
   real(r8), intent(in) :: grpla(plond)       ! sum(n) of lnps(n,m)*P(n,m)*m/a
   real(r8), intent(in) :: grpma(plond)       ! sum(n) of lnps(n,m)*H(n,m)

   real(r8), intent(out) :: dps(plond)
   real(r8), intent(out) :: dpsl(plond)
   real(r8), intent(out) :: dpsm(plond)

   integer, intent(in) :: lat             ! latitude index
   integer, intent(in) :: nlon            ! number of longitudes
!
!---------------------------Local workspace-----------------------------
!
   real(r8) :: duh(plond,plev) ! 
   real(r8) :: dvh(plond,plev) ! 
   real(r8) :: dth(plond,plev) ! 

   real(r8) pmid(plond,plev)   ! pressure at model levels
   real(r8) pint(plond,plevp)  ! pressure at model interfaces
   real(r8) pdel(plond,plev)   ! pdel(k) = pint(k+1) - pint(k)
   real(r8) pdelb(plond,plev)  ! pressure diff bet intfcs (press defined using the "B" part 
                               ! of the hybrid grid only)
   real(r8) qfcst1(plond,plev,pcnst) ! workspace to please lf95 compiler
   real(r8) hcwavaw            ! 0.5*cwava*w(lat)
   real(r8) sum
!
   real(r8) rcoslat            ! 1./cosine(latitude)
   real(r8) dotproda           ! dot product
   real(r8) dotprodb           ! dot product
#if ( ! defined USEFFTLIB )
   real(r8) work((plon+1)*plev)
#else
   real(r8) work((plon+1)*pcray) ! workspace needed by fft991
#endif
   integer isign           ! +1 => transform spectral to grid
   integer inc             ! distance between transform elements
   integer i,k,m           ! longitude, level, constituent indices
   integer klev            ! top level where hybrid coordinates apply
!
!-----------------------------------------------------------------------
!
   qfcst1(1:nlon,:,:) = qfcst(1:nlon,:,:)
!
! Set local integer pointers into part of buffer which is not written to SSD.
!
   inc = 1
   isign = +1
!
! Assemble northern and southern hemisphere grid values from the
! symmetric and antisymmetric fourier coefficients. 
! 1. Determine the fourier coefficients for the northern or southern
!    hemisphere latitude. 
! 2. Transform to gridpoint values
! 3. Clean up
!
   if (lat > plat/2) then                       ! Northern hemisphere
      do k=1,plev
         do i=1,nlon+2
            vort(i,k) = grzs(i,k) + grza(i,k)
            div(i,k) = grds(i,k) + grda(i,k)
            duh(i,k) = gruhs(i,k) + gruha(i,k)
            dvh(i,k) = grvhs(i,k) + grvha(i,k)
            dth(i,k) = grths(i,k) + grtha(i,k)
            u3(i,k) = grus(i,k) + grua(i,k)
            v3(i,k) = grvs(i,k) + grva(i,k)
            t3(i,k) = grts(i,k) + grta(i,k)
         end do
      end do

      do i=1,nlon+2
         dps(i) = grdps(i) + grdpa(i)
         ps(i) = grpss(i) + grpsa(i)
         dpsl(i) = grpls(i) + grpla(i)
         dpsm(i) = grpms(i) + grpma(i)
      end do

   else                                          ! Southern hemisphere

      do k=1,plev
         do i=1,nlon+2
            vort(i,k) = grzs(i,k) - grza(i,k)
            div(i,k) = grds(i,k) - grda(i,k)
            duh(i,k) = gruhs(i,k) - gruha(i,k)
            dvh(i,k) = grvhs(i,k) - grvha(i,k)
            dth(i,k) = grths(i,k) - grtha(i,k)
            u3(i,k) = grus(i,k) - grua(i,k)
            v3(i,k) = grvs(i,k) - grva(i,k)
            t3(i,k) = grts(i,k) - grta(i,k)
         end do
      end do

      do i=1,nlon+2
         dps(i) = grdps(i) - grdpa(i)
         ps(i) = grpss(i) - grpsa(i)
         dpsl(i) = grpls(i) - grpla(i)
         dpsm(i) = grpms(i) - grpma(i)
      end do
   end if
!
! Transform from fourier coefficients to gridpoint values.
! ps,vort,div,duh,dvh,dth,dpsl,dpsm,dps
!
   call fft991 (ps,   work, trig(1,lat), ifax(1,lat), inc, plond, nlon, 1,    isign)
   call fft991 (vort, work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
   call fft991 (div,  work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
   call fft991 (duh,  work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
   call fft991 (dvh,  work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
   call fft991 (dth,  work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
   call fft991 (dpsl, work, trig(1,lat), ifax(1,lat), inc, plond, nlon, 1,    isign)
   call fft991 (dpsm, work, trig(1,lat), ifax(1,lat), inc, plond, nlon, 1,    isign)
   call fft991 (dps,  work, trig(1,lat), ifax(1,lat), inc, plond, nlon, 1,    isign)
!
! u,v,t (SLT) [If you want to do spectral transport, do q as well]
!
   call fft991 (u3, work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
   call fft991 (v3, work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
   call fft991 (t3, work, trig(1,lat), ifax(1,lat), inc, plond, nlon, plev, isign)
!
! Remove cosine(latitude) from momentum variables
!
   rcoslat = 1./cos(clat(lat))
   do k=1,plev
      do i=1,nlon
         u3(i,k) = u3(i,k)*rcoslat
         v3(i,k) = v3(i,k)*rcoslat
         duh(i,k) = duh(i,k)*rcoslat
         dvh(i,k) = dvh(i,k)*rcoslat
      end do
   end do
!
! Copy transformed surface pressure back to in-core array, converting from 
! log(ps) to ps.
!
   do i=1,nlon
      ps(i) = exp(ps(i))
   end do
!
! Diagnose pressure arrays needed by DIFCOR
!
   call plevs0 (nlon, plond, plev, ps, pint, pmid, pdel)
   call pdelb0 (ps, pdelb, nlon)
!
! Accumulate mass integrals
!
   sum = 0.
   do i=1,nlon
      sum = sum + ps(i)
   end do
   tmass(lat) = w(lat)*rga*sum/nlon
!
! Finish horizontal diffusion: add pressure surface correction term to t and
! q diffusions; add kinetic energy dissipation to internal energy (temperature)
!
   klev = max(kmnhd4,nprlev)
   call difcor (klev,        ztodt,  dps,    u3,     v3, &
                q3m1(1,1,1), pdel,   pint,   t3,     dth, &
                duh,         dvh,    nlon)
!
! Calculate SLT moisture and constituent integrals
!
   hcwavaw = 0.5*cwava*w(lat)
   do m=1,pcnst
      hw2al(m) = 0.
      hw2bl(m) = 0.
      hw3al(m) = 0.
      hw3bl(m) = 0.
      hwxal(m,1) = 0.
      hwxal(m,2) = 0.
      hwxal(m,3) = 0.
      hwxal(m,4) = 0.
      hwxbl(m,1) = 0.
      hwxbl(m,2) = 0.
      hwxbl(m,3) = 0.
      hwxbl(m,4) = 0.
      do k=1,plev
         dotproda = 0.
         dotprodb = 0.
         do i=1,nlon
            dotproda = dotproda + qfcst1(i,k,m)*pdela(i,k)
            dotprodb = dotprodb + qfcst1(i,k,m)*pdelb(i,k)
         end do
         hw2al(m) = hw2al(m) + hcwavaw*dotproda
         hw2bl(m) = hw2bl(m) + hcwavaw*dotprodb
      end do
   end do

   call qmassd (cwava, etamid, w(lat), qminus, qfcst1, &
                pdela, hw3al, nlon)

   call qmassd (cwava, etamid, w(lat), qminus, qfcst1, &
                pdelb, hw3bl, nlon)

   if (pcnst.gt.1) then
      call xqmass (cwava, etamid, w(lat), qminus, qfcst1, &
                   qminus, qfcst1, pdela, pdelb, hwxal, &
                   hwxbl, nlon)
   end if

   call outfld ('DTH     ',dth     ,plond   ,lat     )

   return
end subroutine spegrd