c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INIT                       ######
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 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     PURPOSE:
c
c     Initialize the model parameters and variables.
c
c     This is the main driver for model initializations.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/5/92.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
c
c     9/14/1992 (M. Xue)
c     Different surface drag coefficients defined for mementum, heat and 
c     moisture. 
c     Three options included for the Coriolis force calculations. 
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     OUTPUT:
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        z-component of velocity at all time 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       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       vbar     Base state y-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 computational space (m)
c       y        y-coordinate of grid points in computational 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       km       Turbulent mixing coef. Filled with zeros.
c
c       temxy1   2-d temporary array
c       temxy2   2-d temporary array
c       temxy3   2-d temporary array
c       temxy4   2-d temporary array
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#######################################################################
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          ! The 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
                                ! from that of base state atmosphere (Kelvin).
      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 
                                ! w-point on the staggered grid.
      real hterain(nx,ny)       ! The height of terrain.

      real j1    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - del( zp )/del( x ).
      real j2    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! - del( zp )/del( y ).
      real j3    (nx,ny,nz)     ! Coordinate transform Jacobian defined as
                                ! del( zp )/del( z ).

      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for 
                                ! momentum. ( m**2/s )      

      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.
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external initpara 
      external init0
      external inigrd 
      external initvar 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     First, initialize the model run parameters. Most of them are 
c     global parameters passed among subroutines through common blocks.
c
c     These parameters are declared in include files 
c     dims.inc, globcst.inc, bndry.inc debug.inc etc.
c     
c
c#######################################################################
c
      CALL initpara(nx,ny,nz)
c
C#######################################################################
c
c     Set all arrays to zero before starting the model.
c
c     This step is done for safety reasons only. The initialization
c     of model variables should not rely on this.
c
C#######################################################################
c      
      CALL init0(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar, km,
     :           temxy1,temxy2,temxy3,temxy4,
     :           tem1,tem2,tem3,tem4,tem5,tem6,tem7,
     :           tem8,tem9,tem10,tem11,tem12,tem13)
c
c#######################################################################
c
c     INItialization of model GRiD definition arrays. 
c
c#######################################################################
c

      CALL inigrd(nx,ny,nz,
     :            x,y,z,zp,hterain,j1,j2,j3)

c
c#######################################################################
c
c     INITialize model array VARiables.
c
c#######################################################################
c
      CALL initvar(nx,ny,nz, 
     :             x,y,z,zp,hterain,j1,j2,j3,
     :             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,
     :             tem1)
c
C#######################################################################
c
c     End of model initialization.
c
C#######################################################################
c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INITPARA                   ######
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 INITPARA(nx,ny,nz)

c
c#######################################################################
c
c     PURPOSE:
c
c     Iniitialize the model control PARAmeters. Most of them are read in 
c     from an input file. This subroutine also produces a log file which can be used 
c     as the input file for replicating a particular model run.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c
c     3/17/1991. 
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c       Added full documentation.
c
c     5/25/92 (M. Xue)
c       Reworked to provide a friendly user interface for control 
c       parameter input, and write out a log file.
C
c     6/04/92 (M. Xue)
c       Further facelift.
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       Added output control parameter input etc.
c
c#######################################################################
c
c     INPUT:
c     
c       nx       Number of grid points in the x-direction
c       ny       Number of grid points in the y-direction
c       nz       Number of grid points in the z-direction
c
c     OUTPUT:
c
c       Parameters declared in include file 'bccntl.inc':
c
c       ebc      Parameter defining the east   boundary condition type.
c       wbc      Parameter defining the west   boundary condition type.
c       nbc      Parameter defining the north  boundary condition type.
c       sbc      Parameter defining the south  boundary condition type.
c       tbc      Parameter defining the top    boundary condition type.
c       bbc      Parameter defining the bottom boundary condition type.
c
c                When the value of these parameters are:
c
c                  1, boundary condition is the regid wall type;
c                  2, boundary condition is periodic;
c                  3, zero normal gradient boundary condition is specified;
c                  4, radiation (open) boundary condition is specified;
c                  5, boundary conditon is by specified user externally.
c
c       Parameters declared in include file 'globcst.inc':
c
c       runname  The name of this job (a string at least 6 character long).
c
c       initopt  Option control parameter for model initialization.
c
c                initopt = 1, Self initialization (e.g. analytical field);
c                initopt = 2, Restart run, initialize from restart file;
c                initopt = 3, Initialize from external data set.         
c
c       rstinf   File name of restart input data.
c       rstiunt  FORTRAN I/O unit used by restart data read.
c
c       inifile  File name of initial fields data for option 3
c       inifunt  FORTRAN I/O unit number for initial data input.
c
c       restrt   Parameter designating if this is a restart run.
c
c                For a non-restart run restrt = 1;
c                for a restart run restrt = 0. 
c
c       sndfile  Name of the input sounding file.
c       sfunit   FORTRAN I/O unit used by sounding input.
c
c       ptpert0  The magnitude of initial potential temperature 
c                perturbation (thermal bubble) (K)
c       pt0radx  Radius of bubble in x-direction (m) (Default value set)
c       pt0rady  Radius of bubble in y-direction (m) (Default value set)
c       pt0radz  Radius of bubble in z-direction (m) (Default value set)
c       pt0ctrx  x coordinate of bubble center (m) (Default value set)
c       pt0ctry  y coordinate of bubble center (m) (Default value set)
c       pt0ctrz  z coordinate of bubble center (m) (Default value set)
c
c       dx       Grid spacing in x-direction in computational 
c                as well as in physical space (m)
c       dy       Grid spacing in y-direction in computational 
c                as well as in physical space (m)
c       dz       Grid spacing in vertical direction in 
c                transformed computataional space (m).
c       dxinv    Reciprocal of dx (i.e. 1/dx)
c       dyinv    Reciprocal of dy (i.e. 1/dy)
c       dzinv    Reciprocal of dz (i.e. 1/dz)
c
c       xl       Physical scale of the model domain in x-direction (m). 
c                In the uniform spacing case xl = (nx-3)*dx.
c       yl       Physical scale of the model domain in y-direction (m). 
c                In the uniform spacing case yl = (ny-3)*dy.
c       zh       Physical scale of the model domain in z-direction (m)
c
c       csopt    Parameter that determines the definition of sound
c                wave speed used in the pressure equation. 
c
c                csopt = 1, csound = cp/cv * rd * T;
c                csopt = 2, csound = cp/cv * rd * T * csfactr;
c                csopt = 3, csound = specfied constant. 
c
c       csfactr  A multiplication factor for sound speed (m/s).
c       csound   User specified constant sound speed (m/s). 
c
c       dtbig    Model integration time step (s).
c       dtsml    Small time step size for the acoustic wave integrations (s).
c
c       tstart   The time when the model is initialized (s). 
c       tstop    The time when time integration is to be stopped (s).
c       curtim   Current model time (s).
c
c       nstep    The number time steps already taken.
c       nsteps   The total number of time steps to be taken in this run.
c       nsmstp   The number of small acoustic wave time 
c                integration steps per large time step.
c
c       flteps   Coef. of Robert-Asselin time filter (non-dimensional).
c
c       tmixopt  Control parameter defining the turbulent mixing options.
c
c                tmixopt = 0, No turbulent mixing,
c                tmixopt = 1, Constant mixing coefficient km = tmixcst,
c                tmixopt = 2, Smagorinsky formulation of mixing coeff.,
c                tmixopt = 3, Smagorinsky mixing coefficient 
c                             plus a constant coeffcient tmixcst.
c
c       prantl   Turbulent prandtl number, prantl = km/kh.
c       tmixcst  Value of the constant mixing coefficient (m**2/s).
c
c       kcmix2   Control switch for 2nd order comp. mixing. =0 or 1 
c       cfcmh2   Coef. of 2nd order horizontal computational mixing (m**2/s)
c       cfcmv2   Coef. of 2nd order vertical computational mixing (m**2/s)
c
c       kcmix4   Control switch for 4th order comp. mixing. =0 or 1
c       cfcmh4   Coef. of 4th order horizontal computational mixing (m**4/s)
c       cfcmv4   Coef. of 4th order vertical computational mixing (m**4/s)
c
c       raydmp   Control parameter for Rayleigh damping.
c                raydmp = 1, damping is on; = 0, damping is off.
c       zbrdmp   Height of the bottom of Rayleigh damping (sponge) layer (m).
c       ztrdmp   Height of the top of Rayleigh damping (sponge) layer (m).
c       cfrdmp   Rayleigh damping coefficient (1/second).
c
c       divdmp   Control parameter for acoustic wave divergence damping. 
c
c                divdmp = 1, div. damping is on; 
c                       = 0, div. damping is off.
c
c       cdvdmp   Coef. of the divergence damping term. (m**2/s). 
c
c       mphyopt  Parameter specifying the microphysics parameterization 
c                scheme to be chosed.
c
c                mphyopt = 0, No microphysics. 
c                mphyopt = 1, Kessler type warm rain microphysics.
c                mphyopt = 2, More elaborate scheme - GEMS of Straka. 
c                             Not implemented yet.
c
c       moist    Control parameter for moist processes in the  model run. 
c
c                If moist = 0, the run is dry, all processes 
c                related to moisture are switched off, including advection.
c                If moist = 1, moist processes are activated.
c
c       ice      Control parameter for ice microphysical processes 
c                in the model run.
c
c                If ice   = 0, all processes related to ice are switched 
c                              off, including the advection of ice variables.  
c                If moist = 1, ice processes are activated.
c
c       umove    Speed of model domain translation in x direction (m/s).
c       vmove    Speed of model domain translation in y direction (m/s).
c
c       sfcphy   Control parameter specifying the option of
c                surface physics parameterizations. 
c
c                sfcphy = 0, no surface physics parameterization
c                sfcphy = 1, simple bulk aerodynamic formulation. 
c 
c       cdm      Bulk aerodynamic surface drag coefficient for momentum 
c       cdh      Bulk aerodynamic drag coefficient for surface heat flux 
c       cdq      Bulk aerodynamic drag coefficient for surface mositure flux
c                cdm, cdh and cdq are nondimensional.
c
c       coriopt  Options switch for Coriolis term defintions
c                coriopt = 0, Coriolis terms are set to zero,
c                coriopt = 1, Only Coriolis terms with horiziontal winds are 
c                             included. 
c                coriopt = 2, Coriolis terms with vertical wind are also 
c                             included. 
c
c       latitud  The Earth latitude in degree at the center of model domain.
c                fcorio1  Coriolis parameter 2*omega*sin(latitud) (1/s)
c                fcorio2  Coriolis parameter 2*omega*cos(latitud) (1/s)
c
c       nfmtprt  Time steps between formatted print out
c       tfmtprt  Time interval between formatted print out (s)
c
c       nmaxmin  Time steps between max/min statistics calculations.
c       tmaxmin  Time interval between max/min statistics calxulations (s)
c
c       nenergy  Time steps between energy statistics calculations.
c       tenergy  Time interval between energy statistics calxulations (s)
c
c       nrstout  Time steps    between restart data dumps
c       trstout  Time interval between restart data dumps (s)
c
c       nhisdmp  Time steps    between history data dumps
c       thisdmp  Time interval between history data dumps (s)
c
c       basout   Parameter controlling base state field output. 
c                basout = 0 or 1.
c       mstout   Parameter controlling output of moist variables
c                qv, qc, qr, qi, qs and qh.
c       iceout   Parameter controlling output of ice variables
c                qi, qs and qh. iceout = 1 or 0 for yes or no.
c
c       hdmpfmt  Parameter specifying the format of history data dump. 
c
c                hdmpfmt = 1, unformatted binary data dump;
c                hdmpfmt = 2, formatted ascii data dump;
c                hdmpfmt = 3, NCSA HDF format data dump;
c                hdmpfmt = 4, Packed binary data dump; 
c                hdmpfmt = 5, Data dump for Savi3D. 
c
c       nchdmp   FORTRAN I/O channel for history data dump.
c       hdmpfn   History data dump file name.
c
c       rstount  FORTRAN I/O unit used by restart data dump.
c       rstoutf  File name of restart data dump.
c
c
c       nchmax   FORTRAN I/O channel for max/min statistics output.
c
c       Parameters declared in include file 'debug.inc':
c
c       subtrac  Control switch for debug information printing.
c
c                If subtrac=1, information is printed when each 
c                              subroutine is called.
c                If subtrac=0, printing is switched off. 
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations. (Local Variables)
c
c#######################################################################
c
      implicit none 

      integer nx,ny,nz      ! The number of grid points in 3 directions
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c

      integer i

      character logfn*80    ! A string used as the log file filename.
      integer   llogfn      ! The length of the log file filename.
      integer   logfunt     ! FORTRAN unit number for log file output.
      integer   lenstr      ! Length of a string
      integer   istat       ! Flag set by open statement on the status
                            ! of file opening
      real      alpha       ! Non-dimensional coeff. of div. damping.
      logical   iexist      ! Flag set by inquire statement for file
                            ! existence
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c

