      PROGRAM ARPS30   
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######      Advanced Regional Prediction System (ARPS)      ######
c     ######                   Version 3.0                        ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
c
c     For further information, contact:
c
c     Center for Analysis and Prediction of Storms 
c     Sarkeys Energy Center, Rm 1110
c     100 East Boyd
c     Norman, OK  73019  USA
c     Phone:  (405) 325-3041
c     FAX  :  (405) 325-7614
c     E-mail: kdroege@geohub.gcn.uoknor.edu  
c         or: mxue@geohub.gcn.uoknor.edu  
c
c#######################################################################
c
c     PROGRAM ARPS 
c
c     This program is the driver routine for the entire ARPS model
c     system.  The ARPS is a fully three-dimensional, compressible,
c     nonhydrostatic model based on the Arakawa C grid.  The model
c     solves three momentum equations, a thermodynamic energy
c     equation, a pressure equation, six moisture equations (water
c     vapor, cloud water, rainwater, cloud ice, hail, and snow),
c     and also provides for the statistical representation of 
c     unresolvable processes via a turbulence parameterization
c     (user option).  The model is written in an terrain-following 
c     coordinate system, and has provision for several
c     boundary condition options.
c
c#######################################################################
c
c     MODIFICATION HISTORY:  
c
c     10/01/91  (Ming Xue)
c
c       Model dynamic framework.
c
c     11/01/91 (Ming Xue) 
c
c       Smagorinsky's subgrid scale turbulence and options.
c
c     11/20/91 (Ming Xue) 
c
c       Add the option to use three different formulations of the 
c       advection terms. 
c
c     2/28/92 (Ming Xue)
c
c       Major restructuring of the model.
c       work arrays urho, vrho, wrho, ptdif, qvdif, qcdif, 
c       qrdif removed from the main program.
c
c       The dif terms for scalars are calculated where they are used.
c       Diffusion/mixing subroutines are made to take input arrays 
c       at a single time level.
c
c       The integrations of water quantity equations are completely 
c       separated.
c
c     2/29/92 (Ming Xue)
c
c
c       Arrays for water quantities are defined at tow time levels.
c
c       Warm microphysics parameterization included. 
c       Model results for dry case are not affected by this addition.
c
c     3/11/1992 (Ming Xue)
c
c       Main program restructured.
c       Subroutines initout, output and chkstab added. 
c       History dump routines rewritten.
c       Model results are not affect by the structure change.
c
c     3/12/1992 (Ming Xue)
c
c       Arrays uflux,vflux,wflux,udif,vdif,wdif,pflux,wbuoy,psrc
c       removed from the main program.
c       Arrays uforce,vforce,wforce and pforce added.
c       All acoustically inactive forcing terms are now accumulated
c       together and saved in a single forcing array for each variable.
c       Results are not affect by the structure change.
c
c     3/13/1992 (Ming Xue)
c
c       Divergence damping added and results verified.
c
c     3/17/1992 (Ming Xue)
c
c       Base state subroutine basest enhanced to read in arbitry 
c       base state sounding from a file
c
c     4/1/1992  (Ming Xue)
c
c       Model restart data dump and restart run data read added.
c
c     4/7/1992 (Ming Xue)
c
c       A simple surface flux scheme added. These fluxes are calculated
c       in subroutine sfcflx and passed into turbulent mixing
c       subroutines.
c
c     4/17/1992 Version 2.2a (Ming Xue)
c
c       Array name changes made. 
c       uflux, vflux, wflux, pflux, ptflux changed to 
c       uadv, vadv, wadv, padv and ptadv. 
c
c       udif, vdif, wdif, ptdif, qdif, qvdif changed to 
c       umix, vmix, wmix, ptmix, qmix, qvmix.
c
c       variable names ckh2, ckv2, ckh4, ckv4 , kdif, cstkm
c       changed to cfcmh2, cfcmv2, cfcmh4, cfcmv4, tmixopt, tmixcst.
c
c       subroutine names basest, initia changed to inibase, initdvr.
c
c       More documentations added.
c
c       Results verified.
c
c     4/19/1992. (Ming Xue)
c
c       Ice variables qi, qs and qh added.
c       Results verified.
c
c     4/19/1992. (Ming Xue)
c
c       Boundary conditions other than periodic implemented for the 
c       deformation tensor in subroutine bcdefm.
c 
c       Some of the loop bounds in bcdefm corrected, which affects only 
c       km above the top and below the bottom boundary. Previous
c       testing results are still valid since km values outside the 
c       top and bottom boundaries are not required.
c       Results verified.
c
c     4/19/1992 (Ming Xue)
c
c       The definition of Vaisala frequency above top and below bottom 
c       boundary modefied in subroutine cftmix. 
c       Results verified.
c
c     4/20/1992 Version 2.3 (Ming Xue)
c
c       Coordinate arrays x, y, z and zp added. 
c       Transformation Jacobians J1, J2 and J3 are calculated in inigrd. 
c       To be verified.
c
c     4/20/1992 (Ming Xue)
c
c       Coriolis force terms added.
c       To be verified.
c
c     4/21/1992 (Ming Xue)
c
c       Rayleigh damping included.
c       To be verified.
c
c     4/22/1992 (Ming Xue)
c
c       Further facelift per discussion with Kelvin
c       Subroutine stepdt changed to cordintg, which is then split into
c       calls to tinteg and acoust. 
c       Initialization routines re-organized.
c
c     4/26/92 (K. Droegemeier)
c       Added full documentation
c
c     5/03/1992 (Ming Xue)
c       Many changes in varables names, more standard documentations.
c
c     5/21/1992 (Ming Xue)
c       subroutine advq added. several subroutine names changed.
c
c     5/23/1992 (Ming Xue)
c       Three instead of two time levels are defined for water variables,
c       Related structure changes made, the structure changes
c       cause the results to be slightly different from previous
c       version due to the order that time filter is applied to 
c       water variables. When time filter is switched off, 
c       both versions produced the same results.
c
c     5/25/1992 (Ming Xue)
c       Subroutine INIPARA rewritten to provide a more friendly
c       user interface for setting control parameters.
c
c       globcst.inc modified to define more control parameters.
c
c     6/06/1992 (Ming Xue)
c
c       Regrouping of certain subroutines. Order of subroutine
c       arguments changed in INIT, INITDVR, ADVP and SOLVPT.
c
c       SOLVUVW, SOLVP, BCU, BCV, BCW and BCP modified to take in 
c       only fields at tfuture.
c
c     6/11/1992 (MX)
c       uforce,vforce,wforce,pforce and usflx,vsflx,ptsflx and qvsflx
c       removed from the main program and redefined as temporary arrays.
c
c     6/20/1992 (MX)
c       Corrections made in raydmp3d.f in two places.
c 
c     6/22/1992 (MX)
c       Documenation cleaning up after proofreading by users.
c       Results not affected.
c
c     7/8/92 (MX)
c       A bug fix in subroutine inibase in do loop 30.
c
c     7/9/92 (MX)
c       Included in the initial bubble size and location parameters into
c       the read in parameter list. 
c
c     7/10/92 (MX)
c       Modefied the history data dump format, to include more 
c       data information.       
c
c     7/21/92 (MX)
c       Subroutine name strlen changed to strlnth to avoid name conflict 
c       with the same name in C library. 
c
c     7/21/92 (MX)
c       Energy statistics calculations added. Affected files:
c       energy.f, init3d.f, out3d.f and globcst.inc 
c       Input file is also changed.
c
c     7/27/92 (MX)
c       Added the option for initial random pt perturbation, 
c       model run mode control, option for 2 LBC formulations and that 
c       additional history dump format. The number of input parameters 
c       is increased
c
c     7/30/92 (M. Xue and K. Brewster)
c       Implemented new history data dump format (format version 3.0). 
c       Added packed binary history data dump option. 
c
c     8/03/92 (M. Xue)
c       The grid scale in the formula of divergecne damping coefficient
c       is changed from (dx*dy*dz)**(2/3.0) to min(dx,dy,dz)**2.
c       This will affect the results compared with previous runs.
c
c     8/04/92 (M. Xue)
c       History data dumping for Savi3D implemented. A lot of changes 
c       to the input file.
c
c       Time tendencies at the lateral boundaries included in subroutines
c       rstin and rstout.
c
c     8/23/92 (M. Xue)
c       History dump calls streamlined. A single call to dtadump is now
c       performed.
c
c     9/16/92 (M. Xue)
c       Full Coriolis formualtion implemented.
c       Surface flux formulations changed.
c
c#######################################################################
c
c     GRID STRUCTURE AND STENCIL NUMBERING CONVENTION
c
c     Each of the three velocity components is located at a physical
c     boundary in the model.  The numbering convention for the
c     staggered grid is shown below for each of the coordinate
c     directions, where S denotes a scalar variable.
c
c     Arakawa C-grid is used by this model.
c
c
c     X-Direction (East.......West):
c
c           West boundary                       East boundary
c                 |                                   |
c                 <--------   Physical Domain   ------>
c     +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
c     U     S     U     S     U     S     U     S     U     S     U
c
c i=  1     1     2     2     . . .     NX-2  NX-2  NX-1  NX-1   NX
c
c
c     Y-Direction (North......South):
c
c            South boundary                     North boundary
c                 |                                   |
c                 <-------    Physical Domain  ------->
c     +-----+-----+-----+-----+-----+-----+-----+-----+-----+-----+
c     V     S     V     S     V     S     V     S     V     S     V
c
c j=  1     1     2     2     . . .     NY-2  NY-2  NY-1  NY-1   NY
c
c
c     Z-Direction:
c
c                         NZ  + W   
c                             |
c                             |
c                             |
c                        NZ-1 + S
c                             |
c                             |
c                             |
c                        NZ-1 + W --  Physical Top boundary.
c                             |
c                             .
c                             .
c                             .
c
c                             |
c                          2  + S
c                             |
c                             |
c                             |
c                          2  + W --  Model Ground Level
c                             |
c                             |
c                             |
c                          1  + S
c                             |
c                             |
c                             |
c                    k=    1  + W
c                             
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none        ! Force explicit declarations
c
      include 'dims.inc'

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

