dyn.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.
! Make sure to make appropriate coding changes where necessary.
#if ( defined PVP )

subroutine dyn(irow    ,grlps1  ,grt1    ,grz1    ,grd1    ,  &
               grfu1   ,grfv1   ,grut1   ,grvt1   ,grrh1   ,  &
               grlps2  ,grt2    ,grz2    ,grd2    ,grfu2   ,  &
               grfv2   ,grut2   ,grvt2   ,grrh2   )

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Combine undifferentiated and longitudinally differentiated Fourier
! coefficient terms for later use in the Gaussian quadrature
!
! Method: 
! Computational note: Index "2*m-1" refers to the real part of the
! complex coefficient, and "2*m" to the imaginary.
!
! The naming convention is as follows:
!  - t, q, d, z refer to temperature, specific humidity, divergence
!     and vorticity
!  - "1" suffix to an array => symmetric component of current latitude pair
!  - "2" suffix to an array => antisymmetric component
!
! Author: 
! Original version:  J. Rosinski
! Standardized:      J. Rosinski, June 1992
! Reviewed:          D. Williamson, B. Boville, J. Hack, August 1992
! Reviewed:          D. Williamson, March 1996
! Reviewed:          B. Boville, April 1996
!
!-----------------------------------------------------------------------
!
! $Id: dyn.F90,v 1.4 2001/10/19 17:50:31 eaton Exp $
! $Author: eaton $
!
!-----------------------------------------------------------------------
   use precision
   use pmgrid
   use pspect
   use rgrid
   use commap
   use dynconst, only: rearth
   use time_manager, only: get_step_size, is_first_step
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------
!
! Input arguments
!
   integer, intent(in) :: irow                ! latitude pair index
!
! Input/output arguments
!
   real(r8), intent(inout) :: grlps1(2*pmmax)        ! sym. surface pressure equation term
   real(r8), intent(inout) :: grt1(2*pmmax,plev)     ! sym. undifferentiated term in t eqn.
   real(r8), intent(inout) :: grz1(2*pmmax,plev)     ! sym. undifferentiated term in z eqn.
   real(r8), intent(inout) :: grd1(2*pmmax,plev)     ! sym. undifferentiated term in d eqn.
   real(r8), intent(inout) :: grfu1(2*pmmax,plev)    ! sym. nonlinear terms in u eqn.
   real(r8), intent(inout) :: grfv1(2*pmmax,plev)    ! sym. nonlinear terms in v eqn.
   real(r8), intent(inout) :: grut1(2*pmmax,plev)    ! sym. lambda derivative term in t eqn.
   real(r8), intent(inout) :: grvt1(2*pmmax,plev)    ! sym. mu derivative term in t eqn.
   real(r8), intent(inout) :: grrh1(2*pmmax,plev)    ! sym. RHS of divergence eqn (del^2 term)
   real(r8), intent(inout) :: grlps2(2*pmmax)        ! antisym. surface pressure equation term
   real(r8), intent(inout) :: grt2(2*pmmax,plev)     ! antisym. undifferentiated term in t eqn.
   real(r8), intent(inout) :: grz2(2*pmmax,plev)     ! antisym. undifferentiated term in z eqn.
   real(r8), intent(inout) :: grd2(2*pmmax,plev)     ! antisym. undifferentiated term in d eqn.
   real(r8), intent(inout) :: grfu2(2*pmmax,plev)    ! antisym. nonlinear terms in u eqn.
   real(r8), intent(inout) :: grfv2(2*pmmax,plev)    ! antisym. nonlinear terms in v eqn.
   real(r8), intent(inout) :: grut2(2*pmmax,plev)    ! antisym. lambda derivative term in t eqn.
   real(r8), intent(inout) :: grvt2(2*pmmax,plev)    ! antisym. mu derivative term in t eqn.
   real(r8), intent(inout) :: grrh2(2*pmmax,plev)    ! antisym. RHS of divergence eqn (del^2 term)
!
!---------------------------Local workspace-----------------------------
!
   real(r8) tmp1,tmp2              ! temporaries
   real(r8) zxm(pmmax)             ! m*2dt/(a*cos(lat)**2)
   real(r8) zrcsj                  ! 1./(a*cos(lat)**2)
   real(r8) dtime                  ! timestep size [seconds]
   real(r8) ztdtrc                 ! 2dt/(a*cos(lat)**2), 1dt/..... at nstep=0
   integer m                   ! Fourier wavenumber index
   integer k                   ! level index
