c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STEPU                      ######
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 STEPU(nx,ny,nz, dtsml1,
     :           u,uforce,upgrad,rhobar)
c
c#######################################################################
c
c     PURPOSE:
c
c     Advance the u velocity field forward one time step. The time integration 
c     scheme is forward in time relative to the acoustically active pressure 
c     gradient term, but is centered in time relative to the other acoustically 
c     inactive terms that are stored in uforce. This routine is called by the 
c     acoustic wave integration routine ACOUST. The small time step dtsml1 is 
c     used for the acoustic wave integration.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/04/92 (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       dtsml1   Integration time step. 
c                It is the small time step for acoustic wave integration.
c
c       u        x-component of velocity at a given time level (m/s).
c       uforce   all acoustically inactive forcing terms in u-eq.
c                (kg/(m*s)**2)
c       upgrad   Pressure gradient force in u-eq. (kg/(m*s)**2)
c       rhobar   Base state density (kg/m**3)
c
c     OUTPUT:  
c
c       u        x-component of velocity updated in time by dtsml1 (m/s).
c
c#######################################################################
c           

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

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

      real u     (nx,ny,nz)     ! u-velocity at a given time level (m/s).
      real uforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in u-momentum equation (kg/(m*s)**2)
      real upgrad(nx,ny,nz)     ! Pressure gradient force in u-eq.
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3).

c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
      real dtsml1
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Note here that needed average on rhobar in x-direction is not
c     performed since rhobar is assumed to be independent of x.
c
c     When terrain is included, the average must be performed.
c
c#######################################################################
c
      DO 100 k=2,nz-2
      DO 100 j=2,ny-2
      DO 100 i=2,nx-1
        u(i,j,k)=u(i,j,k)
     :       +dtsml1 *(uforce(i,j,k)-upgrad(i,j,k))/rhobar(i,j,k)
100   CONTINUE
    
c
      RETURN
      END


c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STEPV                      ######
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 STEPV(nx,ny,nz, dtsml1,
     :           v,vforce,vpgrad,rhobar)
c
c#######################################################################
c
c     PURPOSE:
c
c     Advance the v-velocity field forward one time step. The time integration 
c     scheme is forward in time relative to the acoustically active pressure 
c     gradient term, but is centered in time relative to the other acoustically 
c     active terms that are stored in vforce. This routine is called by the 
c     acoustic wave integration routine ACOUST. The small time step dtsml1 is 
c     used for the acoustic wave integration.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/04/92 (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       dtsml1   Integration time step.
c                It is the small time step for acoustic wave integration.
c
c       v        y-component of velocity at a given time level (m/s).
c       vforce   all acoustically inactive forcing terms in v-eq.
c                (kg/(m*s)**2)
c       vpgrad   Pressure gradient force in v-eq. (kg/(m*s)**2)
c       rhobar   Base state density (kg/m**3)
c
c     OUTPUT:
c
c       v        y-component of velocity updated in time by dtsml1 (m/s).
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

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

      real v     (nx,ny,nz)     ! v-velocity at a given time level (m/s).
      real vforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in v-momentum equation (kg/(m*s)**2)
      real vpgrad(nx,ny,nz)     ! Pressure gradient force in v-eq.
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3).
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
      real dtsml1
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c 
c#######################################################################
c
c     Note here that needed average on rhobar in y-direction is not
c     performed since rhobar is assumed to be independent of y.
c
c     When terrain is included, the average must be performed.
c
c#######################################################################
c
      DO 100 k=2,nz-2
      DO 100 j=2,ny-1
      DO 100 i=2,nx-2
        v(i,j,k)=v(i,j,k)
     :       +dtsml1 *(vforce(i,j,k)-vpgrad(i,j,k))/rhobar(i,j,k)
100   CONTINUE
c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STEPW                      ######
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 STEPW(nx,ny,nz, dtsml1,
     :      w,wforce,wpgrad,rhobar, tem1)
c
c#######################################################################
c
c     PURPOSE:
c
c     Advance the w-velocity field forward one time step.The time integration 
c     scheme is forward in time relative to the acoustically active pressure 
c     gradient term, but is centered in time relative to the other acoustically 
c     inactive terms stored in wforce. This routine is called by the acoustic 
c     wave integration routine ACOUST. The small time step dtsml1 is used for 
c     the acoustic wave integration.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/04/92 (M. Xue)
c     Added full documentation.
c
c     6/05/92 (M. Xue)
c     Further facelift. 
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       dtsml1   Integration time step.
c                It is the small time step for acoustic wave integration.
c
c       w        z-component of velocity at a given time level (m/s).
c       wforce   all the acoustically inactive forcing terms in w-eq. 
c                (kg/(m*s)**2)
c       wpgrad    Pressure gradient force in w-eq. (kg/(m*s)**2)
c       rhobar   Base state density (kg/m**3)
c
c     OUTPUT:
c
c       w        z-component of velocity updated in time by dtsml1 (m/s).
c
c     WORK ARRAY:
c
c       tem1     work array.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

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

      real w     (nx,ny,nz)     ! w-velocity at a given time level (m/s).
      real wforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in w-momentum equation (kg/(m*s)**2)
      real wpgrad(nx,ny,nz)     ! Pressure gradient force in w-eq.
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3).

      real tem1  (nx,ny,nz)     ! Local temporary array.
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
      integer onvf          ! indicate if output var. is on volume face
      real dtsml1
