c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE INITOUT                    ######
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 INITOUT(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,km,
     :           x,y,z,zp,hterain, j1,j2,j3,
     :           tem1,tem2,tem3,tem4,tem5, svitem)
c
c#######################################################################
c
c     PURPOSE:
c
c     Handle the model data output at the initial time.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/13/91.
c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
c
c     7/13/92 (K. Brewster)
c     Added comment line variables to arg list to be passed to the
c     data dumping routines.
c
c     7/20/92 (M. Xue)
c     Added call to energy and base state dump routines.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        Vertical component of Cartesian velocity (m/s)
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c       qv       Water vapor specific humidity (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
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       km       Turbulent mixing coefficient (m**2/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       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       None.
c
c     WORK ARRAYS:
c
c       tem1     work array.
c       tem2     work array.
c       tem3     work array.
c       tem4     work array.
c       tem5     work array.
c       svitem   work array.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

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

      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
      real v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
      real w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
      real ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
c
      real qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
      real qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)
c
      real 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 )
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 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 svitem(3*nx*ny*nz)   ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer istat
      integer tlevel            ! Time level at which data are printed.
      character basdmpfn*80
      integer lbasdmpf
      integer grdbas, stagrd
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external maxmin
      external energy
      external basprt
      external fmtprt
      external gtbasfn
      external gtdmpfn
      external dtadump 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
c     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      tlevel = tpresent
c
c#######################################################################
c
c     Max./min. statistics calculation and printing of initial fields
c
c#######################################################################
c
      CALL maxmin(nx,ny,nz, tlevel, u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh)
c
c#######################################################################
c
c     Energy statistics calculation and printing at the initial time
c
c#######################################################################
c
      CALL energy(nx,ny,nz,
     :            u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,rhobar,
     :            x,y,z,zp,hterain, j3,
     :            tem1,tem2,tem3,tem4,tem5)
c
c#######################################################################
c
c     Formatted printing of base state variables
c
c#######################################################################
c

      CALL basprt(nx,ny,nz,ubar,vbar,ptbar,pbar,rhobar,qvbar)

      write(6,'(/'' The initial fields:''/)')
c
c#######################################################################
c
c     Formatted printing of time dependent fields
c
c#######################################################################
c

      CALL fmtprt(nx,ny,nz, tlevel,
     :            u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :            x,y,z,zp,hterain, j1,j2,j3, tem1)

c
c#######################################################################
c
c     Check to see if any history data dump is to be produced.
c
c#######################################################################
c

      IF( hdmpfmt.eq.0 ) THEN

        write(6,'(1x,a)')
     :       'History data dump option=0, no data is dumped.'
        GOTO 800

      ENDIF

c
c#######################################################################
c
c     Data dump of the model grid and base state arrays:
c
c#######################################################################
c
      IF( hdmpfmt.eq.5 ) GOTO 700 
c
C#######################################################################
c
c     First find a unique name basdmpfn(1:lbasdmpf) for the grid and
c     base state array dump file
c
C#######################################################################
c
      CALL gtbasfn(runname(1:6),dirname,ldirnam,hdmpfmt,
     :             basdmpfn,lbasdmpf)

      write(6,'(1x,a,a)')
     :     'Data dump of grid and base state arrays into file ',
     :      basdmpfn(1:lbasdmpf)


      grdbas = 1      ! Dump out grd and base state arrays only
      stagrd = 1      ! Arrays are staggered.

      CALL dtadump(nx,ny,nz,
     :             hdmpfmt,nchdmp,basdmpfn(1:lbasdmpf),grdbas,stagrd,
     :             u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :             ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar,
     :             x,y,z,zp,hterain,j1,j2,j3, tem2,tem3,svitem)

700   CONTINUE
c
c#######################################################################
c
c     Data dump of the model time dependent arrays at initial time:
c
c#######################################################################
c

c
C#######################################################################
c
c     Find a unique name hdmpfn(1:ldmpf) for history dump data set 
c     at time 'curtim'.
c
C#######################################################################
c
c
      CALL gtdmpfn(runname(1:6),dirname,ldirnam,curtim,hdmpfmt,
     :             hdmpfn, ldmpf )

      write(6,'(1x,a,a)') 'History data dump in file ',hdmpfn(1:ldmpf)


      grdbas = 0      ! Not just dump out grd and base state arrays only
      stagrd = 1      ! Arrays are staggered.

      CALL dtadump(nx,ny,nz,
     :             hdmpfmt,nchdmp,hdmpfn(1:ldmpf),grdbas,stagrd,
     :             u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :             ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar,
     :             x,y,z,zp,hterain, j1,j2,j3, tem2,tem3,svitem)