c
c#######################################################################
c
c     Global constants and parameters, most of them specify the 
c     model run options. 
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Control parameters defining the boundary condition types. 
c
c#######################################################################
c
      include 'bndry.inc'
c
c#######################################################################
c
c     Parameters controlling debug information printing. 
c
c#######################################################################
c
      include 'debug.inc'
c
c#######################################################################
c
c     Universal physical constants such as gas constants. 
c
c#######################################################################
c
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external strlnth
      external gtlogfn
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      write(6,'(/ 16(/5x,a)//)')

     :'###############################################################',
     :'###############################################################',
     :'#####                                                     #####',
     :'#####                  Welcome to                         #####',
     :'#####                                                     #####',
     :'#####   The Advanced Regional Predication System  (ARPS)  #####',
     :'#####                                                     #####',
     :'#####                  Version 3.0                        #####',
     :'#####                                                     #####',
     :'#####                Copyright (c) 1992                   #####',
     :'#####  Center for the Analysis and Prediction of Storms   #####',
     :'#####    University of Oklahoma.  All rights reserved.    #####',
     :'#####                                                     #####',
     :'###############################################################',
     :'###############################################################'
c
c
c#######################################################################
c
c     Now we begin to read in the values of parameters:
c
c####################################################################### 
c
      Write(6,'(4(/5x,a),/6(/5x,a)/)') 
     :  'You begin running this model by inputing a number of control',
     :  'parameters. These parameters are prompted for if you run the',
     :  'job interactively. In a batch mode, they are read in from an',
     :  'input file.',
     :  'At the end of all parameter input, a log file is produced',
     :  'which can be directly used as the input file when you want',
     :  'to replicate the same job. This file is named runnam.log.nn',
     :  'where runnam is a standard prefix for all output files that',
     :  'are produced by this job and nn is a number appended to the',
     :  'file name when file runnam.log.nn-1 already exists.'
c
c#######################################################################
c
c     First, read in some comments lines on this job and the name of 
c     this run designated by a string at least 6 character long.
c
c#######################################################################
c 

      Write(6,'(/2x,a,/3(/5x,a),a)') 

     :     'COMMENT LINES AND JOB NAME',
     :     'You can provide the model with a few lines of comments',
     :     'which designate the nature of this model run.',
     :     'You need to tell the model first how many lines of ',
     :     'comments you have.'

      Write(6,'(/5x,a/)')
     :     'Input the number of comment lines (0 to 50):'

      read (5,*) nocmnt

      IF( nocmnt.gt.50 ) THEN

        write(6,'(/1x,a,/1x,a,i3,a/)') 
     :  'Too many lines of comments, the number can not be more ',
     :  'than 50, input was ',nocmnt,'. Job is stopped in INITPARA.'

      ENDIF

      write(6,'(/5x,i4,a)')
     :     nocmnt,' comment lines will be read.'

      IF( nocmnt.ne.0 ) THEN

        Write(6,'(/5x,a,/5x,a/)')
     :  'Now input the specified number of comment lines as character',
     :  'strings (in quotes):'
      
        DO 10 i = 1, nocmnt
          read (5,*) cmnt(i)
10      CONTINUE

        write(6,'(/5x,a/)')
     :       'The comment lines for this run are:'

        DO 12 i = 1, nocmnt
          write(6,'(5x,a)') cmnt(i)
12      CONTINUE

      ENDIF

      Write(6,'(/5x,a,/5x,a,/5x,a,/5x,a/))') 
     :     'The model requires a 6 character long string to construct',
     :     'the names of output files. This string is also the name of '
     :     ,'this job.',
     :     'Please input such a string (at least 6 character long):'

      read (5,*) runname

      write(6,'(/5x,a,a)') 'The name of this run is: ', runname 

c
C#######################################################################
c
c     Read in the parameter that controls the model run mode:
c
C#######################################################################
c
      Write(6,'(/8(/5x,a),/4(/5x,a)))') 
     :     'The model can be run in 3D, 2D xz, 2D yz and 1D z mode.',
     :     'For each case the internally determined grid scale is ',
     :     'calculated accordingly. Other parameters should be properly'
     :    ,'set by the user. e.g. for 2D xz case, set ny>=4, dy=dx,',
     :     'and periodic boundary condition in the y direction, then ',
     :     'initially y-independent fields should remain so.',
     :     'Please set the option value for the mode of model run,',
     :     'by selecting a number:',
     :     '  1 for 3-D run;',
     :     '  2 for 2-D xz plane run;',
     :     '  3 for 2-D yz plane run;',
     :     '  4 for vertical 1-D run.'

      read(5,*) runmod 
c
C#######################################################################
c
c     Read in control parameter INITOPT for model initialization 
c
c     INITOPT = 1, Self initialization (e.g. specify perturbation using 
c                  analytical functions),
c             = 2, Restart run, initialize the model using previous
c                  model output,
c             = 3, Initialize the model using external input data file.
c
c     For options 2 and 3, the names of input files need to be provided.
c
C#######################################################################
c

      write(6,'(/2x,a,//5x,a,//5x,a,a,/5x,a,/5x,a,a,/5x,a,//5x,a/)')

     :     'MODEL INITIALIZATION OPTIONS:',
     :     'Three options:',
     :     '  1  Self initialization (e.g. specify initial ',
     :           'perturbations',
     :     '     using analytical functions);',
     :     '  2  Restart run, initialize model using previous ',
     :           'model output;',
     :     '  3  Initialize the model from external data file.',
     :     'Select a number for the model initialization:'

      read (5,*) initopt
      write(6,'(5x,a,i4)') 'Input was ', initopt

      IF( initopt.eq.2 ) THEN

        restrt = 1

      ELSE

        restrt = 0

      ENDIF
c
C#######################################################################
c
c     Specify the types of initial potential temperature perturbations: 
c
C#######################################################################
c
      write(6,'(2(/5x,a),//5x,a,/5x,a/)')
     :    'For initialization option 1, two types of initial potential',
     :    'temperature perturbation can be specified. Type in',
     :    '  1 for bubble-shaped perturbation;',
     :    '  2 for random perturbation (the horizontal average of the',
     :    '    perturbation will be zero):'

      read(5,*) pt0opt
      write(6,'(/5x,a,i4/)') 'Input was ', pt0opt