c
c#######################################################################
c
c     Permanent arrays defined in the model. 
c
c     These arrays must be saved (can not be overwritten) between 
c     time steps.
c
c#######################################################################
c

c
c#######################################################################
c
c     Time dependent variables:
c
c#######################################################################
c
      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
                                ! from that of base state atmosphere (K)
      real pprt  (nx,ny,nz,nt)  ! Perturbation pressure from that
                                ! of base state atmosphere (Pascal)
c
      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 km    (nx,ny,nz)     ! The turbulent mixing coefficient for 
                                ! momentum. ( m**2/s )      

c
c#######################################################################
c
c     Time tendencies of certain time dependent variables at the lateral 
c     boundaries, which are used to set the lateral boundary conditions
c     for options 4 and 5.
c
c#######################################################################
c
      real udteb (ny,nz)        ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)        ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)        ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)        ! T-tendency of u at s-boundary (m/s**2)

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

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

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

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

      real qvdteb(ny,nz)        ! T-tendency of qv at e-boundary (kg/(kg*s))
      real qvdtwb(ny,nz)        ! T-tendency of qv at w-boundary (kg/(kg*s))
      real qvdtnb(nx,nz)        ! T-tendency of qv at n-boundary (kg/(kg*s))
      real qvdtsb(nx,nz)        ! T-tendency of qv at s-boundary (kg/(kg*s))

      real qcdteb(ny,nz)        ! T-tendency of qc at e-boundary (kg/(kg*s))
      real qcdtwb(ny,nz)        ! T-tendency of qc at w-boundary (kg/(kg*s))
      real qcdtnb(nx,nz)        ! T-tendency of qc at n-boundary (kg/(kg*s))
      real qcdtsb(nx,nz)        ! T-tendency of qc at s-boundary (kg/(kg*s))

      real qrdteb(ny,nz)        ! T-tendency of qr at e-boundary (kg/(kg*s))
      real qrdtwb(ny,nz)        ! T-tendency of qr at w-boundary (kg/(kg*s))
      real qrdtnb(nx,nz)        ! T-tendency of qr at n-boundary (kg/(kg*s))
      real qrdtsb(nx,nz)        ! T-tendency of qr at s-boundary (kg/(kg*s))

      real qidteb(ny,nz)        ! T-tendency of qi at e-boundary (kg/(kg*s))
      real qidtwb(ny,nz)        ! T-tendency of qi at w-boundary (kg/(kg*s))
      real qidtnb(nx,nz)        ! T-tendency of qi at n-boundary (kg/(kg*s))
      real qidtsb(nx,nz)        ! T-tendency of qi at s-boundary (kg/(kg*s))

      real qsdteb(ny,nz)        ! T-tendency of qs at e-boundary (kg/(kg*s))
      real qsdtwb(ny,nz)        ! T-tendency of qs at w-boundary (kg/(kg*s))
      real qsdtnb(nx,nz)        ! T-tendency of qs at n-boundary (kg/(kg*s))
      real qsdtsb(nx,nz)        ! T-tendency of qs at s-boundary (kg/(kg*s))

      real qhdteb(ny,nz)        ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)        ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)        ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)        ! T-tendency of qh at s-boundary (kg/(kg*s))
c
c#######################################################################
c
c     Base state variables:
c
c#######################################################################
c
      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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity 
                                ! (kg/kg)
c
c#######################################################################
c
c     Arrays related to model grid definition: 
c
c#######################################################################
c
      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 hterain(nx,ny)       ! The height of the terrain.

      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ). 
c
c#######################################################################
c
c     Pure work arrays that do not carry physical meaning in the code
c
c#######################################################################
c
      real temxy1(nx,ny)        ! 2-D temporary array
      real temxy2(nx,ny)        ! 2-D temporary array
      real temxy3(nx,ny)        ! 2-D temporary array
      real temxy4(nx,ny)        ! 2-D temporary array
c
      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 svitem(3*nxyz)       ! A work space required by Savi3D 
                                ! library routines.
c
c#######################################################################
c
c     Make equivalence of work arrays so that savitem takes the same
c     memory space as tem11, tem12 and tem13 combined.  
c
c#######################################################################
c 
      equivalence (svitem(0*nxyz+1),tem11) 
      equivalence (svitem(1*nxyz+1),tem12)
      equivalence (svitem(2*nxyz+1),tem13)
c
      integer frstep            ! Flag for the initial time step


      real umax,umin,vmax,vmin,wmax,wmin  ! added by G.C.
      integer i,j,k,it     ! added by G. C.

c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'      ! Unchanging physical constants
      include 'globcst.inc'     ! Global constants that control model
      include 'bndry.inc'       ! Boundary condition control parameters
c
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external init
      external initout
      external initbdt 
      external cordintg
      external output 
      external chkstab
      external grafclose
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Initialize all arrays and variables, and set initial
c     condition for thermal disturbance if desired. 
c
c#######################################################################
c

      CALL INIT(nx,ny,nz,
     :     u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :     udteb, udtwb, udtnb, udtsb, vdteb, vdtwb, vdtnb, vdtsb,
     :     wdteb, wdtwb, wdtnb, wdtsb, ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :     pdteb, pdtwb, pdtnb, pdtsb,
     :     qvdteb,qvdtwb,qvdtnb,qvdtsb,qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :     qrdteb,qrdtwb,qrdtnb,qrdtsb,qidteb,qidtwb,qidtnb,qidtsb,
     :     qsdteb,qsdtwb,qsdtnb,qsdtsb,qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :     ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :     x,y,z,zp,hterain, j1,j2,j3, km, 
     :     temxy1,temxy2,temxy3,temxy4,
     :     tem1,tem2,tem3,tem4,tem5,tem6,tem7,
     :     tem8,tem9,tem10,tem11,tem12,tem13)

c
c#######################################################################
c
c     The current model time (initial time): 
c
c#######################################################################
c
      curtim = tstart

      write(6,'(1x,a,f13.3,a)') 
     :     'The initial model time is at ', curtim,' seconds.'
c
c#######################################################################
c
c     Output the initial fields
c
c#######################################################################
c
      CALL INITOUT(nx,ny,nz,
     :             u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :             ubar,vbar,ptbar,pbar,rhobar,qvbar,km,
     :             x,y,z,zp,hterain,j1,j2,j3, 
     :             tem1,tem2,tem3,tem4,tem5, svitem)
c
c#######################################################################
c
c     Total number of time steps to be taken: 
c
c#######################################################################
c

      nsteps = int( (tstop-tstart+0.1*dtbig)/dtbig )
      
c
c#######################################################################
c
c     Time integration loop begins  ----------------------------------->
c
c#######################################################################
c


      DO 900 nstep = 1, nsteps
      
      IF( (nstep.eq.1) .and. (restrt.ne.1) ) THEN 

        frstep=1           ! Indicate that this is the initial step of
                           ! model integration. For this step forward
                           ! time integration scheme will be used.

      ELSE                 ! For non-first step or restart run

        frstep=0           ! When frstep=0, leapfrog scheme is used.

      ENDIF

c
c#######################################################################
c
c     Perform one time step of integration for all equations
c     On exit of this routine, all time dependent fields are 
c     advanced in time by one time step.
c
c#######################################################################
c
      CALL CORDINTG(frstep, nx,ny,nz,
     :     u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :     udteb, udtwb, udtnb, udtsb, vdteb, vdtwb, vdtnb, vdtsb,
     :     wdteb, wdtwb, wdtnb, wdtsb, ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :     pdteb, pdtwb, pdtnb, pdtsb,
     :     qvdteb,qvdtwb,qvdtnb,qvdtsb,qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :     qrdteb,qrdtwb,qrdtnb,qrdtsb,qidteb,qidtwb,qidtnb,qidtsb,
     :     qsdteb,qsdtwb,qsdtnb,qsdtsb,qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :     ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :     x,y,z,zp, j1,j2,j3, km,
     :     temxy1,temxy2,temxy3,temxy4,
     :     tem1,tem2,tem3,tem4,tem5,tem6,tem7,
     :     tem8,tem9,tem10,tem11,tem12,tem13)
