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

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Transfrom variables from spherical harmonic coefficients 
! to grid point values during second gaussian latitude scan (scan2)
! 
! Method: 
! 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
!
! Author: 
! Original version:  J. Rosinski
! Standardized:      J. Rosinski, June 1992
! Reviewed:          B. Boville, J. Hack, August 1992
! Reviewed:          B. Boville, April 1996
! Modified:          P. Worley, September 2002
! 
!-----------------------------------------------------------------------
!


subroutine spegrd_bft (lat     ,nlon_fft, & 2,3
                       grdps   ,grzs    ,grds    ,gruhs   ,grvhs   , &
                       grths   ,grpss   ,grus    ,grvs    ,grts    , &
                       grpls   ,grpms   ,grdpa   ,grza    ,grda    , &
                       gruha   ,grvha   ,grtha   ,grpsa   ,grua    , &
                       grva    ,grta    ,grpla   ,grpma   ,fftbuf  )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Preparation for transform of variables from spherical harmonic 
! coefficients to grid point values during second gaussian latitude scan 
! (scan2)
! 
! Original version:  J. Rosinski
! Standardized:      J. Rosinski, June 1992
! Reviewed:          B. Boville, J. Hack, August 1992
! Reviewed:          B. Boville, April 1996
! Modified:          P. Worley, September 2002
! 
!-----------------------------------------------------------------------
!
! $Id: spegrd.F90,v 1.12.4.9 2004/04/28 23:43:58 eaton Exp $
! $Author: eaton $
!
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
   use comspe, only: maxm, numm
!-----------------------------------------------------------------------
   implicit none
!---------------------------------------------------------------------
!
! Arguments
!
   integer, intent(in) :: lat                  ! latitude index
   integer, intent(in) :: nlon_fft             ! first dimension of FFT work array