800   CONTINUE ! Entry if hdmpfmt=0

      RETURN

      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE OUTPUT                     ######
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 OUTPUT(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           udteb, udtwb, udtnb, udtsb,
     :           vdteb, vdtwb, vdtnb, vdtsb,
     :           wdteb, wdtwb, wdtnb, wdtsb,
     :           ptdteb,ptdtwb,ptdtnb,ptdtsb,
     :           pdteb ,pdtwb ,pdtnb ,pdtsb,
     :           qvdteb,qvdtwb,qvdtnb,qvdtsb,
     :           qcdteb,qcdtwb,qcdtnb,qcdtsb,
     :           qrdteb,qrdtwb,qrdtnb,qrdtsb,
     :           qidteb,qidtwb,qidtnb,qidtsb,
     :           qsdteb,qsdtwb,qsdtnb,qsdtsb,
     :           qhdteb,qhdtwb,qhdtnb,qhdtsb,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,km,
     :           x,y,z,zp,hterain, j1,j2,j3,
     :           tem1,tem2,tem3,tem4,tem5,svitem)
c
c#######################################################################
c
c     PURPOSE:
c
c     Handle the output of model data at a particular time.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/13/91.
c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
c
c     7/13/92 (K. Brewster)
c     Added comment line variables to arg list to be passed to the
c     data dumping routines.
C
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        Vertical component of Cartesian velocity (m/s)
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c       qv       Water vapor specific humidity (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (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       km       Turbulent mixing coefficient (m**2/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       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       None.
c
c     WORK ARRAYS:
c
c       tem1     work array.
c       tem2     work array.
c       tem3     work array.
c       tem4     work array.
c       tem5     work array.
c       svitem   work array.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nt           ! The no. of t-levels of t-dependent arrays.
      integer tpast        ! Index of time level for the past time.
      integer tpresent     ! Index of time level for the present time.
      integer tfuture      ! Index of time level for the future time.

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

      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
      real v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
      real w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
      real ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
c
      real qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
      real qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)
c
      real udteb (ny,nz)        ! T-tendency of u at e-boundary (m/s**2)
      real udtwb (ny,nz)        ! T-tendency of u at w-boundary (m/s**2)
      real udtnb (nx,nz)        ! T-tendency of u at n-boundary (m/s**2)
      real udtsb (nx,nz)        ! T-tendency of u at s-boundary (m/s**2)

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

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

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

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

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

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

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

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

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

      real qhdteb(ny,nz)        ! T-tendency of qh at e-boundary (kg/(kg*s))
      real qhdtwb(ny,nz)        ! T-tendency of qh at w-boundary (kg/(kg*s))
      real qhdtnb(nx,nz)        ! T-tendency of qh at n-boundary (kg/(kg*s))
      real qhdtsb(nx,nz)        ! T-tendency of qh at s-boundary (kg/(kg*s))
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 )
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 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 svitem(3*nx*ny*nz)   ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer istat
      integer tlevel            ! Time level at which data are printed.
      integer grdbas, stagrd
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'phycst.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external rstout
      external maxmin
      external energy
      external fmtprt
      external gtdmpfn
      external dtadump
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
C#######################################################################
c
c     Dump restart data files every nrstout number of time steps.
c
C#######################################################################
c
      IF( mod(nstep,nrstout).eq.0) THEN

         CALL rstout(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)

      ENDIF

c
C#######################################################################
c
c     Printing, diagnostics and history dump of data at time tpresent
c
C#######################################################################
c
      tlevel = tpresent
c
C#######################################################################
c
c     Calculate the  max and min values of variables and write them
c     into a separate file.
c
C#######################################################################
c
      IF( mod(nstep,nmaxmin).eq.0) THEN

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

      ENDIF
c
c#######################################################################
c
c     Energy statistics calculations and printing
c
c#######################################################################
c
      IF( mod(nstep,nenergy).eq.0) THEN

        CALL energy(nx,ny,nz,
     :              u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,rhobar,
     :              x,y,z,zp,hterain, j3,
     :              tem1,tem2,tem3,tem4,tem5)

      ENDIF
c
C#######################################################################
c
c     Print out formatted 2-d slices of 3-d arrays
c
C#######################################################################
c
      IF( mod(nstep,nfmtprt).eq.0) THEN

        write(6,'(1x,a,f13.3,a,i10)')
     :       'Print fields at time=',curtim,'(s) with nstep=',nstep

        CALL fmtprt(nx,ny,nz, tlevel,
     :       u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :       x,y,z,zp,hterain, j1,j2,j3, tem1)

      ENDIF

c
c#######################################################################
c
c     Produce history data dump every nhisdmp number of time steps.
c
c#######################################################################
c
      IF( mod(nstep,nhisdmp).eq.0 .and. hdmpfmt.ne.0 ) THEN

c
C#######################################################################
c
c     Find a unique name hdmpfn(1:ldmpf) for the history dump data 
c     set at time 'curtim'.
c
c     For the savi3D data dump case, the file name is specified only 
c     once in INITOUT.
c
C#######################################################################
c
c
        IF( hdmpfmt.ne.5) THEN

          CALL gtdmpfn(runname(1:6),dirname,ldirnam,curtim,hdmpfmt,
     :                 hdmpfn, ldmpf )

        ENDIF

        write(6,'(1x,a,a)')
     :       'History data dump in file ',hdmpfn(1:ldmpf)

        grdbas = 0      ! Not just dump out grd and base state arrays only
        stagrd = 1      ! Arrays are staggered.

        CALL dtadump(nx,ny,nz,
     :               hdmpfmt,nchdmp,hdmpfn(1:ldmpf),grdbas,stagrd,
     :               u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :               ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar,
     :               x,y,z,zp,hterain, j1,j2,j3, tem2,tem3,svitem)


      ENDIF

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE CHKSTAB                    ######
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 CHKSTAB(nx,ny,nz,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,km,
     :           x,y,z,zp,hterain, j1,j2,j3, tem1,tem2,tem3, svitem)
c
c#######################################################################
c
c     PURPOSE:
c
c     Check the stability of the time integration, and if unstable, dump
c     out model data and stop the model run.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/13/91.
c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
c
c     7/13/92 (K. Brewster)
c     Added comment line variables to arg list to be passed to the
c     data dumping routines.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        Vertical component of Cartesian velocity (m/s)
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c       qv       Water vapor specific humidity (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
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       km       Turbulent mixing coefficient (m**2/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       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       None.
c
c     WORK ARRAYS:
c
c       tem1     work array.
c       tem2     work array.
c       tem3     work array.
c       svitem   work array.
c
c#######################################################################
c

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

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

      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      real u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
      real v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
      real w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
      real ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
c
      real qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
      real qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)
c
      real 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 )
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 )

      real tem1  (nx,ny,nz)     ! Temporary
      real tem2  (nx,ny,nz)     ! Temporary
      real tem3  (nx,ny,nz)     ! Temporary
      real svitem(3*nx*ny*nz)   ! Temporary work array
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer istat
      integer i,j,k
      integer tlevel
      real    absmax,vellmt
      integer grdbas, stagrd
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c

      include 'phycst.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external maxmin
      external fmtprt
      external gtdmpfn
      external dtadump
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      tlevel = tpresent