c
c#######################################################################
c
c     Routine called: 
c
c#######################################################################
c
      external avgz
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
      onvf = 1
      CALL avgz(rhobar, onvf,
     :      nx,ny,nz, 2,nx-2, 2,ny-2, 2,nz-1, tem1)

      DO 100 k=2,nz-1
      DO 100 j=2,ny-2
      DO 100 i=2,nx-2
        w(i,j,k)=w(i,j,k)
     :       +dtsml1 *(wforce(i,j,k)-wpgrad(i,j,k))/tem1(i,j,k)
100   CONTINUE    
c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STEPP                      ######
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 STEPP(nx,ny,nz, dtsml1,
     :           pprt,rhobar,pforce,pdiv)
c
c#######################################################################
c
c     PURPOSE:
c
c     Advance the pressure perturbation forward by small time steps dtsml1. 
c     The scheme used is backward in time relative to the acoustically 
c     active divergence term, but is centered in time (leapfrog) relative 
c     to the other inactive terms stored in pforce.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/04/92 (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       dtsml1   Integration time step.
c                It is the small time step for acoustic wave integration.
c
c       pprt     Perturbation pressure at a given time level (Pascal)
c       rhobar   Base state density (kg/m**3)
c       pforce   Acoustically inactive forcing terms in p-eq (Pascal/s)
c       pdiv     The divergence term in pressure equation (kg*/(m*s**3))
c
c     OUTPUT:  
c
c       pprt     Perturbation pressure updated in time by dtsml1 (Pascal)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

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

      real pprt (nx,ny,nz)      ! Perturbation potential temperature
                                ! (Kelvin).
      real pdiv(nx,ny,nz)       ! The divergence term in the pressure
                                ! Equation (kg/(m*s**3))
      real pforce(nx,ny,nz)     ! Acoustically inactive forcing terms
                                ! in pressure equation (Pascal/s)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3).
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
      real dtsml1
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      DO 100 k=2,nz-2
      DO 100 j=2,ny-2
      DO 100 i=2,nx-2
        pprt(i,j,k)=pprt(i,j,k)
     :             +dtsml1 *(pforce(i,j,k)-pdiv(i,j,k))
100   CONTINUE
c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STEPPT                     ######
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 STEPPT(nx,ny,nz, dtbig1,
     :           ptprt,rhobar,ptadv,ptmix,ptsrc)
c
c#######################################################################
c
c     PURPOSE:
c
c     Advance the potential temperature perturbation forward by one large time 
c     step dtbig1.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/04/92 (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       dtbig1   Integration time step.
c                It is the small time step for acoustic wave integration.
c
c       ptprt    Perturbation potential temperature at tpast and tpresent(K)
c       ptadv    Advection term in the potential temp. eq. (kg*K/(m**3*s))
c       ptmix    Total mixing in the potential temp. eq. (kg*K/(m**3*s))
c       ptsrc    Other source/sink term in potential temperature eq.
c       rhobar   Base state density (kg/m**3)
c
c     OUTPUT:  
c
c       ptprt    Perturbation potential temperature at time future (K)
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx, ny, nz   ! Number of grid points in x, y, z direction

      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)

      real ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K).
      real ptadv (nx,ny,nz)     ! Advection term
      real ptmix (nx,ny,nz)     ! Total mixing
      real ptsrc (nx,ny,nz)     ! Other source/sink term
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3).
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
      real dtbig1, dt2
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
      dt2 = dtbig1*2
      DO 100 k=2,nz-2
      DO 100 j=2,ny-2
      DO 100 i=2,nx-2
        ptprt(i,j,k,tfuture)=ptprt(i,j,k,tpast)
     :       +dt2*(-ptadv(i,j,k)+ptmix(i,j,k)+ptsrc(i,j,k)) 
     :       /rhobar(i,j,k)
100   CONTINUE
c
      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE STEPQ                      ######
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 STEPQ(nx,ny,nz, dtbig1, q, qadv,qmix, rhobar)
c
c#######################################################################
c
c     PURPOSE:
c
c     Advance a water substance quantity q forward by one large time step dtbig1.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/91.
c
c     MODIFICATION HISTORY:
c
c     5/04/92 (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       dtbig1   Integration time step.
c                It is the small time step for acoustic wave integration.
c
c       q        Mixing ratio of a water substance at tpast and tpresent. 
c       qadv     Advection term of q-eq. (kg/(m**3*s)) 
c       qmix     Total mixing on q. (kg/(m**3*s))
c       rhobar   Base state density (kg/m**3)
c
c     OUTPUT:  
c
c       q        Mixing ratio of a water substance at time tfuture (kg/kg)     
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none

      integer nx, ny, nz   ! Number of grid points in x, y, z direction

      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)

      real q     (nx,ny,nz,nt)  ! Mixing ratio of a water substance (kg/kg) 

      real qadv  (nx,ny,nz)     ! Advection of q. 
      real qmix  (nx,ny,nz)     ! Total mixing term.
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3).
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i, j, k
      real dtbig1, dt2
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      dt2 = dtbig1*2
      DO 100 k=2,nz-2
      DO 100 j=2,ny-2
      DO 100 i=2,nx-2
        q(i,j,k,tfuture)=q(i,j,k,tpast)
     :       +dt2*(-qadv(i,j,k)+qmix(i,j,k))/rhobar(i,j,k)
100   CONTINUE
c
      RETURN
      END