c
C#######################################################################
c
c     Magnitude, size and location of the initial bubble-shaped 
c     thermal perturbation: 
c
C#######################################################################
c
c
C#######################################################################
c
c     Magnitude, size and location of the initial bubble-shaped 
c     thermal perturbation: 
c
C#######################################################################
c
      write(6,'(3(/5x,a)/)')
     :   'For initialization option 1, a bubble-shaped perturbation',
     :   'is specified for the potential temperature field.',
     :   'Please give the magnitude of this perturbation (K):'    

      read (5,*) ptpert0
      write(6,'(/5x,a,f10.3,a)')
     :   'The magnitude of the initial perturbation is ', ptpert0,' K.' 

      write(6,'(3(/5x,a)/)')
     :   'Please give the radius of thermal bubble in x direction (m)',
     :   '(If the value given is negative, it''s treated as infinitely',
     :   'large so that the bubble is invariant in this direction):' 
      read (5,*) pt0radx 

      write(6,'(3(/5x,a)/)')
     :   'Please give the radius of thermal bubble in y direction (m)',
     :   '(If the value given is negative, it''s treated as infinitely',
     :   'large so that the bubble is invariant in this direction):' 
      read (5,*) pt0rady 

      write(6,'(3(/5x,a)/)')
     :   'Please give the radius of thermal bubble in z direction (m)',
     :   '(If the value given is negative, it''s treated as infinitely',
     :   'large so that the bubble is invariant in this direction):' 
      read (5,*) pt0radz 

      write(6,'(2(/5x,a)/)')
     :   'Please give the bubble center location in x direction (m)',
     :   '(x coordinate of the bubble center):'
      read (5,*) pt0ctrx 

      write(6,'(2(/5x,a)/)')
     :   'Please give the bubble center location in y direction (m)',
     :   '(y coordinate of the bubble center):'
      read (5,*) pt0ctry 

      write(6,'(2(/5x,a)/)')
     :   'Please give the bubble center location in z direction (m)',
     :   '(z coordinate of the bubble center):'
      read (5,*) pt0ctrz 

      write(6,'(/5x,a,3e10.3,a,/5x,a,a,/5x,a,e10.3,a,e10.3,a,e10.3,a))')
     :   'The input radii of the thermal bubble are ',
     :    pt0radx,pt0rady,pt0radz,' (meter)',
     :   'in x, y and z direction recpectively, and the center is ',
     :   'located at','x=',pt0ctrx,' y=',pt0ctry,' z=',pt0ctrz,
     :   ' (meter).' 
c
C#######################################################################
c
c     Data files for initializations:
c
C#######################################################################
c
 
      write(6,'(//2x,a,//5x,a,/5x,a/)')
     :     'INITIAL DATA FILES:',
     :     'Input the name of the re-start data file used by ',
     :     'initialization option (a character string in quotes):'

      read(5,*) rstinf 

      lenstr = 80
      CALL strlnth( rstinf, lenstr)
      write(6,'(/5x,a,a)')
     :'The restart data file to be read in is ', rstinf(1:lenstr) 

      rstiunt = 30

      write(6,'(/5x,a,i2,a)') 
     :     'FORTRAN unit ',rstiunt,
     :     ' will be used for restart data input.'

      write(6,'(/5x,a,/5x,a/)')
     :     'Input the name of initial data file used by initialization',
     :     'option 3 (a character string in quotes):'

      read(5,*) inifile

      lenstr = 80
      CALL strlnth( inifile, lenstr)
      write(6,'(/5x,a,a)')
     :     'The initial data file to be read is ', inifile(1:lenstr) 

      inifunt = 31

      write(6,'(/5x,a,i2,5x,a)') 
     :     'FORTRAN unit ',inifunt,
     :     ' will be used for initial data input.'

      rstount = 32
c
C#######################################################################
c
c     The environmental sounding.
c
C#######################################################################
c
      Write(6,'(//2x,a,//5x,a, /5x,a,a, /5x,a,a/)')

     :     'THE SOUNDING FILE:',
     :     'For all three options of model initializations, a sounding',
     :     'file should be provided to initialize the model base ',
     :     'state. ',
     :     'Please give the name of the sounding file ',
     :     '(a character string):'

      read (5,*) sndfile

      lenstr = 80
      CALL strlnth( sndfile, lenstr )
      write(6,'(/5x,a,a)') 
     :     'Sounding file to be used is ', sndfile(1:lenstr) 

      sfunit = 33
      write(6,'(/5x,a,i2,a)') 
     :     'FORTRAN unit ',sfunit,
     :     ' will be used for sounding file input.'


c
C#######################################################################
c
c     Grid size and time step size
c
C#######################################################################
c
      write(6,'(//2x,a/)')
     :   'THE MODEL GRID INTERVALS (uniform grid interval at present):'

      write(6,'(/5x,a/)')
     :     'Grid interval in x direction (meter):'

      read (5,*) dx
      write(6,'(/5x,a,f10.3,a)')
     :     'Input was dx=',dx,' meters'

      write(6,'(/5x,a/)')
     :     'Grid interval in y direction (meter):'

      read (5,*) dy
      write(6,'(/5x,a,f10.3,a)')
     :     'Input was dy=',dy,' meters'

      write(6,'(/5x,a/)')
     :     'Grid interval in z direction (meter):'

      read (5,*) dz
      write(6,'(/5x,a,f10.3,a)')
     :     'Input was dz=',dz,' meters'

      write(6,'(//2x,a,//5x,a,16(/5x,a))')

     :   'TIME INTEGRATION STEP SIZE:',
     :   'Mode splitting scheme is used for the time integration.',
     :   'Acoustically active terms in the equations are treated in',
     :   'small time step integrations, so that the major time step',
     :   'size is not limited by the acoustic wave speed. This major',
     :   '(big) time step is limited by the CFL linear stability',
     :   'criterion for advection, gravity wave propagation, and',
     :   'that for diffusive processes.',
     :   ' ',
     :   'A rough estimate of the large time step is given by',
     :   ' ',
     :   'dtbig < 1/[ (Vmax+C)*sqrt(1/dx**2+1/dy**2+1/dz**2) )].',
     :   ' ',
     :   'where Vmax is the maximum wind speed, C the internal gravity',
     :   'wave speed. When dx>>dz, adjustment is needed to the formula.'

      write(6,'(/5x,a/)')
     :     'Input the value of large time step size (seconds):'

      read (5,*) dtbig
      write(6,'(/5x,a,f10.3,a)')
     :     'Input was ',dtbig,' seconds.'

      write(6,'(2(/5x,a),2(//5x,a))')
     :     'The size of small time step for acoustic wave integration',
     :     'is limited by',
     :     'dtsml < 1/[ Cs*sqrt(1/dx**2+1/dy**2+1/dz**2) )].',
     :     'where Cs is the acoustic wave speed.'

      write(6,'(5(/5x,a),2(//5x,a),/5x,a,/5x,a,//5x,a/))')
     :     'In most meteorological applications, super-compressbile',
     :     'approximation can be made with no significant effect on',
     :     'atmospheric flows. With this approximation, the acoustic',
     :     'wave speed is artificially reduced so that a larger time',
     :     'step can be used for acoustic wave integrations.',
     :     'We have three options for choosing the acoustic wave speed:'
     :    ,'  1  calculated internally  as Cs = (Cp/Cv*Rd*T);',
     :     '  2  multiplied by a factor as Cs = (csfactr*Cp/Cv*Rd*T);',
     :     '  3  specified externally as a constant.' ,
     :     'Please choose an option (1 suggested):'

      read (5,*) csopt
      write(6,'(/5x,a,i4)') 'Input was ',csopt

      write(6,'(2(/5x,a)/)') 
     :     'For option 2, please input the value of multiplcation',
     :     'factor (a value between 0.5 and 1.0 is suggested):'

      read (5,*) csfactr
      write(6,'(/5x,a,f10.3)') 'Input was ', csfactr
 
      write(6,'(2(/5x,a)/)') 
     :     'For option 3, please input the speed of acoustic',
     :     '(a value between 150.0 and 350.0 m/s is suggested):'

      read (5,*) csound
      write(6,'(/5x,a,f10.3,a)') 'Input was ', csound,' m/s.'
 

      write(6,'(/5x,a/)')
     :     'Now input the value of small time step size (seconds):'

      read (5,*) dtsml
      write(6,'(/5x,a,f10.3,a)') 
     :     'Input was ',dtsml,' seconds.'

      nsmstp = int( 2*dtbig/dtsml ) 

      IF( mod(2*dtbig,dtsml) .ne.0.0 ) nsmstp =nsmstp +1

      dtsml = (2*dtbig)/nsmstp

      write(6,'(//5x,a,i3,a)')
     :     'There are ',nsmstp,' small steps in every large step.'

      write(6,'(/5x,a,f10.3,a)') 
     :     'The actual small time step size to be used is ',
     :     dtsml,' seconds.'
c
C#######################################################################
c
c     Boundary condition control parameters:
c
C#######################################################################
c
      write(6,'(//2x,a/)') 
     :     'BOUNDARY CONDITION CONTROL PARAMETERS'

      write(6,'(2(/5x,a),//5x,a,4(/5x,a),//5x,a,a)') 
     :     'Five of types lateral boundary conditions are available.',
     :     'They are:',
     :     '  1  rigid wall (mirror) type boundary condition;',
     :     '  2  periodic boundary condition;',
     :     '  3  zero gradient boundary condition;',
     :     '  4  radiation (open) boundary condition;',
     :     '  5  externally spcified boundary condition.',
     :     'Choose the type of condition for each boundary ',
     :     'by selecting a number.'

      write(6,'(/5x,a/)') 
     :     'Select a number for the west boundary (2 is suggested):'
      read (5,*) wbc

      write(6,'(/5x,a/)') 
     :     'Select a number for the east boundary (2 is suggested):'
      read (5,*) ebc

      write(6,'(/5x,a/)') 
     :     'Select a number for the south boundary (2 is suggested):'
      read (5,*) sbc

      write(6,'(/5x,a/)') 
     :     'Select a number for the north boundary (2 is suggested):'
      read (5,*) nbc

      write(6,'(//5x,a, /4(/5x,i3,a)/)') 
     '     'The numbers Selected were:',
     '     wbc,' for  west boundary,',ebc,' for  east boundary,',
     '     sbc,' for south boundary,',nbc,' for north boundary.'

      write(6,'(2(/5x,a),/3(/5x,a),//5x,a,a)') 
     :     'Three types of conditions are available for the top and',
     :     'bottom boundaries. They are:',
     :     '  1  nonpenetrative ground and rigid lid;',
     :     '  2  periodic boundary condition;',
     :     '  3  zero gradient boundary condition;',
     :     'Choose the type of condition for each boundary ',
     :     'by Selecting a number.'


      write(6,'(/5x,a/)') 
     :     'Select a number for the top boundary (1 is suggested):'
      read (5,*) tbc

      write(6,'(/5x,a/)') 
     :     'Select a number for the bottom boundary (1 is suggested):'
      read (5,*) bbc

      write(6,'(/5x,a, /2(/5x,i3,a))') 
     '     'The numbers Selected were:',
     '     tbc,' for top boundary,',bbc,' for bottom boundary.'