c
c#######################################################################
c
c     Update physical time and generate output files
c
c#######################################################################
c
      curtim = curtim + dtbig  
c
c#######################################################################
c
c     Print timestep marker
c
c#######################################################################
c
      write(6,'(a,i5,a,f8.2,a)') 
     :  '  End of integeration time step ', nstep,
     :  ', Model time=',curtim,' (S)'

      CALL OUTPUT(nx,ny,nz,
     :     u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :     udteb, udtwb, udtnb, udtsb, vdteb, vdtwb, vdtnb, vdtsb,
     :     wdteb, wdtwb, wdtnb, wdtsb, ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :     pdteb, pdtwb, pdtnb, pdtsb,
     :     qvdteb,qvdtwb,qvdtnb,qvdtsb,qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :     qrdteb,qrdtwb,qrdtnb,qrdtsb,qidteb,qidtwb,qidtnb,qidtsb,
     :     qsdteb,qsdtwb,qsdtnb,qsdtsb,qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :     ubar,vbar,ptbar,pbar,rhobar,qvbar,km,
     :     x,y,z,zp,hterain, j1,j2,j3, 
     :     tem1,tem2,tem3,tem4,tem5, svitem)

c
c#######################################################################
c
c     Check the stability of the integration.  If the model solution
c     appears to be exceeding the linear stability condition, then
c     perform a data dump for post mortum. 
c
c#######################################################################
c
      CALL CHKSTAB(nx,ny,nz,
     :             u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :             ubar,vbar,ptbar,pbar,rhobar,qvbar,km,
     :             x,y,z,zp,hterain, j1,j2,j3, tem1,tem2,tem3, svitem)
c
c#######################################################################
c
c     Time integration loop ends   <------------------------------------
c
c#######################################################################
c

c *********************************************** added by Gang Cheng
        umax = -10000.
        umin = 1000.0
        vmax = umax
        vmin = umin
        wmax = umax
        wmin = umin

        do 100 it=1,3
        do 100 i = 1,nx
        do 100 j = 1,ny
        do 100 k = 1,nz
          umax = amax1(umax,u(i,j,k,it))
          umin = amin1(umin,u(i,j,k,it))
          vmax = amax1(vmax,v(i,j,k,it))
          vmin = amin1(vmin,v(i,j,k,it))
          wmax = amax1(wmax,w(i,j,k,it))
          wmin = amin1(wmin,w(i,j,k,it))
100     continue

        print *,'umax=',umax,' umin',umin
        print *,'vmax=',vmax,' vmin',vmin
        print *,'wmax=',wmax,' wmin',wmin
c *********************************************** added by Gang Cheng


900   CONTINUE                ! loop over timesteps

c
c#######################################################################
c
c     End of entire model time integration. The program stops.
c
c     If ARPS is dumping Savi3D GRAF format history files,
c     then grafclose is called to close the GRAF file.
c
c#######################################################################
c

      IF (hdmpfmt .eq. 5) THEN

        CALL grafclose(grafh)

      ENDIF
    
      STOP
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CORDINTG                   ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE CORDINTG(frstep, nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           udteb, udtwb, udtnb, udtsb,
     :           vdteb, vdtwb, vdtnb, vdtsb,
     :           wdteb, wdtwb, wdtnb, wdtsb,
     :           ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :           pdteb ,pdtwb ,pdtnb ,pdtsb,
     :           qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :           qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :           qrdteb,qrdtwb,qrdtnb,qrdtsb,
     :           qidteb,qidtwb,qidtnb,qidtsb,
     :           qsdteb,qsdtwb,qsdtnb,qsdtsb,
     :           qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :           x,y,z,zp, j1,j2,j3, km,
     :           usflx,vsflx,ptsflx,qvsflx,
     :           tem1,tem2,tem3,tem4,tem5,tem6,tem7,
     :           tem8,tem9,tem10,tem11,tem12,tem13)
