c
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE DTAREAD                    ######
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 DTAREAD(nx,ny,nz, 
     :           hinfmt ,grdbasfn,lengbf,datafn,lendtf,time, 
     :           x,y,zc,z, uprt ,vprt ,wprt ,ptprt, pprt ,
     :           qvprt, qc, qr, qi, qs, qh, km,
     :           ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :           ireturn, tem1)

c
c#######################################################################
c
c
c     hinfmt  = 1 for machine dependent unformatted binary dump,
c     hinfmt  = 2 for formatted ascii dump,
c     hinfmt  = 3 for NCSA HDF file format.
c     hinfmt  = 4 for machine dependent packed unformatted binary dump
c
c#######################################################################
c
c     AUTHOR: Ming Xue, CAPS/OU. 
c       2/19/1992.     
c
c     MODIFICATION HISTORY:
c
c#######################################################################
c
c     INPUT:
c
c       nx,ny,nz The dimension of data arrays 
c       hinfmt   The format of the history data dump
c       grdbasfn Name of the grid/base state array file
c       lengbf   Length of the grid/base state data file name string
c       datafn   Name of the other time dependent data file
c       lendtf   Length of the data file name string 
c     
c     DATA ARRAYS READ IN:
c
c       x        x coordinate of grid points in physical/comp. space (km)
c       y        y coordinate of grid points in physical/comp. space (km)
c       zc       z coordinate of grid points in computational space (km)
c       z        z coordinate of grid points in physical space (km)
c
c       uprt     x component of perturbation velocity (m/s)
c       vprt     y component of perturbation velocity (m/s)
c       wprt     vertical component of perturbation velocity in Cartesian
c                coordinates (m/s).
c
c       ptprt    perturbation potential temperature (K)
c       pprt     perturbation pressure (Pascal)
c
c       qvprt    perturbation water vapor mixing ratio (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       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       wbar     Base state z velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor mixing ratio (kg/kg)
c
c     OUTPUT: 
c
c       time     The time of the input data (s)
c       x        x coordinate of grid points in physical/comp. space (km)
c       y        y coordinate of grid points in physical/comp. space (km)
c       zc       z coordinate of grid points in computational space (km)
c       z        z coordinate of grid points in physical space (km)
c
c       uprt     x component of perturbation velocity (m/s)
c       vprt     y component of perturbation velocity (m/s)
c       wprt     vertical component of perturbation velocity in Cartesian
c                coordinates (m/s).
c
c       ptprt    perturbation potential temperature (K)
c       pprt     perturbation pressure (Pascal)
c
c       qvprt    perturbation water vapor mixing ratio (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       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       wbar     Base state z velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor mixing ratio (kg/kg)
c
c     WORK ARRAYS:
c
c       tem1     Temporary work array.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 

      integer nx, ny, nz
      real time
      integer hinfmt ,lengbf, lendtf
      character grdbasfn*(*)
      character datafn*(*)
c
      real x    (nx,ny,nz)
      real y    (nx,ny,nz)
      real zc   (nx,ny,nz)
      real z    (nx,ny,nz)

      real uprt (nx,ny,nz)
      real vprt (nx,ny,nz)
      real wprt (nx,ny,nz) 
      real ptprt(nx,ny,nz)
      real pprt (nx,ny,nz)
      real qvprt(nx,ny,nz)
      real qc   (nx,ny,nz)
      real qr   (nx,ny,nz)
      real qi   (nx,ny,nz)
      real qs   (nx,ny,nz)
      real qh   (nx,ny,nz)
      real km    (nx,ny,nz)     
      real ubar  (nx,ny,nz)
      real vbar  (nx,ny,nz)
      real wbar  (nx,ny,nz)
      real ptbar (nx,ny,nz)
      real rhobar(nx,ny,nz)
      real pbar  (nx,ny,nz)
      real qvbar (nx,ny,nz)
c
      real tem1(nx,ny,nz)

      integer grdread
      data grdread /0/
      save grdread

      integer ngchan,nchanl,ireturn,istat
      integer grdbas

      real btime ! The time of the base state data
      real amin, amax
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'indtflg.inc'
      include 'globcst.inc'
c
c#######################################################################
c
c     Subroutines called:
c
c#######################################################################
c
      external binread
      external ascread
      external hdfread
      external pakread
      external a3dmax0 
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Read data fields.
c
c#######################################################################
c
      nchanl = 64
      grdbas = 0

      IF( hinfmt .eq.1 ) THEN 

        open(unit=nchanl,file=datafn(1:lendtf),
     :		status='old',form='unformatted',iostat=istat)

        IF( istat.ne.0 ) GOTO 998

        CALL binread(nx,ny,nz, grdbas, nchanl , time, x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn)

        close(unit=nchanl)

      ELSEIF( hinfmt .eq.2 ) THEN 

        open(unit=nchanl,file=datafn(1:lendtf),
     :		status='old',form='formatted',iostat=istat)

        IF( istat.ne.0 ) GOTO 998

        CALL ascread(nx,ny,nz, grdbas, nchanl , time, x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn)

        close(unit=nchanl)

      ELSEIF( hinfmt .eq.3 ) THEN 

        CALL hdfread(nx,ny,nz, grdbas, datafn(1:lendtf),time,x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn, tem1)

      ELSEIF( hinfmt .eq.4 ) THEN 

        open(unit=nchanl,file=datafn(1:lendtf),
     :		status='old',form='unformatted',iostat=istat)

        IF( istat.ne.0 ) GOTO 998

        CALL pakread(nx,ny,nz, grdbas, nchanl , time, x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn, tem1)

        close(unit=nchanl)

      ELSE
	  
        write(6,'(a,i3,a)') 
     :        ' Data format flag had an invalid value ',
     :         	hinfmt ,' program stopped.'
        STOP
      ENDIF

      CLOSE(unit=nchanl)
c
c#######################################################################
c
c     Open and read grid and base data file depending on the values of 
c     parameters grdin and basin, which are read in from the 
c     time dependent data set. If grdin or basin is zero, the grid and 
c     base state arrays have to be read in from a separate file. 
c
c#######################################################################
c
      IF( grdread.eq.0 .and. grdin.eq.0 ) THEN 

      ngchan= 63
      grdbas = 1
c
c#######################################################################
c
c     Read grid and base fields.
c
c#######################################################################
c
        IF( hinfmt .eq.1 ) THEN 

          open(unit=ngchan,file=grdbasfn(1:lengbf),
     :		status='old',form='unformatted',iostat=istat)

          IF( istat .ne. 0 ) GOTO 998

          CALL binread(nx,ny,nz, grdbas, ngchan , btime, x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn)

          close(unit=ngchan)

        ELSEIF( hinfmt .eq.2 ) THEN 

          open(unit=ngchan,file=grdbasfn(1:lengbf),
     :		status='old',form='formatted',iostat=istat)

          IF( istat .ne. 0 ) GOTO 998

          CALL ascread(nx,ny,nz, grdbas, ngchan , btime, x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn)

          close(unit=ngchan)

        ELSEIF( hinfmt .eq.3 ) THEN 

          CALL hdfread(nx,ny,nz, grdbas, grdbasfn(1:lengbf),
     :               btime,x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn, tem1)

        ELSEIF( hinfmt .eq.4 ) THEN 

          open(unit=ngchan,file=grdbasfn(1:lengbf),
     :		status='old',form='unformatted',iostat=istat)

          IF( istat .ne. 0 ) GOTO 998

          CALL pakread(nx,ny,nz,grdbas,  ngchan , btime, x,y,zc,z,
     :               uprt ,vprt ,wprt ,ptprt, pprt , 
     :               qvprt, qc, qr, qi, qs, qh, km,
     :               ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :               ireturn, tem1)

          close(unit=ngchan)

        ELSE
	  
          write(6,'(a,i3,a)') 
     :        ' Data format flag had an invalid value ',
     :         	hinfmt ,' program stopped.'
          STOP

        ENDIF

        grdread = 1

        CLOSE(unit=ngchan)

      ENDIF

      write(6,'(/1x,a/)') 'Min. and max. of the data arrays read in:'

      CALL a3dmax0(x,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(/1x,2(a,e13.6))') 'xmin    = ', amin,',  xmax    =',amax

      CALL a3dmax0(y,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'ymin    = ', amin,',  ymax    =',amax

      CALL a3dmax0(z,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'zmin    = ', amin,',  zmax    =',amax

      CALL a3dmax0(ubar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,',  ubarmax =',amax

      CALL a3dmax0(vbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,',  vbarmax =',amax

      CALL a3dmax0(wbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'wbarmin = ', amin,',  wbarmax =',amax

      CALL a3dmax0(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,',  ptbarmax=',amax

      CALL a3dmax0(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,',  pbarmax =',amax

      IF( mstin.eq.1) THEN

      CALL a3dmax0(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin,',  qvbarmax=',amax

      ENDIF

      CALL a3dmax0(uprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'uprtmin = ', amin,',  uprtmax =',amax

      CALL a3dmax0(vprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'vprtmin = ', amin,',  vprtmax =',amax

      CALL a3dmax0(wprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'wprtmin = ', amin,',  wprtmax =',amax

      CALL a3dmax0(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,',  ptprtmax=',amax

      CALL a3dmax0(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,',  pprtmax =',amax

      IF( mstin.eq.1) THEN

      CALL a3dmax0(qvprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'qvprtmin= ', amin,',  qvprtmax=',amax

      CALL a3dmax0(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'qcmin   = ', amin,',  qcmax   =',amax

      CALL a3dmax0(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'qrmin   = ', amin,',  qrmax   =',amax

      IF( icein.eq.1) THEN

      CALL a3dmax0(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'qimin   = ', amin,',  qimax   =',amax

      CALL a3dmax0(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'qsmin   = ', amin,',  qsmax   =',amax

      CALL a3dmax0(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.6))') 'qhmin   = ', amin,',  qhmax   =',amax

      ENDIF
      ENDIF

      if (trbin.eq.1) THEN 

      CALL a3dmax0(km,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      write(6,'(1x,2(a,e13.4))') 'kmmin   = ', amin,',  kmmax   =',amax
 
      ENDIF

c
      RETURN

998   CONTINUE
      write(6,'(1x,a,a,/1x,i3,a)')
     :  'Error occured when opening file ',datafn(1:lendtf),
     :  'using FORTRAN unit ',nchanl,' Program stopped in DTAREAD.'
      STOP

999   CONTINUE
      write(6,'(1x,a,a,/1x,i3,a)')
     :  'Error occured when opening file ',grdbasfn(1:lengbf),
     :  'using FORTRAN unit ',ngchan,' Program stopped in DTAREAD.'
      STOP

      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BINREAD                    ######
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 binread(nx,ny,nz, grdbas, inch, time, x,y,z,zp,
     :           uprt, vprt, wprt, ptprt, pprt, 
     :           qvprt, qc, qr, qi, qs, qh, km,
     :           ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :           ireturn)

c
c#######################################################################
c
c     PURPOSE:
c
c     Read in data field from binary dataset created by ARPS 3.0
c
c     All data read in are located at the grid cell centers.  Variables 
c     u, v and w were averaged to the cell center before being dumped. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     2/27/92.    
c
c     MODIFICATION HISTORY:
c
c     6/08/92  Added full documentation (K. Brewster)
c
c     7/14/92 (K. Brewster)
c     Added runname, comment and version number reading
C
c     8/20/92 (M. Xue)
c     Added data reading of computational z coordinate array z. 
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       grdbas   Data read mode flag. 
c                If grdbas=1, only grid and base state arrays will be read 
c                If grdbas=0, all arrays will be read based on data 
c                             parameter setting.  
c       inch     Channel number for binary reading.
c                This channel must be opened for unformatted reading
c                by the calling routine.
c
c     OUTPUT:
c
c       time     Time in seconds of data read from "filename"
c
c       x        x coordinate of grid points in physical/comp. space (km)
c       y        y coordinate of grid points in physical/comp. space (km)
c       z        z coordinate of grid points in computational space (km)
c       zp       z coordinate of grid points in physical space (km)
c
c       uprt     x component of perturbation velocity (m/s)
c       vprt     y component of perturbation velocity (m/s)
c       wprt     Vertical component of perturbation velocity in Cartesian 
c                coordinates (m/s).
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c
c       qvprt    Perturbation water vapor mixing ratio (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       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       wbar     Base state z velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor mixing ratio (kg/kg)
c
c       ireturn  Return status indicator
c                ireturn=0   successful read of all data
c                ireturn=1   error reading data
c                ireturn=2   end-of-file reached during read attempt
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nx,ny,nz
c
      real x     (nx,ny,nz)    
      real y     (nx,ny,nz)   
      real z     (nx,ny,nz)  
      real zp    (nx,ny,nz)     
c
      integer grdbas
      integer inch
      real time
c
      real uprt  (nx,ny,nz)     ! Perturbation u-velocity (m/s)
      real vprt  (nx,ny,nz)     ! Perturbation v-velocity (m/s)
      real wprt  (nx,ny,nz)     ! Perturbation w-velocity (m/s)
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)
c
      real qvprt (nx,ny,nz)     ! Perturbation water vapor mixing ratio (kg/kg)
      real qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
      real km    (nx,ny,nz)     ! Turbulent mixing coefficient
c
      real ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real wbar  (nx,ny,nz)     ! Base state w-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 mixing ratio

      integer ireturn
c
c#######################################################################
c
c     Parameters describing routine that wrote the gridded data
c
c#######################################################################
c
      character*40 fmtver,fmtverin
      character*10 tmunit
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer lchanl
      parameter (lchanl=6)                    ! listing unit number
c
      integer i
      character label*12
      integer knt,kntscal
      integer kntin,kntscin
      integer nxin,nyin,nzin
      integer bgrdin, bbasin, bvarin, bmstin, bicein, btrbin
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'indtflg.inc' 
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external kntary
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Read header info
c
c#######################################################################
c
      read(inch,err=110,end=120) fmtverin
      read(inch,err=110,end=120) runname
      read(inch,err=110,end=120) nocmnt
      IF( nocmnt.gt.0 ) THEN
        DO 10 i=1,nocmnt
          read(inch,err=110,end=120) cmnt(i)
 10     CONTINUE
      ENDIF
      read(inch,err=110,end=120) time,tmunit

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

      IF( nocmnt.gt.0 ) THEN
        DO 12 i=1,nocmnt
          write(6,'(1x,a)') cmnt(i)
 12     CONTINUE
      ENDIF

      IF( grdbas.eq.1 ) THEN ! Read grid and base state arrays


        write(lchanl,'(1x,a,f8.1,a,f8.3,a/)')
     :       'To read grid and base state data at time ', time,
     :       ' secs = ',(time/60.),' mins.'

        read(inch,err=110,end=120) 
     :       bgrdin,bbasin,bvarin,bmstin,bicein,btrbin
        read(inch,err=110,end=120) kntin
        read(inch,err=110,end=120) kntscin

        CALL kntary(fmtver,bgrdin,bbasin,bvarin,bmstin,bicein,btrbin,
     :              knt,kntscal)

      ELSE ! Normal data reading

        write(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:',
     :       time,' secs = ',(time/60.),' mins.'

        read(inch,err=110,end=120) grdin,basin,varin,mstin,icein,trbin
        read(inch,err=110,end=120) kntin
        read(inch,err=110,end=120) kntscin

        CALL kntary(fmtver,grdin,basin,varin,mstin,icein,trbin,
     :              knt,kntscal)
c
c     Check if fmtverin=fmtver, kntin=knt, kntscal=kntscin
c     If not, there is inconsistency in the data.
c

      ENDIF
c
C#######################################################################
c
c     Get dimensions of data in binary file and check against
c     the dimensions passed to BINREAD
c
C#######################################################################
c
      read(inch,err=110,end=120) nxin, nyin, nzin

      IF( nxin.ne.nx .or. nyin.ne.ny .or. nzin.ne.nz ) THEN 
        write(6,'(1x,a)') 
     :       ' Dimensions in BINREAD inconsistent with data.'
        write(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
        write(6,'(1x,a)') 
     :       ' Program aborted in BINREAD.'
        STOP
      ENDIF
c
C#######################################################################
c
c     Read in x,y and z at grid cell centers (scalar points).
c
C######################################################################
c
      IF( grdin.eq.1 .or. grdbas.eq.1 ) THEN
        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) x
        write(lchanl,910) label,' x.' 

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) y
        write(lchanl,910) label,' y.' 

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) z
        write(lchanl,910) label,' z.' 

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) zp
        write(lchanl,910) label,' zp.'
      ENDIF  ! grdin
c
C#######################################################################
c
c     Read in base state fields
c
C######################################################################
c
      IF( basin.eq.1 .or. grdbas.eq.1 ) THEN 

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) ubar
        write(lchanl,910) label,' ubar.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) vbar
        write(lchanl,910) label,' vbar.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) wbar
        write(lchanl,910) label,' wbar.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) ptbar
        write(lchanl,910) label,' ptbar.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) pbar
        write(lchanl,910) label,' pbar.'


        IF( mstin.eq.1) THEN 
          read(inch,err=110,end=120) label
          read(inch,err=110,end=120) qvbar
          write(lchanl,910) label,' qvbar.'
        ENDIF

      ENDIF

      IF( grdbas.eq.1 ) GOTO 930 

      IF( varin.eq.1 ) THEN
c
C#######################################################################
c
c     Read in u,v and w averaged to grid cell centers (scalar points).
c
C######################################################################
c
        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) uprt
        write(lchanl,910) label,' uprt.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) vprt
        write(lchanl,910) label,' vprt.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) wprt
        write(lchanl,910) label,' wprt.'
c
C#######################################################################
c
c     Read in scalars
c
C######################################################################
c
        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) ptprt
        write(lchanl,910) label,' ptprt.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) pprt
        write(lchanl,910) label,' pprt.'

      ENDIF
c
C#######################################################################
c
c     Read in moisture variables
c
C######################################################################
c
      IF( mstin .eq.1 ) THEN

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) qvprt
        write(lchanl,910) label,' qvprt.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) qc 
        write(lchanl,910) label,' qc.'

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) qr 
        write(lchanl,910) label,' qr.'

        IF( icein.eq.1 ) THEN 

          read(inch,err=110,end=120) label
          read(inch,err=110,end=120) qi 
          write(lchanl,910) label,' qi.'

          read(inch,err=110,end=120) label
          read(inch,err=110,end=120) qs 
          write(lchanl,910) label,' qs.'

          read(inch,err=110,end=120) label
          read(inch,err=110,end=120) qh 
          write(lchanl,910) label,' qh.'

        ENDIF
      ENDIF  

      IF( trbin.eq.1 ) THEN

        read(inch,err=110,end=120) label
        read(inch,err=110,end=120) km
        write(lchanl,910) label,' km.'

      ENDIF


910   FORMAT(1x,'Field ',a12,' was read into array',a)

c
C#######################################################################
c
c     Friendly exit message
c
C######################################################################
c
930   CONTINUE

      write(6,'(/a,F8.1,a/)')
     : ' Data at time=', time/60,' (min) were successfully read.'

      ireturn = 0 
      RETURN
c
C#######################################################################
c
c     Destination for error during read
c
C######################################################################
c

110   CONTINUE
      write(6,'(/a/)') ' Error reading data in BINREAD'
      ireturn=1
      RETURN
c
C#######################################################################
c
c     Destination for end-of-file during read
c
C######################################################################
c

120   CONTINUE
      write(6,'(/a/)') ' End of file reached in BINREAD'
      ireturn=2
      RETURN
      END     
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ASCREAD                    ######
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 ascread(nx,ny,nz, grdbas, inch, time, x,y,z,zp,
     :           uprt, vprt, wprt, ptprt, pprt, 
     :           qvprt, qc, qr, qi, qs, qh, km,
     :           ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,
     :           ireturn)

c
c#######################################################################
c
c     PURPOSE:
c
c     Read history data from channel nchanl in formatted ASCII data
c     format.
c
c     All data read in are located at the grid cell centers.  Variables 
c     u, v and w were averaged to the cell center before being dumped. 
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     6/02/92.    
c
c     MODIFICATION HISTORY:
c
c     7/14/92 (K. Brewster)
c     Added runname, comment and version number reading
C
c     8/20/92 (M. Xue)
c     Added data reading of computational z coordinate array z. 
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       grdbas   Data read mode flag. 
c                If grdbas=1, only grid and base state arrays will be read 
c                If grdbas=0, all arrays will be read based on data 
c                             parameter setting.  
c       inch     Channel number for ASCII reading.
c                This channel must be opened for formatted reading
c                by the calling routine.
c
c     OUTPUT:
c
c       time     Time in seconds of data read from "filename"
c
c       x        x coordinate of grid points in physical/comp. space (km)
c       y        y coordinate of grid points in physical/comp. space (km)
c       z        z coordinate of grid points in computational space (km)
c       zp       z coordinate of grid points in physical space (km)
c
c       uprt     x component of perturbation velocity (m/s)
c       vprt     y component of perturbation velocity (m/s)
c       wprt     Vertical component of perturbation velocity in Cartesian 
c                coordinates (m/s).
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c
c       qvprt    Perturbation water vapor mixing ratio (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       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       wbar     Base state z velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor mixing ratio (kg/kg)
c
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nx,ny,nz
c
      integer grdbas
      integer inch
c
      real time
c
      real x     (nx,ny,nz)     
      real y     (nx,ny,nz)    
      real z     (nx,ny,nz)   
      real zp    (nx,ny,nz)     
c
      real uprt  (nx,ny,nz)     ! Perturbation u-velocity (m/s)
      real vprt  (nx,ny,nz)     ! Perturbation v-velocity (m/s)
      real wprt  (nx,ny,nz)     ! Perturbation w-velocity (m/s)
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      real pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)
c
      real qvprt (nx,ny,nz)     ! Water vapor mixing ratio (kg/kg)
      real qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
      real qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
      real qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
      real qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
      real qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
      real km    (nx,ny,nz)     ! Turbulent mixing coefficient
c
      real ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real wbar  (nx,ny,nz)     ! Base state w-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 mixing ratio

      integer ireturn
c
c#######################################################################
c
c     Parameters describing routine that wrote the gridded data
c
c#######################################################################
c
      character*40 fmtver,fmtverin
      character*10 tmunit
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer lchanl
      parameter (lchanl=6)                ! listing unit, 6=std output
c
      integer i
      character label*12
      integer knt,kntscal
      integer kntin,kntscin
      integer nxin,nyin,nzin
      integer bgrdin, bbasin, bvarin, bmstin, bicein, btrbin
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'indtflg.inc' 
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external kntary
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c#######################################################################
c
c     Read header info
c
c#######################################################################
c
      read(inch,'(1x,a40)',err=110,end=120) fmtverin
      read(inch,'(1x,a80)',err=110,end=120) runname
      write(6,'(//''  The name of this run is:  '',a//)') runname
c
      read(inch,'(1x,i4)',err=110,end=120) nocmnt
      IF( nocmnt.gt.0 ) THEN
        DO 10 i=1,nocmnt
          read(inch,'(1x,a80)',err=110,end=120) cmnt(i)
 10     CONTINUE
      ENDIF

      IF( nocmnt.gt.0 ) THEN
        DO 12 i=1,nocmnt
          write(6,'(1x,a)') cmnt(i)
 12     CONTINUE
      ENDIF
c
      read(inch,'(1x,e16.8,1x,a10)',err=110,end=120) time,tmunit

      IF( grdbas.eq.1 ) THEN ! Read grid and base state arrays


        write(lchanl,'(1x,a,f8.1,a,f8.3,a/)')
     :       'To read grid and base state data at time ', time,
     :       ' secs = ',(time/60.),' mins.'

        read(inch,'(6(1x,i4))',err=110,end=120)
     :       bgrdin,bbasin,bvarin,bmstin,bicein,btrbin
        read(inch,'(1x,i12)',err=110,end=120) kntin
        read(inch,'(1x,i12)',err=110,end=120) kntscin

        CALL kntary(fmtver,bgrdin,bbasin,bvarin,bmstin,bicein,btrbin,
     :              knt,kntscal)

      ELSE ! Normal data reading

        write(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:',
     :       time,' secs = ',(time/60.),' mins.'

        read(inch,'(6(1x,i4))',err=110,end=120)
     :       grdin,basin,varin,mstin,icein,trbin
        read(inch,'(1x,i12)',err=110,end=120) kntin
        read(inch,'(1x,i12)',err=110,end=120) kntscin

        CALL kntary(fmtver,grdin,basin,varin,mstin,icein,trbin,
     :              knt,kntscal)
c
c     Check if fmtverin=fmtver, kntin=knt, kntscal=kntscin
c     If not, there is inconsistency in the data.
c

      ENDIF
c
C#######################################################################
c
c     Get dimensions of data in ASCII file and check against
c        the dimensions passed to ASCREAD
c
C#######################################################################
c
      read(inch,'(1x, 3i12)',err=110,end=120) nxin,nyin,nzin
      IF( nxin.ne.nx .or. nyin.ne.ny .or. nzin.ne.nz ) THEN
        write(6,'(1x,a)') 
     :       ' Dimensions in ASCREAD inconsistent with data.'
        write(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
        write(6,'(1x,a)') 
     :       ' Program aborted in ASCREAD.'
        STOP
      ENDIF
c
C#######################################################################
c
c     Read in x,y and z at grid cell centers (scalar points).
c
C######################################################################
c
      IF( grdin.eq.1 .or. grdbas.eq.1 ) THEN
        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) x 
        write(lchanl,920) label,' x.' 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) y 
        write(lchanl,920) label,' y.' 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) z 
        write(lchanl,920) label,' z.' 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) zp
        write(lchanl,920) label,' zp.' 
      ENDIF
c
C#######################################################################
c
c     Read in base state fields
c
C######################################################################
c
      IF( basin.eq.1 .or. grdbas.eq.1 ) THEN 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) ubar 
        write(lchanl,920) label,' ubar.' 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) vbar 
        write(lchanl,920) label,' vbar.' 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) wbar
        write(lchanl,920) label,' wbar.'

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) ptbar
        write(lchanl,920) label,' ptbar.'

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) pbar 
        write(lchanl,920) label,' pbar.' 

        IF( mstin.eq.1 ) THEN 
          read(inch,900,err=110,end=120) label 
          read(inch,910,err=110,end=120) qvbar
          write(lchanl,920) label,' qvbar.'
        ENDIF

      ENDIF

      IF( grdbas.eq.1 ) GOTO 930 

      IF( varin.eq.1 ) THEN
c
C#######################################################################
c
c     Read in u,v and w averaged to grid cell centers (scalar points).
c
C######################################################################
c
        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) uprt
        write(lchanl,920) label,' uprt.'

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) vprt
        write(lchanl,920) label,' vprt.'

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) wprt
        write(lchanl,920) label,' wprt.'
c
C#######################################################################
c
c     Read in scalars
c
C######################################################################
c
        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) ptprt 
        write(lchanl,920) label,' ptprt.'

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) pprt 
        write(lchanl,920) label,' pprt.'
      ENDIF

      IF( mstin.eq.1 ) THEN 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) qvprt
        write(lchanl,920) label,' qvprt.'

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) qc 
        write(lchanl,920) label,' qc.'

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) qr
        write(lchanl,920) label,' qr.' 

        IF( icein.eq.1 ) THEN 

          read(inch,900,err=110,end=120) label 
          read(inch,910,err=110,end=120) qi
          write(lchanl,920) label,' qi.'

          read(inch,900,err=110,end=120) label 
          read(inch,910,err=110,end=120) qs 
          write(lchanl,920) label,' qs.'

          read(inch,900,err=110,end=120) label 
          read(inch,910,err=110,end=120) qh
          write(lchanl,920) label,' qh.'

        ENDIF

      ENDIF

      IF( trbin.eq.1 ) THEN 

        read(inch,900,err=110,end=120) label 
        read(inch,910,err=110,end=120) km
        write(lchanl,920) label,' km.'

      ENDIF
c
C#######################################################################
c
c     Friendly exit message
c
C######################################################################
c
930   CONTINUE

      write(6,'(/a,F8.1,a/)')
     : ' Data at time=', time/60,' (min) were successfully read.'

      ireturn = 0

      RETURN

900   format(1x,a12)
910   format(1x,8e16.10)
920   FORMAT(1x,'Field ',a12,' was read into array ',a)
c
C#######################################################################
c
c     Destination for error during read
c
C######################################################################
c

110   CONTINUE
      write(6,'(/a/)') ' Error reading data in ASCREAD'
      ireturn=1
      RETURN
c
C#######################################################################
c
c     Destination for end-of-file during read
c
C######################################################################
c

120   CONTINUE
      write(6,'(/a/)') ' End of file reached in ASCREAD'
      ireturn=2
      RETURN
      END     