c
C#######################################################################
c
c     Radiation lateral boundary condition options: 
c
C#######################################################################
c
      write(6,'(2(/5x,a),/2(/5x,a),//5x,a/)') 
     :     'Two optiontions of the radiation lateral boundary ',
     :     'conditions are avalaible. They are:',
     :     '  1  one dimensional formualtion;',
     :     '  2  multi-dimensional formulation.',
     :     'Please select an option:'
      read (5,*) rbcopt 

      write(6,'(/5x,a,i4/)') 
     :     'Input was ',rbcopt 
c
C#######################################################################
c
c     A few derived variables
c
C#######################################################################
c

      dxinv = 1.0/dx
      dyinv = 1.0/dy
      dzinv = 1.0/dz 

      xl = (nx-3)*dx
      yl = (ny-3)*dy
      zh = (nz-3)*dz
c
C#######################################################################
c
c     Calculate the grid scale depending on the model run mode: 
c
C#######################################################################
c
      IF( runmod. eq. 1 ) THEN

        grdscl = (dx*dy*dz)**(1.0/3.0)

      ELSEIF( runmod.eq.2 ) THEN

        grdscl = sqrt(dx*dz)

      ELSEIF( runmod.eq.3 ) THEN

        grdscl = sqrt(dy*dz)

      ELSEIF( runmod.eq.4 ) THEN

        grdscl = dz 

      ENDIF
c
C#######################################################################
c
c     Control parameters for microphysics parameterizations 
c
C#######################################################################
c

      write(6,'(//2x,a,2(//5x,a),3(/5x,a),//5x,a/)')

     :     'OPTIONS FOR MICROPHYSICS PARAMETERIZATIONS:',
     :     'Three options:',
     :     '  0  No microphysics parameterization;',
     :     '  1  Keseler warm rain microphysics parameterization;',
     :     '  2  Generalized ice phase microphysics (GEMS)',
     :     '     (not implemented yet).',
     :     'Select a number (1 suggested):'

      read (5,*) mphyopt

      write(6,'(/5x,a,i4)') 'Input was ',mphyopt

      write(6,'(//2x,a,2(//5x,a),(/5x,a),//5x,a/)')
     :     'FURTHER CONTROL ON MICROPHYSICS PROCESSES:',
     :     'Control parameter for moist process:',
     :     '  0  for dry model run;',
     :     '  1  for a moist run either with or without ice process.',
     :     'Select a number (1 suggested):'

      read (5,*) moist
      write(6,'(/5x,a,i4)') 'Input was ',moist

      write(6,'(2(//5x,a),(/5x,a),//5x,a/)')
     :     'Control parameter for ice process:',
     :     '  0  for a model run without ice processes;',
     :     '  1  for a model run with ice process.',
     :     'Select a number (0 suggested):'

      read (5,*) ice
      write(6,'(/5x,a,i4)') 'Input was ',ice

c
C#######################################################################
c
c     Input control parameters for subgrid scale turbulence mixing 
c
C#######################################################################
c

      write(6,'(//2x,a,2(//5x,a),4(/5x,a),//5x,a/)')

     :     'SUBGRID SCALE TURBULENCE MIXING PARAMETERIZATIONS:',
     :     'Four options:',
     :     '  0  no turbulent mixing parameterization;',
     :     '  1  turbulent mixing with a constant coefficient;',
     :     '  2  Smagorinsky formulation of mixing coefficient;',
     :     '  3  Smagorinsky formulation plus a constant coeffcient',
     :     '     as background mixing.',
     :     'Select a number (2 suggested):'

      read (5,*) tmixopt
      write(6,'(/5x,a,i4)')
     :     'Input was ',tmixopt
     :

      write(6,'(/5x,a,/5x,a/)') 
     :     'Input the value of constant mixing coefficient for options',
     :     '1 and 3 (m**2/s) (no effect for other options):'

      read (5,*) tmixcst 
      write(6,'(/5x,a,f10.3)') 'Input was ',tmixcst
c
C#######################################################################
c
c      Turbulent Prandtl number used in the turbulent mixing terms.
c
C#######################################################################
c
      write(6,'(/5x,a,/5x,a/)') 
     :    'Input the value of turbulent Prandtl number used by options',
     :    '2 and 3 of turbulent mixing (1.0 suggested):'

      read (5,*) prantl 
      write(6,'(/5x,a,f10.3)') 'Input was ',prantl
c
C#######################################################################
c
c     Input control parameters for computational mixing 
c
C#######################################################################
c

      write(6,'(//2x,a,/5x,a,2(//5x,a),/5x,a,//5x,a/)')
     :     'SECOND AND FOURTH ORDER COMPUTATIONAL MIXING (both can be',
     :     '     switched on at the same time)',
     :     'Control parameter for second order computational mixing:',
     :     '  0  switch off 2nd order computational mixing;',
     :     '  1  switch on  2nd order computational mixing.',
     :     'Select a number (1 suggested):'

      read (5,*) kcmix2
      write(6,'(/5x,a,i4)')  'Input was ',kcmix2


      write(6,'(/5x,a/)') 
     :     'Input 2nd order horizontal comp. mixing coefficient(m**2/s)' 

      read (5,*) cfcmh2
      write(6,'(/5x,a,f10.3,a)')  'Input was ',cfcmh2,' m**2/s.'

      write(6,'(/5x,a/)') 
     :     'Input 2nd order vertical comp. mixing coefficient (m**2/s):' 

      read (5,*) cfcmv2
      write(6,'(/5x,a,f10.3,a)')  'Input was ',cfcmv2,' m**2/s.'

      write(6,'(//5x,a,//5x,a,/5x,a,//5x,a/)')
     :     'Control parameter for fourth order computational mixing:',
     :     '  0  switch off 4th order computational mixing;',
     :     '  1  switch on  4th order computational mixing.',
     :     'Select a number (0 suggested):'

      read (5,*) kcmix4
      write(6,'(/5x,a,i4)')  'Input was ',kcmix4

      write(6,'(/5x,a/)') 
     :     'Input 4th order horizontal comp. mixing coefficient(m**4/s)' 
      read (5,*) cfcmh4
      write(6,'(/5x,a,e13.3,a)')  'Input was ',cfcmh4,' m**4/s.'

      write(6,'(/5x,a/)') 
     :     'Input 4th order vertical comp. mixing coefficient (m**4/s):' 

      read (5,*) cfcmv4
      write(6,'(/5x,a,e13.3,a)')  'Input was ',cfcmv4,' m**4/s.'

c
C#######################################################################
c
c     Controls for divergence damping on acoustic waves
c
C#######################################################################
c
      write(6,'(//2x,a,2(//5x,a),/5x,a,//5x,a/)') 
     :     'CONTROL PARAMETERS FOR DIVERGENCE DAMPING',
     :     'The switch for divergence damping:',
     :     '  0 switch off divergence damping;',
     :     '  1 switch on  divergence damping.',
     :     'Select a number (1 suggested):'

      read (5,*) divdmp
      write(6,'(/5x,a,i4)')  'Input was ', divdmp
c
      write(6,'(/5x,a,2(/5x,a)/5x,a,a,/5x,a)')
     :     'Input non-dimensional coeff. for div. damping ALPHA.',
     :     'where ALPHA is used to calculate dimensional divergence',
     :     'coefficient CDVDMP = ALPHA* DELTA**2/dtsml.',
     :     'DELTA is a measure of the grid size taken as the minimum ',
     :     'of DX, DY and DZ.',
     :     'Input ALPHA (0.01 suggested):' 

      read (5,*) alpha 
      cdvdmp = alpha * min(dx,dy,dz)**2 / dtsml
      write(6,'(/5x,a,f10.3,/5x,a,f10.3)')
     :     'Input was ',alpha,
     :     'Dimensional damping coefficient is ',cdvdmp

c
C#######################################################################
c 
c     Parameters for upper level Rayleigh damping:
c
C#######################################################################
c
      write(6,'(//2x,a)')
     :     'PARAMETERS FOR UPPER LEVEL RAYLEIGH DAMPING'

      write(6,'(/5x,a,//5x,a,/5x,a,//5x,a/)')
     :     'The switch for Rayleigh damping:',
     :     '  0  switch off Rayleight damping;',
     :     '  1  switch on  Rayleight damping.',
     :     'Select a number (1 is suggested):'

      read (5,*) raydmp
      write(6,'(/5x,a,i4)') 'Input was ',raydmp

      write(6,'(/5x,a,/5x,a/)')
     :     'Input the coefficient of the Rayleigh damping (1/second)',
     :     '( 1.0/(50*dtbig) is suggested):'

      read (5,*) cfrdmp
      write(6,'(/5x,a,f10.3)') 'Input was ', cfrdmp

      write(6,'(/5x,a,/5x,a/)')
     :     'Input the height of the bottom of the damping layer (m) ',
     :     '(12,000.0 meters suggested):'

      read (5,*) zbrdmp
      write(6,'(/5x,a,f10.3)') 'Input was ',zbrdmp

      write(6,'(/5x,a,/5x,a/)')
     :     'Input the height of the top of the damping layer (m) ',
     :     '(20,000.0 meters suggested):'

      read (5,*) ztrdmp
      write(6,'(/5x,a,f10.3)') 'Input was ',ztrdmp

c
C#######################################################################
c
c     Options of the surface flux parameterization.
c
C#######################################################################
c
      write(6,'(//2x,a)')
     :     'PARAMETERS FOR SURFACE PHYSICS PARAMETERIZATION'

      write(6,'(/5x,a,//5x,a,/5x,a,//5x,a/)')
     :     'Options for surface physics calculations:',
     :     '  0  no surface flux;',
     :     '  1  bulk aerodynamic formulation of surface fluxes.',
     :     'Select a number (0 suggested):'
      
      read (5,*) sfcphy
      write(6,'(/5x,a,i4)') 'Input was ', sfcphy
c
C#######################################################################
c
c     Surface drag coefficeints
c
C#######################################################################
c
      write(6,'(/5x,a,/5x,a/)')
     :     'Input the surface momentum drag coefficient',
     :     '(0.0 suggested):'

      read (5,*) cdm
      write(6,'(/5x,a,f10.3)') 'Input was ', cdm