!
   do m=1,2*nmmax(irow)
      tmp1 = 0.5*(grlps2(m) + grlps1(m))
      tmp2 = 0.5*(grlps2(m) - grlps1(m))
      grlps1(m) = tmp1
      grlps2(m) = tmp2
   end do
   do k=1,plev
      do m=1,2*nmmax(irow)
         tmp1 = 0.5*(grt2(m,k) + grt1(m,k))
         tmp2 = 0.5*(grt2(m,k) - grt1(m,k))
         grt1(m,k) = tmp1
         grt2(m,k) = tmp2
!
         tmp1 = 0.5*(grz2(m,k) + grz1(m,k))
         tmp2 = 0.5*(grz2(m,k) - grz1(m,k))
         grz1(m,k) = tmp1
         grz2(m,k) = tmp2
!
         tmp1 = 0.5*(grd2(m,k) + grd1(m,k))
         tmp2 = 0.5*(grd2(m,k) - grd1(m,k))
         grd1(m,k) = tmp1
         grd2(m,k) = tmp2
!
         tmp1 = 0.5*(grfu2(m,k) + grfu1(m,k))
         tmp2 = 0.5*(grfu2(m,k) - grfu1(m,k))
         grfu1(m,k) = tmp1
         grfu2(m,k) = tmp2
!
         tmp1 = 0.5*(grfv2(m,k) + grfv1(m,k))
         tmp2 = 0.5*(grfv2(m,k) - grfv1(m,k))
         grfv1(m,k) = tmp1
         grfv2(m,k) = tmp2
!
         tmp1 = 0.5*(grut2(m,k) + grut1(m,k))
         tmp2 = 0.5*(grut2(m,k) - grut1(m,k))
         grut1(m,k) = tmp1
         grut2(m,k) = tmp2
!
         tmp1 = 0.5*(grvt2(m,k) + grvt1(m,k))
         tmp2 = 0.5*(grvt2(m,k) - grvt1(m,k))
         grvt1(m,k) = tmp1
         grvt2(m,k) = tmp2
!
         tmp1 = 0.5*(grrh2(m,k) + grrh1(m,k))
         tmp2 = 0.5*(grrh2(m,k) - grrh1(m,k))
         grrh1(m,k) = tmp1
         grrh2(m,k) = tmp2
      end do
   end do
!
! Set constants
!
   dtime = get_step_size()
   zrcsj = 1./(cs(irow)*rearth)
   if (is_first_step()) then
      ztdtrc = dtime*zrcsj
   else
      ztdtrc = 2.0*dtime*zrcsj
   end if
!
! Combine constants with Fourier wavenumber m
!
   do m=1,nmmax(irow)
      zxm(m) = ztdtrc*xm(m)
   end do
!
! Combine undifferentiated and longitudinal derivative terms for
! later use in Gaussian quadrature
!
   do k=1,plev
      do m=1,nmmax(irow)
         grt1(2*m-1,k) = grt1(2*m-1,k) + zxm(m)*grut1(2*m,k)
         grt1(2*m,k)   = grt1(2*m,k)   - zxm(m)*grut1(2*m-1,k)
         grd1(2*m-1,k) = grd1(2*m-1,k) - zxm(m)*grfu1(2*m,k)
         grd1(2*m,k)   = grd1(2*m,k)   + zxm(m)*grfu1(2*m-1,k)
         grz1(2*m-1,k) = grz1(2*m-1,k) - zxm(m)*grfv1(2*m,k)
         grz1(2*m,k)   = grz1(2*m,k)   + zxm(m)*grfv1(2*m-1,k)
!
         grt2(2*m-1,k) = grt2(2*m-1,k) + zxm(m)*grut2(2*m,k)
         grt2(2*m,k)   = grt2(2*m,k)   - zxm(m)*grut2(2*m-1,k)
         grd2(2*m-1,k) = grd2(2*m-1,k) - zxm(m)*grfu2(2*m,k)
         grd2(2*m,k)   = grd2(2*m,k)   + zxm(m)*grfu2(2*m-1,k)
         grz2(2*m-1,k) = grz2(2*m-1,k) - zxm(m)*grfv2(2*m,k)
         grz2(2*m,k)   = grz2(2*m,k)   + zxm(m)*grfv2(2*m-1,k)
      end do
   end do
   return