c
c#######################################################################
c
c     PURPOSE:          
c
c     Coordinates the forward integration of all field variables by one time
c     step.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/01/1991.
c
c     MODIFICATION HISTORY:
c
c     4/26/92 (K. Droegemeier)  
c     Added full documentation.
c
c#######################################################################
c
c     INPUT:
c
c       frstep   Flag to determine if this is the initial time step. 
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at times tpast and tpresent (m/s)
c       v        y component of velocity at times tpast and tpresent (m/s)
c       w        Vertical component of Cartesian velocity at times 
c                tpast and tpresent (m/s)
c       wcont    Vertical component of contravariant velocity in
c                computational coordinates (m/s)
c       ptprt    Perturbation potential temperature at times tpast and
c                tpresent (K)
c       pprt     Perturbation pressure at times tpast and tpresent (Pascal)
c       qv       Water vapor specific humidity at times tpast and tpresent (kg/kg)
c       qc       Cloud water mixing ratio at times tpast and tpresent (kg/kg)
c       qr       Rainwater mixing ratio at times tpast and tpresent (kg/kg)
c       qi       Cloud ice mixing ratio at times tpast and tpresent (kg/kg)
c       qs       Snow mixing ratio at times tpast and tpresent (kg/kg)
c       qh       Hail mixing ratio at times tpast and tpresent (kg/kg)
c
c       udteb    Time tendency of u field at east boundary (m/s**2)
c       udtwb    Time tendency of u field at west boundary (m/s**2)
c       udtnb    Time tendency of u field at north boundary (m/s**2)
c       udtsb    Time tendency of u field at south boundary (m/s**2)
c
c       vdteb    Time tendency of v field at east boundary (m/s**2)
c       vdtwb    Time tendency of v field at west boundary (m/s**2)
c       vdtnb    Time tendency of v field at north boundary (m/s**2)
c       vdtsb    Time tendency of v field at south boundary (m/s**2)
c
c       wdteb    Time tendency of w field at east boundary (m/s**2)
c       wdtwb    Time tendency of w field at west boundary (m/s**2)
c       wdtnb    Time tendency of w field at north boundary (m/s**2)
c       wdtsb    Time tendency of w field at south boundary (m/s**2)
c
c       ptdteb   Time tendency of ptprt field at east boundary (K/s)
c       ptdtwb   Time tendency of ptprt field at west boundary (K/s)
c       ptdtnb   Time tendency of ptprt field at north boundary (K/s)
c       ptdtsb   Time tendency of ptprt field at south boundary (K/s)
c
c       pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
c       pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
c       pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
c       pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
c
c       qvdteb   Time tendency of qv field at east boundary (kg/(kg*s))
c       qvdtwb   Time tendency of qv field at west boundary (kg/(kg*s))
c       qvdtnb   Time tendency of qv field at north boundary (kg/(kg*s))
c       qvdtsb   Time tendency of qv field at south boundary (kg/(kg*s))
c
c       qcdteb   Time tendency of qc field at east boundary (kg/(kg*s))
c       qcdtwb   Time tendency of qc field at west boundary (kg/(kg*s))
c       qcdtnb   Time tendency of qc field at north boundary (kg/(kg*s))
c       qcdtsb   Time tendency of qc field at south boundary (kg/(kg*s))
c
c       qrdteb   Time tendency of qr field at east boundary (kg/(kg*s))
c       qrdtwb   Time tendency of qr field at west boundary (kg/(kg*s))
c       qrdtnb   Time tendency of qr field at north boundary (kg/(kg*s))
c       qrdtsb   Time tendency of qr field at south boundary (kg/(kg*s))
c
c       qidteb   Time tendency of qi field at east boundary (kg/(kg*s))
c       qidtwb   Time tendency of qi field at west boundary (kg/(kg*s))
c       qidtnb   Time tendency of qi field at north boundary (kg/(kg*s))
c       qidtsb   Time tendency of qi field at south boundary (kg/(kg*s))
c
c       qsdteb   Time tendency of qs field at east boundary (kg/(kg*s))
c       qsdtwb   Time tendency of qs field at west boundary (kg/(kg*s))
c       qsdtnb   Time tendency of qs field at north boundary (kg/(kg*s))
c       qsdtsb   Time tendency of qs field at south boundary (kg/(kg*s))
c
c       qhdteb   Time tendency of qh field at east boundary (kg/(kg*s))
c       qhdtwb   Time tendency of qh field at west boundary (kg/(kg*s))
c       qhdtnb   Time tendency of qh field at north boundary (kg/(kg*s))
c       qhdtsb   Time tendency of qh field at south boundary (kg/(kg*s))
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       hterain  The height of terrain (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       u        x component of velocity at times tpast and tpresent but 
c                advanced in time by dtbig (m/s).
c       v        y component of velocity at times tpast and tpresent but
c                advanced in time by dtbig (m/s).
c       w        Vertical component of Cartesian velocity at times 
c                tpast and tpresent but advanced in time by dtbig (m/s).
c       ptprt    Perturbation potential temperature at times tpast and
c                tpresent but advanced in time by dtbig (K).
c       pprt     Perturbation pressure at times tpast and tpresent but 
c                advanced in time by dtbig (Pascal).
c       qv       Water vapor specific humidity at times tpast and 
c                tpresent but advanced in time by dtbig (kg/kg).
c       qc       Cloud water mixing ratio at times tpast and tpresent 
c                but advanced in time by dtbig (kg/kg).
c       qr       Rainwater mixing ratio at times tpast and tpresent but 
c                advanced in time by dtbig (kg/kg).
c       qi       Cloud ice mixing ratio at times tpast and tpresent but 
c                advanced in time by dtbig (kg/kg).
c       qs       Snow mixing ratio at times tpast and tpresent but 
c                advanced in time by dtbig (kg/kg).
c       qh       Hail mixing ratio at times tpast and tpresent but 
c                advanced in time by dtbig (kg/kg).
c
c       udteb    Time tendency of u field at east boundary (m/s**2)
c       udtwb    Time tendency of u field at west boundary (m/s**2)
c       udtnb    Time tendency of u field at north boundary (m/s**2)
c       udtsb    Time tendency of u field at south boundary (m/s**2)
c
c       vdteb    Time tendency of v field at east boundary (m/s**2)
c       vdtwb    Time tendency of v field at west boundary (m/s**2)
c       vdtnb    Time tendency of v field at north boundary (m/s**2)
c       vdtsb    Time tendency of v field at south boundary (m/s**2)
c
c       wdteb    Time tendency of w field at east boundary (m/s**2)
c       wdtwb    Time tendency of w field at west boundary (m/s**2)
c       wdtnb    Time tendency of w field at north boundary (m/s**2)
c       wdtsb    Time tendency of w field at south boundary (m/s**2)
c
c       ptdteb   Time tendency of ptprt field at east boundary (K/s)
c       ptdtwb   Time tendency of ptprt field at west boundary (K/s)
c       ptdtnb   Time tendency of ptprt field at north boundary (K/s)
c       ptdtsb   Time tendency of ptprt field at south boundary (K/s)
c
c       pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
c       pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
c       pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
c       pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
c
c       qvdteb   Time tendency of qv field at east boundary (kg/(kg*s))
c       qvdtwb   Time tendency of qv field at west boundary (kg/(kg*s))
c       qvdtnb   Time tendency of qv field at north boundary (kg/(kg*s))
c       qvdtsb   Time tendency of qv field at south boundary (kg/(kg*s))
c
c       qcdteb   Time tendency of qc field at east boundary (kg/(kg*s))
c       qcdtwb   Time tendency of qc field at west boundary (kg/(kg*s))
c       qcdtnb   Time tendency of qc field at north boundary (kg/(kg*s))
c       qcdtsb   Time tendency of qc field at south boundary (kg/(kg*s))
c
c       qrdteb   Time tendency of qr field at east boundary (kg/(kg*s))
c       qrdtwb   Time tendency of qr field at west boundary (kg/(kg*s))
c       qrdtnb   Time tendency of qr field at north boundary (kg/(kg*s))
c       qrdtsb   Time tendency of qr field at south boundary (kg/(kg*s))
c
c       qidteb   Time tendency of qi field at east boundary (kg/(kg*s))
c       qidtwb   Time tendency of qi field at west boundary (kg/(kg*s))
c       qidtnb   Time tendency of qi field at north boundary (kg/(kg*s))
c       qidtsb   Time tendency of qi field at south boundary (kg/(kg*s))
c
c       qsdteb   Time tendency of qs field at east boundary (kg/(kg*s))
c       qsdtwb   Time tendency of qs field at west boundary (kg/(kg*s))
c       qsdtnb   Time tendency of qs field at north boundary (kg/(kg*s))
c       qsdtsb   Time tendency of qs field at south boundary (kg/(kg*s))
c
c       qhdteb   Time tendency of qh field at east boundary (kg/(kg*s))
c       qhdtwb   Time tendency of qh field at west boundary (kg/(kg*s))
c       qhdtnb   Time tendency of qh field at north boundary (kg/(kg*s))
c       qhdtsb   Time tendency of qh field at south boundary (kg/(kg*s))
c
c       km       Turbulent mixing coefficient for momentum.
c
c     WORK ARRAYS:
c
c       usflx    Surface flux of u-momentum (kg/(m*s**2))
c       vsflx    Surface flux of v-momentum (kg/(m*s**2))
c       ptsflx   Surface heat flux (K*kg/(m**2 * s ))
c       qvsflx   Surface moisture flux of (kg/(m**2 * s))
c
c       tem1     Temporary work array.
c       tem2     Temporary work array.
c       tem3     Temporary work array.
c       tem4     Temporary work array.
c       tem5     Temporary work array.
c       tem6     Temporary work array.
c       tem7     Temporary work array.
c       tem8     Temporary work array.
c       tem9     Temporary work array.
c       tem10    Temporary work array.
c       tem11    Temporary work array.
c       tem12    Temporary work array.
c       tem13    Temporary work array.
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work 
c       arrays before you alter the code.)
c
c#######################################################################
c