c
      write(6,'(/5x,a,/5x,a/)')
     :     'Input the drag coefficient for surface heat flux',
     :     '(0.0 suggested):'

      read (5,*) cdh
      write(6,'(/5x,a,f10.3)') 'Input was ', cdh
c
      write(6,'(/5x,a,/5x,a/)')
     :     'Input the drag coefficient for surface moisture flux',
     :     '(0.0 suggested):'

      read (5,*) cdq
      write(6,'(/5x,a,f10.3)') 'Input was ', cdq
c
C#######################################################################
c
c     Coriolis parameters: 
c
C#######################################################################
c
      write(6,'(//5x,a/)')
     :     'Select the option of the Coriolis term definitions (0/1/2):' 

      read (5,*) coriopt
      write(6,'(/5x,a,i3)') 'Input was ', coriopt

      write(6,'(//5x,a/)')
     :     'Specifiy the latitude of the model domain center in degree:'

      read (5,*) latitud 
      write(6,'(/5x,a,f10.3)') 'Input was ', latitud

      IF( coriopt.eq.0 ) THEN 

        fcorio1 = 0.0
        fcorio2 = 0.0

      ELSEIF( coriopt.eq.1 ) THEN

        fcorio1 = 2*omega* sin( atan(1.0)/45.0 * latitud )
        fcorio2 = 0.0 

      ELSEIF( coriopt.eq.2 ) THEN

        fcorio1 = 2*omega* sin( atan(1.0)/45.0 * latitud )
        fcorio2 = 2*omega* cos( atan(1.0)/45.0 * latitud )

      ENDIF

      write(6,'(/5x,a,2e10.4)') 
     :  'Coriolis parameters fcorio1 and fcorio2 =', fcorio1, fcorio2 
c
C#######################################################################
c
c     Robert-Asselin time filter coefficient:
c
C#######################################################################
c
      write(6,'(/2x,a,//5x,a,/5x,a,/5x,a,//5x,a,/5x,a/)')
     :     'ROBERT-ASSELIN TIME FILTER',
     :     'The model uses the leapfrog scheme for time integration.',
     :     'A time filter needs to be applied to variables predicted',
     :     'by this scheme to suppress the computational mode.',
     :     'Please set the coefficient of the (Robert-Asselin) ',
     :     'time filter (0.1 suggested):'

      read (5,*) flteps
      write(6,'(/5x,a,f10.3)') 'Input was ',flteps
c
C#######################################################################
c
c     DOmain translation speed relative to earth:
c
C#######################################################################
c
      write(6,'(/2x,a,//5x,a,/5x,a,/5x,a)')
     :     'THE MODEL DOMAIN TRANSLATION SPEED',
     :     'In some simulations, the model domain (i.e. coordinate',
     :     'framework) is assumed to move at a given speed, so that',
     :     'interested systems remain inside the model domain.'

      write(6,'(/5x,a,/5x,a/)')
     :     'Specify the domain translation speed in x direction (m/s) ',
     :     '(0.0 suggested):'

      read (5,*) umove
      write(6,'(/5x,a,f10.3)') 'Input was ',umove

      write(6,'(/5x,a,/5x,a/)')
     :     'Specify the domain translation speed in y direction (m/s) ',
     :     '(0.0 suggested):'

      read (5,*) vmove
      write(6,'(/5x,a,f10.3)') 'Input was ',vmove

      write(6,'(/5x,a,/5x,a/)')
     :     'Velocity (umove, vmove) will be subtracted from',
     :     'the wind speed in the base state sounding.' 
c
C#######################################################################
c
c     Control parameters for time integration and input/output control 
c
C#######################################################################
c

      write(6,'(/2x,a)')
     :     'TIME INTEGRATION AND MODEL I/O CONTROLS:'
c
C#######################################################################
c
c     Model stop time and number of time steps to be taken. 
c
C#######################################################################
c

      write(6,'(/5x,a,/5x,a/)')
     :     'Give the time when you want to stop the model',
     :     'integration (seconds):'

      read (5,*) tstop

      write(6,'(/5x,a,f10.3)') 'Input was ',tstop

c
C#######################################################################
c
c     Model output parameter:
c    
c     First, give the name of the directory into which output files 
c     will be written: 
c
C#######################################################################
c
      write(6,'(3(/5x,a),/5x,a,a,/5x,a/)')
     :  'You can redirect the output data files into a separate',
     :  'directory rather than the current work directory, by',
     :  'specifying the full path name of such a directory.',
     :  'If you give a blank string, the current work directory ',
     :  'is assumed.',
     :  'The directory for output is:'

      read (5,*) dirname

      ldirnam = 80
      CALL strlnth( dirname, ldirnam)
      
      IF( ldirnam.eq.0 ) THEN
        dirname = '.'
        ldirnam=1
      ENDIF

      IF( dirname(1:ldirnam).ne.' ') THEN 
c
C#######################################################################
c
c     Check if the specified output directory exists, if not,
c     abort the job.
c
C#######################################################################
c
        inquire(file=dirname(1:ldirnam),exist=iexist)

        IF( .not.iexist ) THEN

c          write(6,'(/5x,a,2(/5x,a))')
c     :      'Specified output directory '//dirname(1:ldirnam)//
c     :      ' not found.',
c     :      'Please create it before starting the model.',
c     :      'Job aborted in subroutine INITPARA.'
          STOP

        ENDIF

c        write(6,'(/5x,a)')
c     :    'Output files will be in directory '//dirname(1:ldirnam)//'.'

      ELSE

        write(6,'(/5x,a)')
     :    'Output files will be in the current work directory.'

      ENDIF
c
c
C#######################################################################
c
c     Formatted print out: 
c
C#######################################################################
c

      write(6,'(/5x,a/)')
     :     'Specify the time interval between formatted print out (s):'

      read (5,*) tfmtprt
      write(6,'(/5x,a,f10.3,a)') 'Input was ',tfmtprt,' seconds.'

      nfmtprt =  nint(tfmtprt/dtbig)
      write(6,'(/5x,a,i6,a)')
     :     'Formatted printing is done every ', nfmtprt,' time steps.'
c
C#######################################################################
c
c     History data dump: 
c
C#######################################################################
c

      write(6,'(//5x,a/)')
     :     'Specify the time interval between history data dumps (s):'

      read (5,*) thisdmp

      write(6,'(/5x,a,f10.3,a)') 'Input was ',thisdmp,' seconds '

      write(6,'(2(//5x,a),4(/5x,a),//5x,a/)')
     :     'Now choose a format for the history data dump:',
     :     '  1  for unformatted binary data dump;',
     :     '  2  for formatted ASCII data dump;',
     :     '  3  for NCSA HDF format data dump;',
     :     '  4  for packed binary data dump;',
     :     '  5  for Savi3D data dump.',
     :     'Select a number (1 suggested):'

      read (5,*) hdmpfmt
      write(6,'(/5x,a,i4)') 'Input was ',hdmpfmt

      IF( hdmpfmt.lt.0 .or. hdmpfmt.gt.5) THEN 

        write(6,'(/5x,a,/5x,a/)')
     :       'Invalid input values for history dump format option,',
     :       'Job stopped in INITPARA.'
        STOP

      ENDIF

      nhisdmp =  nint(thisdmp/dtbig)
      nchdmp = 57

      IF( hdmpfmt.eq.1) THEN

        write(6,'(/5x,a)')
     :       'History dump files '//runname(1:6)//'.bintttttt'

      ELSEIF( hdmpfmt.eq.2) THEN

        write(6,'(/5x,a)')
     :       'History dump files '//runname(1:6)//'.asctttttt'

      ELSEIF( hdmpfmt.eq.3) THEN

        write(6,'(/5x,a)')
     :       'History dump files '//runname(1:6)//'.hdftttttt'

      ELSEIF( hdmpfmt.eq.4) THEN

        write(6,'(/5x,a)')
     :       'History dump files '//runname(1:6)//'.paktttttt'

      ELSEIF( hdmpfmt.eq.5) THEN

        write(6,'(/5X,a)')
     :       'History dump files '//runname(1:6)//'.svi'

      ENDIF

      write(6,'(5x,a,i6,a,/5x,a,a,/5x,a,i2,a)')
     :     'will be produced every ', nhisdmp,' time steps',
     :     'where tttttt (if any) stands for the time of the data ',
     :     'in seconds.',
     :     'FORTRAN channel ',nchdmp ,' is used for history data dump.'
c
C#######################################################################
c
c     Set control parameters for select filed output: 
c
C#######################################################################
c
      write(6,'(2(/5x,a)/)')
     :     'Do you want to dump out model grid arrays along with the',
     :     'time dependent fields? (select 0 or 1)'
 
      read(5,*) grdout
      write(6,'(/5x,a,i4/)') 'Input was ', grdout 

      write(6,'(2(/5x,a)/)')
     :     'Do you want to dump out the base state arrays',
     :     'along with the time dependent fields? (select 0 or 1)'
 
      read(5,*) basout
      write(6,'(/5x,a,i4/)') 'Input was ', basout 

      write(6,'(2(/5x,a)/)')
     :     'Do you want to dump out u, v, w, ptprt and pprt arrays?',
     :     '(select 0 or 1)'

      read(5,*) varout
      write(6,'(/5x,a,i4/)') 'Input was ', varout 

      write(6,'(2(/5x,a)/)')
     :     'Do you want to dump out any water and/or ice arrays? ',
     :     '(select 0 or 1)'

      read(5,*) mstout
      write(6,'(/5x,a,i4/)') 'Input was ', mstout 

      write(6,'(/5x,a/)')
     :     'Do you want to dump out any ice arrays? (select 0 or 1)'

      read(5,*) iceout

      write(6,'(/5x,a,i4/)') 'Input was ', iceout 

      write(6,'(2(/5x,a)/)')
     :     'Do you want to dump out turbulent mixing coefficient KM?', 
     :     '(select 0 or 1)'

      read(5,*) trbout
      write(6,'(/5x,a,i4/)') 'Input was ', trbout 
