tstep.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

#include <misc.h>
#include <params.h>
! Note that this routine has 2 complete blocks of code for PVP vs.
! non-PVP.  This is due to the fact that spectral coefficients are
! stored consecutively along diagonals of M-N wavenumber space when the
! target architecture is PVP (optimal for vectorization), and along
! total wavenumber N otherwise (optimal for message-passing).
#if ( defined PVP )
subroutine tstep(n       ,ztdtsq  )
!-----------------------------------------------------------------------
!
! Purpose:
! Solution of the vertically coupled system of equations arising
! from the semi-implicit equations for each spectral element along
! the n(th) diagonal. (Note, n is distinct from the two dimensional
! wavenumber which is also often denoted n.) The inverse matrix depends
! only on two dimensional wavenumber and the reference atmosphere.
! It is precomputed and stored for use during the forecast. The routine
! overwrites the d,T and lnps coefficients with the new values.
!
! Original version:  CCM1
!
!-----------------------------------------------------------------------
!
! $Id: tstep.F90,v 1.3 2000/12/20 18:02:17 rosinski Exp $
! $Author: rosinski $
!
!-----------------------------------------------------------------------

  use pmgrid
  use pspect
  use comspe
  use comslt, only: epssld
  implicit none
#include <comhyb.h>
use commap

!------------------------------Arguments--------------------------------
!
  integer , intent(in)   :: n               ! index of spectral diagonal being calculated