c
c#######################################################################
c
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none        ! Force explicit declarations
c
      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

      integer nx, ny, nz        ! 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 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 udteb (ny,nz)        ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)        ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)        ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)        ! T-tendency of u at s-boundary (m/s**2)

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

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

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

      real pdteb (ny,nz)        ! T-tendency of pprt at e-boundary (PASCAL/s)
      real pdtwb (ny,nz)        ! T-tendency of pprt at w-boundary (PASCAL/s)
      real pdtnb (nx,nz)        ! T-tendency of pprt at n-boundary (PASCAL/s)
      real pdtsb (nx,nz)        ! T-tendency of pprt at s-boundary (PASCAL/s)

      real qvdteb(ny,nz)        ! T-tendency of qv at e-boundary (kg/(kg*s))
      real qvdtwb(ny,nz)        ! T-tendency of qv at w-boundary (kg/(kg*s))
      real qvdtnb(nx,nz)        ! T-tendency of qv at n-boundary (kg/(kg*s))
      real qvdtsb(nx,nz)        ! T-tendency of qv at s-boundary (kg/(kg*s))

      real qcdteb(ny,nz)        ! T-tendency of qc at e-boundary (kg/(kg*s))
      real qcdtwb(ny,nz)        ! T-tendency of qc at w-boundary (kg/(kg*s))
      real qcdtnb(nx,nz)        ! T-tendency of qc at n-boundary (kg/(kg*s))
      real qcdtsb(nx,nz)        ! T-tendency of qc at s-boundary (kg/(kg*s))

      real qrdteb(ny,nz)        ! T-tendency of qr at e-boundary (kg/(kg*s))
      real qrdtwb(ny,nz)        ! T-tendency of qr at w-boundary (kg/(kg*s))
      real qrdtnb(nx,nz)        ! T-tendency of qr at n-boundary (kg/(kg*s))
      real qrdtsb(nx,nz)        ! T-tendency of qr at s-boundary (kg/(kg*s))

      real qidteb(ny,nz)        ! T-tendency of qi at e-boundary (kg/(kg*s))
      real qidtwb(ny,nz)        ! T-tendency of qi at w-boundary (kg/(kg*s))
      real qidtnb(nx,nz)        ! T-tendency of qi at n-boundary (kg/(kg*s))
      real qidtsb(nx,nz)        ! T-tendency of qi at s-boundary (kg/(kg*s))

      real qsdteb(ny,nz)        ! T-tendency of qs at e-boundary (kg/(kg*s))
      real qsdtwb(ny,nz)        ! T-tendency of qs at w-boundary (kg/(kg*s))
      real qsdtnb(nx,nz)        ! T-tendency of qs at n-boundary (kg/(kg*s))
      real qsdtsb(nx,nz)        ! T-tendency of qs at s-boundary (kg/(kg*s))

      real qhdteb(ny,nz)        ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)        ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)        ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)        ! T-tendency of qh at s-boundary (kg/(kg*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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                                ! (kg/kg)

      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

      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian -d(zp)/d(x)
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian -d(zp)/d(y)
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian  d(zp)/d(z)
c
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
c
      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 flux of heat (K*kg/(m*s**2))
      real qvsflx(nx,ny)        ! surface flux of moisture (kg/(m**2*s))
c
      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.
c
c#######################################################################
c
c     Miscellaneous local variables:
c
c#######################################################################
c 
      integer nocalls
      data nocalls / 0 /
      save nocalls
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'      ! Unchanging physical constants
      include 'globcst.inc'     ! Global constants that control model
      include 'bndry.inc'       ! Boundary condition control parameters
      include 'debug.inc'       ! Debug information print control
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external sfcflx
      external tinteg 
      external microph
      external tfilt
      external tflip
      external radbdt 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF( subtrac .eq. 1) THEN

        nocalls = nocalls + 1

        write(6,'(1x,a,i8,a)') 
     :       'Subroutine CORDINTG is called ',nocalls,' times.'

      ENDIF
c
c#######################################################################
c
c     Calculate the surface fluxes of momentum, heat, and water
c     vapor for use in the turbulent mixing subroutine
c
c#######################################################################
c
      CALL SFCFLX(nx,ny,nz, 
     :            u(1,1,1,1),v(1,1,1,1),ptprt(1,1,1,1),qv(1,1,1,1),
     :            rhobar,ptbar,qvbar,
     :            usflx,vsflx,ptsflx,qvsflx)
c
c#######################################################################
c
c     Integrate the conservation equations of momentum, heat, pressure
c     and water variables forward by one timestep
c
c#######################################################################
c
      CALL TINTEG(frstep, nx,ny,nz,
     :            u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :            udteb, udtwb, udtnb, udtsb,
     :            vdteb, vdtwb, vdtnb, vdtsb,
     :            wdteb, wdtwb, wdtnb, wdtsb,
     :            ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :            pdteb, pdtwb, pdtnb, pdtsb,
     :            qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :            qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :            qrdteb,qrdtwb,qrdtnb,qrdtsb,
     :            qidteb,qidtwb,qidtnb,qidtsb,
     :            qsdteb,qsdtwb,qsdtnb,qsdtsb,
     :            qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :            ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :            x,y,z,zp, j1,j2,j3, 
     :            usflx,vsflx,ptsflx,qvsflx, km,
     :            tem1,tem2,tem3,tem4,tem5,tem6,tem7,
     :            tem8,tem9,tem10,tem11,tem12,tem13)
c
c#######################################################################
c
c     Call the microphysics package which adjusts thermodynamic
c     and water variables.
c
c     On exit of this subroutine supersaturation is removed. 
c
c     If the model is being run dry, skip these calcuations 
c     to save computer time.
c
c#######################################################################
c

      IF( mphyopt .eq. 0 ) THEN        ! No microphysics
c
c       No microphysics process! Do nothing.
c 
      ELSEIF( mphyopt .eq. 1 ) THEN    ! Option 1
                                

        CALL MICROPH(nx,ny,nz, ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :               rhobar,pbar,ptbar,
     :               tem1,tem2,tem3,tem4,tem5)

      ELSEIF( mphyopt .eq. 2 ) THEN    ! Option 2

        write(6,'(1x,a,/1x,a)')
     :       'Option 2 of microphysics parameterization scheme not',
     :       'implemented yet. Reset the option. Job stopped here.'
        STOP

      ELSE                             ! Invalid option
                                
        write(6,'(1x,a,/1x,a,i3)')
     :       'Invalid option for microphysics parameterization.',
     :       'Job stopped here. MPHYOPT was ',mphyopt
        STOP
                                
      ENDIF
c
c#######################################################################
c        
c     Apply asselin time filter to all variables
c 
c#######################################################################
c

      IF( flteps.ne.0 ) THEN 

        CALL TFILT(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh)

      ENDIF
c
c#######################################################################
c        
c     Set the time tendencies of variables at the lateral boundaries. 
c     For option 4, radiation boundary condition is used to calculate 
c     the tendencies, for option 5, the tendencies are obtained using 
c     external information.
c
c     If none of the four lateral boundaries has radiation or user 
c     specified condition, these time tendencies need not to be set.
c 
c#######################################################################
c
      IF( ebc.eq.4 .or. wbc.eq.4 .or. nbc.eq.4 .or. sbc.eq.4 ) THEN

        CALL RADBDT(nx,ny,nz,dtbig,
     :              u, v, w, ptprt, pprt, qv, qc, qr, qi, qs, qh,
     :              udteb , udtwb , udtnb , udtsb ,
     :              vdteb , vdtwb , vdtnb , vdtsb ,
     :              wdteb , wdtwb , wdtnb , wdtsb ,
     :              ptdteb, ptdtwb, ptdtnb, ptdtsb,
     :              pdteb , pdtwb , pdtnb , pdtsb ,
     :              qvdteb, qvdtwb, qvdtnb, qvdtsb,
     :              qcdteb, qcdtwb, qcdtnb, qcdtsb,
     :              qrdteb, qrdtwb, qrdtnb, qrdtsb,
     :              qidteb, qidtwb, qidtnb, qidtsb,
     :              qsdteb, qsdtwb, qsdtnb, qsdtsb,
     :              qhdteb, qhdtwb, qhdtnb, qhdtsb)

      ELSEIF( ebc.eq.5 .or. wbc.eq.5 .or. nbc.eq.5 .or. sbc.eq.5 ) THEN

c       CALL EXTBDT( )

        write(6,'(1x,a)') 
     :       'Option 5 of boundary conditions not implemented yet.'

      ENDIF
c
c#######################################################################
c
c     Swap time levels for all dependent variables to prepare for
c     next timestep integration.  
c
c     The swap is as follows:
c
c     Fields at time tpresent is assigned to tpast,
c     fields at time tfuture  is assigned to tpresent,
c     fields at time tfuture  is not changed. 
c
c#######################################################################
c

      CALL TFLIP(nx,ny,nz, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh)

c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TINTEG                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
c
      SUBROUTINE TINTEG(frstep, nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           udteb, udtwb, udtnb, udtsb,
     :           vdteb, vdtwb, vdtnb, vdtsb,
     :           wdteb, wdtwb, wdtnb, wdtsb,
     :           ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :           pdteb ,pdtwb ,pdtnb ,pdtsb,
     :           qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :           qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :           qrdteb,qrdtwb,qrdtnb,qrdtsb,
     :           qidteb,qidtwb,qidtnb,qidtsb,
     :           qsdteb,qsdtwb,qsdtnb,qsdtsb,
     :           qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :           x,y,z,zp, j1,j2,j3, 
     :           usflx,vsflx,ptsflx,qvsflx, km,
     :           uforce,vforce,wforce,pforce, 
     :           tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
c
c#######################################################################
c
c     PURPOSE:
c
c     Orchestrate the time integration of the basic governing equations 
c     for a single timestep. Note that the time tendencies of the 
c     non-conservative processes (e.g., microphysics, radiation, etc.) 
c     are not included in this routine but are applied in their 
c     corresponding routines.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/21/92.
c
c     MODIFICATION HISTORY:
c
c     5/03/92 (K. Droegemeier)
c       Added full documentation.
c
c     5/20/92 (M. Xue)
c       Reformatted certain documentations.
c
c#######################################################################
c
c     INPUT :
c
c       frstep   Flag to determine if this is the initial time step.
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at times tpast and tpresent (m/s)
c       v        y component of velocity at times tpast and tpresent (m/s)
c       w        Vertical component of Cartesian velocity at times
c                tpast and tpresent (m/s)
c       wcont    Vertical component of contravariant velocity in
c                computational coordinates (m/s)
c       ptprt    Perturbation potential temperature at times tpast and
c                tpresent (K)
c       pprt     Perturbation pressure at times tpast and tpresent (Pascal)
c       qv       Water vapor specific humidity at times tpast and tpresent (kg/kg)
c       qc       Cloud water mixing ratio at times tpast and tpresent (kg/kg)
c       qr       Rainwater mixing ratio at times tpast and tpresent (kg/kg)
c       qi       Cloud ice mixing ratio at times tpast and tpresent (kg/kg)
c       qs       Snow mixing ratio at times tpast and tpresent (kg/kg)
c       qh       Hail mixing ratio at times tpast and tpresent (kg/kg)
c
c       udteb    Time tendency of u field at east boundary (m/s**2)
c       udtwb    Time tendency of u field at west boundary (m/s**2)
c       udtnb    Time tendency of u field at north boundary (m/s**2)
c       udtsb    Time tendency of u field at south boundary (m/s**2)
c
c       vdteb    Time tendency of v field at east boundary (m/s**2)
c       vdtwb    Time tendency of v field at west boundary (m/s**2)
c       vdtnb    Time tendency of v field at north boundary (m/s**2)
c       vdtsb    Time tendency of v field at south boundary (m/s**2)
c
c       wdteb    Time tendency of w field at east boundary (m/s**2)
c       wdtwb    Time tendency of w field at west boundary (m/s**2)
c       wdtnb    Time tendency of w field at north boundary (m/s**2)
c       wdtsb    Time tendency of w field at south boundary (m/s**2)
c
c       ptdteb   Time tendency of ptprt field at east boundary (K/s)
c       ptdtwb   Time tendency of ptprt field at west boundary (K/s)
c       ptdtnb   Time tendency of ptprt field at north boundary (K/s)
c       ptdtsb   Time tendency of ptprt field at south boundary (K/s)
c
c       pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
c       pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
c       pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
c       pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
c
c       qvdteb   Time tendency of qv field at east boundary (kg/(kg*s))
c       qvdtwb   Time tendency of qv field at west boundary (kg/(kg*s))
c       qvdtnb   Time tendency of qv field at north boundary (kg/(kg*s))
c       qvdtsb   Time tendency of qv field at south boundary (kg/(kg*s))
c
c       qcdteb   Time tendency of qc field at east boundary (kg/(kg*s))
c       qcdtwb   Time tendency of qc field at west boundary (kg/(kg*s))
c       qcdtnb   Time tendency of qc field at north boundary (kg/(kg*s))
c       qcdtsb   Time tendency of qc field at south boundary (kg/(kg*s))
c
c       qrdteb   Time tendency of qr field at east boundary (kg/(kg*s))
c       qrdtwb   Time tendency of qr field at west boundary (kg/(kg*s))
c       qrdtnb   Time tendency of qr field at north boundary (kg/(kg*s))
c       qrdtsb   Time tendency of qr field at south boundary (kg/(kg*s))
c
c       qidteb   Time tendency of qi field at east boundary (kg/(kg*s))
c       qidtwb   Time tendency of qi field at west boundary (kg/(kg*s))
c       qidtnb   Time tendency of qi field at north boundary (kg/(kg*s))
c       qidtsb   Time tendency of qi field at south boundary (kg/(kg*s))
c
c       qsdteb   Time tendency of qs field at east boundary (kg/(kg*s))
c       qsdtwb   Time tendency of qs field at west boundary (kg/(kg*s))
c       qsdtnb   Time tendency of qs field at north boundary (kg/(kg*s))
c       qsdtsb   Time tendency of qs field at south boundary (kg/(kg*s))
c
c       qhdteb   Time tendency of qh field at east boundary (kg/(kg*s))
c       qhdtwb   Time tendency of qh field at west boundary (kg/(kg*s))
c       qhdtnb   Time tendency of qh field at north boundary (kg/(kg*s))
c       qhdtsb   Time tendency of qh field at south boundary (kg/(kg*s))
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c       usflx    Surface flux of u-momentum (kg/(m*s**2))
c       vsflx    Surface flux of v-momentum (kg/(m*s**2))
c       ptsflx   Surface heat flux (K*kg/(m**2 * s ))
c       qvsflx   Surface moisture flux of (kg/(m**2 * s))
c
c     OUTPUT:
c
c       u        x component of velocity at time tfuture (m/s)
c       v        y component of velocity at time tfuture (m/s)
c       w        Vertical component of Cartesian velocity at tfuture (m/s)
c       ptprt    Perturbation potential temperature at time tfuture (K)
c       pprt     Perturbation pressure at time tfuture (Pascal)
c       qv       Water vapor specific humidity at time tfuture (kg/kg)
c       qc       Cloud water mixing ratio at time tfuture (kg/kg)
c       qr       Rainwater mixing ratio at time tfuture (kg/kg)
c       qi       Cloud ice mixing ratio at time tfuture (kg/kg)
c       qs       Snow mixing ratio at time tfuture (kg/kg)
c       qh       Hail mixing ratio at time tfuture (kg/kg)
c
c       udteb    Time tendency of u field at east boundary (m/s**2)
c       udtwb    Time tendency of u field at west boundary (m/s**2)
c       udtnb    Time tendency of u field at north boundary (m/s**2)
c       udtsb    Time tendency of u field at south boundary (m/s**2)
c
c       vdteb    Time tendency of v field at east boundary (m/s**2)
c       vdtwb    Time tendency of v field at west boundary (m/s**2)
c       vdtnb    Time tendency of v field at north boundary (m/s**2)
c       vdtsb    Time tendency of v field at south boundary (m/s**2)
c
c       wdteb    Time tendency of w field at east boundary (m/s**2)
c       wdtwb    Time tendency of w field at west boundary (m/s**2)
c       wdtnb    Time tendency of w field at north boundary (m/s**2)
c       wdtsb    Time tendency of w field at south boundary (m/s**2)
c
c       ptdteb   Time tendency of ptprt field at east boundary (K/s)
c       ptdtwb   Time tendency of ptprt field at west boundary (K/s)
c       ptdtnb   Time tendency of ptprt field at north boundary (K/s)
c       ptdtsb   Time tendency of ptprt field at south boundary (K/s)
c
c       pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
c       pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
c       pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
c       pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
c
c       qvdteb   Time tendency of qv field at east boundary (kg/(kg*s))
c       qvdtwb   Time tendency of qv field at west boundary (kg/(kg*s))
c       qvdtnb   Time tendency of qv field at north boundary (kg/(kg*s))
c       qvdtsb   Time tendency of qv field at south boundary (kg/(kg*s))
c
c       qcdteb   Time tendency of qc field at east boundary (kg/(kg*s))
c       qcdtwb   Time tendency of qc field at west boundary (kg/(kg*s))
c       qcdtnb   Time tendency of qc field at north boundary (kg/(kg*s))
c       qcdtsb   Time tendency of qc field at south boundary (kg/(kg*s))
c
c       qrdteb   Time tendency of qr field at east boundary (kg/(kg*s))
c       qrdtwb   Time tendency of qr field at west boundary (kg/(kg*s))
c       qrdtnb   Time tendency of qr field at north boundary (kg/(kg*s))
c       qrdtsb   Time tendency of qr field at south boundary (kg/(kg*s))
c
c       qidteb   Time tendency of qi field at east boundary (kg/(kg*s))
c       qidtwb   Time tendency of qi field at west boundary (kg/(kg*s))
c       qidtnb   Time tendency of qi field at north boundary (kg/(kg*s))
c       qidtsb   Time tendency of qi field at south boundary (kg/(kg*s))
c
c       qsdteb   Time tendency of qs field at east boundary (kg/(kg*s))
c       qsdtwb   Time tendency of qs field at west boundary (kg/(kg*s))
c       qsdtnb   Time tendency of qs field at north boundary (kg/(kg*s))
c       qsdtsb   Time tendency of qs field at south boundary (kg/(kg*s))
c
c       qhdteb   Time tendency of qh field at east boundary (kg/(kg*s))
c       qhdtwb   Time tendency of qh field at west boundary (kg/(kg*s))
c       qhdtnb   Time tendency of qh field at north boundary (kg/(kg*s))
c       qhdtsb   Time tendency of qh field at south boundary (kg/(kg*s))
c
c       km       Turbulent mixing coefficient for momentum (m**2/s).
c
c     WORK ARRAYS:
c
c       uforce   Acoustically inactive forcing terms in u-Eq. ((kg/(m*s)**2)
c       vforce   Acoustically inactive forcing terms in v-Eq. ((kg/(m*s)**2)
c       wforce   Acoustically inactive forcing terms in w-Eq. ((kg/(m*s)**2)
c       pforce   Acoustically inactive forcing terms in p-eq. (Pascal/s)
c
c       tem1     Temporary work array
c       tem2     Temporary work array
c       tem3     Temporary work array
c       tem4     Temporary work array
c       tem5     Temporary work array
c       tem6     Temporary work array
c       tem7     Temporary work array
c       tem8     Temporary work array
c       tem9     Temporary work array
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work 
c       arrays before you alter the code.)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none        ! Force explicit declarations
c
      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      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)
c
      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)
c
      real udteb (ny,nz)        ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)        ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)        ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)        ! T-tendency of u at s-boundary (m/s**2)

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

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

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

      real pdteb (ny,nz)        ! T-tendency of pprt at e-boundary (PASCAL/s)
      real pdtwb (ny,nz)        ! T-tendency of pprt at w-boundary (PASCAL/s)
      real pdtnb (nx,nz)        ! T-tendency of pprt at n-boundary (PASCAL/s)
      real pdtsb (nx,nz)        ! T-tendency of pprt at s-boundary (PASCAL/s)

      real qvdteb(ny,nz)        ! T-tendency of qv at e-boundary (kg/(kg*s))
      real qvdtwb(ny,nz)        ! T-tendency of qv at w-boundary (kg/(kg*s))
      real qvdtnb(nx,nz)        ! T-tendency of qv at n-boundary (kg/(kg*s))
      real qvdtsb(nx,nz)        ! T-tendency of qv at s-boundary (kg/(kg*s))

      real qcdteb(ny,nz)        ! T-tendency of qc at e-boundary (kg/(kg*s))
      real qcdtwb(ny,nz)        ! T-tendency of qc at w-boundary (kg/(kg*s))
      real qcdtnb(nx,nz)        ! T-tendency of qc at n-boundary (kg/(kg*s))
      real qcdtsb(nx,nz)        ! T-tendency of qc at s-boundary (kg/(kg*s))

      real qrdteb(ny,nz)        ! T-tendency of qr at e-boundary (kg/(kg*s))
      real qrdtwb(ny,nz)        ! T-tendency of qr at w-boundary (kg/(kg*s))
      real qrdtnb(nx,nz)        ! T-tendency of qr at n-boundary (kg/(kg*s))
      real qrdtsb(nx,nz)        ! T-tendency of qr at s-boundary (kg/(kg*s))

      real qidteb(ny,nz)        ! T-tendency of qi at e-boundary (kg/(kg*s))
      real qidtwb(ny,nz)        ! T-tendency of qi at w-boundary (kg/(kg*s))
      real qidtnb(nx,nz)        ! T-tendency of qi at n-boundary (kg/(kg*s))
      real qidtsb(nx,nz)        ! T-tendency of qi at s-boundary (kg/(kg*s))

      real qsdteb(ny,nz)        ! T-tendency of qs at e-boundary (kg/(kg*s))
      real qsdtwb(ny,nz)        ! T-tendency of qs at w-boundary (kg/(kg*s))
      real qsdtnb(nx,nz)        ! T-tendency of qs at n-boundary (kg/(kg*s))
      real qsdtsb(nx,nz)        ! T-tendency of qs at s-boundary (kg/(kg*s))

      real qhdteb(ny,nz)        ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)        ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)        ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)        ! T-tendency of qh at s-boundary (kg/(kg*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 rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity 
                                ! (kg/kg) 

      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 j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).
c
      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 flux of heat (K*kg/(m*s**2))
      real qvsflx(nx,ny)        ! Surface flux of moisture (kg/(m**2*s))
c
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
c
      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)