c
C#######################################################################
c
c     Restart data dump: 
c
C#######################################################################
c
      write(6,'(//5x,a/)')
     :     'Specify the time interval between restart data dumps (s):'

      read (5,*) trstout
      write(6,'(/5x,a,f10.3,a)') 'Input was ',trstout,' seconds '

      nrstout =  nint(trstout/dtbig)
      write(6,'(/5x,a,/5x,a,i6,a,/5x,a,/5x,a,i2,a)')
     :     'Restart data files '//runname(1:6)//'.rsttttttt',
     :     'will be produced every ', nrstout,' time steps',
     :     'where tttttt stands for the time of the data in seconds.',
     :     'FORTRAN I/O unit ',rstount,
     :     ' will be used for restart data output.'

c
C#######################################################################
c
c     Maximum and minimum statistics calculations: 
c
C#######################################################################
c

      write(6,'(//5x,a,a/)')
     :     'Specify the time interval between max./min. ',
     :     'statistics calculations (s):'

      read (5,*) tmaxmin
      write(6,'(/5x,a,f10.3,a)') 'Input was ',tmaxmin,' seconds '

      nmaxmin = nint(tmaxmin/dtbig)
      nchmax  = 36

      write(6,'(/5x,a,i6,a,/5x,a,/5x,a,i4/)')
     :     'Max. min. statistics are done every ', nmaxmin,
     :     ' time steps.',
     :     'and results written into file '//runname(1:6)//'.maxmin ',
     :     'using FORTRAN I/O channel ',nchmax

c
C#######################################################################
c
c     Energy/ptprt variance statistics calculations:
c
C#######################################################################
c

      write(6,'(//5x,a,a/)')
     :     'Specify the time interval between energy ',
     :     'statistics calculations (s):'

      read (5,*) tenergy 
      write(6,'(/5x,a,f10.3,a)') 'Input was ',tenergy,' seconds '

      nenergy = nint(tenergy/dtbig)
      ncheng  = 37

      write(6,'(/5x,a,i6,a,/5x,a,/5x,a,i4/)')
     :     'Energy statistics are done every ', nenergy,' time steps.',
     :     'and results written into file '//runname(1:6)//'.eng ',
     :     'using FORTRAN I/O channel ',ncheng

C#######################################################################
c
c     Debug information print controls: 
c
C#######################################################################
c
      write(6,'(3(/5x,a)/)')
     :     'You can print out some information when each subroutine',
     :     'is called. If you wish to do so, input value 1, otherwise',
     :     'input 0 (suggested value):'
      
      read (5,*) subtrac
      write(6,'(/5x,a,i4)') 'Input was ', subtrac
c
C#######################################################################
c
c     Print out the input parameters:
c
C#######################################################################
c

      write(6,'(///2x,a//)')
     :    'PRINT OUT OF MODEL PARAMETERS FOR THIS RUN:'

      write(6,'(/5x,a/)')
     :     'Comment lines for this run:'

      DO 22 i = 1, nocmnt
        write(6,'(5x,a)') cmnt(i)
22    CONTINUE
      write(6,*) ' '

      write(6,102) 'nx      =',nx    
      write(6,102) 'ny      =',ny    
      write(6,102) 'nz      =',nz    
      write(6,101) 'dx      =',dx    
      write(6,101) 'dy      =',dy    
      write(6,101) 'dz      =',dz    
      write(6,101) 'xl      =',xl    
      write(6,101) 'yl      =',yl    
      write(6,101) 'zh      =',zh    
      write(6,101) 'dtbig   =',dtbig 
      write(6,102) 'csopt   =',csopt 
      write(6,101) 'csfactr =',csfactr
      write(6,101) 'csound  =',csound
      write(6,101) 'dtsml   =',dtsml 

      write(6,101) 'tstop   =',tstop 

      write(6,102) 'initopt =',initopt
   
      write(6,101) 'ptpert0 =',ptpert0
      write(6,101) 'pt0radx =',pt0radx
      write(6,101) 'pt0rady =',pt0rady
      write(6,101) 'pt0radz =',pt0radz
      write(6,101) 'pt0ctrx =',pt0ctrx
      write(6,101) 'pt0ctry =',pt0ctry
      write(6,101) 'pt0ctrz =',pt0ctrz

      write(6,102) 'ebc     =',ebc 
      write(6,102) 'wbc     =',wbc 
      write(6,102) 'nbc     =',nbc 
      write(6,102) 'sbc     =',sbc 
      write(6,102) 'tbc     =',tbc 
      write(6,102) 'bbc     =',bbc 
      write(6,102) 'rbcopt  =',rbcopt

      write(6,102) 'mphyopt =',mphyopt 
      write(6,102) 'moist   =',moist 
      write(6,102) 'ice     =',ice 

      write(6,102) 'tmixopt =',tmixopt
      write(6,101) 'tmixcst =',tmixcst 
      write(6,101) 'prantl  =',prantl  

      write(6,102) 'kcmix2  =',kcmix2
      write(6,101) 'cfcmh2  =',cfcmh2
      write(6,101) 'cfcmv2  =',cfcmv2

      write(6,102) 'kcmix4  =',kcmix4
      write(6,103) 'cfcmh4  =',cfcmh4
      write(6,103) 'cfcmv4  =',cfcmv4

      write(6,102) 'divdmp  =',divdmp
      write(6,101) 'cdvdmp  =',cdvdmp

      write(6,102) 'raydmp  =',raydmp  
      write(6,101) 'zbrdmp  =',zbrdmp  
      write(6,101) 'ztrdmp  =',ztrdmp  
      write(6,101) 'cfrdmp  =',cfrdmp  

      write(6,102) 'sfcphy  =',sfcphy
      write(6,101) 'cdm     =',cdm   
      write(6,101) 'cdh     =',cdh   
      write(6,101) 'cdq     =',cdq   

      write(6,101) 'flteps  =',flteps  

      write(6,102) 'coriopt =',coriopt 
      write(6,101) 'latitud =',latitud 

      write(6,101) 'umove   =',umove   
      write(6,101) 'vmove   =',vmove   

      write(6,101) 'tfmtprt =',tfmtprt  
      write(6,101) 'thisdmp =',thisdmp
      write(6,102) 'hdmpfmt =',hdmpfmt 
      write(6,101) 'trstdmp =',trstout
      write(6,101) 'tmaxmin =',tmaxmin
      write(6,101) 'tenergy =',tenergy

      write(6,102) 'subtrac =',subtrac

      write(6,102) 'grdout  =',grdout 
      write(6,102) 'basout  =',basout 
      write(6,102) 'varout  =',varout 
      write(6,102) 'mstout  =',mstout 
      write(6,102) 'iceout  =',iceout 
      write(6,102) 'trbout  =',trbout 

101   format(5x,a,f12.3)
102   format(5x,a,i8)
103   format(5x,a,e13.3)

      write(6,*) ' '

c
C#######################################################################
c
c     Now let's write out a log file of this job which can be used 
c     as the input file for replicating this run. 
c
c     First get a name for the log file:
c
C#######################################################################
c
      CALL gtlogfn(runname, logfn, llogfn )

      logfunt = 40

      open (unit=logfunt, file=logfn(1:llogfn),status='new',
     :            iostat=istat)

      IF(istat.ne.0) THEN

        write(6,'(/5x,a)') 'Error in opening log file ', logfn
        write(6,'(5x,a/)') 'Job stopped in subroutine INITPARA.'
        STOP

      ELSE

        write(6,'(/5x,a,i3,a)')
     :       'FORTRAN unit ',logfunt,
     :       ' will be used for log file output.'

      ENDIF

      write(logfunt,'( i12 ,a)') nocmnt,
     :      '   ! nocmnt, number of comment lines' 

      DO 32 i = 1, nocmnt
        write(logfunt,'(a,a,a)') '''',cmnt(i),''''
