linemsdyn.f90

CCSM Research Tools: Community Atmosphere Model (CAM)

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

subroutine linemsdyn(lat     ,psm1    ,psm2    ,u3m1    ,          &
                     u3m2    ,v3m1    ,v3m2    ,t3m1    ,t3m2    , &
                     q3m1    ,etadot  ,etamid  ,                   &
                     ztodt   ,grlps1  ,grt1    ,grz1    ,grd1    , &
                     grfu1   ,grfv1   ,grut1   ,grvt1   ,grrh1   , &
                     grlps2  ,grt2    ,grz2    ,grd2    ,grfu2   , &
                     grfv2   ,grut2   ,grvt2   ,grrh2   ,          &
                     vcour   ,vmax    ,vmaxt   ,                   &
                     detam   ,t2      ,fu      ,fv      ,          &
                     divm1   ,vortm2  ,divm2   ,vortm1  ,phis    , &
                     dpsl    ,dpsm    ,omga    ,nlon    ,nmmax   )
!----------------------------------------------------------------------- 
! 
! Purpose: 
! Control non-linear dynamical terms, FFT and combine terms
! in preparation for Fourier -> spectral quadrature.
! 
! Method: 
! The naming convention is as follows:
!  - prefix gr contains grid point values before FFT and Fourier
!     coefficients after
!  - t, q, d, z and ps refer to temperature, specific humidity,
!     divergence, vorticity and surface pressure
!  - "1" suffix to an array => symmetric component current latitude pair
!  - "2" suffix to an array => antisymmetric component.
!
! Note: cpp token "PVP" has an "ifdef" associated with it in this routine
! because the Fourier coefficient arrays (gr...) require a different storage
! order to optimize performance on vector vs. message-passing architectures.
! 
! Author: 
! Original version:  CCM3
!
!-----------------------------------------------------------------------
!
! $Id: linemsdyn.F90,v 1.12 2001/10/19 17:50:32 eaton Exp $
! $Author: eaton $

   use precision
   use pmgrid
   use constituents, only: pcnst, pnats
   use pspect
   use commap
   use history, only: outfld
   use time_manager, only: get_step_size

   implicit none

#include <comctl.h>
#include <comfft.h>
#include <comhyb.h>
#include <comlun.h>
!
! Input arguments
!     
   integer lat               ! latitude index for S->N storage
   integer nlon
   integer nmmax

   real(r8), intent(in) :: psm1(plond)        ! surface pressure (time n)
   real(r8), intent(in) :: psm2(plond)      ! surface pressure (time n-1)
   real(r8), intent(in) :: u3m1(plond,plev)   ! u-wind (time n)
   real(r8), intent(in) :: u3m2(plond,plev) ! u-wind (time n-1)
   real(r8), intent(in) :: v3m1(plond,plev)   ! v-wind (time n)
   real(r8), intent(in) :: v3m2(plond,plev) ! v-wind (time n-1)
   real(r8), intent(in) :: t3m1(plond,plev)   ! temperature (time n)
   real(r8), intent(in) :: q3m1(plond,plev,pcnst+pnats)   ! constituent conc(time n: h2o first)
   real(r8), intent(inout) :: etadot(plon,plevp) ! vertical motion (3-d used by slt)
   real(r8), intent(in) :: etamid(plev)     ! midpoint values of eta (a+b)
   real(r8), intent(in) :: ztodt            ! 2*timestep unless nstep = 0
   real(r8), intent(in) :: detam(plev)      ! maximum Courant number in vert.
!     
! Input/Output arguments
!     
   real(r8), intent(inout) :: t2(plond,plev)   ! t tend
   real(r8), intent(inout) :: fu(plond,plev)   ! nonlinear term - u momentum eqn.
   real(r8), intent(inout) :: fv(plond,plev)   ! nonlinear term - v momentum eqn.
   real(r8), intent(inout) :: divm1(plond,plev)
   real(r8), intent(inout) :: vortm2(plond,plev)
   real(r8), intent(inout) :: divm2(plond,plev)
   real(r8), intent(inout) :: vortm1(plond,plev)
   real(r8), intent(inout) :: phis(plond)
   real(r8), intent(inout) :: dpsl(plond)
   real(r8), intent(inout) :: dpsm(plond)
   real(r8), intent(inout) :: omga(plond,plev)
   real(r8), intent(inout) :: t3m2(plond,plev) ! temperature (time n-1)