#else

   subroutine dyn(irow    ,grlps1  ,grt1    ,grz1    ,grd1    ,  &
                  grfu1   ,grfv1   ,grut1   ,grvt1   ,grrh1   ,  &
                  grlps2  ,grt2    ,grz2    ,grd2    ,grfu2   ,  &
                  grfv2   ,grut2   ,grvt2   ,grrh2   )
!-----------------------------------------------------------------------
!
! Combine undifferentiated and longitudinally differentiated Fourier
! coefficient terms for later use in the Gaussian quadrature
!
! Computational note: Index "2*m-1" refers to the real part of the
! complex coefficient, and "2*m" to the imaginary.
!
! The naming convention is as follows:
!  - t, q, d, z refer to temperature, specific humidity, divergence
!     and vorticity
!  - "1" suffix to an array => symmetric component of current latitude pair
!  - "2" suffix to an array => antisymmetric component
!
!---------------------------Code history--------------------------------
!
! Original version:  J. Rosinski
! Standardized:      J. Rosinski, June 1992
! Reviewed:          D. Williamson, B. Boville, J. Hack, August 1992
! Reviewed:          D. Williamson, March 1996
!
!-----------------------------------------------------------------------
!
! $Id: dyn.F90,v 1.4 2001/10/19 17:50:31 eaton Exp $
! $Author: eaton $
!
!-----------------------------------------------------------------------
      use precision
      use pmgrid
      use pspect
      use rgrid
      use commap
      use dynconst, only: rearth
      use time_manager, only: get_step_size, is_first_step

      implicit none

!
! Input arguments
!
      integer irow                ! latitude pair index
!
! Input/output arguments
!
      real(r8) grlps1(2*pmmax)        ! sym. surface pressure equation term
      real(r8) grt1(plev,2*pmmax)     ! sym. undifferentiated term in t eqn.
      real(r8) grz1(plev,2*pmmax)     ! sym. undifferentiated term in z eqn.
      real(r8) grd1(plev,2*pmmax)     ! sym. undifferentiated term in d eqn.
      real(r8) grfu1(plev,2*pmmax)    ! sym. nonlinear terms in u eqn.
      real(r8) grfv1(plev,2*pmmax)    ! sym. nonlinear terms in v eqn.
      real(r8) grut1(plev,2*pmmax)    ! sym. lambda derivative term in t eqn.
      real(r8) grvt1(plev,2*pmmax)    ! sym. mu derivative term in t eqn.
      real(r8) grrh1(plev,2*pmmax)    ! sym. RHS of divergence eqn (del^2 term)
      real(r8) grlps2(2*pmmax)        ! antisym. surface pressure equation term
      real(r8) grt2(plev,2*pmmax)     ! antisym. undifferentiated term in t eqn.
      real(r8) grz2(plev,2*pmmax)     ! antisym. undifferentiated term in z eqn.
      real(r8) grd2(plev,2*pmmax)     ! antisym. undifferentiated term in d eqn.
      real(r8) grfu2(plev,2*pmmax)    ! antisym. nonlinear terms in u eqn.
      real(r8) grfv2(plev,2*pmmax)    ! antisym. nonlinear terms in v eqn.
      real(r8) grut2(plev,2*pmmax)    ! antisym. lambda derivative term in t eqn.
      real(r8) grvt2(plev,2*pmmax)    ! antisym. mu derivative term in t eqn.
      real(r8) grrh2(plev,2*pmmax)    ! antisym. RHS of divergence eqn (del^2 term)
!
!---------------------------Local workspace-----------------------------
!
      real(r8) tmp1,tmp2              ! temporaries
      real(r8) zxm(pmmax)             ! m*2dt/(a*cos(lat)**2)
      real(r8) zrcsj                  ! 1./(a*cos(lat)**2)
      real(r8) dtime                  ! timestep size [seconds]
      real(r8) ztdtrc                 ! 2dt/(a*cos(lat)**2)  1dt/..... at nstep=0
      integer m                   ! Fourier wavenumber index
      integer k                   ! level index