32    CONTINUE

      write(logfunt,'(a,a,a)') '''',runname,'''   ! Name of this run'

      write(logfunt,'( i12 ,a)') runmod ,' ! runmod ' 
      write(logfunt,'( i12 ,a)') initopt,' ! initopt' 
      write(logfunt,'( i12 ,a)') pt0opt ,' ! pt0opt ' 
      write(logfunt,'(f12.3,a)') ptpert0,' ! ptpert0' 
      write(logfunt,'(e13.3,a)') pt0radx,' ! pt0radx' 
      write(logfunt,'(e13.3,a)') pt0rady,' ! pt0rady' 
      write(logfunt,'(e13.3,a)') pt0radz,' ! pt0radz' 
      write(logfunt,'(e13.3,a)') pt0ctrx,' ! pt0ctrx' 
      write(logfunt,'(e13.3,a)') pt0ctry,' ! pt0ctry' 
      write(logfunt,'(e13.3,a)') pt0ctrz,' ! pt0ctrz' 

      lenstr = 80
      CALL strlnth( rstinf, lenstr )
      write(logfunt,'(a,a,a)') '''',rstinf (1:lenstr),
     :     '''   ! Restart file'

      lenstr = 80
      CALL strlnth( inifile,lenstr )
      write(logfunt,'(a,a,a)') '''',inifile(1:lenstr),
     :     '''   ! Input file of initial data'

      lenstr = 80
      CALL strlnth( sndfile,lenstr )
      write(logfunt,'(a,a,a)') '''',sndfile(1:lenstr),
     :     '''   ! Sounding file'

      write(logfunt,'(f12.3,a)') dx      ,'  ! dx     '
      write(logfunt,'(f12.3,a)') dy      ,'  ! dy     '
      write(logfunt,'(f12.3,a)') dz      ,'  ! dz     '
      write(logfunt,'(f12.3,a)') dtbig   ,'  ! dtbig  '
      write(logfunt,'( i12 ,a)') csopt   ,'  ! csopt  '
      write(logfunt,'(f12.3,a)') csfactr ,'  ! csfactr'
      write(logfunt,'(f12.3,a)') csound  ,'  ! csound '
      write(logfunt,'(f12.3,a)') dtsml   ,'  ! dtsml  '
      write(logfunt,'( i12 ,a)') wbc     ,'  ! wbc    '
      write(logfunt,'( i12 ,a)') ebc     ,'  ! ebc    '
      write(logfunt,'( i12 ,a)') sbc     ,'  ! sbc    '
      write(logfunt,'( i12 ,a)') nbc     ,'  ! nbc    '
      write(logfunt,'( i12 ,a)') tbc     ,'  ! tbc    '
      write(logfunt,'( i12 ,a)') bbc     ,'  ! bbc    '
      write(logfunt,'( i12 ,a)') rbcopt  ,'  ! rbcopt '
      write(logfunt,'( i12 ,a)') mphyopt ,'  ! mphyopt'
      write(logfunt,'( i12 ,a)') moist   ,'  ! moist  '
      write(logfunt,'( i12 ,a)') ice     ,'  ! ice    '
      write(logfunt,'( i12 ,a)') tmixopt ,'  ! tmixopt'
      write(logfunt,'(f12.3,a)') tmixcst ,'  ! tmixcst'
      write(logfunt,'(f12.3,a)') prantl  ,'  ! prantl '
      write(logfunt,'( i12 ,a)') kcmix2  ,'  ! kcmix2 '
      write(logfunt,'(e13.3,a)') cfcmh2  , ' ! cfcmh2 '
      write(logfunt,'(e13.3,a)') cfcmv2  , ' ! cfcmv2 '
      write(logfunt,'( i12 ,a)') kcmix4  ,'  ! kcmix4 '
      write(logfunt,'(e13.3,a)') cfcmh4  , ' ! cfcmh4 '
      write(logfunt,'(e13.3,a)') cfcmv4  , ' ! cfcmv4 '
      write(logfunt,'( i12 ,a)') divdmp  ,'  ! divdmp '
      write(logfunt,'(f12.4,a)') alpha   ,'  ! alpha  '
      write(logfunt,'( i12 ,a)') raydmp  ,'  ! raydmp '
      write(logfunt,'(f12.4,a)') cfrdmp  ,'  ! cfrdmp '
      write(logfunt,'(f12.3,a)') zbrdmp  ,'  ! zbrdmp '
      write(logfunt,'(f12.3,a)') ztrdmp  ,'  ! ztrdmp '
      write(logfunt,'( i12 ,a)') sfcphy  ,'  ! sfcphy '
      write(logfunt,'(f12.3,a)') cdm     ,'  ! cdm    '
      write(logfunt,'(f12.3,a)') cdh     ,'  ! cdh    '
      write(logfunt,'(f12.3,a)') cdq     ,'  ! cdq    '
      write(logfunt,'( i12 ,a)') coriopt ,'  ! coriopt'
      write(logfunt,'(f12.3,a)') latitud ,'  ! latitud'
      write(logfunt,'(f12.4,a)') flteps  ,'  ! flteps '
      write(logfunt,'(f12.4,a)') umove   ,'  ! umove  '
      write(logfunt,'(f12.4,a)') vmove   ,'  ! vmove  '
      write(logfunt,'(f12.2,a)') tstop   ,'  ! tstop  '
c      write(logfunt,'(a,a)') ''''//dirname(1:ldirnam)//'''',
c     :      '  ! output directory name'
      write(logfunt,'(f12.3,a)') tfmtprt ,'  ! tfmtprt'
      write(logfunt,'(f12.3,a)') thisdmp ,'  ! thisdmp'
      write(logfunt,'( i12 ,a)') hdmpfmt ,'  ! hdmpfmt'
      write(logfunt,'( i12 ,a)') grdout  ,'  ! grdout '
      write(logfunt,'( i12 ,a)') basout  ,'  ! basout '
      write(logfunt,'( i12 ,a)') varout  ,'  ! varout '
      write(logfunt,'( i12 ,a)') mstout  ,'  ! mstout '
      write(logfunt,'( i12 ,a)') iceout  ,'  ! iceout '
      write(logfunt,'( i12 ,a)') trbout  ,'  ! trbout '
      write(logfunt,'(f12.3,a)') trstout ,'  ! trstout'
      write(logfunt,'(f12.3,a)') tmaxmin ,'  ! tmaxmin'
      write(logfunt,'(f12.3,a)') tenergy ,'  ! tenergy'
      write(logfunt,'( i12 ,a)') subtrac ,'  ! subtrac'

      close(unit=logfunt)

      write(6,'(/5x,a,a,a/)')
     :     'Log file ',logfn(1:llogfn),' was produced for this job.'
      
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INIGRD                     ######
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 INIGRD(nx,ny,nz, 
     :                  x,y,z,zp,hterain,j1,j2,j3)

c
c#######################################################################
c
c     PURPOSE:
c
c     Initialize of the model GRiD variables.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further 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     OUTPUT:
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

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz     ! The number of grid points in 3 directions
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 J1(nx,ny,nz)    ! coordinate transform Jacobian defined as
                           ! - del( zp )/del( x ).
      real J2(nx,ny,nz)    ! coordinate transform Jacobian defined as
                           ! - del( zp )/del( y ).
      real J3(nx,ny,nz)    ! coordinate transform Jacobian defined as
                           ! del( zp )/del( z ) = sqrt( G ).

      real hterain(nx,ny)  ! the height of terrain.
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c 
      integer i,j,k
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc' 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Define the model grid. For now, define a uniform grid.       
c
c#######################################################################
c
      DO 100 i=1,nx
        x(i) = (i-2) * dx
100   CONTINUE

      DO 200 j=1,ny
        y(j) = (j-2) * dy
200   CONTINUE

      DO 300 k=1,nz
        z(k) = (k-2) * dz
300   CONTINUE

      DO 400 j=1,ny
      DO 400 i=1,nx
        hterain(i,j) = 0.0 
400   CONTINUE
c
c#######################################################################
c
c     Physical height of computational grid defined as
c
c     Zp=Zt*(z-hterain)/(Zt-hterain) where Zt the height of the model top.
c
c#######################################################################
c
      DO 500 k=1,nz
      DO 500 j=1,ny
      DO 500 i=1,nx
        zp(i,j,k) = z(nz-1) * (z(k)-hterain(i,j))/(z(nz-1)-hterain(i,j)) 
500   CONTINUE
c
c#######################################################################
c
c     Calculate transformation Jacobian J1 defined as -del(zp)/del(x)
c
c#######################################################################
c
      DO 600 k=1,nz
      DO 600 j=1,ny-1
      DO 600 i=2,nx-1
        j1(i,j,k) = -2 * (zp(i,j,k)-zp(i-1,j,k)) / (x(i+1)-x(i-1)) 
600   CONTINUE
c
c#######################################################################
c
c     X - boundary conditions of j1
c
c#######################################################################
c
      DO 610 k=1,nz
      DO 610 j=1,ny-1
        j1( 1,j,k)=j1( 2  ,j,k)
        j1(nx,j,k)=j1(nx-1,j,k)
610   CONTINUE

c
c#######################################################################
c
c     Calculate transformation Jacobian J2 defined as -del(zp)/del(y)
c
c#######################################################################
c
      DO 700 k=1,nz
      DO 700 i=1,nx-1
      DO 700 j=2,ny-1
        j2(i,j,k) = -2 * (zp(i,j,k)-zp(i,j-1,k)) / (y(j+1)-y(j-1)) 
700   CONTINUE
c
c#######################################################################
c
c     Y - boundary conditions of j2
c
c#######################################################################
c
      DO 710 k=1,nz
      DO 710 i=1,nx-1
        j2(i, 1,k)=j2(i, 2  ,k)
        j2(i,ny,k)=j2(i,ny-1,k)
710   CONTINUE

c
c#######################################################################
c
c     Calculate transformation Jacobian J3 defined as del(zp)/del(z)
c
c#######################################################################
c
      DO 800 k=1,nz-1
      DO 800 j=1,ny-1
      DO 800 i=1,nx-1
        j3(i,j,k) = (zp(i,j,k+1)-zp(i,j,k))/(z(k+1)-z(k)) 
800   CONTINUE

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INITVAR                    ######
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 INITVAR(nx,ny,nz, 
     :           x,y,z,zp,hterain,j1,j2,j3,
     :           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,
     :           tem1)

c
c#######################################################################
c
c     PURPOSE:
c
c     INITialize the model array VARiables.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/17/1992.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
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       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     OUTPUT:
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 
c                at all time 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       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     WORK ARRAYS:
c
c       tem1     Temporary work array.
c
c      (These are arrays that are defined in this routine and possibly
c       passed into routines called by this one. They share the same
c       memory locations as certain temporary arrays in the calling
c       routine, their memory location could therefore be reused in
c       the calling routine for other purposes on exit of this one.)
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)
c
      integer nx,ny,nz          ! The number of grid points in 3 directions
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
      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)

      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
      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
      real tem1  (nx,ny,nz)     ! Temporary array.
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c

c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'bndry.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external inibase
      external initdvr 
      external initbdt 
      external rstin 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     If initopt = 1, initialize the model fields using intialization 
c                     routines. Typically from an analytical 
c                     definition of initial perturbations. 
c
c     If initopt = 2, initialize the model fields from a restart file 
c                     produced by a previous model run.
c
c     If initopt = 3, initialize the model fields by reading in an 
c                     external data file.  
c 
C#######################################################################
c

      IF( initopt .eq. 1 ) THEN 

c
C#######################################################################
c
c     Set up the base state of atmosphere.
c
C#######################################################################
c      
        CALL inibase(nx,ny,nz, ubar,vbar,ptbar,pbar,rhobar,qvbar)
c
C#######################################################################
c
c     Initialize time dependent model variables.
c
C#######################################################################
c      
        CALL initdvr(nx,ny,nz,
     :       ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :       x,y,z,zp,hterain, j1,j2,j3, 
     :       u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :       tem1)

        tstart = 0.0 
c
c
c#######################################################################
c
c     Set the time tendencies of variables at the lateral boundaries,
c     for options 4 and 5 of lateral boundary condition, for the first 
c     time step of integration.
c
c     For the later time steps, these tendencies are set in the previous
c     time steps.
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) .or.
     :      (ebc.eq.5 .or. wbc.eq.5 .or. nbc.eq.5 .or. sbc.eq.5) ) THEN

          CALL initbdt(nx,ny,nz,
     :                 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)


        ENDIF

      ELSEIF( initopt .eq. 2) THEN   ! restart run
