#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))!