!                                           ! this call (not two dimensional wavenumber)
  real(r8), intent(in)   :: ztdtsq(2*pnmax) ! 2*dt*(n(n+1)/a^2 where n is 2-d wavenumber
!
!---------------------------Local workspace-----------------------------
!
  real(r8) hhref           ! href/2 (reference hydrostatic matrix / 2)
  real(r8) hbps            ! bps/2 (ref. coeff. for lnps term in div. eq. / 2)
  real(r8) onepeps         ! decentering coefficient
  integer ne               ! index into ztdtsq
  integer m                ! diagonal element (index) of complex array
  integer k,kk             ! level indices
  integer irh              ! index into levels of spectral arrays
  integer irhr,irhi        ! index into real, imaginary coefficients
  integer isp              ! index into spectral arrays
!
!-----------------------------------------------------------------------
!
! Complete rhs of helmholtz eq.
! Set offsets for beginning of diagonal being calculated this call
!
  isp = nco2(n) - 2
  ne = 2*(n-1)
  onepeps = 1. + epssld
!
  do k=1,plev
!
! Coefficients for diagonal terms
!
     hhref = onepeps*0.5*href(k,k)
     hbps  = onepeps*0.5*bps (k)
!
! Loop along current diagonal (in spectral space)
! Add lnps and diagonal (vertical space) T terms to initialize hs
!
     do m=1,2*nm(n)
        hs(isp+m,k) = ztdtsq(ne+m)*(hhref*t(isp+m,k) + hbps*alps(isp+m))
     end do
     if (k.lt.plev) then
        do kk=k+1,plev
!
! Add off-diagonal (vertical space) T terms to hs
!
           hhref = onepeps*0.5*href(kk,k)
           do m=1,2*nm(n)
              hs(isp+m,k) = hs(isp+m,k) + ztdtsq(ne+m)*hhref*t(isp+m,kk)
           end do
        end do
     end if
  end do                    ! k=1,plev (calculation level)
!
! Transform semi-implicit vectors to vertical normal mode space
!
  do k = 1,plev
     do m = 1,2*nm(n)
        dsnm(isp+m,k) = 0.
        hsnm(isp+m,k) = 0.
        vznm(isp+m,k) = 0.
        do kk = 1,plev
           dsnm(isp+m,k) = dsnm(isp+m,k) + bmi(kk,k)*d (isp+m,kk)
           hsnm(isp+m,k) = hsnm(isp+m,k) + bmi(kk,k)*hs(isp+m,kk)
           vznm(isp+m,k) = vznm(isp+m,k) + bmi(kk,k)*vz(isp+m,kk)
        end do
     end do
  end do
!
  return
#else
  subroutine tstep(m       ,ztdtsq  )
!-----------------------------------------------------------------------
!
! Purpose:
! Solution of the vertically coupled system of equations arising
! from the semi-implicit equations for each spectral element along
! two dimensional wavenumber n.  The inverse matrix depends
! only on two dimensional wavenumber and the reference atmosphere.
! It is precomputed and stored for use during the forecast. The routine
! overwrites the d,T and lnps coefficients with the new values.
!
! Original version:  CCM1
!
!-----------------------------------------------------------------------
!
! $Id: tstep.F90,v 1.3 2000/12/20 18:02:17 rosinski Exp $
! $Author: rosinski $
!
!-----------------------------------------------------------------------

    use precision
    use pmgrid
    use pspect
    use comspe
    use comslt, only: epssld
    use commap

    implicit none

#include <comhyb.h>

!------------------------------Arguments--------------------------------
!
! Input arguments
!
    integer , intent(in)   :: m             ! Fourier wavenumber               
    real(r8), intent(in)   :: ztdtsq(pnmax) ! 2*dt*(n(n+1)/a^2 where n is 2-d wavenumber
!
!---------------------------Local workspace-----------------------------
!
    real(r8) hhref           ! href/2 (reference hydrostatic matrix / 2)
    real(r8) hbps            ! bps/2 (ref. coeff. for lnps term in div. eq. / 2)
    real(r8) onepeps         ! decentering coefficient
    integer n                ! 2-d wavenumber index
    integer k,kk             ! level indices
    integer mr,mc            ! real and imaginary spectral indices
    integer ir,ii            ! real and imaginary spectral indices
!
!-----------------------------------------------------------------------
!
! Complete rhs of helmholtz eq.
!
    mr = nstart(m)
    mc = 2*mr
    onepeps = 1. + epssld
!
    do k=1,plev
!
! Coefficients for diagonal terms
!
       hhref = onepeps*0.5*href(k,k)
       hbps = onepeps*0.5*bps(k)
!
! Loop along total wavenumber index (in spectral space)
! Add lnps and diagonal (vertical space) T terms to initialize hs
!
       do n=1,nlen(m)
          ir = mc + 2*n - 1
          ii = ir + 1
          hs(ir,k) = ztdtsq(n+m-1)*(hhref*t(ir,k) + hbps*alps(ir))
          hs(ii,k) = ztdtsq(n+m-1)*(hhref*t(ii,k) + hbps*alps(ii))
       end do
       if (k.lt.plev) then
          do kk=k+1,plev
!
! Add off-diagonal (vertical space) T terms to hs
!
             hhref = onepeps*0.5*href(kk,k)
             do n=1,nlen(m)
                ir = mc + 2*n - 1
                ii = ir + 1
                hs(ir,k) = hs(ir,k) + ztdtsq(n+m-1)*hhref*t(ir,kk)
                hs(ii,k) = hs(ii,k) + ztdtsq(n+m-1)*hhref*t(ii,kk)
             end do
          end do
       end if
    end do                    ! k=1,plev (calculation level)
!
! Transform semi-implicit vectors to vertical normal mode space
!
    do k = 1,plev
       do n=1,nlen(m)
          ir = mc + 2*n - 1
          ii = ir + 1
          dsnm(ir,k) = 0.
          dsnm(ii,k) = 0.
          hsnm(ir,k) = 0.
          hsnm(ii,k) = 0.
          vznm(ir,k) = 0.
          vznm(ii,k) = 0.
          do kk = 1,plev
             dsnm(ir,k) = dsnm(ir,k) + bmi(kk,k)*d (ir,kk)
             dsnm(ii,k) = dsnm(ii,k) + bmi(kk,k)*d (ii,kk)
             hsnm(ir,k) = hsnm(ir,k) + bmi(kk,k)*hs(ir,kk)
             hsnm(ii,k) = hsnm(ii,k) + bmi(kk,k)*hs(ii,kk)
             vznm(ir,k) = vznm(ir,k) + bmi(kk,k)*vz(ir,kk)
             vznm(ii,k) = vznm(ii,k) + bmi(kk,k)*vz(ii,kk)
          end do
       end do
    end do
!
    return
#endif
  end subroutine tstep