c
C#######################################################################
c
c     To read in restart data from a restart file to initialize fields
c     u,v,w,prprt,pprt,qv,qc,qr,qi,qs,qh at time tpast and tpresent,
c     the base state variables ubar,vbar,ptbar,pbar,rhobar,qvbar,
c     and the time tendencies of variables at the lateral boundries.
c
c     Fields at tfuture are set to equal to those at tpresent.
c
c     This subroutine also sets the value of tstart.
c
C#######################################################################
c
c
        CALL rstin(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)

      ELSEIF( initopt .eq. 3) THEN  ! External data input.

        write(6,'(1x,a,/1x,a)')
     :       'Third option of initialization (external input) not ',
     :       'implemented yet, please respecify the option.
     :        Job stopped.'

        STOP

      ELSE

        write(6,'(1x,a,i3,a)')
     :       'Invalid option for intialization, INITOPT was ',initopt,
     :       ' Job stopped.'

        STOP

      ENDIF
c
C#######################################################################
c
c     End of initialization stage and output the initial fields
c
C#######################################################################
c
      RETURN 
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INITDVR                    ######
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 INITDVR(nx,ny,nz,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar, 
     :           x,y,z,zp,hterain, j1,j2,j3,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           tem1)

c
c#######################################################################
c
c     PURPOSE:
c
C     INItialize the model Time Dependent VaRiables.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/1991.
c
c     MODIFICATION HISTORY:
c
c     5/25/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further 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       ubar     Base state x-velocity component (m/s)
c       vbar     Base state y-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 computational space (m)
c       y        y-coordinate of grid points in computational 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     OUTPUT:
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        z-component of velocity at all time 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     WORK ARRAYS:
c
c       tem1     Temporary work array.
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)
c
      integer nx,ny,nz          ! The number of grid points in 3 directions
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
      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
      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 tem1(nx,ny,nz)       ! Temporary array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      real xs, ys, zs
      real beta, pi2
      integer i,j,k, n
      integer iseed

c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
      include 'phycst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external ranary
      external satmrpt
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Specify the initial potential temperature perturbation.
c
c#######################################################################
c

      IF( pt0opt.eq.1 ) THEN  ! Bubble shaped perturbation

c
c#######################################################################
c
c     Set the potential temperature perturbation inside a bubble volume.
c
c#######################################################################
c

        pi2=atan(1.0)*4.0/2.0
        pi2 = 1.5707963267949

        IF( ptpert0.eq. 0.0) THEN 

          DO 55 k= 1,nz-1
          DO 55 j= 1,ny-1
          DO 55 i= 1,nx-1
            ptprt(i,j,k,1) = 0.0
55        CONTINUE

        ELSE

          DO 50 k= 1,nz-1
          DO 50 j= 1,ny-1
          DO 50 i= 1,nx-1
 
            xs = (x(i)+x(i+1))*0.5 
            ys = (y(j)+y(j+1))*0.5
            zs = (zp(i,j,k)+zp(i,j,k+1))*0.5

            IF( pt0rady.lt.0.0 ) THEN   ! 2-d bubble in x-z plane.
 
              beta = sqrt( ((xs-pt0ctrx)/pt0radx)**2 
     :                   + ((zs-pt0ctrz)/pt0radz)**2 )
 
            ELSEIF( pt0radx.lt.0.0 ) THEN ! 2-d bubble in y-z plane.
 
              beta = sqrt( ((ys-pt0ctry)/pt0rady)**2 
     :                   + ((zs-pt0ctrz)/pt0radz)**2 )
 
            ELSE                        ! 3-d bubble

              beta = sqrt( ((xs-pt0ctrx)/pt0radx)**2 +
     :                     ((ys-pt0ctry)/pt0rady)**2 +
     :                     ((zs-pt0ctrz)/pt0radz)**2 )
            ENDIF

            IF(beta.ge.1.0) THEN
              ptprt(i,j,k,1) = 0.0
            ELSE
              ptprt(i,j,k,1) = ptpert0*(cos(pi2*beta)**2)
            ENDIF

50        CONTINUE

        ENDIF

      ELSE
c
c#######################################################################
c
c     Set the potential temperature perturbation as a random function,
c     which ensures the horizontal average of the perturbation in a 
c     specified domain to be zero.
c
c#######################################################################
c
        iseed = -100 

        DO 60 k= 1,nz-1

          CALL ranary(nx,ny,1,nx-1,1,ny-1,iseed,ptpert0,ptprt(1,1,k,1)) 

60      CONTINUE

      ENDIF

      DO 100 n=1,3
      DO 100 k=1,nz-1
      DO 100 j=1,ny-1
      DO 100 i=1,nx-1
        ptprt(i,j,k,n)=ptprt(i,j,k,1)
        pprt(i,j,k,n)=0.0
100   CONTINUE

      DO 200 n=1,3
      DO 200 k=1,nz-1
      DO 200 j=1,ny-1
      DO 200 i=1,nx
        u(i,j,k,n)=ubar(i,j,k)
200   CONTINUE

      DO 300 n=1,3
      DO 300 k=1,nz-1
      DO 300 j=1,ny
      DO 300 i=1,nx-1
        v(i,j,k,n)=vbar(i,j,k)
300   CONTINUE

      DO 400 n=1,3
      DO 400 k=1,nz
      DO 400 j=1,ny-1
      DO 400 i=1,nx-1
        w(i,j,k,n)=0.0
400   CONTINUE


c      CALL satmrpt(nx,ny,nz,pbar, ptbar, tem1)


      DO 500 n=1,3
      DO 500 k=1,nz-1
      DO 500 j=1,ny-1
      DO 500 i=1,nx-1
        qv(i,j,k,n)= qvbar(i,j,k)
500   CONTINUE

      DO 600 n=1,3
      DO 600 k=1,nz-1
      DO 600 j=1,ny-1
      DO 600 i=1,nx-1
        qc(i,j,k,n)=0.0
600   CONTINUE

      DO 700 n=1,3
      DO 700 k=1,nz-1
      DO 700 j=1,ny-1
      DO 700 i=1,nx-1
        qr(i,j,k,n)= 0.0
700   CONTINUE

      DO 710 n=1,3
      DO 710 k=1,nz-1
      DO 710 j=1,ny-1
      DO 710 i=1,nx-1
        qi(i,j,k,n)= 0.0
710   CONTINUE

      DO 720 n=1,3
      DO 720 k=1,nz-1
      DO 720 j=1,ny-1
      DO 720 i=1,nx-1
        qs(i,j,k,n)= 0.0
720   CONTINUE

      DO 730 n=1,3
      DO 730 k=1,nz-1
      DO 730 j=1,ny-1
      DO 730 i=1,nx-1
        qh(i,j,k,n)= 0.0
730   CONTINUE

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INIT0                      ######
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 INIT0(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar, km,
     :           temxy1,temxy2,temxy3,temxy4,
     :           tem1,tem2,tem3,tem4,tem5,tem6,tem7,
     :           tem8,tem9,tem10,tem11,tem12,tem13)

c
c#######################################################################
c
c     PURPOSE:
c
C     Initializes all defined arrays with zeros (0).
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further 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     OUTPUT:
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        z-component of velocity at all time 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       ubar     Base state x-velocity component (m/s)
c       vbar     Base state y-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 computational space (m)
c       y        y-coordinate of grid points in computational 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       km       Turbulent mixing coef. Filled with zeros.
c
c       temxy1   2-d temporary array
c       temxy2   2-d temporary array
c       temxy3   2-d temporary array
c       temxy4   2-d temporary array
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#######################################################################
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)
c
      integer nx,ny,nz          ! The 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 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 km    (nx,ny,nz)     ! The turbulent mixing coefficient for 
                                ! momentum. ( m**2/s )      
  
      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.
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer nxyz
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external flzero
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      nxyz = nx*ny*nz

      CALL flzero(u,nxyz*3)
      CALL flzero(v,nxyz*3)
      CALL flzero(w,nxyz*3)
      CALL flzero(ptprt,nxyz*3)
      CALL flzero(pprt,nxyz*3)
      CALL flzero(qv,nxyz*3)
      CALL flzero(qc,nxyz*3)
      CALL flzero(qr,nxyz*3)
      CALL flzero(qi,nxyz*3)
      CALL flzero(qs,nxyz*3)
      CALL flzero(qh,nxyz*3)

      CALL flzero(ubar,nxyz)
      CALL flzero(vbar,nxyz)
      CALL flzero(ptbar,nxyz)
      CALL flzero(pbar,nxyz)
      CALL flzero(rhobar,nxyz)
      CALL flzero(qvbar,nxyz)
      
      CALL flzero(km,nxyz)

      CALL flzero(temxy1, nx*ny)
      CALL flzero(temxy2, nx*ny)
      CALL flzero(temxy3, nx*ny)
      CALL flzero(temxy4, nx*ny)

      CALL flzero(tem1 ,nxyz)
      CALL flzero(tem2 ,nxyz)
      CALL flzero(tem3 ,nxyz)
      CALL flzero(tem4 ,nxyz)
      CALL flzero(tem5 ,nxyz)
      CALL flzero(tem6 ,nxyz)
      CALL flzero(tem7 ,nxyz)
      CALL flzero(tem8 ,nxyz)
      CALL flzero(tem9 ,nxyz)
      CALL flzero(tem10,nxyz)
      CALL flzero(tem11,nxyz)
      CALL flzero(tem12,nxyz)
      CALL flzero(tem13,nxyz)

      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE FLZERO                     ######
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 FLZERO(a,n)

c
c#######################################################################
c
c     PURPOSE:
c
c     Fill an entire array with zeros. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     10/10/1991.
c
c     MODIFICATION HISTORY:
c
c     5/02/92 (M. Xue)
c     Added full documentation.
C
c     5/03/92 (M. Xue)
c     Further documentation. 
C
c#######################################################################
c
c     INPUT :
c
c       n        The dimension of 1-D array a.
c
c     OUTPUT:
c
c       a        1-D array to be filled with zeros.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, n

      real a(n)

c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      DO 100 i=1,n
        a(i)=0.0
100   CONTINUE

      RETURN
      END