c
      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
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      real dtbig1               ! local value of big timestep
      integer frstep            ! Flag for the initial time step
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'
      include 'globcst.inc'
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external frcuvw
      external frcp
      external acoust
      external solvpt
      external solvqv
      external solvq
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c
c#######################################################################
c
c     To initialize the three-time level scheme, we assume for
c     the first step that the values of all variables at time 
c     past equal the values at time present.  We then perform a
c     forward-in-time integration for the first step only, and
c     centered-in-time integrations for all subsequent steps
c     using the same algorithm.  The size of the timestep is
c     adjusted accordingly, such that the leapfrog step size is
c     twice that of the first forward-in-time step.
c
c#######################################################################
c
      IF(frstep.eq.1) THEN      ! frstep is an interation counter
                                ! that equals 1 on first step only
                                ! to indicate forward-in-time integ
        dtbig1 = dtbig/2

      ELSE

        dtbig1 = dtbig

      ENDIF
c
c
c#######################################################################
c
c     Compute the acoustically inactive terms in the momentum and
c     pressure equations that are held fixed during the small time
c     step computations.  This includes advection, buoyancy, mixing
c     (both physical and computational), and the Coriolis force.
c     These forcing terms are accumulated into arrays for each
c     of the momentum equations, e.g., uforce for the u-equation,
c     vforce for the v-equation, and wforce for the w-equation.
c
c#######################################################################
c
      CALL FRCUVW(nx,ny,nz,
     :            u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :            ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :            usflx,vsflx, x,y,z,zp, j1,j2,j3,
     :            uforce,vforce,wforce, km,
     :            tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
c
c#######################################################################
c
c     Now, do the pressure equation, with the forcing stored in
c     array pforce.
c
c#######################################################################
c
      CALL FRCP(nx,ny,nz,
     :          u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :          ptbar,pbar,rhobar,qvbar,
     :          pforce,
     :          tem1,tem2,tem3,tem4)

c
c#######################################################################
c
c     Integrate the momentum and pressure equations (i.e., the
c     acoustically-active equations) in time using a mode-splitting
c     approach.  The mometum components are advanced in time using 
c     a forward scheme (relative to the pressure gradient force terms) 
c     and the pressure is advanced using a backward scheme (relative to 
c     the divergence term, which is evaluated using the newly updated 
c     velocities hence the backward scheme). These small time steps
c     bring the momentum and pressure from time tpast to tfuture, and
c     during these steps the acoustically-inactive forcing terms (e.g., 
c     uforce, pforce, etc) are held fixed at time tpresent (i.e. they
c     are evaluated only once at time tpresent each large time step
c     integration cycle), therefore the time integration is leapfrog 
c     in time relative to these forcing terms. 
c
c#######################################################################
c

      CALL ACOUST(frstep, nx,ny,nz,
     :            u,v,w,pprt, 
     :            udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,
     :            wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,
     :            pbar,rhobar, 
     :            uforce,vforce,wforce,pforce, 
     :            x,y,z,zp, j1,j2,j3,
     :            tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
c
c
c#######################################################################
c 
c     Integrate the potential temperature equation forward one
c     timestep.  Note that, since this equation contains no
c     acoustically-active terms, it is integrated on the large
c     timestep only.  Sources and sinks due to both phase changes and
c     radiation processes are handled separately.
c
c#######################################################################
c

      CALL SOLVPT(nx,ny,nz, dtbig1,
     :            ptprt,u,v,w, ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :            ptbar,rhobar,km,ptsflx,
     :            x,y,z,zp, j1,j2,j3,
     :            tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9)
c
c
c#######################################################################
c 
c     Integrate the liquid water substance continuity equations 
c     forward one timestep.
c
c     Sources and sinks due to phase changes and radiative processes 
c     are handled separately.
c
c     If the simulation is for dry dynamics, then we skip over
c     the moisture computations to save computer resources.  
c     The flag for this choice is "moist", i.e.,
c
c     Moist = 0 for dry run
c     Moist = 1 for moist run
c
c#######################################################################
c

      IF( moist.eq.1) THEN   ! Test if this run is dry or moist.

c
c#######################################################################
c
c     Water vapor equation
c
c#######################################################################
c

        CALL SOLVQV(nx,ny,nz, dtbig1,
     :              qv,u,v,w, qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :              rhobar,qvbar, km, qvsflx,
     :              x,y,z,zp, j1,j2,j3,
     :              tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

c
c#######################################################################
c
c     Cloud water equation
c
c#######################################################################
c
        CALL SOLVQ(nx,ny,nz, dtbig1,
     :             qc,u,v,w, qcdteb,qcdtwb,qcdtnb,qcdtsb, rhobar,km,
     :             x,y,z,zp, j1,j2,j3,
     :             tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

c
c#######################################################################
c
c     Rain water equation
c
c#######################################################################
c
        CALL SOLVQ(nx,ny,nz, dtbig1,
     :             qr,u,v,w, qrdteb,qrdtwb,qrdtnb,qrdtsb, rhobar,km,
     :             x,y,z,zp, j1,j2,j3,
     :             tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
c
      ENDIF

c
c#######################################################################
c
c     Integrate the frozen water substance continuity equations 
c     forward one timestep.
c
c     Sources and sinks due to phase changes and radiative processes 
c     are handled separately.
c
c     If the simulation does not consider ice processes, we skip over
c     the computations for ice variables to save computer resources.  
c     The flags for these choices are "moist" and "ice", i.e.,
c
c     Moist = 0 for dry run
c     Moist = 1 for moisture with or without ice processes. 
c     and moist =1, Ice = 1 for a moist run with ice. 
c
c#######################################################################
c

      IF( moist.eq.1 .and. ice.eq.1) THEN    ! If this run includes ice. 

c
c#######################################################################
c
c     Cloud ice equation
c
c#######################################################################
c
        CALL SOLVQ(nx,ny,nz, dtbig1,
     :             qi,u,v,w, qidteb,qidtwb,qidtnb,qidtsb, rhobar,km,
     :             x,y,z,zp, j1,j2,j3,
     :             tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

c
c#######################################################################
c
c     Snow equation
c
c#######################################################################
c
        CALL SOLVQ(nx,ny,nz, dtbig1,
     :             qs,u,v,w, qsdteb,qsdtwb,qsdtnb,qsdtsb, rhobar,km,
     :             x,y,z,zp, j1,j2,j3,
     :             tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

c
c#######################################################################
c
c     Hail equation
c
c#######################################################################
c
        CALL SOLVQ(nx,ny,nz, dtbig1,
     :             qh,u,v,w, qhdteb,qhdtwb,qhdtnb,qhdtsb, rhobar,km,
     :             x,y,z,zp, j1,j2,j3,
     :             tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)
c
      ENDIF
c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ACOUST                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE ACOUST(frstep, nx,ny,nz,
     :           u,v,w,pprt, 
     :           udteb,udtwb,udtnb,udtsb,vdteb,vdtwb,vdtnb,vdtsb,
     :           wdteb,wdtwb,wdtnb,wdtsb,pdteb,pdtwb,pdtnb,pdtsb,
     :           pbar,rhobar, 
     :           uforce,vforce,wforce,pforce, 
     :           x,y,z,zp, j1,j2,j3,
     :           tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

C#######################################################################
c
c     PURPOSE:
c
c     Coordinate the integration of the acoustically active parts of 
c     the dynamic equations. 
c
C#######################################################################
c
c     AUTHOR: Ming Xue
c     3/21/92.
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (K. Droegemeier and M. Xue)
c     Added full documentation.
c
c#######################################################################
c
c     INPUT:
c
c       frstep   Flag to determine if this is the initial time step.
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity at times tpast and tpresent (m/s)
c       v        y component of velocity at times tpast and tpresent (m/s)
c       w        Vertical component of Cartesian velocity at times
c                tpast and tpresent (m/s)
c       pprt     Perturbation pressure at times tpast and tpresent (Pascal)
c
c       udteb    Time tendency of u field at east boundary (m/s**2)
c       udtwb    Time tendency of u field at west boundary (m/s**2)
c       udtnb    Time tendency of u field at north boundary (m/s**2)
c       udtsb    Time tendency of u field at south boundary (m/s**2)
c
c       vdteb    Time tendency of v field at east boundary (m/s**2)
c       vdtwb    Time tendency of v field at west boundary (m/s**2)
c       vdtnb    Time tendency of v field at north boundary (m/s**2)
c       vdtsb    Time tendency of v field at south boundary (m/s**2)
c
c       wdteb    Time tendency of w field at east boundary (m/s**2)
c       wdtwb    Time tendency of w field at west boundary (m/s**2)
c       wdtnb    Time tendency of w field at north boundary (m/s**2)
c       wdtsb    Time tendency of w field at south boundary (m/s**2)
c
c       pdteb    Time tendency of pprt field at east boundary (PASCAL/s)
c       pdtwb    Time tendency of pprt field at west boundary (PASCAL/s)
c       pdtnb    Time tendency of pprt field at north boundary (PASCAL/s)
c       pdtsb    Time tendency of pprt field at south boundary (PASCAL/s)
c
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state density (kg/m**3)
c
c       uforce   Acoustically inactive forcing terms in u-eq. (kg/(m*s)**2)
c       vforce   Acoustically inactive forcing terms in v-eq. (kg/(m*s)**2)
c       wforce   Acoustically inactive forcing terms in w-eq. (kg/(m*s)**2)
c       pforce   Acoustically inactive forcing terms in p-eq. (Pascal/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       u        x component of velocity at time tfuture (m/s) 
c       v        y component of velocity at time tfuture (m/s)
c       w        Vertical component of Cartesian velocity at tfuture (m/s) 
c       pprt     Perturbation pressure at time tfuture (Pascal)
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array
c       tem2     Temporary work array
c       tem3     Temporary work array
c       tem4     Temporary work array
c       tem5     Temporary work array
c       tem6     Temporary work array
c       tem7     Temporary work array
c       tem8     Temporary work array
c
c      (These arrays are defined and used locally (i.e. inside this
c       subroutine), they may also be passed into routines called by
c       this one. Exiting the call to this subroutine, these temporary
c       work arrays may be used for other purposes, and therefore their
c       contents may be overwritten. Please examine the usage of work
c       arrays before you alter the code.)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none        ! Force explicit declarations
c
      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

      integer nx, ny, nz        ! Number of grid points in 3 directions
c
      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 pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
c
      real udteb (ny,nz)        ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)        ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)        ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)        ! T-tendency of u at s-boundary (m/s**2)

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

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

      real pdteb (ny,nz)        ! T-tendency of pprt at e-boundary (PASCAL/s)
      real pdtwb (ny,nz)        ! T-tendency of pprt at w-boundary (PASCAL/s)
      real pdtnb (nx,nz)        ! T-tendency of pprt at n-boundary (PASCAL/s)
      real pdtsb (nx,nz)        ! T-tendency of pprt at s-boundary (PASCAL/s)
c
      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)

      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 j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - d( zp )/d( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! d( zp )/d( z ).

      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 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
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer frstep
      real    dtsml1            ! local value of small timestep
      integer i, j, k
      integer ismstp
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'
      include 'globcst.inc'
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external solvuvw
      external solvp
c
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     To initialize the three-time level scheme, we assume for
c     the first step that the values of all variables at time
c     past equal the values at time present.  We then perform a
c     forward-in-time integration for the first step only, and
c     centered-in-time integrations for all subsequent steps
c     using the same algorithm.  The timestep size for the first
c     time step is therefore half that of the following steps.
c
c#######################################################################
c
      IF(frstep.eq.1) THEN 

        dtsml1 = dtsml/2

      ELSE

        dtsml1 = dtsml

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

      DO 900 ismstp =1, nsmstp

c
c#######################################################################
c
c     On the first time through the small timestepping loop, set
c     the values at time tfuture equal to those at tpast.  By doing
c     this, we simplify the forward-in-time integrations because the
c     equations take the form:
c
c     u(i,j,k,future) = u(i,j,k,future) 
c    :                + dtsml*( acoustic +  non-acoustic forcing terms )
c
c     that is, we need only use a single time level in this two-time
c     level scheme, and we finish with u at time tfuture automatically.
c
c#######################################################################
c

        IF(ismstp.eq.1) THEN

          DO 950 k=1,nz
          DO 950 j=1,ny
          DO 950 i=1,nx
            pprt(i,j,k,tfuture)=pprt(i,j,k,tpast)
            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)
950       CONTINUE

        ENDIF
c
c#######################################################################
c
c     Advance the three momentum equations in time one small
c     timestep. (Putting arguments like u(1,1,1,tfuture) means
c     passing the array values at time tfuture only into the subroutine.
c     Inside, these arrays are only defined at a single time level,
c     which match the array values at tfuture outside).
c
c#######################################################################
c
        CALL SOLVUVW(nx,ny,nz, dtsml1,
     :               u(1,1,1,tfuture),v   (1,1,1,tfuture),
     :               w(1,1,1,tfuture),pprt(1,1,1,tfuture),
     :               udteb,udtwb,udtnb,udtsb,
     :               vdteb,vdtwb,vdtnb,vdtsb,
     :               wdteb,wdtwb,wdtnb,wdtsb,
     :               rhobar,uforce,vforce,wforce,
     :               x,y,z,zp, j1,j2,j3,
     :               tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8)

c
c#######################################################################
c
c     Advance the pressure equation in time one small timestep.
c
c#######################################################################
c
        CALL SOLVP(nx,ny,nz,dtsml1, 
     :             pprt(1,1,1,tfuture),u(1,1,1,tfuture),
     :             v   (1,1,1,tfuture),w(1,1,1,tfuture), 
     :             pdteb,pdtwb,pdtnb,pdtsb,
     :             rhobar,pbar, pforce,
     :             x,y,z,zp, j1,j2,j3,
     :             tem1,tem2,tem3,tem4)

900   CONTINUE        ! End small timestep loop

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TFILT                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE TFILT(nx,ny,nz, 
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh)
c
c#######################################################################
c
c     PURPOSE:
c
c     Apply the Asselin time filter to all time dependent variables. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91. 
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (K. Droegemeier and M. Xue)
c     Added full documentation.
c
c#######################################################################
c            

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

      integer nt                ! The no. of t-levels of t-dependent arrays.
      parameter (nt=3)

      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)
c
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'debug.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external aselin
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF( subtrac.eq.1 ) THEN 

        write(6,'(1x,a,f13.3)') 'TFILT was called with flteps =', flteps

      ENDIF

      
      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)

      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( moist.ne.0 ) 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.ne.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)

        ENDIF

      ENDIF

      RETURN
      END
      
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ASELIN                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE ASELIN(a, nx,ny,nz, nx1,nx2,ny1,ny2,nz1,nz2, flteps)

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

c
c#######################################################################
c
c     INPUT:
c
c       a        A variable at three time levels.
c
c       nx1,nx2  The bounds of 1st index within which filter is applied
c       ny1,ny2  The bounds of 2nd index within which filter is applied
c       nz1,nz2  The bounds of 3rd index within which filter is applied
c
c       flteps   The asselin time filter coefficient. 
c                Defined in globcst.inc
c
c     OUTPUT: 
c
c       a        The new array values at time tpresent after the time 
c                filter is applied. 
c
c#######################################################################
c
c
c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

      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
c
c#######################################################################
c
c     Misc. local variable:
c
c#######################################################################
c
      integer i, j, k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      DO 100 k=nz1,nz2
      DO 100 j=ny1,ny2
      DO 100 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))
