!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE CORDINTG                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!EMK BMJ

SUBROUTINE cordintg(mptr,frstep, nx,ny,nz,nxndg,nyndg,nzndg,            & 3,69
           rbufsz,nstyps,exbcbufsz,                                     &
           u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth,        &
           udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,             &
           wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,             &
           sdteb,sdtwb,sdtnb,sdtsb,                                     &
           ubar,vbar,ptbar,pbar,ptbari,pbari,                           &
           rhostr,rhostri,qvbar,ppi,csndsq,                             &
           x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv,             &
           trigs1,trigs2,ifax1,ifax2,                                   &
           wsave1,wsave2,vwork1,vwork2,                                 &
           sinlat,kmh,kmv,rprntl,                                       &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,ptsfc,qvsfc,        &
           ptcumsrc,qcumsrc,raing,rainc,prcrate,w0avg,nca,kfraincv,     &
           cldefi,xland,bmjraincv,                                      &
           radfrc,radsw,rnflx,rad2d,radbuf,exbcbuf,bcrlx,               &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           uincr,vincr,wincr,pincr,ptincr,qvincr,                       &
           qcincr,qrincr,qiincr,qsincr,qhincr,                          &
           phydro,tem1,tem2,tem3,tem4,tem5,tem6,tem7,                   &
           tem8,tem9,tem10,tem11,tem12,tem13,                           &
           tem14,tem15,tem16,tem17,tem18,tem19,                         &
           tem20,tem21,tem22,tem23,tem24,tem25,tem26,                   &
           tem1_0,tem2_0,tem3_0)
!EMK END
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Coordinate the forward integration of all model time-dependent
!  variables.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/01/1991.
!
!  MODIFICATION HISTORY:
!
!  4/26/92 (K. Droegemeier)
!  Added full documentation.
!
!  4/10/93 (M. Xue & Hao Jin)
!  Add the terrain.
!
!  01/24/94 (Yuhe Liu)
!  Added array tsfc, qvsfc, tsoil, wetsfc, wetdp, wetcanp, and some of
!  temporary arrays to the arguments of subroutine SFCFLX, for
!  prediction of the surface temperature and specific humidity
!
!  9/10/94 (D. Weber & Y. Lu)
!  Cleaned up documentation.
!
!  02/07/1995 (Yuhe Liu)
!  Added a new 2-D array, veg(nx,ny), to the argument list
!
!  7/6/1995 (M. Xue)
!  Pressure detrending is applied to the base grid (as opposed to
!  nested grids) only.
!
!  10/31/95 (D. Weber)
!  Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
!  radiation condition.
!
!  01/25/96 (Donghai Wang & Yuhe Liu)
!  Added the map projection factor to ARPS governing equations.
!
!  01/31/1996 (V. Wong and X. Song)
!  Added a parameter, qpfgfrq, that controls the frequency of calling
!  the subroutine QPFGRID.
!
!  08/14/1996 (Yuhe Liu)
!  Moved the definition of radiation buffer to the ARPS main program
!  and passed it through argument list.
!
!  07/22/97 (D. Weber)
!  Added  wsave1,wsave2,vwork1,vwork2 for use in the upper w-p
!  radiation condition (fftopt=2).
!
!  08/01/97 (Zonghui Huo)
!  Added Kain-fritsch cumulus parameterization scheme.
!
!  10/21/97 (Donghai Wang)
!  Using total density (rho) in the calculation of the pressure
!  gradient force terms, and added the second order terms
!  in the linerized buoyancy terms.
!
!  11/05/97 (D. Weber)
!  Added phydro array for use in the bottom boundary condition for
!  perturbation pressure (hydrostatic).
!
!  11/06/97 (D. Weber)
!  Added three additional levels to the mapfct array.  The three
!  levels (4,5,6) represent the inverse of the first three in order.
!  The inverse map factors are computed to improve efficiency.
!
!  2/11/1998 (M. Xue)
!  Corrections made to ensure no moist process is activatived when
!  moist=0.
!
!  4/15/1998 (Donghai Wang)
!  Added the source terms to the right hand terms of the qc,qr,qi,qs
!  equations due to the K-F cumulus parameterization.
!
!  4/15/1998 (Donghai Wang)
!  Added the running average vertical velocity (array w0avg)
!  for the K-F cumulus parameterization scheme.
!
!  9/18/98 (D. Weber)
!  Added precomputed averages of j3 in the x,y, and z directions
!  to improve code efficiency.
!
!  8/31/1998 (K. Brewster)
!  Added call to NUDGEALL for nudging to observation increments.
!
!  11/18/98 (Keith Brewster)
!  Changed pibar to ppi (full pi).
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  12/05/2001 (Yunheng Wang)
!  Replaced sfcflx call with sfcphysics call which was a structure
!  change by M. Xue.
!
!  07/10/2001 (K. Brewster)
!  Added increment arrays to argument list and to call to nudgeall
!
!  07/23/2001 (K. Brewster)
!  Added mptr to call to tinteg, needed for printing of diagnostic noise
!  parameter.
!
!  13 March 2002 (Eric Kemp)
!  Added arrays for WRF BMJ cumulus scheme.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    mptr     Grid identifier.
!    frstep   Flag to determine if this is the initial time step.
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at times tpast and tpresent (m/s)
!    v        y component of velocity at times tpast and tpresent (m/s)
!    w        Vertical component of Cartesian velocity at times
!             tpast and tpresent (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!    ptprt    Perturbation potential temperature at times tpast and
!             tpresent (K)
!    pprt     Perturbation pressure at times tpast and tpresent (Pascal)
!    qv       Water vapor specific humidity at times tpast and tpresent (kg/kg)
!    qc       Cloud water mixing ratio at times tpast and tpresent (kg/kg)
!    qr       Rainwater mixing ratio at times tpast and tpresent (kg/kg)
!    qi       Cloud ice mixing ratio at times tpast and tpresent (kg/kg)
!    qs       Snow mixing ratio at times tpast and tpresent (kg/kg)
!    qh       Hail mixing ratio at times tpast and tpresent (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!    pbldpth  Planetary boundary layer depth (m)
!
!    udteb    Time tendency of u field at east boundary (m/s**2)
!    udtwb    Time tendency of u field at west boundary (m/s**2)
!
!    vdtnb    Time tendency of v field at north boundary (m/s**2)
!    vdtsb    Time tendency of v field at south boundary (m/s**2)
!
!    pdteb    Time tendency of pprt field at east boundary (Pascal/s)
!    pdtwb    Time tendency of pprt field at west boundary (Pascal/s)
!    pdtnb    Time tendency of pprt field at north boundary (Pascal/s)
!    pdtsb    Time tendency of pprt field at south boundary (Pascal/s)
!
!    phydro   Big time step forcing term for use in computing the
!             hydrostatic pressure at k=1.
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    ptbari   Inverse Base state potential temperature (K)
!    pbari    Inverse Base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    rhostri  Inverse base state density rhobar times j3 (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!    ppi      Exner function.
!    csndsq   Sound wave speed squared.
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!    mapfct   Map factors at scalar, u and v points
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    aj3x     Avgx of the coordinate transformation Jacobian  d(zp)/dz
!    aj3y     Avgy of the coordinate transformation Jacobian  d(zp)/dz
!    aj3z     Avgz of the coordinate transformation Jacobian  d(zp)/dz
!    j3inv    Inverse of the coordinate transformation Jacobian  d(zp)/dz
!
!    trigs1   Array containing pre-computed trig function for
!             fftopt=1.
!    trigs2   Array containing pre-computed trig function for
!             fftopt=1.
!    ifax1    Array containing the factors of nx for fftopt=1.
!    ifax2    Array containing the factors of ny for fftopt=1.
!
!    vwork1   2-D work array for fftopt=2.
!    vwork2   2-D work array for fftopt=2.
!    wsave1   Work array for fftopt=2.
!    wsave2   Work array for fftopt=2.
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Ground surface temperature (K)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    qvsfc    Effective qv at sfc.
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    ptsfc    Ground surface potential temperature (K)
!    qvsfc    Effective S.H. at sfc.
!
!  OUTPUT:
!    u        Updated (by dtbig) x velocity component at times tpast
!             and tpresent (m/s).
!    v        Updated (by dtbig) y velocity component at times tpast
!             and tpresent (m/s).
!    w        Updated (by dtbig) z velocity component at times tpast
!             and tpresent (m/s).
!    ptprt    Updated (by dtbig) perturbation potential temperature
!             at times tpast and tpresent (K).
!    pprt     Updated (by dtbig) perturbation pressure at times tpast
!             and tpresent (Pascal).
!    qv       Updated (by dtbig) water vapor specific humidity at times
!             tpast and tpresent (kg/kg).
!    qc       Updated (by dtbig) cloud water mixing ratio at times
!             tpast and tpresent (kg/kg).
!    qr       Updated (by dtbig) rainwater mixing ratio at times tpast
!             and tpresent (kg/kg).
!    qi       Updated (by dtbig) cloud ice mixing ratio at times tpast
!             and tpresent (kg/kg).
!    qs       Updated (by dtbig) snow mixing ratio at times tpast and
!             tpresent (kg/kg).
!    qh       Updated (by dtbig) hail mixing ratio at times tpast and
!             tpresent (kg/kg).
!
!    udteb    Time tendency of u field at east boundary (m/s**2)
!    udtwb    Time tendency of u field at west boundary (m/s**2)
!
!    vdtnb    Time tendency of v field at north boundary (m/s**2)
!    vdtsb    Time tendency of v field at south boundary (m/s**2)
!
!    pdteb    Time tendency of pprt field at east boundary (Pascal/s)
!    pdtwb    Time tendency of pprt field at west boundary (Pascal/s)
!    pdtnb    Time tendency of pprt field at north boundary (Pascal/s)
!    pdtsb    Time tendency of pprt field at south boundary (Pascal/s)
!
!    sinlat   Sin of latitude at each grid point
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!    rprntl   Reciprocal of Prandtl number
!
!    tsfc     Ground surface temperature (K)
!    qvsfc    Effective qv at sfc.
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    ptsfc    Ground surface potential temperature (K)
!    qvsfc    Effective S.H. at sfc.
!
!    ptcumsrc Source term in pt-equation due to cumulus parameterization
!    qcumsrc Source term in water equations due to cumulus parameterization
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    kfraincv   K-F convective rainfall (cm)
!    nca      K-F counter for CAPE release
!    cldefi   BMJ cloud efficiency
!    xland    BMJ land/sea mask
!    bmjraincv   BMJ convective rainfall (cm)
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the ground surface (W/m**2)
!    rnflx    Net radiation flux at the ground surface (W/m**2)
!    rad2d    2-D arrays for radiation calculation
!    radbuf   Buffer array to carry temporary working arrays for
!             radiation calculation
!
!  WORK ARRAYS:
!
!    udtnb    Time tendency of u field at north boundary (m/s**2)
!    udtsb    Time tendency of u field at south boundary (m/s**2)
!
!    vdteb    Time tendency of v field at east boundary (m/s**2)
!    vdtwb    Time tendency of v field at west boundary (m/s**2)
!
!    wdteb    Time tendency of w field at east boundary (m/s**2)
!    wdtwb    Time tendency of w field at west boundary (m/s**2)
!    wdtnb    Time tendency of w field at north boundary (m/s**2)
!    wdtsb    Time tendency of w field at south boundary (m/s**2)
!
!    ptdteb   Time tendency of ptprt field at east boundary (K/s)
!    ptdtwb   Time tendency of ptprt field at west boundary (K/s)
!    ptdtnb   Time tendency of ptprt field at north boundary(K/s)
!    ptdtsb   Time tendency of ptprt field at south boundary(K/s)
!
!    sdteb    Time tendency of a scalar at e-boundary
!    sdtwb    Time tendency of a scalar at w-boundary
!    sdtnb    Time tendency of a scalar at n-boundary
!    sdtsb    Time tendency of a scalar at s-boundary
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!    tem6     Temporary work array.
!    tem7     Temporary work array.
!    tem8     Temporary work array.
!    tem9     Temporary work array.
!    tem10    Temporary work array.
!    tem11    Temporary work array.
!    tem12    Temporary work array.
!    tem13    Temporary work array.
!    tem14    Temporary work array.
!    tem15    Temporary work array.
!    tem16    Temporary work array.
!    tem17    Temporary work array.
!    tem18    Temporary work array.
!    tem19    Temporary work array.
!    tem20    Temporary work array.
!    tem21    Temporary work array.
!    tem22    Temporary work array.
!    tem23    Temporary work array.
!    tem24    Temporary work array.
!    tem25    Temporary work array.
!    tem26    Temporary work array.
!
!    tem1_0   Temporary work array.
!    tem2_0   Temporary work array.
!    tem3_0   Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE             ! Force explicit declarations

  INTEGER :: mptr              ! Grid identifier.

  INCLUDE 'timelvls.inc'

  INTEGER :: nx, ny, nz        ! Number of grid points in 3 directions
  INTEGER :: nxndg,nyndg,nzndg ! Number of grid points in 3 directions

  INTEGER :: frstep            ! Indicator of the first time step call

  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
  REAL :: wcont (nx,ny,nz)     ! Contravariant vertical velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature
                               ! from that of base state atmosphere (K)
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure from that
                               ! of base state atmosphere (Pascal)
  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)
  REAL :: tke   (nx,ny,nz,nt)  ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: pbldpth(nx,ny,nt)    ! Planetary boundary layer depth (m)

  REAL :: udteb (ny,nz)        ! Time tendency of u at e-boundary (m/s**2)
  REAL :: udtwb (ny,nz)        ! Time tendency of u at w-boundary (m/s**2)
  REAL :: udtnb (nx,nz)        ! Time tendency of u at n-boundary (m/s**2)
  REAL :: udtsb (nx,nz)        ! Time tendency of u at s-boundary (m/s**2)

  REAL :: vdteb (ny,nz)        ! Time tendency of v at e-boundary (m/s**2)
  REAL :: vdtwb (ny,nz)        ! Time tendency of v at w-boundary (m/s**2)
  REAL :: vdtnb (nx,nz)        ! Time tendency of v at n-boundary (m/s**2)
  REAL :: vdtsb (nx,nz)        ! Time tendency of v at s-boundary (m/s**2)

  REAL :: wdteb (ny,nz)        ! Time tendency of w at e-boundary (m/s**2)
  REAL :: wdtwb (ny,nz)        ! Time tendency of w at w-boundary (m/s**2)
  REAL :: wdtnb (nx,nz)        ! Time tendency of w at n-boundary (m/s**2)
  REAL :: wdtsb (nx,nz)        ! Time tendency of w at s-boundary (m/s**2)

  REAL :: pdteb (ny,nz)        ! Time tendency of pprt at e-boundary (Pascal/s)
  REAL :: pdtwb (ny,nz)        ! Time tendency of pprt at w-boundary (Pascal/s)
  REAL :: pdtnb (nx,nz)        ! Time tendency of pprt at n-boundary (Pascal/s)
  REAL :: pdtsb (nx,nz)        ! Time tendency of pprt at s-boundary (Pascal/s)

  REAL :: phydro(nx,ny)        ! Big time step forcing for computing
                               ! hydrostatic pprt at k=1.

  REAL :: sdteb (ny,nz)        ! Time tendency of a scalar at e-boundary
  REAL :: sdtwb (ny,nz)        ! Time tendency of a scalar at w-boundary
  REAL :: sdtnb (nx,nz)        ! Time tendency of a scalar at n-boundary
  REAL :: sdtsb (nx,nz)        ! Time tendency of a scalar at s-boundary

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: ptbari(nx,ny,nz)     ! Inverse Base state pot. temperature (K)
  REAL :: pbari (nx,ny,nz)     ! Inverse Base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: rhostri(nx,ny,nz)    ! Inverse base state density rhobar times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific
                               ! humidity (kg/kg)
  REAL :: ppi   (nx,ny,nz)     ! Exner function.
  REAL :: csndsq(nx,ny,nz)     ! Sound wave speed squared.

  REAL :: x     (nx)           ! The x-coord. of the physical and
                               ! computational grid. Defined at u-point.
  REAL :: y     (ny)           ! The y-coord. of the physical and
                               ! computational grid. Defined at v-point.
  REAL :: z     (nz)           ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of staggered grid.

  REAL :: mapfct(nx,ny,8)      ! Map factors at scalar, u and v points

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian -d(zp)/d(x)
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian -d(zp)/d(y)
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian  d(zp)/d(z)
  REAL :: aj3x  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
  REAL :: aj3y  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
  REAL :: aj3z  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.

  REAL :: j3inv (nx,ny,nz)     ! Coordinate transformation Jacobian  d(zp)/d(z)

  REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.

  INTEGER :: ifax1(13)         ! Array containing the factors of nx for
                               ! fftopt=1.
  INTEGER :: ifax2(13)         ! Array containing the factors of ny for
                               ! fftopt=1.

  REAL :: vwork1 (nx+1,ny+1)   ! 2-D work array for fftopt=2.
  REAL :: vwork2 (ny,nx+1)     ! 2-D work array for fftopt=2.
  REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt=2.
  REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt=2.

  REAL :: sinlat(nx,ny)        ! Sin of latitude at each grid point

  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: rprntl(nx,ny,nz)     ! Reciprocal of Prandtl number

  INTEGER :: nstyps               ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps)! Soil type at each point
  REAL :: stypfrct(nx,ny,nstyps)  ! Fraction of soil types
  INTEGER :: vegtyp (nx,ny)       ! Vegetation type at each point

  REAL :: lai    (nx,ny)          ! Leaf Area Index
  REAL :: roufns (nx,ny)          ! Surface roughness
  REAL :: veg    (nx,ny)          ! Vegetation fraction

  REAL :: tsfc   (nx,ny,0:nstyps) ! Temperature at ground(K)(in top 1 cm layer)
  REAL :: tsoil  (nx,ny,0:nstyps) ! Deep soil temperature(K)(in deep 1 m layer)
  REAL :: wetsfc (nx,ny,0:nstyps) ! Surface soil moisture in the top 1 cm layer
  REAL :: wetdp  (nx,ny,0:nstyps) ! Deep soil moisture in the deep 1 m layer
  REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount
  REAL :: snowdpth(nx,ny)         ! Snow depth (m)

  REAL :: ptsfc  (nx,ny)          ! Potential temperature at the ground level (K)
  REAL :: qvsfc  (nx,ny,0:nstyps) ! Effective qv at the surface (kg/kg)


  REAL :: ptcumsrc(nx,ny,nz)   ! Source term in pt-equation due
                               ! to cumulus parameterization
  REAL :: qcumsrc(nx,ny,nz,5)  ! Source term in water equations due
                               ! to cumulus parameterization:
                               ! qcumsrc(1,1,1,1) for qv equation
                               ! qcumsrc(1,1,1,2) for qc equation
                               ! qcumsrc(1,1,1,3) for qr equation
                               ! qcumsrc(1,1,1,4) for qi equation
                               ! qcumsrc(1,1,1,5) for qs equation
  REAL :: w0avg(nx,ny,nz)      ! a closing running average vertical
                               ! velocity in 10min for K-F scheme
  REAL :: kfraincv (nx,ny)     ! K-F convective rainfall (cm)
  INTEGER :: nca (nx,ny)       ! K-F counter for CAPE release

!EMK BMJ
  REAL,INTENT(INOUT) :: cldefi(nx,ny) ! BMJ cloud efficiency
  REAL,INTENT(IN) :: xland(nx,ny)     ! BMJ land mask
                                      ! (1.0 = land, 2.0 = sea)
  REAL,INTENT(INOUT) :: bmjraincv (nx,ny) ! BMJ convective rainfall
                                          ! (cm)
!EMK END

  REAL :: raing  (nx,ny)       ! Gridscale rainfall (mm)
  REAL :: rainc  (nx,ny)       ! Cumulus rainfall (mm)
  REAL :: prcrate(nx,ny,4)     ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precipitation rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate
!
!-----------------------------------------------------------------------
!
!  Arrays for radiation
!
!-----------------------------------------------------------------------
!
  INCLUDE 'radcst.inc'      ! Radiation constants and parameters

  INTEGER :: rbufsz            ! Radiation buffer size
  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: rnflx (nx,ny)        ! Net absorbed radiation flux at surface

  REAL :: rad2d(nx,ny,nrad2d)
           ! Buffur array to carry the variables calculated or used in
! radiation calculation. The last index defines the variables:
! 1  = nrsirbm,  Solar IR surface albedo for beam
! 2  = nrsirdf,  Solar IR surface albedo for diffuse
! 3  = nrsuvbm,  Solar UV surface albedo for beam
! 4  = nrsuvdf,  Solar UV surface albedo for diffuse
! 5  = ncosz,    Cosine of zenith
! 6  = ncosss,   Cosine of angle between sun light and slope
! 7  = nfdirir,  all-sky direct downward IR flux (0.7-10 micron)
!                at the surface
! 8  = nfdifir,  all-sky diffuse downward IR flux
!                at the surface
           ! 9  = nfdirpar, all-sky direct downward and par flux
!                (0.4-0.7 micron) at the surface
! 10 = nfdifpar, all-sky diffuse downward and par flux
!                at thediffuse downward and par flux
!                at the surface

  REAL :: radbuf(rbufsz)       ! Buffer array storing temporary working
                               ! arrays for radiation computing

  INTEGER :: exbcbufsz         ! EXBC buffer size
  REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
  REAL :: bcrlx(nx,ny)         ! EXBC relaxation coefficients
!
  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))