!     
! Output arguments
!     
   real(r8), intent(out) :: grlps1(2*pmmax)  ! sym. undiff. term in lnps eqn.
   real(r8), intent(out) :: grlps2(2*pmmax)  ! antisym undiff. term in lnps eqn.
#if ( defined PVP )
   real(r8), intent(out) :: grt1(2*pmmax,plev) ! sym. undiff. term in t eqn.
   real(r8), intent(out) :: grt2(2*pmmax,plev) ! antisym. undiff. term in t eqn.
   real(r8), intent(out) :: grz1(2*pmmax,plev) ! sym. undiff. term in z eqn.
   real(r8), intent(out) :: grz2(2*pmmax,plev) ! antisym. undiff. term in z eqn.
   real(r8), intent(out) :: grd1(2*pmmax,plev) ! sym. undiff. term in d eqn.
   real(r8), intent(out) :: grd2(2*pmmax,plev) ! antisym. undiff. term in d eqn.
   real(r8), intent(out) :: grfu1(2*pmmax,plev) ! sym. nonlinear terms in u eqn.
   real(r8), intent(out) :: grfu2(2*pmmax,plev) ! antisym. nonlinear terms in u eqn.
   real(r8), intent(out) :: grfv1(2*pmmax,plev) ! sym. nonlinear terms in v eqn.
   real(r8), intent(out) :: grfv2(2*pmmax,plev) ! antisym. nonlinear terms in v eqn.
   real(r8), intent(out) :: grut1(2*pmmax,plev) ! sym. lambda deriv. term in t eqn.
   real(r8), intent(out) :: grut2(2*pmmax,plev) ! antisym. lambda deriv. term in t eqn.
   real(r8), intent(out) :: grvt1(2*pmmax,plev) ! sym. mu derivative term in t eqn.
   real(r8), intent(out) :: grvt2(2*pmmax,plev) ! antisym. mu deriv. term in t eqn.
   real(r8), intent(out) :: grrh1(2*pmmax,plev) ! sym. del**2 term in d eqn.
   real(r8), intent(out) :: grrh2(2*pmmax,plev) ! antisym. del**2 term in d eqn.
#else 
   real(r8), intent(out) :: grt1(plev,2*pmmax) ! sym. undiff. term in t eqn.
   real(r8), intent(out) :: grt2(plev,2*pmmax) ! antisym. undiff. term in t eqn.
   real(r8), intent(out) :: grz1(plev,2*pmmax) ! sym. undiff. term in z eqn.
   real(r8), intent(out) :: grz2(plev,2*pmmax) ! antisym. undiff. term in z eqn.
   real(r8), intent(out) :: grd1(plev,2*pmmax) ! sym. undiff. term in d eqn.
   real(r8), intent(out) :: grd2(plev,2*pmmax) ! antisym. undiff. term in d eqn.
   real(r8), intent(out) :: grfu1(plev,2*pmmax) ! sym. nonlinear terms in u eqn.
   real(r8), intent(out) :: grfu2(plev,2*pmmax) ! antisym. nonlinear terms in u eqn.
   real(r8), intent(out) :: grfv1(plev,2*pmmax) ! sym. nonlinear terms in v eqn.
   real(r8), intent(out) :: grfv2(plev,2*pmmax) ! antisym. nonlinear terms in v eqn.
   real(r8), intent(out) :: grut1(plev,2*pmmax) ! sym. lambda deriv. term in t eqn.
   real(r8), intent(out) :: grut2(plev,2*pmmax) ! antisym. lambda deriv. term in t eqn.
   real(r8), intent(out) :: grvt1(plev,2*pmmax) ! sym. mu derivative term in t eqn.
   real(r8), intent(out) :: grvt2(plev,2*pmmax) ! antisym. mu deriv. term in t eqn.
   real(r8), intent(out) :: grrh1(plev,2*pmmax) ! sym. del**2 term in d eqn.
   real(r8), intent(out) :: grrh2(plev,2*pmmax) ! antisym. del**2 term in d eqn.