100   CONTINUE

      RETURN
      END

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

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none

      integer nt                ! The no. of t-levels of t-dependent arrays.
      parameter (nt=3)

      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) 
c
      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)
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external tswap
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      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( moist.ne.0 ) THEN 

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

        IF( ice.ne.0 ) THEN 

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

        ENDIF

      ENDIF

      RETURN
      END

c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TSWAP                      ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
 
      SUBROUTINE TSWAP(nx,ny,nz, a)
c
c#######################################################################
c
c     PURPOSE:
c
c     Update an array in time.
c
c     a(*,*,*,tpast   ) = a(*,*,*,tpresent)
c     a(*,*,*,tpresent) = a(*,*,*,tfuture )
c     a(*,*,*,tfuture ) is not changed. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91. 
c
c     MODIFICATION HISTORY:
c
c     5/20/92 (K. Droegemeier and M. Xue)
c     Added full documentation.
c
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       a        A 3-D array at three time levels that will be updated
c                in time.
c
c     OUPTUT:
c
c       a        Its new values at time tpast and tfuture.
c
c#######################################################################

c
c#######################################################################
c
c     Variable Declarations:
c
c#######################################################################
c
      implicit none

      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

      parameter (nt=3, tpast=1, tpresent=2, tfuture=3)

      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
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
c
c#######################################################################
c
c     Routines called: None.
c
c#######################################################################
c
 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      DO 100 k=1,nz
      DO 100 j=1,ny
      DO 100 i=1,nx
        a(i,j,k,tpast   )=a(i,j,k,tpresent)
        a(i,j,k,tpresent)=a(i,j,k,tfuture )
100   CONTINUE

      RETURN
      END