c
C#######################################################################
c
c     Check for instabilty of model integration.
c
c     When instability is detected, stop the model run and dump out the
c     model data.
c
C#######################################################################
c
      absmax = abs(u(1,1,1,tlevel))
c
      DO 110 k=1,nz-1
      DO 110 j=1,ny-1
      DO 110 i=1,nx
        absmax = max(absmax,abs(u(i,j,k,tlevel)))
110   CONTINUE
c
      absmax = max(absmax,abs(v(1,1,1,tlevel)))

      DO 120 k=1,nz-1
      DO 120 i=1,nx-1
      DO 120 j=1,ny
        absmax = max(absmax,abs(v(i,j,k,tlevel)))
120   CONTINUE

      absmax = max(absmax,abs(w(1,1,1,tlevel)))

      DO 130 j=1,ny-1
      DO 130 i=1,nx-1
      DO 130 k=1,nz
        absmax = max(absmax,abs(w(i,j,k,tlevel)))
130   CONTINUE
c
c#######################################################################
c
c     Set an upper limit on maximum velocity component, above which
c     the integeration is regarded as unstable.
c
c#######################################################################
c
      vellmt = 100.0

      IF( absmax.gt.vellmt ) THEN

        write(6,'(/1x,a,f5.1,a,/1x,a,i6,/1x,a)')
     :  'Maximum velocity component exceeded ',vellmt,' m/s,',
     :  'Time integration is unstable, Job aborted at the end of step',
     :   nstep,' Fields at this time dumped out.'