!
      do m=1,2*nmmax(irow)
         tmp1 = 0.5*(grlps2(m) + grlps1(m))
         tmp2 = 0.5*(grlps2(m) - grlps1(m))
         grlps1(m) = tmp1
         grlps2(m) = tmp2
      end do
      do m=1,2*nmmax(irow)
         do k=1,plev
            tmp1 = 0.5*(grt2(k,m) + grt1(k,m))
            tmp2 = 0.5*(grt2(k,m) - grt1(k,m))
            grt1(k,m) = tmp1
            grt2(k,m) = tmp2
!
            tmp1 = 0.5*(grz2(k,m) + grz1(k,m))
            tmp2 = 0.5*(grz2(k,m) - grz1(k,m))
            grz1(k,m) = tmp1
            grz2(k,m) = tmp2
!
            tmp1 = 0.5*(grd2(k,m) + grd1(k,m))
            tmp2 = 0.5*(grd2(k,m) - grd1(k,m))
            grd1(k,m) = tmp1
            grd2(k,m) = tmp2
!
            tmp1 = 0.5*(grfu2(k,m) + grfu1(k,m))
            tmp2 = 0.5*(grfu2(k,m) - grfu1(k,m))
            grfu1(k,m) = tmp1
            grfu2(k,m) = tmp2
!
            tmp1 = 0.5*(grfv2(k,m) + grfv1(k,m))
            tmp2 = 0.5*(grfv2(k,m) - grfv1(k,m))
            grfv1(k,m) = tmp1
            grfv2(k,m) = tmp2
!
            tmp1 = 0.5*(grut2(k,m) + grut1(k,m))
            tmp2 = 0.5*(grut2(k,m) - grut1(k,m))
            grut1(k,m) = tmp1
            grut2(k,m) = tmp2
!
            tmp1 = 0.5*(grvt2(k,m) + grvt1(k,m))
            tmp2 = 0.5*(grvt2(k,m) - grvt1(k,m))
            grvt1(k,m) = tmp1
            grvt2(k,m) = tmp2
!
            tmp1 = 0.5*(grrh2(k,m) + grrh1(k,m))
            tmp2 = 0.5*(grrh2(k,m) - grrh1(k,m))
            grrh1(k,m) = tmp1
            grrh2(k,m) = tmp2
         end do
      end do
!
! Set constants
!
      dtime = get_step_size()
      zrcsj = 1./(cs(irow)*rearth)
      if (is_first_step()) then
         ztdtrc = dtime*zrcsj
      else
         ztdtrc = 2.0*dtime*zrcsj
      end if
!
! Combine constants with Fourier wavenumber m
!
      do m=1,nmmax(irow)
         zxm(m) = ztdtrc*xm(m)
      end do
!
! Combine undifferentiated and longitudinal derivative terms for
! later use in Gaussian quadrature
!
      do m=1,nmmax(irow)
         do k=1,plev
            grt1(k,2*m-1) = grt1(k,2*m-1) + zxm(m)*grut1(k,2*m)
            grt1(k,2*m)   = grt1(k,2*m)   - zxm(m)*grut1(k,2*m-1)
            grd1(k,2*m-1) = grd1(k,2*m-1) - zxm(m)*grfu1(k,2*m)
            grd1(k,2*m)   = grd1(k,2*m)   + zxm(m)*grfu1(k,2*m-1)
            grz1(k,2*m-1) = grz1(k,2*m-1) - zxm(m)*grfv1(k,2*m)
            grz1(k,2*m)   = grz1(k,2*m)   + zxm(m)*grfv1(k,2*m-1)
!
            grt2(k,2*m-1) = grt2(k,2*m-1) + zxm(m)*grut2(k,2*m)
            grt2(k,2*m)   = grt2(k,2*m)   - zxm(m)*grut2(k,2*m-1)
            grd2(k,2*m-1) = grd2(k,2*m-1) - zxm(m)*grfu2(k,2*m)
            grd2(k,2*m)   = grd2(k,2*m)   + zxm(m)*grfu2(k,2*m-1)
            grz2(k,2*m-1) = grz2(k,2*m-1) - zxm(m)*grfv2(k,2*m)
            grz2(k,2*m)   = grz2(k,2*m)   + zxm(m)*grfv2(k,2*m-1)
         end do
      end do

      return
#endif
   end subroutine dyn