!
! Symmetric fourier coefficient arrays for all variables transformed
! from spherical harmonics (see grcalc)
!
   real(r8), intent(in) :: grdps(2*maxm)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grzs(2*maxm,plev)   ! sum(n) of z(n,m)*P(n,m) 
   real(r8), intent(in) :: grds(2*maxm,plev)   ! sum(n) of d(n,m)*P(n,m)
   real(r8), intent(in) :: gruhs(2*maxm,plev)  ! sum(n) of K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grvhs(2*maxm,plev)  ! sum(n) of K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grths(2*maxm,plev)  ! sum(n) of K(2i)*t(n,m)*P(n,m)
   real(r8), intent(in) :: grpss(2*maxm)       ! sum(n) of lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grus(2*maxm,plev)   ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grvs(2*maxm,plev)   ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grts(2*maxm,plev)   ! sum(n) of t(n,m)*P(n,m)
   real(r8), intent(in) :: grpls(2*maxm)       ! sum(n) of lnps(n,m)*P(n,m)*m/a
   real(r8), intent(in) :: grpms(2*maxm)       ! 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(2*maxm)       ! sum(n) of K(4)*(n(n+1)/a**2)**2*2dt*lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grza(2*maxm,plev)   ! sum(n) of z(n,m)*P(n,m)
   real(r8), intent(in) :: grda(2*maxm,plev)   ! sum(n) of d(n,m)*P(n,m)
   real(r8), intent(in) :: gruha(2*maxm,plev)  ! sum(n)K(2i)*z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grvha(2*maxm,plev)  ! sum(n)K(2i)*d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grtha(2*maxm,plev)  ! sum(n) of K(2i)*t(n,m)*P(n,m)
   real(r8), intent(in) :: grpsa(2*maxm)       ! sum(n) of lnps(n,m)*P(n,m)
   real(r8), intent(in) :: grua(2*maxm,plev)   ! sum(n) of z(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grva(2*maxm,plev)   ! sum(n) of d(n,m)*H(n,m)*a/(n(n+1))
   real(r8), intent(in) :: grta(2*maxm,plev)   ! sum(n) of t(n,m)*P(n,m)
   real(r8), intent(in) :: grpla(2*maxm)       ! sum(n) of lnps(n,m)*P(n,m)*m/a
   real(r8), intent(in) :: grpma(2*maxm)       ! sum(n) of lnps(n,m)*H(n,m)

   real(r8), intent(out) :: fftbuf(nlon_fft,8,plevp) ! buffer used for in-place FFTs

!
!---------------------------Local workspace-----------------------------
!
   integer i,k                           ! longitude, level indices
   integer rmlength                      ! twice number of local wavenumbers
   integer, parameter :: vortdex = 1     ! indices into fftbuf 
   integer, parameter :: divdex = 2
   integer, parameter :: duhdex = 3
   integer, parameter :: dvhdex = 4
   integer, parameter :: dthdex = 5
   integer, parameter :: u3dex = 6
   integer, parameter :: v3dex = 7
   integer, parameter :: t3dex = 8
   integer, parameter :: dpsdex = 1
   integer, parameter :: psdex = 2
   integer, parameter :: dpsldex = 3
   integer, parameter :: dpsmdex = 4
!
!-----------------------------------------------------------------------
!
! Assemble northern and southern hemisphere grid values from the
! symmetric and antisymmetric fourier coefficients: pre-FFT
!
!DIR$ NOSTREAM
   rmlength = 2*numm(iam)
   if (lat > plat/2) then                       ! Northern hemisphere
!DIR$ PREFERVECTOR
      do k=1,plev
!cdir loopchg
         do i=1,rmlength
            fftbuf(i,vortdex,k) = grzs(i,k) + grza(i,k)
            fftbuf(i,divdex,k) = grds(i,k) + grda(i,k)
            fftbuf(i,duhdex,k) = gruhs(i,k) + gruha(i,k)
            fftbuf(i,dvhdex,k) = grvhs(i,k) + grvha(i,k)
            fftbuf(i,dthdex,k) = grths(i,k) + grtha(i,k)
            fftbuf(i,u3dex,k) = grus(i,k) + grua(i,k)
            fftbuf(i,v3dex,k) = grvs(i,k) + grva(i,k)
            fftbuf(i,t3dex,k) = grts(i,k) + grta(i,k)
         end do
      end do
!
!cdir altcode=(loopcnt)
      do i=1,rmlength
         fftbuf(i,dpsdex,plevp) = grdps(i) + grdpa(i)
         fftbuf(i,psdex,plevp) = grpss(i) + grpsa(i)
         fftbuf(i,dpsldex,plevp) = grpls(i) + grpla(i)
         fftbuf(i,dpsmdex,plevp) = grpms(i) + grpma(i)
      end do

   else                                          ! Southern hemisphere

!DIR$ PREFERVECTOR
      do k=1,plev
!cdir loopchg
         do i=1,rmlength
            fftbuf(i,vortdex,k) = grzs(i,k) - grza(i,k)
            fftbuf(i,divdex,k) = grds(i,k) - grda(i,k)
            fftbuf(i,duhdex,k) = gruhs(i,k) - gruha(i,k)
            fftbuf(i,dvhdex,k) = grvhs(i,k) - grvha(i,k)
            fftbuf(i,dthdex,k) = grths(i,k) - grtha(i,k)
            fftbuf(i,u3dex,k) = grus(i,k) - grua(i,k)
            fftbuf(i,v3dex,k) = grvs(i,k) - grva(i,k)
            fftbuf(i,t3dex,k) = grts(i,k) - grta(i,k)
         end do
      end do

!cdir altcode=(loopcnt)
      do i=1,rmlength
         fftbuf(i,dpsdex,plevp) = grdps(i) - grdpa(i)
         fftbuf(i,psdex,plevp) = grpss(i) - grpsa(i)
         fftbuf(i,dpsldex,plevp) = grpls(i) - grpla(i)
         fftbuf(i,dpsmdex,plevp) = grpms(i) - grpma(i)
      end do

   end if

   return
end subroutine spegrd_bft


subroutine spegrd_ift (nlon_fft_in, nlon_fft_out, fftbuf_in, fftbuf_out) 1,6

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Inverse Fourier transform of variables from spherical harmonic 
! coefficients to grid point values during second gaussian latitude scan 
! (scan2)
! 
! Author:             P. Worley, September 2002
! 
!-----------------------------------------------------------------------
!
! $Id: spegrd.F90,v 1.12.4.9 2004/04/28 23:43:58 eaton Exp $
! $Author: eaton $
!
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
   use rgrid
   use comspe, only: maxm
#if ( defined SPMD ) && ( defined TIMING_BARRIERS )
   use mpishorthand
#endif
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
#include <comfft.h>
!---------------------------------------------------------------------
!
! Arguments
!
!     
! Input arguments
!     
   integer, intent(in) :: nlon_fft_in      ! first dimension of first FFT work array
   integer, intent(in) :: nlon_fft_out     ! first dimension of second FFT work array
#if (defined SPMD)
   real(r8), intent(in) :: fftbuf_in(nlon_fft_in,8,plevp,plat) 
                            ! buffer containing fields dcomposed over wavenumbers
#else
   real(r8), intent(in) :: fftbuf_in(1,1,1,1) 
                            ! buffer unused
#endif
!     
! Input/Output arguments
!     
   real(r8), intent(inout) :: fftbuf_out(nlon_fft_out,8,plevp,beglat:endlat) 
                            ! buffer used for in-place FFTs
!
!---------------------------Local workspace-----------------------------
!
#if ( ! defined USEFFTLIB )
   real(r8) work((plon+1)*8*plevp)
#else
   real(r8) work((plon+1)*pcray) ! workspace needed by fft991
#endif
   integer lat               ! latitude index
   integer isign           ! +1 => transform spectral to grid
   integer ntr             ! number of transforms to perform
   integer inc             ! distance between transform elements
   integer begtrm          ! (real) location of first truncated wavenumber
!
!-----------------------------------------------------------------------
!
!
#if ( defined SPMD )
!
!  reorder Fourier coefficients
!
#if ( defined TIMING_BARRIERS )
   call t_startf ('sync_realloc4b')
   call mpibarrier (mpicom)
   call t_stopf ('sync_realloc4b')
#endif
   call t_startf('realloc4b')
   call realloc4b(nlon_fft_in, nlon_fft_out, fftbuf_in, fftbuf_out)
   call t_stopf('realloc4b')
#endif
!
! Zero elements corresponding to truncated wavenumbers, then
! transform from fourier coefficients to gridpoint values.
! ps,vort,div,duh,dvh,dth,dpsl,dpsm,dps,
! u,v,t (SLT) [If you want to do spectral transport, do q as well]
!
   begtrm = 2*pmmax+1
   inc = 1
   isign = +1
   ntr = 8*plev + 4
   fftbuf_out(begtrm:plond,:,:,:) = 0.0
!$OMP PARALLEL DO PRIVATE (LAT,WORK)
   do lat=beglat,endlat
      call fft991 (fftbuf_out(1,1,1,lat), work, trig(1,lat), ifax(1,lat), inc, &
                   nlon_fft_out, nlon(lat), ntr, isign)
   enddo
!
   return
end subroutine spegrd_ift


subroutine spegrd_aft (ztodt   ,lat     ,nlon    ,nlon_fft, & 1,19
                       cwava   ,qfcst   , &
                       etamid  ,ps      ,u3      ,v3      ,t3      , &
                       qminus  ,vort    ,div     ,hw2al   ,hw2bl   , &
                       hw3al   ,hw3bl   ,hwxal   ,hwxbl   ,q3m1    , &
                       dps     ,dpsl    ,dpsm    ,t3m2    ,engy2alat, &
                       engy2blat,difftalat, difftblat,phis,fftbuf  )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Completion of transformation of variables from spherical harmonic 
! coefficients to grid point values during second gaussian latitude scan 
! (scan2)
! 
! Method: 
! 
! Original version:  J. Rosinski
! Standardized:      J. Rosinski, June 1992
! Reviewed:          B. Boville, J. Hack, August 1992
! Reviewed:          B. Boville, April 1996
! Modified:          P. Worley, September 2002
! 
!-----------------------------------------------------------------------
!
! $Id: spegrd.F90,v 1.12.4.9 2004/04/28 23:43:58 eaton Exp $
! $Author: eaton $
!
   use shr_kind_mod, only: r8 => shr_kind_r8
   use pmgrid
   use pspect
   use commap
   use history, only: outfld
   use physconst, only: rga
   use constituents, only: pcnst, pnats
   use comhd
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
#include <comctl.h>
!-----------------------------------------------------------------------
#include <comfft.h>
!---------------------------------------------------------------------
#include <comhyb.h>
!---------------------------------------------------------------------
#include <comlun.h>
!---------------------------------------------------------------------
#include <comqfl.h>
!---------------------------------------------------------------------
!
! Arguments
!
   integer, intent(in) :: lat                 ! latitude index
   integer, intent(in) :: nlon                ! number of longitudes
   integer, intent(in) :: nlon_fft            ! first dimension of FFT work arrays

   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(out) :: dps(plond)
   real(r8), intent(out) :: dpsl(plond)
   real(r8), intent(out) :: dpsm(plond)
   real(r8), intent(in)  :: t3m2(plond,plev) ! temperature
   real(r8), intent(out) :: engy2alat
   real(r8), intent(out) :: engy2blat
   real(r8), intent(out) :: difftalat
   real(r8), intent(out) :: difftblat
   real(r8), intent(in)  :: phis(plond)
   real(r8), intent(in) :: fftbuf(nlon_fft,8,plevp) ! buffer used for in-place FFTs
!
!---------------------------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
   integer i,k,m               ! longitude, level, constituent indices
   integer klev                ! top level where hybrid coordinates apply
   integer, parameter :: vortdex = 1     ! indices into fftbuf 
   integer, parameter :: divdex = 2
   integer, parameter :: duhdex = 3
   integer, parameter :: dvhdex = 4
   integer, parameter :: dthdex = 5
   integer, parameter :: u3dex = 6
   integer, parameter :: v3dex = 7
   integer, parameter :: t3dex = 8
   integer, parameter :: dpsdex = 1
   integer, parameter :: psdex = 2
   integer, parameter :: dpsldex = 3
   integer, parameter :: dpsmdex = 4
!
!-----------------------------------------------------------------------
!
   qfcst1(1:nlon,:,:) = qfcst(i1:nlon+i1-1,:,:)
!
! Copy 3D fields out of FFT buffer, removing cosine(latitude) from momentum variables
!
   rcoslat = 1./cos(clat(lat))
   do k=1,plev
      do i=1,nlon
         vort(i,k) = fftbuf(i,vortdex,k)
         div(i,k) = fftbuf(i,divdex,k)
         duh(i,k) = fftbuf(i,duhdex,k)*rcoslat
         dvh(i,k) = fftbuf(i,dvhdex,k)*rcoslat
         dth(i,k) = fftbuf(i,dthdex,k)
         u3(i,k) = fftbuf(i,u3dex,k)*rcoslat
         v3(i,k) = fftbuf(i,v3dex,k)*rcoslat
         t3(i,k) = fftbuf(i,t3dex,k)
      end do
   end do
!
! Copy 2D fields out of FFT buffer, converting
! log(ps) to ps.
!
   do i=1,nlon
      dps(i) = fftbuf(i,dpsdex,plevp)
      dpsl(i) = fftbuf(i,dpsldex,plevp)
      dpsm(i) = fftbuf(i,dpsmdex,plevp)
      ps(i) = exp(fftbuf(i,psdex,plevp))
   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, constituent, energy, and temperature integrals
!
   hcwavaw   = 0.5*cwava*w(lat)
   engy2alat = 0.
   engy2blat = 0.
   difftalat = 0.
   difftblat = 0.
   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 engy_te  (cwava ,w(lat) ,t3  ,u3  ,v3 ,phis    ,pdela, engy2alat ,nlon)
   call engy_te  (cwava ,w(lat) ,t3  ,u3  ,v3 ,phis    ,pdelb, engy2blat ,nlon)
   call engy_tdif(cwava ,w(lat) ,t3  ,t3m2             ,pdela, difftalat ,nlon)
   call engy_tdif(cwava ,w(lat) ,t3  ,t3m2             ,pdelb, difftblat ,nlon)

   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_aft