c
c#######################################################################
c
c     Max./min. statistics calculation and printing of initial fields
c
c#######################################################################
c
        CALL maxmin(nx,ny,nz,tlevel,u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh)
c
c#######################################################################
c
c     Formatted printing of time dependent variables
c
c#######################################################################
c
        CALL fmtprt(nx,ny,nz, tlevel,
     :              u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :              x,y,z,zp,hterain, j1,j2,j3, tem1)
c
c
c#######################################################################
c
c     Check to see if any history data dump is to be produced.
c
c#######################################################################
c
        IF( hdmpfmt.eq.0 ) THEN

          write(6,'(1x,a)')
     :         'History data dump option was 0, no data is dumped.'
          GOTO 800

        ENDIF

c
C#######################################################################
c
c     Find a unique name hdmpfn(1:ldmpf) for the history dump data set 
c     at time 'curtim'.
c
c     For the savi3D data dump case, the file name is specified only 
c     once in INITOUT.
c
C#######################################################################
c
c
        IF( hdmpfmt.ne.5) THEN

          CALL gtdmpfn(runname(1:6),dirname,ldirnam,curtim,hdmpfmt,
     :                 hdmpfn, ldmpf )

        ENDIF

        write(6,'(1x,a,a)')
     :       'History data dump in file ',hdmpfn(1:ldmpf)

        grdbas = 0      ! Not just dump out grd and base state arrays only
        stagrd = 1      ! Arrays are staggered.

        CALL dtadump(nx,ny,nz,
     :               hdmpfmt,nchdmp,hdmpfn(1:ldmpf),grdbas,stagrd,
     :               u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :               ubar,vbar,tem1,ptbar,pbar,rhobar,qvbar,
     :               x,y,z,zp,hterain, j1,j2,j3, tem2,tem3,svitem)

c
c#######################################################################
c
c     If ARPS is dumping Savi3D GRAF format history files,
c     call grafclose to close the GRAF file before program stop.
c
c#######################################################################
c

        IF (hdmpfmt .eq. 5) CALL grafclose(grafh)