!
  REAL :: uincr(nxndg,nyndg,nzndg)      ! Analysis increment for u
  REAL :: vincr(nxndg,nyndg,nzndg)      ! Analysis increment for v
  REAL :: wincr(nxndg,nyndg,nzndg)      ! Analysis increment for w
  REAL :: pincr(nxndg,nyndg,nzndg)      ! Analysis increment for p
  REAL :: ptincr(nxndg,nyndg,nzndg)     ! Analysis increment for pt
  REAL :: qvincr(nxndg,nyndg,nzndg)     ! Analysis increment for qv
  REAL :: qcincr(nxndg,nyndg,nzndg)     ! Analysis increment for qc
  REAL :: qrincr(nxndg,nyndg,nzndg)     ! Analysis increment for qr
  REAL :: qiincr(nxndg,nyndg,nzndg)     ! Analysis increment for qi
  REAL :: qsincr(nxndg,nyndg,nzndg)     ! Analysis increment for qs
  REAL :: qhincr(nxndg,nyndg,nzndg)     ! Analysis increment for qh
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: cdm   (nx,ny)        ! Drag coefficient for surface momentum flux
  REAL :: cdh   (nx,ny)        ! Drag coefficient for surface heat flux
  REAL :: cdq   (nx,ny)        ! Drag coefficient for surface moisture flux

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem2  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem3  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem4  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem5  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem6  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem7  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem8  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem9  (nx,ny,nz)     ! Temporary work array.
  REAL :: tem10 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem11 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem12 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem13 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem14 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem15 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem16 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem17 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem18 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem19 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem20 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem21 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem22 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem23 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem24 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem25 (nx,ny,nz)     ! Temporary work array.
  REAL :: tem26 (nx,ny,nz)     ! Temporary work array.

  REAL :: tem1_0(0:nx,0:ny,0:nz)     ! Temporary work array.
  REAL :: tem2_0(0:nx,0:ny,0:nz)     ! Temporary work array.
  REAL :: tem3_0(0:nx,0:ny,0:nz)     ! Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'     ! Global constants that control model execution
  INCLUDE 'bndry.inc'       ! Boundary condition control parameters
  INCLUDE 'phycst.inc'      ! Physical constants
  INCLUDE 'exbc.inc'        ! EXBC constants and parameters
  INCLUDE 'nudging.inc'     ! variables for nudging assimilation
  INCLUDE 'mp.inc'            ! Message passing parameters.
!
!-----------------------------------------------------------------------
!
!  Miscellaneous local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: nocalls
  SAVE nocalls
  DATA nocalls / 0 /
  INTEGER :: i,j,k
  INTEGER :: ireturn           ! Return status indicator
  REAL :: dtbig1
  REAL :: pisum
  INTEGER :: tstrtlvl
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  mgrid = mptr
!
!-----------------------------------------------------------------------
!
!  Calculate the radiation forcing and surface fluxes if desired.
!
!-----------------------------------------------------------------------
!
  CALL set_acct(rad_acct)

  IF ( radopt /= 0 .AND. MOD(nstep-1, nradstp) == 0 ) THEN

    IF (myproc == 0) &
    WRITE(6,'(a,i8,a,f10.2,a)')                                         &
        ' Compute radiation at time step,', nstep,                      &
        ', model time=',curtim,' (s)'

    CALL radiation(nx,ny,nz,rbufsz,                                     &
                   ptprt(1,1,1,1),pprt(1,1,1,1),                        &
                   qv(1,1,1,1),qc(1,1,1,1),qr(1,1,1,1),                 &
                   qi(1,1,1,1),qs(1,1,1,1),qh(1,1,1,1),                 &
                   ptbar,pbar,ppi,rhostr,                               &
                   x,y,z,zp, j3inv,                                     &
                   soiltyp(1,1,1),tsfc(1,1,0),wetsfc(1,1,0),            &
                   snowdpth,radfrc,radsw,rnflx,                         &
                   rad2d(1,1,nrsirbm), rad2d(1,1,nrsirdf),              &
                   rad2d(1,1,nrsuvbm), rad2d(1,1,nrsuvdf),              &
                   rad2d(1,1,ncosz),   rad2d(1,1,ncosss),               &
                   rad2d(1,1,nfdirir), rad2d(1,1,nfdifir),              &
                   rad2d(1,1,nfdirpar),rad2d(1,1,nfdifpar),             &
                   tem1, tem2, tem3, tem4, tem5,                        &
                   tem6, tem7, tem8, tem9, tem10,                       &
                   tem11,tem12,tem13,tem14,tem15,tem16,                 &
                   radbuf)

    IF( lvldbg >= 2 ) THEN
      CALL checkshx(radfrc, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2,              &
                    'radfrc', tem1)
      CALL checkshy(radfrc, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2,              &
                    'radfrc', tem1)
      CALL checkshx(radsw,nx,ny,1,1,nx-1,1,ny-1,1,1,                    &
                    'radsw', tem1)
      CALL checkshy(radsw,nx,ny,1,1,nx-1,1,ny-1,1,1,                    &
                    'radsw', tem1)
      CALL checkshx(rnflx,nx,ny,1,1,nx-1,1,ny-1,1,1,                    &
                    'rnflx', tem1)
      CALL checkshy(rnflx,nx,ny,1,1,nx-1,1,ny-1,1,1,                    &
                    'rnflx', tem1)
    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate the surface fluxes of momentum, heat, and water
!  vapor for use in the turbulent mixing subroutine
!
!-----------------------------------------------------------------------
!
  CALL set_acct(misc_acct)

  IF(sadvopt /= 4) THEN                   ! Leapfrop scheme
    tstrtlvl = tpast
  ELSE                                    ! Forward-based scheme
    tstrtlvl = tpresent
  END IF

  CALL sfcphysics(nx,ny,nz,nstyps,                                      &
                  u(1,1,1,tstrtlvl),v(1,1,1,tstrtlvl),                  &
                  w(1,1,1,tstrtlvl),ptprt(1,1,1,tstrtlvl),              &
                  pprt(1,1,1,tstrtlvl),qv(1,1,1,tstrtlvl),              &
                  qr(1,1,1,tstrtlvl),                                   &
                  rhostr,ptbar,pbar,qvbar,                              &
                  x,y,z, zp, j1,j2,j3,j3inv, prcrate(1,1,1),            &
                  soiltyp,stypfrct,vegtyp,lai,roufns,veg,               &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,       &
                  radsw,rnflx,                                          &
                  cdm,cdh,cdq,usflx,vsflx,ptsflx,qvsflx,pbldpth )


