#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     )   returnend subroutine spegrd