800     CONTINUE

        STOP

      ENDIF


      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE MAXMIN                     ######
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 MAXMIN(nx,ny,nz,tlevel,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate the maximum and minimum of the time dependent fields and
c     write them into a file named runname(1:6).maxmin.
c
c#######################################################################
c
c     AUTHOR: Ming Xue

c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       tlevel   The time level at which the data are printed.
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        Vertical component of Cartesian velocity (m/s)
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c       qv       Water vapor specific humidity (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
c
c     OUTPUT:
c
c       None.
c
c#######################################################################
c

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

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

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

      integer nx,ny,nz     ! Number of grid points in 3 directions.
      integer tlevel       ! Time level at which data are processed.
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
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer ncalls
      data ncalls /0/
      save ncalls

      character maxfn*80    ! Name of max./min. file
      integer   lmaxfn
      integer istat
      integer imax,jmax,kmax,imin,jmin,kmin
      real umax,vmax,wmax,ptmax,pmax,qvmax,qcmax,qrmax,qimax,qsmax,qhmax
      real umin,vmin,wmin,ptmin,pmin,qvmin,qcmin,qrmin,qimin,qsmin,qhmin
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external a3dmax
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      write(6,*) ' '

      CALL a3dmax
     :     (u(1,1,1,tlevel),1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,
     :     umax,umin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(/2(1x,a,f13.7,3(a,i3)))')
     :     'umin =',umin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'umax =',umax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (v(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,
     :     vmax,vmin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'vmin =',vmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'vmax =',vmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (w(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,
     :     wmax,wmin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'wmin =',wmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'wmax =',wmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (ptprt(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     ptmax,ptmin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'ptmin=',ptmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'ptmax=',ptmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (pprt(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     pmax,pmin, imax,jmax,kmax, imin,jmin,kmin)


      write(6,'(2(1x,a,f13.5,3(a,i3)))')
     :     'pmin =',pmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'pmax =',pmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (qv(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     qvmax,qvmin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'qvmin=',qvmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'qvmax=',qvmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (qc(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     qcmax,qcmin, imax,jmax,kmax, imin,jmin,kmin)


      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'qcmin=',qcmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'qcmax=',qcmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (qr(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     qrmax,qrmin, imax,jmax,kmax, imin,jmin,kmin)


      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'qrmin=',qrmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'qrmax=',qrmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (qi(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     qimax,qimin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'qimin=',qimin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'qimax=',qimax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (qs(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     qsmax,qsmin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'qsmin=',qsmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'qsmax=',qsmax,' at i=',imax,', j=',jmax,', k=',kmax

      CALL a3dmax
     :     (qh(1,1,1,tlevel),1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :     qhmax,qhmin, imax,jmax,kmax, imin,jmin,kmin)

      write(6,'(2(1x,a,f13.7,3(a,i3)))')
     :     'qhmin=',qhmin,' at i=',imin,', j=',jmin,', k=',kmin,
     :     'qhmax=',qhmax,' at i=',imax,', j=',jmax,', k=',kmax

      write(6,*) ' '

      IF(ncalls.eq.0) THEN

        IF( dirname.ne.' ' ) THEN

          maxfn  = dirname(1:ldirnam)//'/'//runname(1:6)//'.maxmin'
          lmaxfn = 13 + ldirnam + 1

        ELSE

          maxfn  = runname(1:6)//'.maxmin'
          lmaxfn = 13

        ENDIF

        OPEN(nchmax,form='formatted',status='unknown',
     :              file=maxfn(1:lmaxfn),iostat=istat)

        IF( istat.ne.0) THEN

          write(6,'(/a,i2,/a/)')
     :    ' Error occured when opening file '//runname(1:6)//'.maxmin'//
     :    ' using FORTRAN unit ',nchmax,' Program stopped in MAXMIN.'
          STOP

        ENDIF

        write(nchmax,'(a)') ''''//runname//''''

        write(nchmax,'(t4,a,t15,a,t25,a,t35,a,t45,a,t55,a,t65,a,t75,a,
     :       t85,a,t95,a,t105,a,t115,a,t125,a,t135,a)')
     :       '''time','UMIN','UMAX','VMIN','VMAX','WMIN','WMAX','PTMIN',
     :       'PTMAX','PMIN','PMAX','QVMAX','QCMAX','QRMAX'//''''

      ENDIF

c     write(nchmax,'(f9.3,7f10.5,/9x,f10.5,2f10.3,3f10.5)')

      write(nchmax,'(f9.3,7f10.5, f10.5,2f10.2,3f9.5)')
     :      curtim,umin,umax,vmin,vmax,wmin,wmax,ptmin,ptmax,
     :      pmin,pmax,qvmax*1000,qcmax*1000.0,qrmax*1000.0
c    :      ,qimax*1000,qsmax*1000.0,qhmax*1000.0

      ncalls = 1


      RETURN
      END

c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BASPRT                     ######
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 BASPRT(nx,ny,nz,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar)
c
c#######################################################################
c
c     PURPOSE:
c
c     Print out base state fields in FORTRAN unit nch=6.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     11/13/91.
c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
C
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       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     OUTPUT:
c
c       None.
c
c#######################################################################
c

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

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

      real ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                                ! (kg/kg)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,j,k,mode
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external wrigar
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      write(6,'(/1x,a/)') 'Print out of base state fields.'
c
c#######################################################################
c
c     Formatted printing of x-y slices through the mid-level
c
c#######################################################################
c
      mode = 1

      k = (nz-1)/2+1

      CALL wrigar(ubar,1,nx,1,ny,1,nz,1,nx,1,ny-1,k,k,
     :    'ubar xy',0.0, mode )

      CALL wrigar(vbar,1,nx,1,ny,1,nz,1,nx-1,1,ny,k,k,
     :    'vbar xy',0.0, mode )

      CALL wrigar(ptbar,1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :    'ptbar xy',0.00, mode )

      CALL wrigar(pbar,1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :    'pbar xy',0.0, mode )

      CALL wrigar(rhobar,1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :    'rhobar xy',0.00, mode )

      CALL wrigar(qvbar,1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :    'qvbar xy',0.0, mode )

c
c#######################################################################
c
c     Formatted printing of x-z slices through the domain center
c
c#######################################################################
c

      j = (ny-1)/2+1

      CALL wrigar(ubar,1,nx,1,ny,1,nz,1,nx,j,j,1,nz-1
     :   ,'ubar x-z',0.0, 2 )

      CALL wrigar(vbar,1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :   ,'vbar x-z',0.0, 2 )

      CALL wrigar(ptbar,1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :   ,'ptbar xz',0.00, 2 )

      CALL wrigar(pbar,1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :   ,'pbar xz',0.0, 2 )

      CALL wrigar(rhobar,1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :   ,'rhobar xz',0.0, 2 )

      CALL wrigar(qvbar,1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :   ,'qvbar xz',0.0, 2 )

c
c#######################################################################
c
c     Formatted printing of y-z slices through the domain center
c
c#######################################################################
c

      i = (nx-1)/2+1

      CALL wrigar(ubar,1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'ubar y-z',0.0, 3 )

      CALL wrigar(vbar,1,nx,1,ny,1,nz,i,i,1,ny  ,1,nz-1
     :   ,'vbar y-z',0.0, 3 )


      CALL wrigar(ptbar,1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'ptbar yz',0.00, 3 )

      CALL wrigar(pbar,1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'pbar yz',0.0, 3 )

      CALL wrigar(rhobar,1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'rhobar yz',0.0, 3 )

      CALL wrigar(qvbar,1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'qvbar yz',0.0, 3 )

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE FMTPRT                     ######
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 FMTPRT(nx,ny,nz, tlevel,
     :           u, v, w, ptprt, pprt, qv,qc,qr,qi,qs,qh, km,
     :           x,y,z,zp,hterain, j1,j2,j3, tem1)
c
c#######################################################################
c
c     PURPOSE:
c
c     Produce formatted print out of array tables in FORTRAN unit nch=6.
c     By default, the x-y, x-z and y-z slices through the domain center 
c     are printed.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     5/29/91.
c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
c
c#######################################################################
c
c     INPUT:
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       tlevel   The time level at which the data are printed.
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        Vertical component of Cartesian velocity (m/s)
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c       qv       Water vapor specific humidity (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
c
c       km       Turbulent mixing coefficient (m**2/s)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       hterain  The height of terrain (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c
c     OUTPUT:
c
c       None.
c
c     WORK ARRAYS:
c
c       tem1     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     ! Number of grid points in 3 directions

      integer tlevel       ! Time level at which data are printed.

      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 km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )

      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 tem1(nx,ny,nz)       ! Temporary array.
c
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      integer iplot(3)
      integer mode      ! Control the printing of individual slice
             		! mode = 1: print x-y slice
             		!      = 2: print x-z slice
              		!      = 3: print y-z slice
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external wrigar
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     The following parameters are passed into this subroutine through
c     a common block in globcst.inc, and they controls the output of
c     water and ice variables.
c
c     mstout =0 or 1. If mstout=0, water variables are not printed.
c     iceout =0 or 1. If iceout=0, printing of qi, qs and qh is skipped.
c     basout =0 or 1. If basout=0, base state variables are not dumped.
c
c#######################################################################
c

c
c#######################################################################
c
c     Formatted printing of x-y slices through the mid-level
c
c#######################################################################
c
      mode = 1
      k = (nz-1)/2+1

      CALL wrigar(u(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx,1,ny-1,k,k,
     :    'u x-y',0.0, mode )

      CALL wrigar(v(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny,k,k,
     :    'v x-y',0.0, mode )

      CALL wrigar(w(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k+1,k+1,
     :    'w x-y',0.0, mode )

      CALL wrigar(ptprt(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :    'ptprt xy',0.00, mode )

      CALL wrigar(pprt(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :    'pprt xy',0.0, mode )

      IF( mstout.eq.1) THEN

        CALL wrigar(qv(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :      'qv x-y',0.0, mode )

        CALL wrigar(qc(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :      'qc x-y',0.0, mode )

        CALL wrigar(qr(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :      'qr x-y',0.0, mode )

        IF( iceout.eq.1) THEN

          CALL wrigar(qi(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :        'qi x-y',0.0, mode )

          CALL wrigar(qs(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :        'qs x-y',0.0, mode )

          CALL wrigar(qh(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,
     :        'qh x-y',0.0, mode )

        ENDIF

      ENDIF

      CALL wrigar(km,1,nx,1,ny,1,nz,1,nx-1,1,ny-1,k,k,'km x-y',0.0,mode)
c
c#######################################################################
c
c     Formatted printing of x-z slices through the domain center
c
c#######################################################################
c
      mode = 2         ! Print x-z slices
c
      j = (ny-1)/2+1

      CALL wrigar(u(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx,j,j,1,nz-1
     :     ,'u x-z',0.0, mode )

      CALL wrigar(v(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :     ,'v x-z',0.0, mode )

      CALL wrigar(w(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz
     :     ,'w x-z',0.0, mode )

      CALL wrigar(ptprt(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :     ,'ptprt xz',0.00, mode )

      CALL wrigar(pprt(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :     ,'pprt xz',0.0, mode )
c
      IF( mstout.eq.1) THEN

        CALL wrigar(qv(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :      ,'qv xz',0.0, mode )

        CALL wrigar(qc(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :      ,'qc xz',0.0, mode )

        CALL wrigar(qr(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :      ,'qr xz',0.0, mode )

        IF( iceout.eq.1) THEN

          CALL wrigar(qi(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :        ,'qi xz',0.0, mode )

          CALL wrigar(qs(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :        ,'qs xz',0.0, mode )

          CALL wrigar(qh(1,1,1,tlevel),1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :        ,'qh xz',0.0, mode )

        ENDIF

      ENDIF

      CALL wrigar(km ,1,nx,1,ny,1,nz,1,nx-1,j,j,1,nz-1
     :     ,'km xz',0.0, 2 )
c
c
c#######################################################################
c
c     Formatted printing of y-z slices through the domain center
c
c#######################################################################
c
      mode = 3         ! Print y-z slices
c
      i = (nx-1)/2+1

      CALL wrigar(u(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'u y-z',0.0, mode )

      CALL wrigar(v(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny  ,1,nz-1
     :   ,'v y-z',0.0, mode )

      CALL wrigar(w(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz
     :   ,'w y-z',0.0, mode )

      CALL wrigar(ptprt(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'ptprt yz',0.00, mode )

      CALL wrigar(pprt(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'pprt yz',0.0, mode )

      IF( mstout.eq.1) THEN

        CALL wrigar(qv(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :    ,'qv yz',0.0, mode )

        CALL wrigar(qc(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :    ,'qc yz',0.0, mode )

        CALL wrigar(qr(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :    ,'qr yz',0.0, mode )

        IF( iceout.eq.1) THEN

          CALL wrigar(qi(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :      ,'qi yz',0.0, mode )

          CALL wrigar(qs(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :      ,'qs yz',0.0, mode )

          CALL wrigar(qh(1,1,1,tlevel),1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :      ,'qh yz',0.0, mode )

        ENDIF

      ENDIF

      CALL wrigar(km ,1,nx,1,ny,1,nz,i,i,1,ny-1,1,nz-1
     :   ,'km yz',0.0, mode )

c
c#######################################################################
c
c     Contour plotting of 2-d slices through the domain center using
c     ncargraphic routine CONREC. Not supported.
c
c#######################################################################
c

      iplot(1)=0   ! x-y slices
      iplot(2)=0   ! x-z slices
      iplot(3)=0   ! y-z slices

      DO 200 mode = 1,3

        IF( iplot(mode).eq.0) GOTO 200

        CALL pltary(u(1,1,1,tlevel),tem1,nx,nx,ny,ny-1,nz-1,
     :             'u', 0.0 ,curtim,mode)

        CALL pltary(v(1,1,1,tlevel),tem1,nx,nx-1,ny,ny,nz-1,
     :             'v', 0.0 ,curtim,mode)

        CALL pltary(w(1,1,2,tlevel),tem1,nx,nx-1,ny,ny-1,nz-2,
     :             'w', 0.0 ,curtim,mode)

        CALL pltary(ptprt(1,1,1,tlevel),tem1,nx,nx-1,ny,ny-1,nz-1,
     :             'ptprt', 0.0 ,curtim,mode)

        CALL pltary(pprt(1,1,1,tlevel),tem1,nx,nx-1,ny,ny-1,nz-1,
     :             'pprt', 0.0 ,curtim,mode)

200   CONTINUE

      RETURN
      END