#endif
   real(r8), intent(out) :: vcour(plev)      ! maximum Courant number in vert.
   real(r8), intent(out) :: vmax(plev)       ! maximum wind speed squared (m^2/s^2)
   real(r8), intent(out) :: vmaxt(plev)      ! maximum truncated wind speed (m^2/s^2)
!     
!---------------------------Local workspace-----------------------------
!     
   real(r8) :: dtime          ! timestep size
   real(r8) :: bpstr(plond)   ! 
   real(r8) pmid(plond,plev)  ! pressure at model levels (time n)
   real(r8) rpmid(plond,plev) ! 1./pmid
   real(r8) pint(plond,plevp) ! pressure at model interfaces (n  )
   real(r8) pdel(plond,plev)  ! pdel(k)   = pint  (k+1)-pint  (k)
   real(r8) rpdel(plond,plev) ! 1./pdel
   real(r8) tdyn(plond,plev)   ! temperature for dynamics
   real(r8) vortdyn(plond,plev)
   real(r8) divdyn(plond,plev)
!     
! The "work" array has a different size requirement depending upon whether
! the proprietary Cray assembly language version of the FFT library
! routines, or the all-Fortran version, is being used.
!     
#if ( ! defined USEFFTLIB )
   real(r8) work((plon+1)*7*plev)
#else 
   real(r8) work((plon+1)*pcray) ! workspace array for fft991
#endif
   real(r8) logpsm1(plond)     ! log(psm1)
   real(r8) logpsm2(plond)   ! log(psm2)
   real(r8) engy(plond,plev) ! kinetic energy
!     region for following arrays to 
!     optimize fft performance)
   real(r8) ut(plond,plev)   ! (u*T) - heat flux - zonal
   real(r8) vt(plond,plev)   ! (v*T) - heat flux - meridional
   real(r8) drhs(plond,plev) ! RHS of divergence eqn. (del^2 term)
   real(r8) lvcour           ! local vertical courant number
   real(r8) dtdz             ! dt/detam(k)
   real(r8) ddivdt(plond,plev) ! temporary workspace
   real(r8) ddpn(plond)      ! complete sum of d*delta p
   real(r8) vpdsn(plond)     ! complete sum V dot grad(ln(ps)) delta b
   real(r8) dpslat(plond,plev) ! Pressure gradient term 
   real(r8) dpslon(plond,plev) ! Pressure gradient term 
   real(r8) coslat           ! cosine(latitude)
   real(r8) rcoslat          ! 1./cosine(latitude)
   real(r8) rhypi            ! 1./hypi(plevp)

   real(r8) wind             ! u**2 + v**2 (m/s)
   real(r8) utfac            ! asymmetric truncation factor for courant calculation
   real(r8) vtfac            ! asymmetric truncation factor for courant calculation

   integer i,k,kk            ! longitude,level,constituent indices
   integer inc               ! increment for fft991
   integer isign             ! flag indicates transform direction
   integer ntr               ! number of transforms to perform
!
! This group of arrays are glued together via equivalence to exbuf for
! communication from LINEMSBC.
!
!
!-----------------------------------------------------------------------
!
!
! Compute maximum wind speed this latitude (used in Courant number estimate)
!
   if (ptrm .lt. ptrn) then
      utfac = float(ptrm)/float(ptrn)
      vtfac = 1.
   else if (ptrn .lt. ptrm) then
      utfac = 1.
      vtfac = float(ptrn)/float(ptrm) 
   else if (ptrn .eq. ptrm) then
      utfac = 1.
      vtfac = 1.
   end if
   do k=1,plev
      vmax(k) = 0.
      vmaxt(k) = 0.
      do i=1,nlon
         wind = u3m2(i,k)**2 + v3m2(i,k)**2
         vmax(k) = max(wind,vmax(k))
!