!call test_dump(usflx,"XXXAsfcflx_usflx",nx,ny,1,1,1)
!call test_dump(vsflx,"XXXAsfcflx_vsflx",nx,ny,1,2,1)
!call test_dump(ptsflx,"XXXAsfcflx_ptsflx",nx,ny,1,0,1)
!call test_dump(qvsflx,"XXXAsfcflx_qvsflx",nx,ny,1,0,1)
!call test_dump(pbldpth,"XXXAsfcflx_pbldpth",nx,ny,3,0,1)
!
  CALL set_acct(bc_acct)
  IF( lbcopt == 2 .AND. mgrid == 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  Read in variables from an external boundary file, and
!  calculate the linear time tendencies of the external data.
!
!-----------------------------------------------------------------------
!
    CALL extbdt(nx,ny,nz, ptbar,pbar, ireturn,                          &
                exbcbuf(nu0exb), exbcbuf(nv0exb),                       &
                exbcbuf(nw0exb), exbcbuf(npt0exb),                      &
                exbcbuf(npr0exb),exbcbuf(nqv0exb),                      &
                exbcbuf(nqc0exb),exbcbuf(nqr0exb),                      &
                exbcbuf(nqi0exb),exbcbuf(nqs0exb),                      &
                exbcbuf(nqh0exb),exbcbuf(nudtexb),                      &
                exbcbuf(nvdtexb),exbcbuf(nwdtexb),                      &
                exbcbuf(nptdtexb),exbcbuf(nprdtexb),                    &
                exbcbuf(nqvdtexb),exbcbuf(nqcdtexb),                    &
                exbcbuf(nqrdtexb),exbcbuf(nqidtexb),                    &
                exbcbuf(nqsdtexb),exbcbuf(nqhdtexb),                    &
                tem1,tem2,tem3,tem4,tem5,tem6,                          &
                tem7,tem8,tem9,tem10,tem11)

    IF( ireturn == 0 ) THEN
      GO TO 112
    ELSE IF( ireturn == 1 ) THEN
      WRITE (6,'(a/a)')                                                 &
          'Can not find the external boundary data. Dump the',          &
          'history file and restart file and then STOP the model.'
      GO TO 111
    ELSE IF( ireturn == 2 ) THEN
      WRITE (6,'(a/a)')                                                 &
          'Can not open the external boundary data. Dump the history',  &
          'file and restart file and then STOP the model.'
      GO TO 111
    ELSE IF( ireturn == 3 ) THEN
      WRITE (6,'(a/a)')                                                 &
          'Read errors in the external boundary data file. Dump the',   &
          'history file and restart file and then STOP the model.'
      GO TO 111
    ELSE
      WRITE (6,'(a/a)')                                                 &
          'Other errors in getting the external boundary data. Dump the', &
          'history file and restart file and then STOP the model.'
      GO TO 111
    END IF

    111     CONTINUE

  CALL set_acct(output_acct)
!
!-----------------------------------------------------------------------
!
!  Write out external data arrays and stop the program if data read
!  failed.
!
!-----------------------------------------------------------------------
!
    IF (mp_opt > 0) THEN
      CALL arpsstop("arpstop called from CORDINTG  ",1)
    END IF

    CALL abortdmp(mptr,nx,ny,nz,nstyps,                                 &
                  u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,         &
                  ubar,vbar,ptbar,pbar,rhostr,qvbar,kmh,kmv,            &
                  x,y,z,zp,zp(1,1,2),mapfct, j1,j2,j3,                  &
                  soiltyp,stypfrct,vegtyp,lai,roufns,veg,               &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,             &
                  raing,rainc,prcrate,                                  &
                  radfrc,radsw,rnflx,                                   &
                  usflx,vsflx,ptsflx,qvsflx,                            &
                  tem1,tem2,tem3, tem4)

!  blocking inserted for ordering i/o for message passing
    DO i=0,nprocs-1,max_fopen
      IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN
!EMK BMJ
        CALL rstout(nx,ny,nz,nstyps,exbcbufsz,                          &
                u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                 &
                udteb, udtwb, vdtnb, vdtsb,                             &
                pdteb ,pdtwb ,pdtnb ,pdtsb,                             &
                ubar,vbar,ptbar,pbar,rhostr,qvbar,                      &
                x,y,z,zp,zp(1,1,2),mapfct,                              &
                soiltyp,stypfrct,vegtyp,lai,roufns,veg,                 &
                tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,qvsfc,         &
                ptcumsrc,qcumsrc,w0avg,nca,kfraincv,                    &
                cldefi,xland,bmjraincv,                                 &
                radfrc,radsw,rnflx,                                     &
                raing,rainc,prcrate,exbcbuf,tem1)
!EMK END

      END IF
      IF (mp_opt > 0) CALL mpbarrier
    END DO

    CALL arpsstop("arpsstop called from CORDINTG run unstable, aborted",1)

    112     CONTINUE

  END IF
!
!-----------------------------------------------------------------------
!
!  Adjust potential temperature and Qv fields with grid condensation
!  or warm rain microphysics
!
!-----------------------------------------------------------------------
!
  CALL set_acct(qpfgrd_acct)

  IF ( (moist /= 0) .AND. (cnvctopt == 1) .AND.                         &
         (MOD(curtim+0.001,qpfgfrq) < (0.5*dtbig)) ) THEN

    CALL qpfgrid(nx,ny,nz,prcrate(1,1,2),                               &
         pprt(1,1,1,2),ptprt(1,1,1,2),qv(1,1,1,2),pbar,ptbar,           &
         rhostr,zp,j3,j3inv,raing, tem1(1,1,1),tem1(1,1,2),             &
         tem1(1,1,3),tem1(1,1,4))

    IF ( lbcopt == 2 .AND.  mgrid == 1 ) THEN
      CALL acct_interrupt(bc_acct)
      CALL exbcpt( nx,ny,nz, curtim, ptprt(1,1,1,2),                    &
                   exbcbuf(npt0exb),exbcbuf(nptdtexb) )
      CALL exbcq( nx,ny,nz, 1,curtim, qv(1,1,1,2),                      &
                   exbcbuf(nqv0exb),exbcbuf(nqvdtexb) )
      CALL acct_stop_inter
    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Integrate the conservation equations of momentum, heat, pressure
!  and water variables forward by one timestep
!
!-----------------------------------------------------------------------
!
  CALL set_acct(tinteg_acct)
!EMK BMJ
  CALL tinteg(mptr, frstep, nx,ny,nz,exbcbufsz,                         &
              u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth,     &
              udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,          &
              wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,          &
              phydro,sdteb,sdtwb,sdtnb,sdtsb,                           &
              ubar,vbar,ptbar,pbar,ptbari,pbari,                        &
              rhostr,rhostri,qvbar,ppi,csndsq,                          &
              x,y,z,zp,mapfct,                                          &
              j1,j2,j3,aj3x,aj3y,aj3z,j3inv,                            &
              trigs1,trigs2,ifax1,ifax2,                                &
              wsave1,wsave2,vwork1,vwork2,sinlat,                       &
              ptsfc,qvsfc(1,1,0),prcrate, radfrc,                       &
              usflx,vsflx,ptsflx,qvsflx,kmh,kmv,rprntl,                 &
              ptcumsrc,qcumsrc,raing,rainc,w0avg,nca,kfraincv,          &
              cldefi,xland,bmjraincv,                                   &
              exbcbuf, bcrlx,                                           &
              tem1,tem2,tem3,tem4,tem5,tem6,tem7,                       &
              tem8,tem9,tem10,tem11,tem12,tem13,                        &
              tem14,tem15,tem16,tem17,tem18,tem19,                      &
              tem20,tem21,tem22,tem23,tem24,tem25,tem26,                &
              tem1_0,tem2_0,tem3_0)
!EMK END

!call test_dump (pprt(1,1,1,2),"XXXAtinteg_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,3),"XXXAtinteg_pprt3",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpast),"XXXAtinteg_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAtinteg_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAtinteg_qr3",nx,ny,nz,0,1)
  CALL set_acct(misc_acct)
  IF( pdetrnd == 1 .AND. mptr == 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  Detrend pressure. This should be used only when unwanted trend
!  develops in the pressure field. The precedure forces the
!  domain averaged Exner function to zero.
!
!  Pressure detrending is sometimes necessary when open lateral
!  boundary condition is used. Pressure detrending is applied
!  to the base grid (as opposed to nested grids) only.
!
!-----------------------------------------------------------------------
!
    pisum = 0.0
    DO k=2,nz-2
      DO j=2,ny-2
        DO i=2,nx-2
          pisum = pisum + ppi(i,j,k)/(pbar(i,j,k)+pprt(i,j,k,3))        &
                                    * pprt(i,j,k,3)
        END DO
      END DO
    END DO

    IF (mp_opt > 0) THEN
      CALL mptotal(pisum)
      pisum = pisum/(nproc_x*nproc_y)
    END IF

    pisum = pisum/((nx-3)*(ny-3)*(nz-3))
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          pprt(i,j,k,3) = pprt(i,j,k,3)                                 &
                        - pisum*(pbar(i,j,k)+pprt(i,j,k,3))/ppi(i,j,k)
        END DO
      END DO
    END DO
!call test_dump (pprt(1,1,1,3),"XXXApdtrnmd_pprt3",nx,ny,nz,0,1)

  END IF
!
!-----------------------------------------------------------------------
!
!  Call the microphysics routines which adjust the thermodynamic
!  and water variables.
!
!  Upon exiting this subroutine, supersaturation is removed.
!
!-----------------------------------------------------------------------
!
  IF( moist /= 0 ) THEN

    IF( mphyopt == 1 .OR. (mphyopt == 0.AND.moist == 1) ) THEN

      IF ( cnvctopt /= 1 ) THEN

!-----------------------------------------------------------------------
!
!  Kessler warm rain microphysics or saturation adjustment only.
!
!-----------------------------------------------------------------------

        IF(frstep == 1) THEN
          dtbig1 = dtbig/2
        ELSE
          dtbig1 = dtbig
        END IF

        CALL set_acct(warmrain_acct)

!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_qr3",nx,ny,nz,0,1)
        CALL microph(nx,ny,nz,dtbig1,                                   &
                     ptprt,pprt,qv,qc,qr,qi,qs,qh,                      &
                     raing,prcrate(1,1,4),                              &
                     rhostr,pbar,ptbar,ppi,j3,j3inv,                    &
                     tem1,tem2,tem3,tem4,tem5,tem6)
!call test_dump (qr(1,1,1,tpast),"XXXAmicroph_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAmicroph_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAmicroph_qr3",nx,ny,nz,0,1)

      END IF
!
!-----------------------------------------------------------------------
!
!  The microphysics may change the boundary values of both the
!  potential temperature and specific humidity. For the choice of
!  external boundary conditions, we must adjust them to the external
!  BC.
!
!-----------------------------------------------------------------------
!
      IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN

        CALL acct_interrupt(bc_acct)

        CALL exbcpt(nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture),       &
                    exbcbuf(npt0exb), exbcbuf(nptdtexb))
        CALL exbcq( nx,ny,nz, 1,curtim+dtbig, qv(1,1,1,tfuture),        &
                    exbcbuf(nqv0exb),exbcbuf(nqvdtexb))
        CALL exbcq( nx,ny,nz, 2,curtim+dtbig, qc(1,1,1,tfuture),        &
                    exbcbuf(nqc0exb),exbcbuf(nqcdtexb) )
        CALL exbcq( nx,ny,nz, 3,curtim+dtbig, qr(1,1,1,tfuture),        &
                    exbcbuf(nqr0exb),exbcbuf(nqrdtexb) )
        CALL acct_stop_inter

      END IF

    ELSE IF( mphyopt == 2 ) THEN    ! Ice microphysics.

      IF(frstep == 1) THEN
        dtbig1 = dtbig/2
      ELSE
        dtbig1 = dtbig
      END IF
!
!-----------------------------------------------------------------------
!
!  Lin ice microphysics scheme with modifications by W.G. Tao.
!  (Lin, Y.-L., R. D. Farley, and H. D. Orville, 1983:
!  Bulk parameterization of the snow field in a cloud model.
!  J. Climate Appl. Meteor., 22, 1065-1092.
!  Tao, W.-K., and J. Simpson, 1993: Goddard cumulus ensemble model.
!  Part I: Model description. Terres. Atmos. Ocean Sci., 4, 35-72.)
!
!-----------------------------------------------------------------------
!
      CALL set_acct(linice_acct)

!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_ice_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_ice_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_ice_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXBmicroph_ice_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXBmicroph_ice_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXBmicroph_ice_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
      CALL microph_ice(nx,ny,nz,dtbig1,                                 &
                       ptprt,pprt,qv,qc,qr,qi,qs,qh,                    &
                       raing,prcrate(1,1,4),                            &
                       rhostr,pbar,ptbar,qvbar,ppi,j3,j3inv,            &
                       tem1,tem2,tem3,tem4,tem5,tem6,tem7,              &
                       tem8,tem9,tem10,tem11,tem12,tem13,tem14,tem15)
!call test_dump (qr(1,1,1,tpast),"XXXAmicroph_ice_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAmicroph_ice_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAmicroph_ice_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXAmicroph_ice_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAmicroph_ice_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAmicroph_ice_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!
!-----------------------------------------------------------------------
!
!  The microphysics may change the boundary values of both the
!  potential temperature and specific humidity. For the choice of
!  external boundary conditions, we must adjust them to the external
!  BC.
!
!-----------------------------------------------------------------------
!
      IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN

        CALL acct_interrupt(bc_acct)

        CALL exbcpt( nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture),      &
                     exbcbuf(npt0exb),exbcbuf(nptdtexb) )
        CALL exbcq( nx,ny,nz, 1,curtim+dtbig, qv(1,1,1,tfuture),        &
                     exbcbuf(nqv0exb),exbcbuf(nqvdtexb) )
        CALL exbcq( nx,ny,nz, 2,curtim+dtbig, qc(1,1,1,tfuture),        &
                     exbcbuf(nqc0exb),exbcbuf(nqcdtexb) )
        CALL exbcq( nx,ny,nz, 3,curtim+dtbig, qr(1,1,1,tfuture),        &
                     exbcbuf(nqr0exb),exbcbuf(nqrdtexb) )
        CALL exbcq( nx,ny,nz, 4,curtim+dtbig, qi(1,1,1,tfuture),        &
                     exbcbuf(nqi0exb),exbcbuf(nqidtexb) )
        CALL exbcq( nx,ny,nz, 5,curtim+dtbig, qs(1,1,1,tfuture),        &
                     exbcbuf(nqs0exb),exbcbuf(nqsdtexb) )
        CALL exbcq( nx,ny,nz, 6,curtim+dtbig, qh(1,1,1,tfuture),        &
                     exbcbuf(nqh0exb),exbcbuf(nqhdtexb) )
        CALL acct_stop_inter

      END IF
!call test_dump (qr(1,1,1,tpast),"XXXAexbcq_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAexbcq_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAexbcq_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXAexbcq_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAexbcq_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAexbcq_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)

    ELSE IF ( mphyopt == 3 ) THEN

      IF(frstep == 1) THEN
        dtbig1 = dtbig/2
      ELSE
        dtbig1 = dtbig
      END IF
!
!-----------------------------------------------------------------------
!
!  Schultz ice microphysics scheme (Schultz, P., 1995: An explicit
!  cloud physics parameterization for operational numerical
!  weather prediction. Mon. Wea. Rev., 123, 3331-3343).
!
!-----------------------------------------------------------------------
!
      CALL set_acct(nemice_acct)

!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_nem_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_nem_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_nem_qr3",nx,ny,nz,0,1)
      CALL microph_nem(nx,ny,nz,dtbig1,ptprt,pprt,qv,qc,                &
                       qr,qi,qs,qh,raing,prcrate(1,1,4),rhostr,         &
                       pbar,ptbar,ppi,j3,j3inv,                         &
                       tem1,tem2,tem3,tem4,tem5,tem6,tem7,              &
                       tem8,tem9,tem10,tem11,tem12,tem13,               &
                       tem14,tem15,tem16)
!call test_dump2(qr(1,1,1,tpast),"XXXAmicroph_nem_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAmicroph_nem_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAmicroph_nem_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!
!-----------------------------------------------------------------------
!
!  The microphysics may change the boundary values of both the
!  potential temperature and specific humidity. For the choice of
!  external boundary conditions, we must adjust them to the external
!  BC.
!
!-----------------------------------------------------------------------
!
      IF ( lbcopt == 2 .AND. mgrid == 1 ) THEN
!call test_dump (qr(1,1,1,tpast),"XXXBmicroph_nembc_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBmicroph_nembc_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBmicroph_nembc_qr3",nx,ny,nz,0,1)

        CALL acct_interrupt(bc_acct)

        CALL exbcpt( nx,ny,nz, curtim+dtbig, ptprt(1,1,1,tfuture),      &
                     exbcbuf(npt0exb),exbcbuf(nptdtexb) )
        CALL exbcq( nx,ny,nz, 1,curtim+dtbig, qv(1,1,1,tfuture),        &
                     exbcbuf(nqv0exb),exbcbuf(nqvdtexb) )
        CALL exbcq( nx,ny,nz, 2,curtim+dtbig, qc(1,1,1,tfuture),        &
                     exbcbuf(nqc0exb),exbcbuf(nqcdtexb) )
        CALL exbcq( nx,ny,nz, 3,curtim+dtbig, qr(1,1,1,tfuture),        &
                     exbcbuf(nqr0exb),exbcbuf(nqrdtexb) )
        CALL exbcq( nx,ny,nz, 4,curtim+dtbig, qi(1,1,1,tfuture),        &
                     exbcbuf(nqi0exb),exbcbuf(nqidtexb) )
        CALL exbcq( nx,ny,nz, 5,curtim+dtbig, qs(1,1,1,tfuture),        &
                     exbcbuf(nqs0exb),exbcbuf(nqsdtexb) )
        CALL exbcq( nx,ny,nz, 6,curtim+dtbig, qh(1,1,1,tfuture),        &
                     exbcbuf(nqh0exb),exbcbuf(nqhdtexb) )
        CALL acct_stop_inter

      END IF
!call test_dump (qr(1,1,1,tpast),"XXXAmicroph_nembc_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXAmicroph_nembc_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXAmicroph_nembc_qr3",nx,ny,nz,0,1)
!call test_dump2(qr(1,1,1,tpast),"XXXAmicroph_nembc_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAmicroph_nembc_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAmicroph_nembc_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)

    ELSE                             ! Invalid option
      WRITE(6,'(1x,a,/1x,a,i3)')                                        &
           'Invalid option for microphysics parameterization.',         &
           'Job stopped here. MPHYOPT was ',mphyopt
      CALL arpsstop("arpsstop called from CORDINTG improper mphyopt",1)

    END IF

    DO j=1,ny-1
      DO i=1,nx-1
        prcrate(i,j,1)=prcrate(i,j,2)+prcrate(i,j,3)+prcrate(i,j,4)
      END DO
    END DO

  END IF ! moist.ne.0
!
!-----------------------------------------------------------------------
!
!  Apply nudging
!
!-----------------------------------------------------------------------
!
  CALL set_acct(misc_acct)
  IF(nudgopt > 0 .AND. curtim < ndstop .AND. curtim >= ndstart .AND.    &
        MOD(nstep, nudgstp) == 0 ) THEN
    WRITE(6,'(a,f10.1)')                                                &
        ' Nudging forecast at ',(curtim+dtbig)

    CALL nudgeall(nx,ny,nz, nxndg,nyndg,nzndg,                          &
        u(1,1,1,tpresent),v(1,1,1,tpresent),                            &
        w(1,1,1,tpresent),pprt(1,1,1,tpresent),ptprt(1,1,1,tpresent),   &
        qv(1,1,1,tpresent),qc(1,1,1,tpresent),qr(1,1,1,tpresent),       &
        qi(1,1,1,tpresent),qs(1,1,1,tpresent),qh(1,1,1,tpresent),       &
        uincr,vincr,wincr,pincr,ptincr,qvincr,                          &
        qcincr,qrincr,qiincr,qsincr,qhincr)


    CALL nudgeall(nx,ny,nz, nxndg,nyndg,nzndg,                          &
        u(1,1,1,tfuture),v(1,1,1,tfuture),                              &
        w(1,1,1,tfuture),pprt(1,1,1,tfuture),ptprt(1,1,1,tfuture),      &
        qv(1,1,1,tfuture),qc(1,1,1,tfuture),qr(1,1,1,tfuture),          &
        qi(1,1,1,tfuture),qs(1,1,1,tfuture),qh(1,1,1,tfuture),          &
        uincr,vincr,wincr,pincr,ptincr,qvincr,                          &
        qcincr,qrincr,qiincr,qsincr,qhincr)

  END IF
!
!-----------------------------------------------------------------------
!
!  Apply the Asselin time filter to all variables
!
!-----------------------------------------------------------------------
!

  IF( flteps /= 0 .AND. frstep /= 1 ) THEN

    CALL tfilt(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke)
!call test_dump2(qr(1,1,1,tpast),"XXXAtfilt_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAtfilt_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAtfilt_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)

  END IF
!
!-----------------------------------------------------------------------
!
!
!  Call the Klemp-Lilly (1978)/Durran (1983) version of the radiation
!
!
!-----------------------------------------------------------------------
!

  IF(rbcopt == 2.OR.rbcopt == 3.OR.rbcopt == 4)THEN

    CALL set_acct(bc_acct)

    IF(frstep == 1) THEN
      dtbig1 = dtbig/2
    ELSE
      dtbig1 = dtbig
    END IF

    CALL bdtu(nx,ny,nz,dtbig1, u,ubar,udteb,udtwb)

    CALL bdtv(nx,ny,nz,dtbig1, v,vbar,vdtnb,vdtsb)

  END IF
!
!-----------------------------------------------------------------------
!
!  To prepare for the next timestep integration, swap time levels
!  for all time-dependent variables.
!
!  The swap is as follows:
!
!  Fields at time tpresent are assigned to tpast,
!  Fields at time tfuture  are assigned to tpresent,
!  Fields at time tfuture  are not changed.
!
!-----------------------------------------------------------------------
!
  CALL set_acct(misc_acct)

!call test_dump (u(1,1,1,tpast),"XXXBflip_u1",nx,ny,nz,0,1)
!call test_dump (u(1,1,1,tpresent),"XXXBflip_u2",nx,ny,nz,0,1)
!call test_dump (u(1,1,1,tfuture),"XXXBflip_u3",nx,ny,nz,0,1)
!call test_dump (v(1,1,1,tpast),"XXXBflip_v1",nx,ny,nz,0,1)
!call test_dump (v(1,1,1,tpresent),"XXXBflip_v2",nx,ny,nz,0,1)
!call test_dump (v(1,1,1,tfuture),"XXXBflip_v3",nx,ny,nz,0,1)
!call test_dump (w(1,1,1,tpast),"XXXBflip_w1",nx,ny,nz,0,1)
!call test_dump (w(1,1,1,tpresent),"XXXBflip_w2",nx,ny,nz,0,1)
!call test_dump (w(1,1,1,tfuture),"XXXBflip_w3",nx,ny,nz,0,1)
!call test_dump (ptprt(1,1,1,tpast),"XXXBflip_ptprt1",nx,ny,nz,0,1)
!call test_dump (ptprt(1,1,1,tpresent),"XXXBflip_ptprt2",nx,ny,nz,0,1)
!call test_dump (ptprt(1,1,1,tfuture),"XXXBflip_ptprt3",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,tpast),"XXXBflip_pprt1",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,tpresent),"XXXBflip_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,tfuture),"XXXBflip_pprt3",nx,ny,nz,0,1)
!call test_dump (qv(1,1,1,tpast),"XXXBflip_qv1",nx,ny,nz,0,1)
!call test_dump (qv(1,1,1,tpresent),"XXXBflip_qv2",nx,ny,nz,0,1)
!call test_dump (qv(1,1,1,tfuture),"XXXBflip_qv3",nx,ny,nz,0,1)
!call test_dump (qc(1,1,1,tpast),"XXXBflip_qc1",nx,ny,nz,0,1)
!call test_dump (qc(1,1,1,tpresent),"XXXBflip_qc2",nx,ny,nz,0,1)
!call test_dump (qc(1,1,1,tfuture),"XXXBflip_qc3",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpast),"XXXBflip_qr1",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tpresent),"XXXBflip_qr2",nx,ny,nz,0,1)
!call test_dump (qr(1,1,1,tfuture),"XXXBflip_qr3",nx,ny,nz,0,1)
!call test_dump (qi(1,1,1,tpast),"XXXBflip_qi1",nx,ny,nz,0,1)
!call test_dump (qi(1,1,1,tpresent),"XXXBflip_qi2",nx,ny,nz,0,1)
!call test_dump (qi(1,1,1,tfuture),"XXXBflip_qi3",nx,ny,nz,0,1)
!call test_dump (qs(1,1,1,tpast),"XXXBflip_qs1",nx,ny,nz,0,1)
!call test_dump (qs(1,1,1,tpresent),"XXXBflip_qs2",nx,ny,nz,0,1)
!call test_dump (qs(1,1,1,tfuture),"XXXBflip_qs3",nx,ny,nz,0,1)
!call test_dump (qh(1,1,1,tpast),"XXXBflip_qh1",nx,ny,nz,0,1)
!call test_dump (qh(1,1,1,tpresent),"XXXBflip_qh2",nx,ny,nz,0,1)
!call test_dump (qh(1,1,1,tfuture),"XXXBflip_qh3",nx,ny,nz,0,1)
!call test_dump (tke(1,1,1,tpast),"XXXBflip_tke1",nx,ny,nz,0,0)
!call test_dump (tke(1,1,1,tpresent),"XXXBflip_tke2",nx,ny,nz,0,0)
!call test_dump (tke(1,1,1,tfuture),"XXXBflip_tke3",nx,ny,nz,0,0)
  CALL tflip(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke)
!
!-----------------------------------------------------------------------
!
!  Print max./min. for the nested grid case.
!
!-----------------------------------------------------------------------

!  IF( mgrid.eq.1 .and. nestgrd.eq.1 ) THEN

  curtim = curtim + dtbig

!    CALL maxmin(mptr,nx,ny,nz,tlevel,rhobar,
!    :              u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv,
!    :              x,y,z,zp,mapfct,
!    :              tsfc,tsoil,wetsfc,wetdp,wetcanp,
!    :              tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
!  Check the stability of the integration.  If the model solution
!  appears to be exceeding the linear stability condition, then
!  perform a data dump for post-mortem.
!
!-----------------------------------------------------------------------
!
  CALL chkstab(mgrid,nx,ny,nz,nstyps,                                   &
               u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,            &
               ubar,vbar,ptbar,pbar,rhostr,qvbar,kmh,kmv,               &
               x,y,z,zp,zp(1,1,2),mapfct,j1,j2,j3,                      &
               soiltyp,stypfrct,vegtyp,lai,roufns,veg,                  &
               tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                &
               raing,rainc,prcrate,                                     &
               radfrc,radsw,rnflx,                                      &
               usflx,vsflx,ptsflx,qvsflx,                               &
               tem1,tem2,tem3, tem4)

  curtim = curtim - dtbig

!  ENDIF

  RETURN
END SUBROUTINE cordintg
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE TINTEG                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!
!EMK BMJ

SUBROUTINE tinteg(mptr, frstep, nx,ny,nz,exbcbufsz,                     & 1,66
           u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth,        &
           udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,             &
           wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,             &
           phydro,sdteb,sdtwb,sdtnb,sdtsb,                              &
           ubar,vbar,ptbar,pbar,ptbari,pbari,                           &
           rhostr,rhostri,qvbar,ppi,csndsq,                             &
           x,y,z,zp, mapfct,                                            &
           j1,j2,j3,aj3x,aj3y,aj3z,j3inv,                               &
           trigs1,trigs2,ifax1,ifax2,                                   &
           wsave1,wsave2,vwork1,vwork2,sinlat,                          &
           ptsfc,qvsfc,prcrate, radfrc,                                 &
           usflx,vsflx,ptsflx,qvsflx,kmh,kmv,rprntl,                    &
           ptcumsrc,qcumsrc,raing,rainc,w0avg,nca,kfraincv,             &
           cldefi,xland,bmjraincv,                                      &
           exbcbuf,bcrlx,                                               &
           rhofct,defsq,lenscl,                                         &
           uforce,vforce,wforce,pforce,ptforce,                         &
           tem1,tem2,tem3,tem4,tem5,tem6,                               &
           tem7,tem8,tem9,tem10,tem11,tem12,tem13,tem14,tem15,          &
           tem16,tem17,tem18,                                           &
           tem1_0,tem2_0,tem3_0)
!EMK END
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Orchestrate the time integration of the dynamics of the basic governing
!  equations for a single time step. Note that the time tendencies of the
!  non-conservative processes (e.g., microphysics, radiation, etc.)
!  are not included in this routine but are applied in their
!  corresponding routines.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/21/92.
!
!  MODIFICATION HISTORY:
!
!  5/03/92 (K. Droegemeier)
!    Added full documentation.
!
!  5/20/92 (M. Xue)
!    Reformatted certain documentations.
!
!  4/10/93 (M. Xue & Hao Jin)
!    Add the terrain.
!
!  5/26/93 (M. Xue)
!    Added w into the solvpt argument list.
!
!  9/7/94 (M.Xue)
!    Call to ADVP modified.
!
!  9/10/94 (D. Weber & Y. Lu)
!  Cleaned up documentation.
!
!  8/20/95 (M. Xue)
!  Bug fix with the use of ptcumsrc.
!
!  10/31/95 (D. Weber)
!  Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
!  radiation condition.
!
!  01/23/1996 (Donghai Wang, Ming Xue and Yuhe Liu)
!  Added the map projection factor to ARPS governing equations.
!
!  03/11/96 (M. Xue)
!  Modified the call to SOLVTKE, so that all three time levels of
!  ptprt, pprt etc are passed into the routine.
!
!  03/19/96 (Yuhe Liu)
!  Added radiation forcing
!
!  07/10/1997 (Fanyou Kong - CMRP)
!  Fixed a bug in 'solvtke' with corrected temporary arga
!
!  07/22/97 (D. Weber)
!  Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version
!  of the upper w-p radiation condition (fftopt=2).
!
!  08/01/97 (Zonghui Huo)
!  Added Kain-fritsch cumulus parameterization scheme.
!
!  10/21/97 (Donghai Wang)
!  Using total density (rho) in the calculation of the pressure
!  gradient force terms, and added the second order terms
!  in the linerized buoyancy terms.
!
!  11/05/97 (D. Weber)
!  Added phydro array for use in the bottom boundary condition for
!  perturbation pressure (hydrostatic).
!
!  11/06/97 (D. Weber)
!  Added three additional levels to the mapfct array.  The three
!  levels (4,5,6) represent the inverse of the first three in order.
!  The inverse map factors are computed to improve efficiency.
!
!  4/15/1998 (Donghai Wang)
!  Added the source terms to the right hand terms of the qc,qr,qi,qs
!  equations due to the K-F cumulus parameterization.
!
!  4/15/1998 (Donghai Wang)
!  Added the running average vertical velocity (array w0avg)
!  for the K-F cumulus parameterization scheme.
!
!  9/18/98 (D. Weber)
!  Added precomputed averages of j3 in the x, y, and z directions
!  to improve the code efficiency.
!
!  07/23/2001 (K. Brewster)
!  Added mptr to argument list, needed for printing of diagnostic 
!  noise parameter.  Also added mptr to calls to smlstep, solvuw, 
!  solvwpex, and solvwpim.
!
!  13 March 2002 (Eric Kemp)
!  Added arrays for WRF BMJ cumulus scheme.
!
!  April 2002 (Fanyou Kong)
!  Added cnvctopt=5 option for the new WRF K-F (KF_ETA) scheme
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    frstep   Initial time step flag.
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at times tpast and tpresent (m/s)
!    v        y component of velocity at times tpast and tpresent (m/s)
!    w        Vertical component of Cartesian velocity at times
!             tpast and tpresent (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!             computational coordinates (m/s)
!    ptprt    Perturbation potential temperature at times tpast and
!             tpresent (K)
!    pprt     Perturbation pressure at times tpast and tpresent (Pascal)
!    qv       Water vapor specific humidity at times tpast and tpresent (kg/kg)
!    qc       Cloud water mixing ratio at times tpast and tpresent (kg/kg)
!    qr       Rainwater mixing ratio at times tpast and tpresent (kg/kg)
!    qi       Cloud ice mixing ratio at times tpast and tpresent (kg/kg)
!    qs       Snow mixing ratio at times tpast and tpresent (kg/kg)
!    qh       Hail mixing ratio at times tpast and tpresent (kg/kg)
!
!    udteb    Time tendency of u field at east boundary (m/s**2)
!    udtwb    Time tendency of u field at west boundary (m/s**2)
!
!    vdtnb    Time tendency of v field at north boundary (m/s**2)
!    vdtsb    Time tendency of v field at south boundary (m/s**2)
!
!    pdteb    Time tendency of pprt field at east boundary (Pascal/s)
!    pdtwb    Time tendency of pprt field at west boundary (Pascal/s)
!    pdtnb    Time tendency of pprt field at north boundary (Pascal/s)
!    pdtsb    Time tendency of pprt field at south boundary (Pascal/s)
!
!    phydro   Big time step forcing term for use in computing the
!             hydrostatic pressure at k=1.
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    ptbari   Inverse Base state potential temperature (K)
!    pbari    Inverse Base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    rhostri  Inverse base state density rhobar times j3 (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!    cnsd     Sound wave speed.
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!    mapfct   Map factors at scalar, u and v points
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    aj3x     Avgx of the coordinate transformation Jacobian  d(zp)/dz
!    aj3y     Avgy of the coordinate transformation Jacobian  d(zp)/dz
!    aj3z     Avgz of the coordinate transformation Jacobian  d(zp)/dz
!    j3inv    Inverse of the coordinate transformation j3
!    trigs1   Array containing pre-computed trig function for
!             fftopt=1.
!    trigs2   Array containing pre-computed trig function for
!             fftopt=1.
!    ifax1    Array containing the factors of nx for fftopt=1.
!    ifax2    Array containing the factors of ny for fftopt=1.
!
!    vwork1   2-D work array for fftopt=2.
!    vwork2   2-D work array for fftopt=2.
!    wsave1   Work array for fftopt=2.
!    wsave2   Work array for fftopt=2.
!
!    sinlat   Sin of latitude at each grid point
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux (kg/(m**2 * s))
!
!    ptcumsrc Source term in pt-equation due to cumulus parameterization
!    qcumsrc Source term in water equations due to cumulus parameterization
!    kfraincv   K-F convective rainfall (cm)
!    nca      K-F counter for CAPE release
!    cldefi   BMJ cloud efficiency
!    xland    BMJ land/sea mask
!    bmjraincv   BMJ convective rainfall (cm)
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!    rprntl   Reciprocal of Prandtl number
!
!    radfrc   Radiation forcing (K/s)
!
!  OUTPUT:
!
!    u        x component of velocity at time tfuture (m/s)
!    v        y component of velocity at time tfuture (m/s)
!    w        Vertical component of Cartesian velocity at tfuture (m/s)
!    ptprt    Perturbation potential temperature at time tfuture (K)
!    pprt     Perturbation pressure at time tfuture (Pascal)
!    qv       Water vapor specific humidity at time tfuture (kg/kg)
!    qc       Cloud water mixing ratio at time tfuture (kg/kg)
!    qr       Rainwater mixing ratio at time tfuture (kg/kg)
!    qi       Cloud ice mixing ratio at time tfuture (kg/kg)
!    qs       Snow mixing ratio at time tfuture (kg/kg)
!    qh       Hail mixing ratio at time tfuture (kg/kg)
!
!    udteb    Time tendency of u field at east boundary (m/s**2)
!    udtwb    Time tendency of u field at west boundary (m/s**2)
!
!    vdtnb    Time tendency of v field at north boundary (m/s**2)
!    vdtsb    Time tendency of v field at south boundary (m/s**2)
!
!    pdteb    Time tendency of pprt field at east boundary (Pascal/s)
!    pdtwb    Time tendency of pprt field at west boundary (Pascal/s)
!    pdtnb    Time tendency of pprt field at north boundary (Pascal/s)
!    pdtsb    Time tendency of pprt field at south boundary (Pascal/s)
!
!    phydro   Big time step forcing term for use in computing the
!             hydrostatic pressure at k=1.
!
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!    pbldpth  Planetary boundary layer depth (m)
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!    rprntl   Reciprocal of Prandtl number
!
!  WORK ARRAYS:
!
!    defsq    Deformation squared (1/s**2)
!    lenscl   Turbulent mixing length scale (m)
!    uforce   Acoustically inactive forcing terms in u-Eq. ((kg/(m*s)**2)
!    vforce   Acoustically inactive forcing terms in v-Eq. ((kg/(m*s)**2)
!    wforce   Acoustically inactive forcing terms in w-Eq. ((kg/(m*s)**2)
!    pforce   Acoustically inactive forcing terms in p-eq. (Pascal/s)
!    ptforce  Gravity wave inactive forcing terms in pt-eq. (K*kg/(m**3*s))
!
!    rhofct   rho-factor: rhobar/rho
!
!    udtnb    Time tendency of u field at north boundary (m/s**2)
!    udtsb    Time tendency of u field at south boundary (m/s**2)
!
!    vdteb    Time tendency of v field at east boundary (m/s**2)
!    vdtwb    Time tendency of v field at west boundary (m/s**2)
!
!    wdteb    Time tendency of w field at east boundary (m/s**2)
!    wdtwb    Time tendency of w field at west boundary (m/s**2)
!    wdtnb    Time tendency of w field at north boundary (m/s**2)
!    wdtsb    Time tendency of w field at south boundary (m/s**2)
!
!    ptdteb   Time tendency of ptprt field at east boundary (K/s)
!    ptdtwb   Time tendency of ptprt field at west boundary (K/s)
!    ptdtnb   Time tendency of ptprt field at north boundary(K/s)
!    ptdtsb   Time tendency of ptprt field at south boundary(K/s)
!
!    sdteb    Time tendency of a scalar at e-boundary
!    sdtwb    Time tendency of a scalar at w-boundary
!    sdtnb    Time tendency of a scalar at n-boundary
!    sdtsb    Time tendency of a scalar at s-boundary
!
!    tem1     Temporary work array
!    tem2     Temporary work array
!    tem3     Temporary work array
!    tem4     Temporary work array
!    tem5     Temporary work array
!    tem6     Temporary work array
!    tem7     Temporary work array
!    tem8     Temporary work array
!    tem9     Temporary work array
!    tem10    Temporary work array
!    tem11    Temporary work array
!    tem12    Temporary work array
!    tem13    Temporary work array
!    tem14    Temporary work array
!    tem15    Temporary work array
!    tem16    Temporary work array
!    tem17    Temporary work array
!    tem18    Temporary work array
!
!    tem1_0   Temporary work array.
!    tem2_0   Temporary work array.
!    tem3_0   Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE        ! Force explicit declarations

  INCLUDE 'timelvls.inc'

  INTEGER :: mptr

  INTEGER :: nx, ny, nz        ! Number of grid points in 3 directions

  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
  REAL :: wcont (nx,ny,nz)     ! Contravariant vertical velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
  REAL :: tke   (nx,ny,nz,nt)  ! Turbulent Kinetic Energy ((m/s)**2)
  REAL :: pbldpth(nx,ny,nt)    ! Planetary boundary layer depth (m)

  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)

  REAL :: udteb (ny,nz)        ! Time tendency of u at e-boundary (m/s**2)
  REAL :: udtwb (ny,nz)        ! Time tendency of u at w-boundary (m/s**2)
  REAL :: udtnb (nx,nz)        ! Time tendency of u at n-boundary (m/s**2)
  REAL :: udtsb (nx,nz)        ! Time tendency of u at s-boundary (m/s**2)

  REAL :: vdteb (ny,nz)        ! Time tendency of v at e-boundary (m/s**2)
  REAL :: vdtwb (ny,nz)        ! Time tendency of v at w-boundary (m/s**2)
  REAL :: vdtnb (nx,nz)        ! Time tendency of v at n-boundary (m/s**2)
  REAL :: vdtsb (nx,nz)        ! Time tendency of v at s-boundary (m/s**2)

  REAL :: wdteb (ny,nz)        ! Time tendency of w at e-boundary (m/s**2)
  REAL :: wdtwb (ny,nz)        ! Time tendency of w at w-boundary (m/s**2)
  REAL :: wdtnb (nx,nz)        ! Time tendency of w at n-boundary (m/s**2)
  REAL :: wdtsb (nx,nz)        ! Time tendency of w at s-boundary (m/s**2)

  REAL :: pdteb (ny,nz)        ! Time tendency of pprt at e-boundary (Pascal/s)
  REAL :: pdtwb (ny,nz)        ! Time tendency of pprt at w-boundary (Pascal/s)
  REAL :: pdtnb (nx,nz)        ! Time tendency of pprt at n-boundary (Pascal/s)
  REAL :: pdtsb (nx,nz)        ! Time tendency of pprt at s-boundary (Pascal/s)

  REAL :: phydro(nx,ny)        ! Pressure at k=1 computed using pert.
                               ! hydrostatic relation.

  REAL :: sdteb (ny,nz)        ! Time tendency of a variable at e-boundary
  REAL :: sdtwb (ny,nz)        ! Time tendency of a variable at w-boundary
  REAL :: sdtnb (nx,nz)        ! Time tendency of a variable at n-boundary
  REAL :: sdtsb (nx,nz)        ! Time tendency of a variable at s-boundary

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: ptbari(nx,ny,nz)     ! Inverse Base state pot. temperature (K)
  REAL :: pbari (nx,ny,nz)     ! Inverse Base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: rhostri(nx,ny,nz)    ! Inverse base state density rhobar times j3.
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)
  REAL :: ppi   (nx,ny,nz)     ! Exner function
  REAL :: csndsq(nx,ny,nz)     ! Sound wave speed squared.

  REAL :: x     (nx)           ! The x-coord. of the physical and
                               ! computational grid. Defined at u-point.
  REAL :: y     (ny)           ! The y-coord. of the physical and
                               ! computational grid. Defined at v-point.
  REAL :: z     (nz)           ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of the staggered grid.

  REAL :: mapfct(nx,ny,8)      ! Map factors at scalar, u and v points

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! - d( zp )/d( x ).
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! - d( zp )/d( y ).
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).
  REAL :: aj3x  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
  REAL :: aj3y  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
  REAL :: aj3z  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
  REAL :: j3inv (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).
  REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  INTEGER :: ifax1(13)         ! Array containing the factors of nx for
                               ! fftopt=1.
  INTEGER :: ifax2(13)         ! Array containing the factors of ny for
                               ! fftopt=1.

  REAL :: vwork1 (nx+1,ny+1)   ! 2-D work array for fftopt=2.
  REAL :: vwork2 (ny,nx+1)     ! 2-D work array for fftopt=2.
  REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt =2.
  REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt =2.

  REAL :: sinlat(nx,ny)        ! Sin of latitude at each grid point

  REAL :: ptsfc  (nx,ny)       ! Potential temperature at the ground level (K)
  REAL :: qvsfc  (nx,ny)       ! Effective qv at the surface (kg/kg)

  REAL :: prcrate(nx,ny,4)     ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precipitation rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: rprntl(nx,ny,nz)     ! Reciprocal of Prandtl number

  REAL :: defsq (nx,ny,nz)     ! Deformation squared (1/s**2)
  REAL :: lenscl(nx,ny,nz)     ! Turbulent mixing length scale (m)

  REAL :: uforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in u-momentum equation (kg/(m*s)**2)
  REAL :: vforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in v-momentum equation (kg/(m*s)**2)
  REAL :: wforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in w-momentum equation (kg/(m*s)**2)
  REAL :: pforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in pressure equation (Pascal/s)
  REAL :: ptforce(nx,ny,nz)    ! Gravity wave inactive forcing terms
                               ! in pressure equation (Pascal/s)

  REAL :: ptcumsrc(nx,ny,nz)   ! Source term in pt-equation due
                               ! to cumulus parameterization
  REAL :: qcumsrc(nx,ny,nz,5)  ! Source term in water equations due
                               ! to cumulus parameterization:
                               ! qcumsrc(1,1,1,1) for qv equation
                               ! qcumsrc(1,1,1,2) for qc equation
                               ! qcumsrc(1,1,1,3) for qr equation
                               ! qcumsrc(1,1,1,4) for qi equation
                               ! qcumsrc(1,1,1,5) for qs equation
  REAL :: w0avg(nx,ny,nz)      ! a closing running average vertical
                               ! velocity in 10min for K-F scheme
  REAL :: kfraincv(nx,ny)      ! K-F convective rainfall (cm)
  INTEGER :: nca(nx,ny)        ! K-F counter for CAPE release

!EMK BMJ
  REAL,INTENT(INOUT) :: cldefi(nx,ny) ! BMJ cloud efficiency
  REAL,INTENT(IN) :: xland(nx,ny)     ! BMJ land mask
                                      ! (1.0 = land, 2.0 = sea)
  REAL,INTENT(INOUT) :: bmjraincv(nx,ny) ! BMJ convective rainfall
                                         ! (cm)
!EMK END

  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain

  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)

  INTEGER :: exbcbufsz         ! EXBC buffer size
  REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array
  REAL :: bcrlx (nx,ny)        ! EXBC relaxation coefficients

  REAL :: rhofct(nx,ny,nz)     ! rho-factor: rhobar/rho

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array
  REAL :: tem2  (nx,ny,nz)     ! Temporary work array
  REAL :: tem3  (nx,ny,nz)     ! Temporary work array
  REAL :: tem4  (nx,ny,nz)     ! Temporary work array
  REAL :: tem5  (nx,ny,nz)     ! Temporary work array
  REAL :: tem6  (nx,ny,nz)     ! Temporary work array
  REAL :: tem7  (nx,ny,nz)     ! Temporary work array
  REAL :: tem8  (nx,ny,nz)     ! Temporary work array
  REAL :: tem9  (nx,ny,nz)     ! Temporary work array
  REAL :: tem10 (nx,ny,nz)     ! Temporary work array
  REAL :: tem11 (nx,ny,nz)     ! Temporary work array
  REAL :: tem12 (nx,ny,nz)     ! Temporary work array
  REAL :: tem13 (nx,ny,nz)     ! Temporary work array
  REAL :: tem14 (nx,ny,nz)     ! Temporary work array
  REAL :: tem15 (nx,ny,nz)     ! Temporary work array
  REAL :: tem16 (nx,ny,nz)     ! Temporary work array
  REAL :: tem17 (nx,ny,nz)     ! Temporary work array
  REAL :: tem18 (nx,ny,nz)     ! Temporary work array

  REAL :: tem1_0(0:nx,0:ny,0:nz)     ! Temporary work array.
  REAL :: tem2_0(0:nx,0:ny,0:nz)     ! Temporary work array.
  REAL :: tem3_0(0:nx,0:ny,0:nz)     ! Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  REAL :: dtbig1               ! Local value of big time step size
  REAL :: dtsml1               ! Local value of small time step size
  INTEGER :: frstep            ! Flag for the initial time step
  INTEGER :: i,j,k
  INTEGER :: qflag             ! Indicator for the water/ice type
                               ! when calling SOLVQ.
  INTEGER :: ntst
  REAL :: tst
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'bndry.inc'
  INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  To initialize the three-time level scheme, we assume for
!  the first step that the values of all variables at time
!  past equal the values at time present.  We then perform a
!  forward-in-time integration for the first step only and
!  then, using the same algorithm, we perform a centered-in-time
!  integration for all subsequent steps. The size of the
!  timestep is adjusted accordingly, such that the leapfrog step
!  size is twice that of the first forward-in-time step.
!
!-----------------------------------------------------------------------
!
  IF(frstep == 1) THEN      ! frstep=1 on the first time step
                            ! only - indicating a forward-in-
                            ! time integration.
    dtbig1 = dtbig/2
    dtsml1 = dtsml/2

  ELSE

    dtbig1 = dtbig
    dtsml1 = dtsml

  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate wcont at time tpresent. wcont will be used in FRCUVW
!  and FRCP. rhostr averaged to u, v, w points. They are stored
!  in tem1, tem2 and tem3.
!
!-----------------------------------------------------------------------
!
  CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)

!call test_dump (w(1,1,1,tpresent),"XXX1Bwc_w1",nx,ny,nz,1,0)
  CALL wcontra(nx,ny,nz,                                                &
               u(1,1,1,tpresent),v(1,1,1,tpresent),                     &
               w(1,1,1,tpresent),mapfct,j1,j2,j3,aj3z,                  &
               rhostr,tem1,tem2,tem3,wcont,tem4,tem5)
!
!-----------------------------------------------------------------------
!
!  Compute the acoustically inactive terms in the momentum and
!  pressure equations that are held fixed during the small time
!  step computations.  This includes advection, buoyancy, mixing
!  (both physical and computational), and the Coriolis terms.
!  These forcing terms are accumulated into arrays for each
!  of the momentum equations, e.g., uforce for the u-equation,
!  vforce for the v-equation, wforce for the w-equation and
!  pforce for the p-equation.
!
!  Note: Arrays, pforce and ptforce, are used as work arrays in
!        subroutine frcuvw.
!
!-----------------------------------------------------------------------
!
!call test_dump2(rhofct,"XXXBfrcuvw_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
  CALL frcuvw(nx,ny,nz,exbcbufsz,dtbig1,                                &
              u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,pbldpth,     &
              ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,qvbar,           &
              usflx,vsflx, x,y,z,zp, mapfct,                            &
              j1,j2,j3,aj3x,aj3y,aj3z,j3inv, sinlat, ptsfc,             &
              uforce,vforce,wforce,kmh,kmv,rprntl,lenscl,defsq,         &
              exbcbuf, bcrlx,rhofct,phydro,                             &
              pforce,ptforce,                                           &
              tem1,tem2,tem3,tem4,tem5,tem6,tem7,                       &
              tem8,tem9,tem10,tem11,tem12,tem13,tem14)
!call test_dump (uforce,"XXX2uforce",nx,ny,nz,1,0)
!call test_dump2(rhofct,"XXXAfrcuvw_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)

  IF (tmixopt == 4) THEN
!
!-----------------------------------------------------------------------
!
!  Integrate the TKE equation.
!  NOTE: Arrays lenscl and defsq should NOT be changed between
!  calls to frcuvw and solvtke.
!  Note: ptforce is used as a work array in the following call.
!
!-----------------------------------------------------------------------
!
    CALL solvtke(nx,ny,nz,dtbig1,u,v,wcont,ptprt,pprt,                  &
                 qv,qc,qr,qi,qs,qh,tke,                                 &
                 ubar,vbar,ptbar,pbar,ptbari,rhostr,rhostri,qvbar,      &
                 x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,j3inv,            &
                 kmh,kmv,rprntl,lenscl,defsq,                           &
                 ptsflx,qvsflx,                                         &
                 tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,               &
                 tem9,tem10,tem11,ptforce, & ! ptforce used as a tem array
             tem1_0,tem2_0,tem3_0)

  END IF

  CALL frcp(nx,ny,nz,exbcbufsz, dtbig1,                                 &
            u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,                   &
            ptbar,pbar,rhostr,qvbar,mapfct,j1,j2,j3,aj3x,aj3y,aj3z,     &
            pforce,                                                     &
            exbcbuf, bcrlx,                                             &
            tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
!call test_dump (pprt(1,1,1,2),"XXXAfrcp_pprt2",nx,ny,nz,0,0)
!call test_dump (pprt(1,1,1,3),"XXXAfrcp_pprt3",nx,ny,nz,0,0)
!call test_dump (pforce,"XXXAfrcp_pforce",nx,ny,nz,0,0)
!
!-----------------------------------------------------------------------
!
!  Calculate gravity wave or acoustic wave inactive terms in the
!  potential temperature equation.
!
!  The force terms are stored in ptforce.
!
!-----------------------------------------------------------------------
!
  CALL frcpt(nx,ny,nz, exbcbufsz, dtbig1,ptprt,u,v,w,wcont,             &
             ptbar,rhostr,rhostri,kmh,kmv,rprntl,                       &
             usflx,vsflx,ptsflx,pbldpth,                                &
             x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,j3inv,ptsfc,           &
             ptforce,                                                   &
             exbcbuf, bcrlx,                                            &
             tem1,tem2,tem3,tem4,tem5,tem6,tem7,                        &
             tem8,tem9,tem10,tem11,                                     &
             tem1_0,tem2_0,tem3_0,tem12)

  CALL set_acct(tinteg_acct)
  IF ( radopt == 2 ) THEN
    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          ptforce(i,j,k) = ptforce(i,j,k)                               &
                         + rhostr(i,j,k)*radfrc(i,j,k)
        END DO
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate source and sink terms in temperature (ptprt) and
!  moisture (qv) equations that are due to subgrid scale cumulus
!  convection.
!
!-----------------------------------------------------------------------
!
!EMK BMJ
!  IF (cnvctopt == 1 .OR. cnvctopt == 2 .OR. cnvctopt == 3) THEN
  IF (cnvctopt == 1 .OR. cnvctopt == 2 .OR. cnvctopt == 3 .OR. &
      cnvctopt == 4 .OR. cnvctopt == 5) THEN
!END EMK

    IF (cnvctopt == 3 .OR. cnvctopt == 5) THEN
      CALL set_acct(kfcum_acct)
!EMK BMJ
    ELSE IF (cnvctopt == 4 ) THEN
      CALL set_acct(bmjcum_acct)
!EMK END
    ELSE 
      CALL set_acct(kuocum_acct)
    ENDIF

!-----------------------------------------------------------------------
!
! Calculate w0avg, is close to a running mean vertical velocity,
! tst is the number of time steps in 10min for K-F scheme
!
!-----------------------------------------------------------------------
!
    IF (cnvctopt == 3 .OR. cnvctopt == 5) THEN

      ntst=nint( 600.0/dtbig )
      tst=FLOAT(ntst)

      IF( (curtim-tstart) <= 300.0 .AND. initopt /= 2 ) THEN
        DO k=1,nz
          DO j=1,ny-1
            DO i=1,nx-1
              w0avg(i,j,k)= (w(i,j,k,1)+w(i,j,k,2))*0.5
            END DO
          END DO
        END DO

      ELSE IF ( (curtim-tstart) > 300.0 .OR. initopt == 2) THEN

        DO k=1,nz
          DO j=1,ny-1
            DO i=1,nx-1
              w0avg(i,j,k)=(w0avg(i,j,k)*(tst-1.)+w(i,j,k,2))/tst
            END DO
          END DO
        END DO
      END IF
    END IF
!
!-----------------------------------------------------------------------
!
!  Calculate vertical velocity for KUO scheme, or running-average
!  vertical velocity (m/s) for Kain Fritsch scheme
!
!-----------------------------------------------------------------------
!
!    IF( cnvctopt == 3 .AND. MOD(curtim+0.001,confrq) <= (0.5*dtbig)     &
!             .AND. ((curtim-tstart) > 300.0                             &
    IF( (cnvctopt == 3 .OR. cnvctopt == 5) .AND. MOD(curtim+0.001,confrq) &
             <= (0.5*dtbig) .AND. ((curtim-tstart) > 300.0                             &
                    .OR. initopt == 2)) THEN

      DO k=1,nz
        DO j=1,ny-1
          DO i=1,nx-1
!        tem1(i,j,k) = (w(i,j,k,1)+w(i,j,k,2))*0.5
            tem1(i,j,k) = w0avg(i,j,k)
          END DO
        END DO
      END DO


    ELSE IF( MOD(curtim+0.001,confrq) <= (0.5*dtbig).OR.                &
              nstep == 1 ) THEN
      DO k=1,nz
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k) = wcont(i,j,k)*j3(i,j,k)
          END DO
        END DO
      END DO
    END IF
!
!-----------------------------------------------------------------------
!
!  Call cumulus parameterization schemes
!  Make sure to reset kfraincv, ptcumsrc and qcumsrc to 0 once nca<=0
!
!-----------------------------------------------------------------------
!
    IF( (moist /= 0) .AND. ((curtim-tstart) > 300.0 .OR. initopt == 2) .AND. &
          (MOD(curtim+0.001,confrq) <= (0.5*dtbig))) THEN

      DO j=1,ny-1
        DO i=1,nx-1
          IF (nca(i,j) <= 0) THEN
            kfraincv(i,j) = 0.0
            prcrate(i,j,3) = 0.0
          END IF
        END DO
      END DO

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            IF (nca(i,j) <= 0) THEN
              ptcumsrc(i,j,k) = 0.0
              qcumsrc(i,j,k,1) = 0.0
              qcumsrc(i,j,k,2) = 0.0
              qcumsrc(i,j,k,3) = 0.0
              qcumsrc(i,j,k,4) = 0.0
              qcumsrc(i,j,k,5) = 0.0
            END IF
          END DO
        END DO
      END DO

!EMK BMJ
      CALL qpfcums(nx,ny,nz,prcrate(1,1,3),qvsflx,                      &
                   u(1,1,1,2),v(1,1,1,2),tem1,                          &
                   pprt(1,1,1,2),ptprt(1,1,1,2),qv(1,1,1,2),            &
                   pbar,ptbar,qvbar,rhostr,zp,j3,                       &
                   ptcumsrc,qcumsrc,rainc,nca,kfraincv,                 &
                   cldefi,xland,bmjraincv,                              &
                   tem2,tem3,tem4,tem5,tem6,                            &
                   tem7,tem8,tem9,tem10,tem11)
!EMK END
    END IF
!
!  Accumulate rainc and update nca. Note: raincv is in cm.
!
    DO j=1,ny-1
      DO i=1,nx-1
        rainc(i,j) = rainc(i,j) + 10.0*kfraincv(i,j) + 10.0*bmjraincv(i,j)
        IF ( nca(i,j) >= 1 ) nca(i,j) = nca(i,j) - 1
      END DO
    END DO

    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          ptforce(i,j,k) = ptforce(i,j,k)                               &
                         + rhostr(i,j,k)*ptcumsrc(i,j,k)
        END DO
      END DO
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  End of cumulus parameterization
!
!-----------------------------------------------------------------------
!
  CALL set_acct(misc_acct)

  IF( lvldbg >= 2 ) THEN
    CALL checkuhx(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2,                &
                  'ufrcx', tem1)
    CALL checkuhy(uforce, nx,ny,nz,2,nx-1,1,ny-1,2,nz-2,                &
                  'ufrcy', tem1)
    CALL checkvhx(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2,                &
                  'vfrcx', tem1)
    CALL checkvhy(vforce, nx,ny,nz,1,nx-1,2,ny-1,2,nz-2,                &
                  'vfrcy', tem1)
    CALL checkwhx(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1,                &
                  'wfrcx', tem1)
    CALL checkwhy(wforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1,                &
                  'wfrcy', tem1)
    CALL checkshx(pforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2,                &
                  'pfrcx', tem1)
    CALL checkshy(pforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2,                &
                  'pfrcy', tem1)
    CALL checkshx(ptforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2,               &
                  'ptfrcx', tem1)
    CALL checkshy(ptforce, nx,ny,nz,1,nx-1,1,ny-1,2,nz-2,               &
                  'ptfrcy', tem1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Integrate the momentum and pressure equations (i.e., the
!  acoustically-active equations) in time using a mode-splitting
!  approach.  The momentum components are advanced in time using
!  a forward scheme (relative to the pressure gradient force terms)
!  and the pressure is advanced using a backward scheme (relative to
!  the divergence term, which is evaluated using the newly updated
!  velocities hence the backward scheme). These small time steps
!  bring the momentum and pressure from time tpast to tfuture.
!  During these steps, the acoustically-inactive forcing terms (e.g.,
!  uforce, pforce, etc.) are held fixed at time tpresent (i.e. they
!  are evaluated only once, at time tpresent, for each large time
!  step integration), hence the time integration is leapfrog-
!  in-time relative to these forcing terms.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  When this solver is called for a nested grid (not the base grid),
!  the information stored in u, v, w and pprt is used to set the
!  time tendency of these variables at the lateral boundaries.
!
!-----------------------------------------------------------------------
!
  IF( mgrid /= 1 .AND. nestgrd == 1 ) THEN

    CALL nestbdt(nx,ny,nz,u, 1,nx,1,ny-1,1,nz-1,dtbig,                  &
                 udteb,udtwb,udtnb,udtsb)

    CALL nestbdt(nx,ny,nz,v, 1,nx-1,1,ny,1,nz-1,dtbig,                  &
                 vdteb,vdtwb,vdtnb,vdtsb)

    CALL nestbdt(nx,ny,nz,w, 1,nx-1,1,ny-1,1,nz,dtbig,                  &
                 wdteb,wdtwb,wdtnb,wdtsb)

    CALL nestbdt(nx,ny,nz,pprt, 1,nx-1,1,ny-1,1,nz-1,dtbig,             &
                 pdteb,pdtwb,pdtnb,pdtsb)

    CALL nestbdt(nx,ny,nz,ptprt,1,nx-1,1,ny-1,1,nz-1,dtbig,             &
                 sdteb,sdtwb,sdtnb,sdtsb)

!
!-----------------------------------------------------------------------
!
!    For the first forward time step, make sure that the fields at
!    tpast equal those at tpresent. This can only be done after
!    the time tendencies on the boundaries are calculated and stored
!    in the time tendency arrays.
!
!-----------------------------------------------------------------------
!
    IF( frstep == 1 ) THEN

      DO k=1,nz
        DO j=1,ny
          DO i=1,nx
            u    (i,j,k,tpast)=u    (i,j,k,tpresent)
            v    (i,j,k,tpast)=v    (i,j,k,tpresent)
            w    (i,j,k,tpast)=w    (i,j,k,tpresent)
            pprt (i,j,k,tpast)=pprt (i,j,k,tpresent)
            ptprt(i,j,k,tpast)=ptprt(i,j,k,tpresent)
          END DO
        END DO
      END DO

    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Note here defsq is used as a work array by ACOUST.
!
!-----------------------------------------------------------------------
!
  CALL set_acct(smlstp_acct)

  CALL smlstep(mptr, nx,ny,nz, exbcbufsz, dtbig1,dtsml1,                &
               u,v,w,wcont,pprt,ptprt,                                  &
               udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,         &
               wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,         &
               phydro,sdteb,sdtwb,sdtnb,sdtsb,                          &
               ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,rhostri,        &
               csndsq,                                                  &
               uforce,vforce,wforce,pforce,ptforce,                     &
               x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv,          &
               trigs1,trigs2,ifax1,ifax2,                               &
               wsave1,wsave2,vwork1,vwork2,                             &
               exbcbuf,rhofct,                                          &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,defsq,      &
               tem10,tem11,tem12,tem13,tem14,tem15,tem16,tem17,         &
               tem18,tem1_0,tem2_0,tem3_0)
!call test_dump (pprt(1,1,1,2),"XXXAsmlstep_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,3),"XXXAsmlstep_pprt3",nx,ny,nz,0,1)

  CALL set_acct(misc_acct)

  IF( lvldbg >= 1 ) THEN
    CALL checkuhx(                  nx,ny,nz,1,nx,1,ny-1,1,nz-1, 'uxbig', tem1)
    CALL checkuhy(                  nx,ny,nz,1,nx,1,ny-1,1,nz-1, 'uybig', tem1)
    CALL checkvhx(                  nx,ny,nz,1,nx-1,1,ny,1,nz-1, 'vxbig', tem1)
    CALL checkvhy(                  nx,ny,nz,1,nx-1,1,ny,1,nz-1, 'vybig', tem1)
    CALL checkwhx(                  nx,ny,nz,1,nx-1,1,ny-1,1,nz, 'wxbig', tem1)
    CALL checkwhy(                  nx,ny,nz,1,nx-1,1,ny-1,1,nz, 'wybig', tem1)
    CALL checkshx(pprt(1,1,1,tfuture),                                  &
                  nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pxbig', tem1)
    CALL checkshy(pprt(1,1,1,tfuture),                                  &
                  nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pybig', tem1)
    CALL checkshx(ptprt(1,1,1,tfuture),                                 &
                  nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptxbig', tem1)
    CALL checkshy(ptprt(1,1,1,tfuture),                                 &
                  nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptybig', tem1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Since wcont was reset in SMLSTP. Now need to re-calculate its
!  value at tpresent. Which will be used in SOLVQV and SOLVQ.
!
!-----------------------------------------------------------------------
!
  CALL set_acct(tinteg_acct)

  CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)

!call test_dump (w(1,1,1,tpresent),"XXX2Bwc_w1",nx,ny,nz,1,0)
  CALL wcontra(nx,ny,nz,                                                &
               u(1,1,1,tpresent),v(1,1,1,tpresent),                     &
               w(1,1,1,tpresent),mapfct,j1,j2,j3,aj3z,                  &
               rhostr,tem1,tem2,tem3,wcont,tem4,tem5)
!
!-----------------------------------------------------------------------
!
!  Since pressure was reset in SMLSTP. Now need to re-calculate its
!  value at tfuture which will be used in microphysics.
!
!-----------------------------------------------------------------------
!
  CALL setppi(nx,ny,nz,nt,tfuture,pprt,pbar,ppi)
!call test_dump (pprt(1,1,1,2),"XXXAsetppi_pprt2",nx,ny,nz,0,1)
!call test_dump (pprt(1,1,1,3),"XXXAsetppi_pprt3",nx,ny,nz,0,1)

  IF( lvldbg >= 1 ) THEN
    CALL acct_interrupt(misc_acct)

    CALL checkwhx(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1,                 &
                  'wcxbig', tem1)
    CALL checkwhy(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1,                 &
                  'wcybig', tem1)
    CALL acct_stop_inter
  END IF
!
!-----------------------------------------------------------------------
!
!  Integrate the liquid water substance continuity equations
!  forward one timestep.
!
!  Sources and sinks due to phase changes and radiative processes
!  are handled separately.
!
!  If the simulation is for dry dynamics, then we skip over
!  the moisture computations to save computer resources.
!  The flag for this choice is "moist", i.e.,
!
!  Moist = 0 for dry run
!  Moist = 1 for moist run
!
!-----------------------------------------------------------------------
!

  IF( moist == 1) THEN   ! Determine if this run is dry or moist.

!
!-----------------------------------------------------------------------
!
!  Water vapor equation. Array uforce and vforce are used as work
!  arrays.
!  Before we call the solve routine for the water varaibles, we
!  need to compute ustr,vstr,wstr and store them into tem9,10,and 11
!  These arrays cannot be altered and will be used in the
!  following six moisture solve calls.
!
!
!-----------------------------------------------------------------------
!
    IF (sadvopt == 1 .OR. sadvopt == 2 .OR. sadvopt == 3 ) THEN
                             ! 2nd or 4th order advection

      CALL rhouvw(nx,ny,nz,rhostr,tem1,tem2,tem3)

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            tem9(i,j,k)=u(i,j,k,2)*tem1(i,j,k)
          END DO
        END DO
      END DO

      DO k=1,nz-1
        DO j=1,ny
          DO i=1,nx-1
            tem10(i,j,k)=v(i,j,k,2)*tem2(i,j,k)
          END DO
        END DO
      END DO

      DO k=1,nz
        DO j=1,ny-1
          DO i=1,nx-1
            tem11(i,j,k)=wcont(i,j,k)*tem3(i,j,k)
          END DO
        END DO
      END DO

    ELSE IF( sadvopt == 4.OR.sadvopt == 5) THEN  ! FCT advection

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem3_0(i,j,k)=rhostr(i,j,k)
          END DO
        END DO
      END DO

      CALL acct_interrupt(bc_acct)

      CALL extndsbc(tem3_0,nx,ny,nz,0,ebc,wbc,nbc,sbc,tbc,bbc)

      CALL acct_stop_inter

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            tem9(i,j,k)=u(i,j,k,2)*(tem3_0(i-1,j,k)+tem3_0(i,j,k))      &
                      *mapfct(i,j,5)*0.5
          END DO
        END DO
      END DO

      DO k=1,nz-1
        DO j=1,ny
          DO i=1,nx-1
            tem10(i,j,k)=v(i,j,k,2)*(tem3_0(i,j-1,k)+tem3_0(i,j,k))     &
                      *mapfct(i,j,6)*0.5
          END DO
        END DO
      END DO

      DO k=1,nz
        DO j=1,ny-1
          DO i=1,nx-1
            tem11(i,j,k)=wcont(i,j,k)                                   &
                      *(tem3_0(i,j,k-1)+tem3_0(i,j,k))*0.5
          END DO
        END DO
      END DO

    END IF

!-----------------------------------------------------------------------
!
!  Water vapor equation.  Array uforce and vforce are used as work
!  arrays. Tem9,10,11 cannot be altered until after the last
!  call to a solvq routine.  They contain the appropriate
!  ustr,vstr, and wstr.
!  Water vapor equation. Array uforce and vforce are used as work
!  arrays.
!
!-----------------------------------------------------------------------

    CALL solvqv(nx,ny,nz, exbcbufsz, dtbig1,                            &
                 qv,u,v,wcont, tem9,tem10,tem11,                        &
                 sdteb,sdtwb,sdtnb,sdtsb,                               &
                 rhostr,rhostri,qvbar,kmh,kmv,rprntl,qvsflx,pbldpth,    &
                 x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,j3inv,             &
                 qcumsrc(1,1,1,1),                                      &
                 usflx,vsflx,ptsflx,ptsfc,qvsfc,ptbar,ptprt,            &
                 exbcbuf, bcrlx,                                        &
                 uforce,vforce,                                         &
                 tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,               &
                 tem1_0,tem2_0,tem3_0)

!-----------------------------------------------------------------------
!
!  Cloud water equation.  Array uforce and vforce are used as work
!  arrays. Tem9,10,11 cannot be altered until after the last
!  call to a solvq routine.  They contain the appropriate
!  ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------

    qflag = 2
    CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag,                      &
               qc,u,v,wcont, tem9,tem10,tem11,                          &
               sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri,                 &
               kmh,kmv,rprntl,x,y,z,zp,mapfct,                          &
               j1,j2,j3,aj3x,aj3y,j3inv,                                &
               qcumsrc(1,1,1,2),qcumsrc(1,1,1,3),                       &
               qcumsrc(1,1,1,4),qcumsrc(1,1,1,5),                       &
               exbcbuf, bcrlx,                                          &
               uforce,vforce,                                           &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,                 &
               tem1_0,tem2_0,tem3_0)

    IF( mphyopt /= 0) THEN
!
!-----------------------------------------------------------------------
!
!  Rain water equation.  Array uforce and vforce are used as work
!  arrays. Tem9,10,11 cannot be altered until after the last
!  call to a solvq routine.  They contain the appropriate
!  ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
      qflag = 3
      CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag,                    &
                 qr,u,v,wcont, tem9,tem10,tem11,                        &
                 sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri,               &
                 kmh,kmv,rprntl,x,y,z,zp,mapfct,                        &
                 j1,j2,j3,aj3x,aj3y,j3inv,                              &
                 qcumsrc(1,1,1,2),qcumsrc(1,1,1,3),                     &
                 qcumsrc(1,1,1,4),qcumsrc(1,1,1,5),                     &
                 exbcbuf, bcrlx,                                        &
                 uforce,vforce,                                         &
                 tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,               &
                 tem1_0,tem2_0,tem3_0)

    END IF

    IF( lvldbg >= 1 ) THEN
      CALL acct_interrupt(misc_acct)
      CALL checkshx(qv(1,1,1,tfuture),                                  &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qvx', tem1)
      CALL checkshy(qv(1,1,1,tfuture),                                  &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qvy', tem1)
      CALL checkshx(qc(1,1,1,tfuture),                                  &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qcx', tem1)
      CALL checkshy(qc(1,1,1,tfuture),                                  &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qcy', tem1)
      CALL checkshx(qr(1,1,1,tfuture),                                  &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qrx', tem1)
      CALL checkshy(qr(1,1,1,tfuture),                                  &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'qry', tem1)
      CALL acct_stop_inter
    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Integrate the frozen water substance continuity equations
!  forward one timestep.
!
!  Sources and sinks due to phase changes and radiative processes
!  are handled separately.
!
!  If the simulation does not include ice processes then, to save
!  computer resources, we skip ice variable computations.
!  The flags for these options are "moist" and "ice", i.e.,
!
!  Moist = 0 dry run
!  Moist = 1, Ice = 0, moist run without ice processes
!  Moist = 1, Ice = 1, moist run with ice processes.
!
!-----------------------------------------------------------------------
!

  IF( moist == 1 .AND. ice == 1 ) THEN

!
!-----------------------------------------------------------------------
!
!  Cloud ice equation.  Array uforce and vforce are used as work
!  arrays. Tem9,10,11 cannot be altered until after the last
!  call to a solvq routine.  They contain the appropriate
!  ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
    qflag = 4
    CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag,                      &
               qi,u,v,wcont, tem9,tem10,tem11,                          &
               sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri,                 &
               kmh,kmv,rprntl,x,y,z,zp,mapfct,                          &
               j1,j2,j3,aj3x,aj3y,j3inv,                                &
               qcumsrc(1,1,1,2),qcumsrc(1,1,1,3),                       &
               qcumsrc(1,1,1,4),qcumsrc(1,1,1,5),                       &
               exbcbuf, bcrlx,                                          &
               uforce,vforce,                                           &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,                 &
               tem1_0,tem2_0,tem3_0)

!
!-----------------------------------------------------------------------
!
!  Snow equation.  Array uforce and vforce are used as work
!  arrays. Tem9,10,11 cannot be altered until after the last
!  call to a solvq routine.  They contain the appropriate
!  ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
    qflag = 5
    CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag,                      &
               qs,u,v,wcont, tem9,tem10,tem11,                          &
               sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri,                 &
               kmh,kmv,rprntl,x,y,z,zp,mapfct,                          &
               j1,j2,j3,aj3x,aj3y,j3inv,                                &
               qcumsrc(1,1,1,2),qcumsrc(1,1,1,3),                       &
               qcumsrc(1,1,1,4),qcumsrc(1,1,1,5),                       &
               exbcbuf, bcrlx,                                          &
               uforce,vforce,                                           &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,                 &
               tem1_0,tem2_0,tem3_0)

!
!-----------------------------------------------------------------------
!
!  Hail equation.  Array uforce and vforce are used as work
!  arrays. Tem9,10,11 cannot be altered until after the last
!  call to a solvq routine.  They contain the appropriate
!  ustr,vstr, and wstr.
!
!-----------------------------------------------------------------------
!
    qflag = 6
    CALL solvq(nx,ny,nz, exbcbufsz, dtbig1, qflag,                      &
               qh,u,v,wcont, tem9,tem10,tem11,                          &
               sdteb,sdtwb,sdtnb,sdtsb, rhostr,rhostri,                 &
               kmh,kmv,rprntl,x,y,z,zp,mapfct,                          &
               j1,j2,j3,aj3x,aj3y,j3inv,                                &
               qcumsrc(1,1,1,2),qcumsrc(1,1,1,3),                       &
               qcumsrc(1,1,1,4),qcumsrc(1,1,1,5),                       &
               exbcbuf, bcrlx,                                          &
               uforce,vforce,                                           &
               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,                 &
               tem1_0,tem2_0,tem3_0)

  END IF

  RETURN
END SUBROUTINE tinteg
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SMLSTEP                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE smlstep(mptr, nx,ny,nz, exbcbufsz, dtbig1,dtsml1,            & 1,21
           u,v,w,wcont,pprt,ptprt,                                      &
           udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,             &
           wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,             &
           phydro,ptdteb,ptdtwb,ptdtnb,ptdtsb,                          &
           ubar,vbar,ptbar,pbar,ptbari,pbari,rhostr,rhostri,            &
           csndsq,                                                      &
           uforce,vforce,wforce,pforce,ptforce,                         &
           x,y,z,zp,mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv,              &
           trigs1,trigs2,ifax1,ifax2,                                   &
           wsave1,wsave2,vwork1,vwork2,                                 &
           exbcbuf,rhofct,                                              &
           rhostru,rhostrv,rhostrw,                                     &
           tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,tem10,          &
           tem11,tem12,tem13,tem14,tem15,tem16,tem17,tem18,tem19)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Coordinate the integration of the acoustically active parts of
!  the dynamic equations.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/21/92.
!
!  MODIFICATION HISTORY:
!
!  5/20/92 (K. Droegemeier and M. Xue)
!  Added full documentation.
!
!  4/10/93 (M. Xue & Hao Jin)
!  Add the terrain.
!
!  9/10/94 (D. Weber & Y. Lu)
!  Cleaned up documentation.
!
!  10/31/95 (D. Weber)
!  Added trigs1,trigs2,ifax1,ifax2 for use in the upper w-p
!  radiation condition.
!
!  1/25/96 (Donghai Wang & Yuhe Liu)
!  Added the map projection factor to ARPS governing equations.
!
!  07/22/97 (D. Weber)
!  Added wsave1,wsave2,vwork1,vwork2 for use in the even fft version
!  of the upper w-p radiation condition (fftopt=2).
!
!  10/21/97 (Donghai Wang)
!  Using total density (rho) in the calculation of the pressure
!  gradient force terms, and added the second order terms
!  in the linerized buoyancy terms.
!
!  11/05/97 (D. Weber)
!  Added phydro array for use in the bottom boundary condition for
!  perturbation pressure (hydrostatic).
!
!  11/06/97 (D. Weber)
!  Added three additional levels to the mapfct array.  The three
!  levels (4,5,6) represent the inverse of the first three in order.
!  The inverse map factors are computed to improve efficiency.
!
!  9/18/98 (D. Weber)
!  Added arrays aj3x,y,z, rhostri,ptbari,pbari, and tem9-12 for
!  use in optimizing this routine.
!  Tem9-19 is used to store commonly used groups of constants.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    dtbig1   Large time step size for this call.
!    dtsml1   Small time step size for this call.
!
!    u        x component of velocity at times tpast and tpresent (m/s)
!    v        y component of velocity at times tpast and tpresent (m/s)
!    w        Vertical component of Cartesian velocity at times
!             tpast and tpresent (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!    pprt     Perturbation pressure at times tpast and tpresent (Pascal)
!    ptprt    Perturbation potential temperature at times tpast and
!             tpresent (K)
!
!    udteb    Time tendency of u field at east boundary (m/s**2)
!    udtwb    Time tendency of u field at west boundary (m/s**2)
!    udtnb    Time tendency of u field at north boundary (m/s**2)
!    udtsb    Time tendency of u field at south boundary (m/s**2)
!
!    vdteb    Time tendency of v field at east boundary (m/s**2)
!    vdtwb    Time tendency of v field at west boundary (m/s**2)
!    vdtnb    Time tendency of v field at north boundary (m/s**2)
!    vdtsb    Time tendency of v field at south boundary (m/s**2)
!
!    wdteb    Time tendency of w field at east boundary (m/s**2)
!    wdtwb    Time tendency of w field at west boundary (m/s**2)
!    wdtnb    Time tendency of w field at north boundary (m/s**2)
!    wdtsb    Time tendency of w field at south boundary (m/s**2)
!
!    pdteb    Time tendency of pprt field at east boundary (Pascal/s)
!    pdtwb    Time tendency of pprt field at west boundary (Pascal/s)
!    pdtnb    Time tendency of pprt field at north boundary (Pascal/s)
!    pdtsb    Time tendency of pprt field at south boundary (Pascal/s)
!
!    ptdteb   Time tendency of ptprt field at east boundary (K/s)
!    ptdtwb   Time tendency of ptprt field at west boundary (K/s)
!    ptdtnb   Time tendency of ptprt field at north boundary(K/s)
!    ptdtsb   Time tendency of ptprt field at south boundary(K/s)
!
!    phydro   Big time step forcing term for use in computing the
!             hydrostatic pressure at k=1.
!
!    ubar     Base state zonal velocity component (m/s)
!    vbar     Base state meridional velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    ptbari   Inverse base state potential temperature (K)
!    pbari    Inverse base state pressure (Pascal)
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    rhostri  Inverse base state density rhobar times j3 (kg/m**3)
!    csndsq   Sound wave speed squared.
!
!    uforce   Acoustically inactive forcing terms in u-eq. (kg/(m*s)**2)
!    vforce   Acoustically inactive forcing terms in v-eq. (kg/(m*s)**2)
!    wforce   Acoustically inactive forcing terms in w-eq. (kg/(m*s)**2)
!    pforce   Acoustically inactive forcing terms in p-eq. (Pascal/s)
!    ptforce  Gravity wave inactive forcing terms in pt-eq. (K*kg/(m**3*s))
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!
!    mapfct   Map factors at scalar, u and v points
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    aj3x     Avgx of the coordinate transformation Jacobian  d(zp)/dz
!    aj3y     Avgy of the coordinate transformation Jacobian  d(zp)/dz
!    aj3z     Avgz of the coordinate transformation Jacobian  d(zp)/dz
!    j3inv    Inverse of the coordinate transformation j3
!    trigs1   Array containing pre-computed trig function for
!             fftopt=1.
!    trigs2   Array containing pre-computed trig function for
!             fftopt=1.
!    ifax1    Array containing the factors of nx for fftopt=1.
!    ifax2    Array containing the factors of ny for fftopt=1.
!
!    vwork1   2-D work array for fftopt = 2.
!    vwork2   2-D work array for fftopt = 2.
!    wsave1   Work array for fftopt = 2.
!    wsave2   Work array for fftopt = 2.
!
!    rhostru  Averaged rhostr at u points (kg/m**3).
!    rhostrv  Averaged rhostr at v points (kg/m**3).
!    rhostrw  Averaged rhostr at w points (kg/m**3).
!
!  OUTPUT:
!
!    u        x component of velocity at time tfuture (m/s)
!    v        y component of velocity at time tfuture (m/s)
!    w        Vertical component of Cartesian velocity at tfuture (m/s)
!    pprt     Perturbation pressure at time tfuture (Pascal)
!    ptprt    Perturbation potential temperature at time tfuture (K)
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array
!    tem2     Temporary work array
!    tem3     Temporary work array
!    tem4     Temporary work array
!    tem5     Temporary work array
!    tem6     Temporary work array
!    tem7     Temporary work array
!    tem8     Temporary work array
!    tem9     Temporary work array
!    tem10    Temporary work array
!    tem11    Temporary work array
!    tem12    Temporary work array
!    tem13    Temporary work array
!    tem14    Temporary work array
!    tem15    Temporary work array
!    tem16    Temporary work array
!    tem17    Temporary work array
!    tem18    Temporary work array
!    tem19    Temporary work array
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE        ! Force explicit declarations

  INCLUDE 'timelvls.inc'

  INTEGER :: mptr 

  INTEGER :: nx, ny, nz        ! Number of grid points in 3 directions

  REAL :: dtbig1               ! Large time step size for this call.
  REAL :: dtsml1               ! Small time step size for this call.

  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
  REAL :: wcont (nx,ny,nz)     ! Contravariant vertical velocity (m/s)
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)

  REAL :: udteb (ny,nz)        ! Time tendency of u at e-boundary (m/s**2)
  REAL :: udtwb (ny,nz)        ! Time tendency of u at w-boundary (m/s**2)
  REAL :: udtnb (nx,nz)        ! Time tendency of u at n-boundary (m/s**2)
  REAL :: udtsb (nx,nz)        ! Time tendency of u at s-boundary (m/s**2)

  REAL :: vdteb (ny,nz)        ! Time tendency of v at e-boundary (m/s**2)
  REAL :: vdtwb (ny,nz)        ! Time tendency of v at w-boundary (m/s**2)
  REAL :: vdtnb (nx,nz)        ! Time tendency of v at n-boundary (m/s**2)
  REAL :: vdtsb (nx,nz)        ! Time tendency of v at s-boundary (m/s**2)

  REAL :: wdteb (ny,nz)        ! Time tendency of w at e-boundary (m/s**2)
  REAL :: wdtwb (ny,nz)        ! Time tendency of w at w-boundary (m/s**2)
  REAL :: wdtnb (nx,nz)        ! Time tendency of w at n-boundary (m/s**2)
  REAL :: wdtsb (nx,nz)        ! Time tendency of w at s-boundary (m/s**2)

  REAL :: pdteb (ny,nz)        ! Time tendency of pprt at e-boundary (Pascal/s)
  REAL :: pdtwb (ny,nz)        ! Time tendency of pprt at w-boundary (Pascal/s)
  REAL :: pdtnb (nx,nz)        ! Time tendency of pprt at n-boundary (Pascal/s)
  REAL :: pdtsb (nx,nz)        ! Time tendency of pprt at s-boundary (Pascal/s)

  REAL :: phydro(nx,ny)        ! Big time step forcing for computing
                               ! hydrostatic pprt at k=1.

  REAL :: ptdteb(ny,nz)        ! Time tendency of ptprt at e-boundary (K/s)
  REAL :: ptdtwb(ny,nz)        ! Time tendency of ptprt at w-boundary (K/s)
  REAL :: ptdtnb(nx,nz)        ! Time tendency of ptprt at n-boundary (K/s)
  REAL :: ptdtsb(nx,nz)        ! Time tendency of ptprt at s-boundary (K/s)

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: ptbari(nx,ny,nz)     ! Inverse base state pot. temperature (K)
  REAL :: pbari (nx,ny,nz)     ! Inverse base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: rhostri(nx,ny,nz)    ! Inverse base state density rhobar times j3.
  REAL :: csndsq(nx,ny,nz)     ! Sound wave speed squared.

  REAL :: x     (nx)           ! The x-coord. of the physical and
                               ! computational grid. Defined at u-point.
  REAL :: y     (ny)           ! The y-coord. of the physical and
                               ! computational grid. Defined at v-point.
  REAL :: z     (nz)           ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of the staggered grid.

  REAL :: mapfct(nx,ny,8)      ! Map factors at scalar, u and v points

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! - d( zp )/d( x ).
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! - d( zp )/d( y ).
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).
  REAL :: aj3x  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE X-DIR.
  REAL :: aj3y  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Y-DIR.
  REAL :: aj3z  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.
  REAL :: j3inv (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).

  REAL :: trigs1(3*(nx-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  REAL :: trigs2(3*(ny-1)/2+1) ! Array containing pre-computed trig
                               ! function for fftopt=1.
  INTEGER :: ifax1(13)         ! Array containing the factors of nx for
                               ! fftopt=1.
  INTEGER :: ifax2(13)         ! Array containing the factors of ny for
                               ! fftopt=1.

  REAL :: vwork1 (nx+1,ny+1)   ! 2-D work array for fftopt = 2.
  REAL :: vwork2 (ny,nx+1)     ! 2-D work array for fftopt = 2.
  REAL :: wsave1 (3*(ny-1)+15) ! Work array for fftopt = 2.
  REAL :: wsave2 (3*(nx-1)+15) ! Work array for fftopt = 2.

  REAL :: uforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in u-momentum equation (kg/(m*s)**2)
  REAL :: vforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in v-momentum equation (kg/(m*s)**2)
  REAL :: wforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in w-momentum equation (kg/(m*s)**2)
  REAL :: pforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                               ! in pressure equation (Pascal/s)
  REAL :: ptforce(nx,ny,nz)    ! Gravity wave inactive forcing terms
                               ! in potential temperature eq. (K*kg/(m**3*s))

  REAL :: rhostru(nx,ny,nz)    ! Averaged rhostr at u points (kg/m**3).
  REAL :: rhostrv(nx,ny,nz)    ! Averaged rhostr at v points (kg/m**3).
  REAL :: rhostrw(nx,ny,nz)    ! Averaged rhostr at w points (kg/m**3).

  INTEGER :: exbcbufsz         ! EXBC buffer size
  REAL :: exbcbuf( exbcbufsz ) ! EXBC buffer array

  REAL :: rhofct(nx,ny,nz)     ! rho-factor: rhobar/rho

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array
  REAL :: tem2  (nx,ny,nz)     ! Temporary work array
  REAL :: tem3  (nx,ny,nz)     ! Temporary work array
  REAL :: tem4  (nx,ny,nz)     ! Temporary work array
  REAL :: tem5  (nx,ny,nz)     ! Temporary work array
  REAL :: tem6  (nx,ny,nz)     ! Temporary work array
  REAL :: tem7  (nx,ny,nz)     ! Temporary work array
  REAL :: tem8  (nx,ny,nz)     ! Temporary work array
  REAL :: tem9  (nx,ny,nz)     ! Temporary work array
  REAL :: tem10 (nx,ny,nz)     ! Temporary work array
  REAL :: tem11 (nx,ny,nz)     ! Temporary work array
  REAL :: tem12 (nx,ny,nz)     ! Temporary work array
  REAL :: tem13 (nx,ny,nz)     ! Temporary work array
  REAL :: tem14 (nx,ny,nz)     ! Temporary work array
  REAL :: tem15 (nx,ny,nz)     ! Temporary work array
  REAL :: tem16 (nx,ny,nz)     ! Temporary work array
  REAL :: tem17 (nx,ny,nz)     ! Temporary work array
  REAL :: tem18 (nx,ny,nz)     ! Temporary work array
  REAL :: tem19 (nx,ny,nz)     ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
  INTEGER :: ismstp

  REAL :: curtsml               ! Current time during small time step
                                ! integration.
  REAL :: tem,tema,temb
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'bndry.inc'
  INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  Integrate the momentum equations using a forward-in-time (relative
!  to the pressure gradient terms) scheme and the pressure equation
!  using a backward-in-time (relative to the mass divergence term)
!  scheme.
!
!  dtsml = (2*dtbig)/nsmstp is used so that after nsmstp small time
!  step integrations, the fields are brought from time t-dtbig
!  to t+dtbig, i.e. from time tpast to time future.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  The first time through the small timestep loop, set the values
!  at time tfuture equal to those at tpast.  By doing this, we
!  simplify the forward-in-time integration because the equations
!  take the form:
!
!  u(i,j,k,future) = u(i,j,k,future)
!    :                + dtsml*( acoustic +  non-acoustic forcing terms )
!
!  i.e., we only use a single time level in this two-time level
!  scheme and automatically finish with u at time tfuture.
!
!-----------------------------------------------------------------------
!
!call test_dump(w(1,1,tfuture),"XXXsml_wfut",nx,ny,nz,3,1)
!call test_dump(w(1,1,tpast),"XXXsml_wpast",nx,ny,nz,3,1)
!call test_dump(w(1,1,tpresent),"XXXsml_wpres",nx,ny,nz,3,1)
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        u    (i,j,k,tfuture)=u    (i,j,k,tpast)
        v    (i,j,k,tfuture)=v    (i,j,k,tpast)
        w    (i,j,k,tfuture)=w    (i,j,k,tpast)
        pprt (i,j,k,tfuture)=pprt (i,j,k,tpast)
      END DO
    END DO
  END DO

  IF(ptsmlstp == 1) THEN

    IF( sadvopt == 4 ) THEN

!-----------------------------------------------------------------------
!
!  When FCT is used to advect ptprt, all (e.g., advection,
!  mixing) except for ptbar advection terms were evaluated at tpresent,
!  and the large time step integration is forward in time.
!  We have to do somthing special here, so that the same ptprt(tfuture)
!  results had the ptbar advection not been included when ptprt is
!  is updated in small steps. Therefore, we require
!
!         (ptprt(tfuture)-ptprt(tpast))/(2*dtbig1)
!       = (ptprt(tfuture)-ptprt(tpresent))/dtbig
!
!-----------------------------------------------------------------------

      tem = 2*dtbig1/dtbig

      IF(nestgrd == 1.AND.mgrid /= 1) THEN
                ! Don't touch the boundary values at tpast for nested grids.

        DO k=1,nz-1
          DO j=2,ny-2
            DO i=2,nx-2
              ptprt(i,j,k,tpast)=ptprt(i,j,k,tfuture)                   &
                   -(ptprt(i,j,k,tfuture)-ptprt(i,j,k,tpresent))*tem
            END DO
          END DO
        END DO

      ELSE

        DO k=1,nz-1
          DO j=1,ny-1
            DO i=1,nx-1
              ptprt(i,j,k,tpast)=ptprt(i,j,k,tfuture)                   &
                   -(ptprt(i,j,k,tfuture)-ptprt(i,j,k,tpresent))*tem
            END DO
          END DO
        END DO

      END IF

    END IF

    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          ptprt(i,j,k,tfuture)=ptprt(i,j,k,tpast)
        END DO
      END DO
    END DO

  END IF
!
!-----------------------------------------------------------------------
!
!  Calculate wcont at future time level.
!
!-----------------------------------------------------------------------
!
  CALL rhouvw(nx,ny,nz,rhostr,rhostru,rhostrv,rhostrw)

!call test_dump (w(1,1,1,tfuture),"XXX3Bwc_w1",nx,ny,nz,1,0)
  CALL wcontra(nx,ny,nz,                                                &
             u(1,1,1,tfuture),v(1,1,1,tfuture),                         &
             w(1,1,1,tfuture),mapfct,j1,j2,j3,aj3z,                     &
             rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)

  IF( lvldbg >= 3 ) THEN
    CALL checkwhx(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1,                 &
                  'wcxsml', tem1)
    CALL checkwhy(wcont, nx,ny,nz,1,nx-1,1,ny-1,2,nz-1,                 &
                  'wcysml', tem1)
  END IF

!-----------------------------------------------------------------------
!
!  Compute a number of static variables (wrt small time step) and
!  pass into the appropriate subroutines.
!  NOTE:  ALL OF THE COMPUTATIONS ASSUME THAT THE VARIABLES USED
!         ARE CONSTANT DURING THE SMALL TIME STEP!!!!!!!
!         THEY CANNOT BE OVERWRITTEN UNTIL AFTER THE W-P SOLVER.
!
!-----------------------------------------------------------------------

  temb = 0.5*dtsml1
  DO k=2,nz-2
    DO j=1,ny-1
      DO i=2,nx-1
        tema = 1.0/rhostru(i,j,k)
        tem13(i,j,k) =temb*(rhofct(i,j,k)+rhofct(i-1,j,k))              &
                          *mapfct(i,j,2)*tema
        uforce(i,j,k) = dtsml1*uforce(i,j,k)*tema
      END DO
    END DO
  END DO
!call test_dump (uforce,"XXX1uforce",nx,ny,nz,1,0)
!call test_dump (rhostru,"XXXrhostru",nx,ny,nz,1,0)
!call test_dump (tem13,"XXXdtmrxrsu_tem13",nx,ny,nz,1,0)
!call test_dump2(rhofct,"XXXdtmrxrsu_rhofct",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)

  DO k=2,nz-2
    DO j=2,ny-1
      DO i=1,nx-1
        tema = 1.0/rhostrv(i,j,k)
        tem14(i,j,k) =temb*(rhofct(i,j,k)+rhofct(i,j-1,k))              &
                           *mapfct(i,j,3)*tema
        vforce(i,j,k) = dtsml1*vforce(i,j,k)*tema
      END DO
    END DO
  END DO

  DO k=2,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem17(i,j,k)=1.0/(pbar(i,j,k)+pbar(i,j,k-1)) ! used in solvwp
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem12(i,j,k)=1.0/rhostrw(i,j,k)              ! used in solvwp
        tem9 (i,j,k)=rhostr(i,j,k)*j3inv(i,j,k)      ! used in solvwp
        tem10(i,j,k)=csndsq(i,j,k)*j3inv(i,j,k)*tem9(i,j,k)  ! solvwp
        tem11(i,j,k)=rhostr(i,j,k)*pbari(i,j,k)      ! used in solvwp
        tem18(i,j,k)=rhostr(i,j,k)*ptbari(i,j,k)     ! used in solvwp
      END DO
    END DO
  END DO

  IF(vimplct == 0)THEN  ! pre-compute wforce, etc.. for wp-ex.
    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem19(i,j,k) =temb*(rhofct(i,j,k)+rhofct(i,j,k-1))            &
                             *tem12(i,j,k)
          wforce(i,j,k) = dtsml1*wforce(i,j,k)*tem12(i,j,k)
        END DO
      END DO
    END DO
  END IF

  tema = dtsml1*dxinv
  DO k=2,nz-2
    DO j=1,ny-1
      DO i=1,nx
        tem15(i,j,k)=tema*aj3x(i,j,k)*mapfct(i,j,5) ! used in solvwp's
      END DO
    END DO
  END DO

  temb = dtsml1*dyinv
  DO k=2,nz-2
    DO j=1,ny
      DO i=1,nx-1
        tem16(i,j,k)=temb*aj3y(i,j,k)*mapfct(i,j,6) ! used in solvwp's
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!    For peqopt=2, set
!
!    pforce=pforce (containging exbc damping only)
!      + rhostr*cs**2/ptbar*d(pt)/dt
!       =pforce (containging exbc damping only)
!      + cs**2/ptbar * (ptforce + buoyancy if not already included).
!
!    Note that adiabet heating/colling is not included in ptforce,
!    therefore is not accounted for in p equation.
!
!-----------------------------------------------------------------------
!
  IF(peqopt == 2) THEN

    IF( ptsmlstp == 0 ) THEN
      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            pforce(i,j,k)=pforce(i,j,k)                                 &
                         +csndsq(i,j,k)*ptforce(i,j,k)*ptbari(i,j,k)
          END DO
        END DO
      END DO
    ELSE
      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            pforce(i,j,k)=pforce(i,j,k)+csndsq(i,j,k)*ptbari(i,j,k)     &
                     *(ptforce(i,j,k)-tem9(i,j,k)*                      &
                       0.25*(ptbar(i,j,k+1)-ptbar(i,j,k-1))*dzinv*      &
                       (w(i,j,k+1,2)+w(i,j,k,2)) )
          END DO
        END DO
      END DO
    END IF

  END IF

  DO k=2,nz-2   ! preparing pforce...for all solvwp codes...
    DO j=1,ny-1
      DO i=1,nx-1
        pforce(i,j,k) = dtsml1*pforce(i,j,k)*j3inv(i,j,k)
      END DO
    END DO
  END DO


  DO ismstp =1, nsmstp
!
!-----------------------------------------------------------------------
!
!  Calculate the current time within small time step. It will be
!  used to calculate the external boundary fields.
!
!-----------------------------------------------------------------------
!
    curtsml = curtim + dtbig - 2*dtbig1 + ismstp*dtsml1
!
!-----------------------------------------------------------------------
!
!  Advance the three momentum equations one small timestep.
!  (Arguments such as u(1,1,1,tfuture) indicate the passing of
!  array values at time tfuture into the subroutine. Inside the
!  subroutine, these arrays are defined at a single time level.
!  u, v, w at tfuture and wcont are updated in time exiting
!  this routine.
!
!  Advance the pressure (p at tfuture) one small timestep.
!
!  tem7 is used to transfer wpgrad between subroutine SOLVUV and
!  subroutine SOLVWPIM or SOLVWPEX.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Calculate the boundary time tendencies of u and v using Klemp &
!  Wilhelmson type open boundary condition (rbcopt=1).
!
!-----------------------------------------------------------------------
!
    IF( rbcopt == 1 .AND. mgrid == 1 ) THEN

      CALL bkwsmldt(nx,ny,nz, u(1,1,1,tfuture),v(1,1,1,tfuture),        &
                    ubar,vbar,udteb,udtwb,vdtnb,vdtsb)

    END IF

    CALL solvuv(mptr, nx,ny,nz, exbcbufsz, dtsml1, curtsml,             &
                u(1,1,1,tfuture),v(1,1,1,tfuture),                      &
                wcont,pprt(1,1,1,tfuture),                              &
                udteb,udtwb,udtnb,udtsb,                                &
                vdteb,vdtwb,vdtnb,vdtsb,                                &
                rhostr, uforce,vforce,ubar,vbar,                        &
                x,y,z,zp,mapfct, j1,j2,j3,j3inv,                        &
                exbcbuf,rhofct,                                         &
                rhostru,rhostrv,rhostrw,tem13,tem14,tem7,               &
                tem1,tem2,tem3,tem4,tem5)
!
!-----------------------------------------------------------------------
!
!  Warning: tem7 is carrying wpgrad can not be overwritten before
!  the subroutine for w and p-equations are called.
!
!-----------------------------------------------------------------------
!

    IF( lvldbg >= 3 ) THEN
      CALL checkuhx(                    nx,ny,nz,1,nx,  1,ny-1,1,nz-1, 'uxsml', tem1)
      CALL checkuhy(                    nx,ny,nz,1,nx,  1,ny-1,1,nz-1, 'uysml', tem1)
      CALL checkvhx(                    nx,ny,nz,1,nx-1,1,ny,  1,nz-1, 'vxsml', tem1)
      CALL checkvhy(                    nx,ny,nz,1,nx-1,1,ny,  1,nz-1, 'vysml', tem1)
    END IF
!
!-----------------------------------------------------------------------
!
!  Advance the w and pressure (at tfuture) one small timestep.
!
!-----------------------------------------------------------------------
!
    IF( vimplct == 0 ) THEN   ! Not implicit
!
!-----------------------------------------------------------------------
!
!  Integrate the w and pressure equations one small time step
!  using explicit scheme.
!
!  NOTE: tem11,tem18,tem12,tem15,tem16,tem9 are assumed to be
!        constant during the small time step!!!
!
!-----------------------------------------------------------------------
!
      CALL solvwpex(mptr, nx,ny,nz, exbcbufsz, dtsml1, curtsml,         &
                    u(1,1,1,tfuture),v(1,1,1,tfuture),                  &
                    w(1,1,1,tfuture),wcont,                             &
                    ptprt(1,1,1,tfuture),pprt(1,1,1,tfuture),phydro,    &
                    wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,    &
                    rhostr,ptbar,ptbari,pbari,csndsq,                   &
                    wforce,tem7,pforce,                                 &
                    x,y,z,zp, mapfct, j1,j2,j3,aj3x,aj3y,aj3z,j3inv,    &
                    rhostru,rhostrv,rhostrw,                            &
                    exbcbuf,rhofct,                                     &
                    tem1,tem2,tem3,tem4,tem5,                           &
                    tem11,tem18,tem19,tem15,tem16,tem9)
!call test_dump(w(1,1,tfuture),"XXXAsolvwpex_w",nx,ny,nz,3,1)
    ELSE   ! Implicit
!
!-----------------------------------------------------------------------
!
!  Integrate the w and pressure equations one small time step
!  using vertically implicit scheme.
!
!  NOTE: tem9,tem10,tem11,tem12,tem15,tem16,tem17,tem18
!        are assumed to be constant during the small time step!!!
!
!-----------------------------------------------------------------------
!
      IF(peqopt == 1) THEN

        CALL solvwpim(mptr,nx,ny,nz,exbcbufsz, dtsml1, curtsml,         &
                    u(1,1,1,tfuture),v(1,1,1,tfuture),                  &
                    w(1,1,1,tfuture),wcont,                             &
                    ptprt(1,1,1,tfuture),pprt(1,1,1,tfuture),phydro,    &
                    wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,    &
                    rhostr,ptbar,ptbari,pbari,csndsq,                   &
                    wforce,tem7,pforce,                                 &
                    x,y,z,zp,mapfct,                                    &
                    j1,j2,j3,aj3x,aj3y,aj3z,j3inv,                      &
                    trigs1,trigs2,ifax1,ifax2,                          &
                    wsave1,wsave2,vwork1,vwork2,                        &
                    rhostru,rhostrv,rhostrw,                            &
                    exbcbuf,rhofct,                                     &
                    tem1,tem2,tem3,tem4,tem5,tem6,tem8,                 &
                    tem9,tem10,tem11,tem12,tem15,tem16,tem17)
!call test_dump(w(1,1,tfuture),"XXXAsolvwpim_w",nx,ny,nz,3,1)

      ELSE

        CALL solvwpim1(nx,ny,nz,exbcbufsz, dtsml1, curtsml,             &
                    u(1,1,1,tfuture),v(1,1,1,tfuture),                  &
                    w(1,1,1,tfuture),wcont,                             &
                    ptprt(1,1,1,tfuture),pprt(1,1,1,tfuture),phydro,    &
                    wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,    &
                    rhostr,ptbar,ptbari,pbari,csndsq,                   &
                    wforce,tem7,pforce,                                 &
                    x,y,z,zp,mapfct,                                    &
                    j1,j2,j3,aj3z,j3inv,                                &
                    trigs1,trigs2,ifax1,ifax2,                          &
                    wsave1,wsave2,vwork1,vwork2,                        &
                    rhostru,rhostrv,rhostrw,                            &
                    exbcbuf,rhofct,                                     &
                    tem1,tem2,tem3,tem4,tem5,tem6,tem8,                 &
                    tem11,tem18,tem12,tem17,tem9)
!call test_dump(w(1,1,tfuture),"XXXAsolvwpim_w",nx,ny,nz,3,1)
      END IF

    END IF

    IF( lvldbg >= 3 ) THEN
      CALL checkwhx(                    nx,ny,nz,1,nx-1,1,ny-1,1,nz,   'wxsml', tem1)
      CALL checkwhy(                    nx,ny,nz,1,nx-1,1,ny-1,1,nz,   'wysml', tem1)
      CALL checkshx(pprt(1,1,1,tfuture),                                &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pxsml', tem1)
      CALL checkshy(pprt(1,1,1,tfuture),                                &
                    nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'pysml', tem1)
    END IF

!-----------------------------------------------------------------------
!
!  Integrate the potential temperature equation one small time step
!  when ptsmlstp=1.
!
!-----------------------------------------------------------------------
!

    IF( ptsmlstp == 1 ) THEN

      CALL solvpt_sml(nx,ny,nz, exbcbufsz, dtbig1,dtsml1,curtsml,       &
                      ptprt,w, ptdteb,ptdtwb,ptdtnb,ptdtsb,             &
                      ptbar,rhostr,rhostri,rhostrw,j3,aj3z,             &
                      ptforce, exbcbuf,                                 &
                      tem1,tem2)

      IF( lvldbg >= 3) THEN
        CALL checkshx(ptprt(1,1,1,tfuture),                             &
                      nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptsmlx',tem1)
        CALL checkshy(ptprt(1,1,1,tfuture),                             &
                      nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, 'ptsmly',tem1)
      END IF

    END IF


  END DO

!
!-----------------------------------------------------------------------
!
!  Integrate the potential temperature equation one large time step
!  when ptsmlstp=0.
!
!-----------------------------------------------------------------------
!
  IF( ptsmlstp == 0 ) THEN

    CALL solvpt_lrg(nx,ny,nz, exbcbufsz, dtbig1, ptprt,                 &
                    ptdteb,ptdtwb,ptdtnb,ptdtsb,rhostr,rhostri,         &
                    ptforce,exbcbuf,tem1)

  END IF

  RETURN
END SUBROUTINE smlstep
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE TFILT                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE tfilt(nx,ny,nz,                                              & 1,12
           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Apply the Asselin time filter to all time-dependent variables.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91.
!
!  MODIFICATION HISTORY:
!
!  5/20/92 (K. Droegemeier and M. Xue)
!  Added full documentation.
!
!  9/10/94 (D. Weber & Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at all time levels (m/s)
!    v        y component of velocity at all time levels (m/s)
!    w        Vertical component of Cartesian velocity at all time
!             levels (m/s)
!    ptprt    Perturbation potential temperature at all time levels (K)
!    pprt     Perturbation pressure at all time levels (Pascal)
!    qv       Water vapor specific humidity at all time levels (kg/kg)
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!    qi       Cloud ice mixing ratio at all time levels (kg/kg)
!    qs       Snow mixing ratio at all time levels (kg/kg)
!    qh       Hail mixing ratio at all time levels (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!  OUTPUT:
!
!    u        x component of velocity at time tpresent (m/s)
!    v        y component of velocity at time tpresent(m/s)
!    w        Vertical component of Cartesian velocity at tpresent (m/s)
!    ptprt    Perturbation potential temperature at time tpresent (K)
!    pprt     Perturbation pressure at time tpresent (Pascal)
!    qv       Water vapor specific humidity at time tpresent (kg/kg)
!    qc       Cloud water mixing ratio at time tpresent (kg/kg)
!    qr       Rainwater mixing ratio at time tpresent (kg/kg)
!    qi       Cloud ice mixing ratio at time tpresent (kg/kg)
!    qs       Snow mixing ratio at time tpresent (kg/kg)
!    qh       Hail mixing ratio at time tpresent (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'timelvls.inc'

  INTEGER :: nx, ny, nz        ! Number of grid points in 3 directions

  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)
  REAL :: tke   (nx,ny,nz,nt)  ! Turbulent kinetic energy
!
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  CALL aselin(u, nx,ny,nz, 1,nx,1,ny-1,1,nz-1, flteps)

  CALL aselin(v, nx,ny,nz, 1,nx-1,1,ny,1,nz-1, flteps)

  CALL aselin(w, nx,ny,nz, 1,nx-1,1,ny-1,1,nz, flteps)

  IF(sadvopt /= 4) CALL aselin(ptprt, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)

  CALL aselin(pprt, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)

  IF( tmixopt == 4) THEN
    IF(sadvopt /= 4) CALL aselin(tke, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
  END IF

  IF( moist /= 0 .AND. sadvopt /= 4 ) THEN

    CALL aselin(qv, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
    CALL aselin(qc, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
    CALL aselin(qr, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)

    IF( ice /= 0 ) THEN

      CALL aselin(qi, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
      CALL aselin(qs, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)
      CALL aselin(qh, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1, flteps)

    END IF

  END IF

  RETURN
END SUBROUTINE tfilt

!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ASELIN                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE aselin(a, nx,ny,nz, nx1,nx2,ny1,ny2,nz1,nz2, flteps) 12

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Apply Asselin time filter.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91.
!
!  MODIFICATION HISTORY:
!
!  5/20/92 (K. Droegemeier and M. Xue)
!  Added full documentation.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    a        A variable at three time levels.
!
!    nx1,nx2  i-index array bounds where the time filter is applied
!    ny1,ny2  j-index array bounds where the time filter is applied
!    nz1,nz2  k-index array bounds where the time filter is applied
!
!    flteps   The asselin time filter coefficient.
!             Defined in globcst.inc
!
!  OUTPUT:
!
!    a        The new array values at time tpresent after the time
!             filter is applied.
!
!-----------------------------------------------------------------------
!
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'timelvls.inc'

  INTEGER :: nx, ny, nz   ! Number of grid points in 3 directions

  REAL :: a(nx,ny,nz,nt)

  INTEGER :: nx1, nx2
  INTEGER :: ny1, ny2
  INTEGER :: nz1, nz2
  REAL :: flteps       ! Coefficient of the Asselin time filter
!
!-----------------------------------------------------------------------
!
!  Misc. local variable:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=nz1,nz2
    DO j=ny1,ny2
      DO i=nx1,nx2
        a(i,j,k,tpresent)=a(i,j,k,tpresent)+flteps*                     &
            (a(i,j,k,tfuture)-2*a(i,j,k,tpresent)+a(i,j,k,tpast))
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE aselin

!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE TFLIP                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE tflip(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke) 1,12

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Update the fields u, v, w, ptprt ,pprt ,qv ,qc ,qr ,qi ,qs ,qh and
!  tke in time. The fields at tpast are replaced by those at tpresent.
!  The fields at tpresent are replaced by those at tfuture. The fields
!  at tfuture are not changed.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91.
!
!  MODIFICATION HISTORY:
!
!  5/20/92 (K. Droegemeier and M. Xue)
!  Added full documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    u        x component of velocity at all time levels (m/s)
!    v        y component of velocity at all time levels (m/s)
!    w        Vertical component of Cartesian velocity at all time
!             levels (m/s)
!    ptprt    Perturbation potential temperature at all time levels (K)
!    pprt     Perturbation pressure at all time levels (Pascal)
!    qv       Water vapor specific humidity at all time levels (kg/kg)
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!    qi       Cloud ice mixing ratio at all time levels (kg/kg)
!    qs       Snow mixing ratio at all time levels (kg/kg)
!    qh       Hail mixing ratio at all time levels (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!  OUTPUT:
!
!    u        x-velocity at tpast and tpresent updated in time (m/s)
!    v        y-velocity at tpast and tpresent updated in time (m/s)
!    w        Vertical component of Cartesian velocity at tpast
!             and tpresent updated in time (m/s)
!    ptprt    Perturbation potential temperature at tpast and tpresent
!             updated in time (K)
!    pprt     Perturbation pressure at tpast and tpresent
!             updated in time (Pascal)
!    qv       Water vapor specific humidity at tpast and tpresent
!             updated in time (kg/kg)
!    qc       Cloud water mixing ratio at tpast and tpresent
!             updated in time (kg/kg)
!    qr       Rainwater mixing ratio at tpast and tpresent
!             updated in time (kg/kg)
!    qi       Cloud ice mixing ratio at tpast and tpresent
!             updated in time (kg/kg)
!    qs       Snow mixing ratio at tpast and tpresent
!             updated in time (kg/kg)
!    qh       Hail mixing ratio at tpast and tpresent
!             updated in time (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'timelvls.inc'

  INTEGER :: nx, ny, nz        ! Number of grid points in 3 directions

  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
!
  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)
  REAL :: tke   (nx,ny,nz,nt)  ! Turbulent kinetic energy
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  CALL tswap(nx,ny,nz, u)
  CALL tswap(nx,ny,nz, v)
  CALL tswap(nx,ny,nz, w)
  CALL tswap(nx,ny,nz, ptprt)
  CALL tswap(nx,ny,nz, pprt)

  IF( tmixopt == 4 ) THEN
    CALL tswap(nx,ny,nz,tke )
  END IF

  IF( moist /= 0 ) THEN

    CALL tswap(nx,ny,nz, qv)
    CALL tswap(nx,ny,nz, qc)
    CALL tswap(nx,ny,nz, qr)

    IF( ice /= 0 ) THEN

      CALL tswap(nx,ny,nz, qi)
      CALL tswap(nx,ny,nz, qs)
      CALL tswap(nx,ny,nz, qh)

    END IF

  END IF

  RETURN
END SUBROUTINE tflip

!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE TSWAP                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE tswap(nx,ny,nz, a) 12
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Update the time levels of a time-dependent array.
!
!  a(*,*,*,tpast   ) = a(*,*,*,tpresent)
!  a(*,*,*,tpresent) = a(*,*,*,tfuture )
!  a(*,*,*,tfuture ) is not changed.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/91.
!
!  MODIFICATION HISTORY:
!
!  5/20/92 (K. Droegemeier and M. Xue)
!  Added full documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    a        A 3-D array at three time levels that will be updated
!             in time.
!
!  OUPTUT:
!
!    a        Its new values at time tpast and tfuture.
!
!-----------------------------------------------------------------------

!
!-----------------------------------------------------------------------
!
!  Variable Declarations:
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INCLUDE 'timelvls.inc'

  INTEGER :: nx, ny, nz   ! Number of grid points in 3 directions

  REAL :: a(nx,ny,nz,nt)  ! Array whose values will be updated in time
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
!
!-----------------------------------------------------------------------
!
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        a(i,j,k,tpast   )=a(i,j,k,tpresent)
        a(i,j,k,tpresent)=a(i,j,k,tfuture )
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE tswap