!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE DTAREAD                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE dtaread(nx,ny,nz,nstyps,                                     & 23,130
           hinfmt, nchanl,grdbasfn,lengbf,datafn,lendtf, time,          &
           x,y,z,zp, uprt ,vprt ,wprt ,ptprt, pprt ,                    &
           qvprt, qc, qr, qi, qs, qh, tke, kmh,kmv,                     &
           ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,                &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                    &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,                                          &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           ireturn, tem1, tem2, tem3)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Coordinate the reading of history data of various formats.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!    2/19/1992.
!
!  MODIFICATION HISTORY:
!
!  10/06/1992  (K. Brewster)
!  Added ENTRY SETGBRD to allow the reading of
!  more than one grid-base file in a program.
!
!  4/23/93 (M. Xue)
!  New data format.
!
!  3/30/94 (M.Xue)
!  GrADS data format and surface fields.
!
!  12/01/95 (Yuhe Liu)
!  Fixed a bug in netread call for dumping grid and base data file.
!  The file name was previously given to datafn, instead of grdbasfn.
!
!  12/01/95 (Yuhe Liu)
!  Changed the order of reading grid and base file and history files.
!  Previously the code read history file first and then checked if
!  need to read the grid and base file. In order to cooperate with
!  GRIB data format which no longer contains perturbation variables,
!  the base fields have to ready so that the perturbation fields can
!  be derived from total fields substracting the base fields.
!
!  12/07/95 (Yuhe Liu)
!  Added a new data dumping format, GRIB. The parameter definitions
!  is mostly in consistance with NMC standard, despite those
!  variables which are not in the WMO and NMC's tables (Table 2)
!
!  01/22/1996 (Yuhe Liu)
!  Changed the ARPS history dumpings to dump the total values of
!  ARPS variables, instead of the perturbated values.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx,ny,nz The dimension of data arrays
!
!    hinfmt   The format of the history data dump
!             =1, machine dependent unformatted binary dump,
!             =2, formatted ascii dump,
!             =3, NCSA HDF file format.
!             =4, machine dependent packed unformatted binary dump
!
!    nchanl   FORTRAN I/O channel number for history data output.
!
!    grdbasfn Name of the grid/base state array file
!    lengbf   Length of the grid/base state data file name string
!    datafn   Name of the other time dependent data file
!    lendtf   Length of the data file name string
!
!  DATA ARRAYS READ IN:
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       z coordinate of grid points in physical space (m)
!
!    uprt     x component of perturbation velocity (m/s)
!    vprt     y component of perturbation velocity (m/s)
!    wprt     vertical component of perturbation velocity
!             in Cartesian coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!    qvprt    perturbation water vapor mixing ratio (kg/kg)
!
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state air density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!    soiltyp  Soil type
!    stypfrct  Soil type fraction
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  OUTPUT:
!
!    time     The time of the input data (s)
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       z coordinate of grid points in physical space (m)
!
!    uprt     x component of perturbation velocity (m/s)
!    vprt     y component of perturbation velocity (m/s)
!    wprt     vertical component of perturbation velocity in
!             Cartesian coordinates (m/s).
!
!    ptprt    perturbation potential temperature (K)
!    pprt     perturbation pressure (Pascal)
!
!    qvprt    perturbation water vapor mixing ratio (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    wbar     Base state z velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state air density (kg/m**3)
!    qvbar    Base state water vapor mixing ratio (kg/kg)
!
!    soiltyp  Soil type
!    stypfrct  Soil type fraction
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at surface (K)
!    tsoil    Deep soil temperature (K)
!    wetsfc   Surface soil moisture
!    wetdp    Deep soil moisture
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    ireturn  Return status indicator
!
!  WORK ARRAYS:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  REAL :: time                 ! The time of the input data (s)
  INTEGER :: hinfmt            ! The format of the history data dump
  INTEGER :: lengbf            ! Length of the grid/base state data
                               ! file name string
  INTEGER :: lendtf            ! Length of the data file name string
  CHARACTER(LEN=*) :: grdbasfn ! Name of the grid/base state array file
  CHARACTER(LEN=*) :: datafn   ! Name of the other time dependent
                               ! data file
  REAL :: x     (nx)           ! x-coord. of the physical and compu
                               ! -tational grid. Defined at u-point(m).
  REAL :: y     (ny)           ! y-coord. of the physical and compu
                               ! -tational grid. Defined at v-point(m).
  REAL :: z     (nz)           ! z-coord. of the computational grid.
                               ! Defined at w-point on the staggered
                               ! grid(m).
  REAL :: zp    (nx,ny,nz)     ! Physical height coordinate defined at
                               ! w-point of the staggered grid(m).


  REAL :: uprt (nx,ny,nz)      ! x component of perturbation velocity
                               ! (m/s)
  REAL :: vprt (nx,ny,nz)      ! y component of perturbation velocity
                               ! (m/s)
  REAL :: wprt (nx,ny,nz)      ! vertical component of perturbation
                               ! velocity in Cartesian coordinates
                               ! (m/s).
  REAL :: ptprt(nx,ny,nz)      ! perturbation potential temperature (K)
  REAL :: pprt (nx,ny,nz)      ! perturbation pressure (Pascal)
  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)      ! Rainwater 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 :: tke  (nx,ny,nz)      ! Turbulent Kinetic Energy ((m/s)**2)

  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: ubar  (nx,ny,nz)     ! Base state x velocity component (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state y velocity component (m/s)
  REAL :: wbar  (nx,ny,nz)     ! Base state z velocity component (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor mixing ratio
                               ! (kg/kg)

  INTEGER :: nstyps             ! Number of soil type

  INTEGER :: soiltyp (nx,ny,nstyps)    ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)    ! Soil type
  INTEGER :: vegtyp (nx,ny)            ! Vegetation type
  REAL :: lai    (nx,ny)            ! Leaf Area Index
  REAL :: roufns (nx,ny)            ! Surface roughness
  REAL :: veg    (nx,ny)            ! Vegetation fraction

  REAL :: tsfc   (nx,ny,0:nstyps)      ! Temperature at surface (K)
  REAL :: tsoil  (nx,ny,0:nstyps)      ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyps)      ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyps)      ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps)      ! Canopy water amount
  REAL :: snowdpth(nx,ny)              ! Snow depth (m)

  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! precipitation rate (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: rnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: tem1(nx,ny,nz)       ! Temporary work array
  REAL :: tem2(nx,ny,nz)       ! Temporary work array
  REAL :: tem3(nx,ny,nz)       ! Temporary work array

  INTEGER :: grdread,iread
  SAVE grdread

  INTEGER :: nchanl            ! FORTRAN I/O channel number for output
  INTEGER :: istat
  INTEGER :: ireturn           ! Return status indicator
  INTEGER :: grdbas            ! Wether this is a grid/base state
                               ! array dump
  REAL :: btime                ! The time of the base state data
  REAL :: p0inv, amin, amax
  INTEGER :: i,j,k
  INTEGER :: ierr
  LOGICAL :: fexist
  INTEGER :: packed
  CHARACTER (LEN=80) :: rname
  INTEGER :: is
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'indtflg.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'mp.inc'            ! mpi parameters.
  INCLUDE 'phycst.inc'

  DATA grdread /0/
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ireturn = 0
  p0inv=1./p0

  IF (lengbf < 1) lengbf = len_trim(grdbasfn)
  IF (lendtf < 1) lendtf = len_trim(datafn)
!
!-----------------------------------------------------------------------
!
!  Open and read grid and base state data file depending on the
!  values of parameters grdin and basin, which are read in from the
!  time dependent data set. If grdin or basin is zero, the grid and
!  base state arrays have to be read in from a separate file.
!
!-----------------------------------------------------------------------
!
  IF ( hinfmt == 9 ) GO TO 500

  IF( grdread == 0 ) THEN

    grdbas = 1
    rname = runname

    INQUIRE(FILE=trim(grdbasfn(1:lengbf)), EXIST = fexist )
    IF( fexist ) GO TO 200

    INQUIRE(FILE=trim(grdbasfn(1:lengbf))//'.Z', EXIST = fexist )
    IF( fexist ) THEN
      CALL uncmprs( trim(grdbasfn(1:lengbf))//'.Z' )
      GO TO 200
    END IF

    INQUIRE(FILE=trim(grdbasfn(1:lengbf))//'.gz', EXIST = fexist )
    IF( fexist ) THEN
      CALL uncmprs( trim(grdbasfn(1:lengbf))//'.gz' )
      GO TO 200
    END IF

    WRITE(6,'(/1x,a,/1x,a/)')                                           &
        'File '//trim(grdbasfn(1:lengbf))//                             &
        ' or its compressed version not found.',                        &
        'Program stopped in DTAREAD.'
    CALL arpsstop('arpsstop called from dtaread during base state read',1)

    200     CONTINUE

!
!-----------------------------------------------------------------------
!
!  Read grid and base state fields.
!
!-----------------------------------------------------------------------
!
    IF( hinfmt == 1 ) THEN

      CALL getunit( nchanl )
!
!-----------------------------------------------------------------------
!
!  Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!

      CALL asnctl ('NEWLOCAL', 1, ierr)
      CALL asnfile(trim(grdbasfn(1:lengbf)), '-F f77 -N ieee', ierr)

      OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))),                   &
           STATUS='old',FORM='unformatted',IOSTAT=istat)

      IF( istat /= 0 ) GO TO 999

      CALL binread(nx,ny,nz,nstyps, grdbas, nchanl, btime, x,y,z,zp,    &
                   uprt ,vprt ,wprt ,ptprt, pprt,                       &
                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
                   ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,        &
                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
                   raing,rainc,prcrate,                                 &
                   radfrc,radsw,rnflx,                                  &
                   usflx,vsflx,ptsflx,qvsflx,                           &
                   ireturn)

      CLOSE(UNIT=nchanl)
      CALL retunit( nchanl )

    ELSE IF( hinfmt == 2 ) THEN

      CALL getunit( nchanl )
      OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))),                   &
           STATUS='old',FORM='formatted',IOSTAT=istat)

      IF( istat /= 0 ) GO TO 999

      CALL ascread(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp,        &
                   uprt ,vprt ,wprt ,ptprt, pprt,                       &
                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
                   ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,        &
                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
                   raing,rainc,prcrate,                                 &
                   radfrc,radsw,rnflx,                                  &
                   usflx,vsflx,ptsflx,qvsflx,                           &
                   ireturn)

      CLOSE(UNIT=nchanl)
      CALL retunit( nchanl )

    ELSE IF( hinfmt == 3 ) THEN

      CALL hdfread(nx,ny,nz,nstyps,grdbas, trim(grdbasfn(1:lengbf)),          &
                   btime,x,y,z,zp,                                      &
                   uprt ,vprt ,wprt ,ptprt, pprt,                       &
                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
                   ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,        &
                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
                   raing,rainc,prcrate,                                 &
                   radfrc,radsw,rnflx,                                  &
                   usflx,vsflx,ptsflx,qvsflx,                           &
                   ireturn, tem1)

    ELSE IF( hinfmt == 4 ) THEN

      CALL getunit( nchanl )
      OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))),                   &
           STATUS='old',FORM='unformatted',IOSTAT=istat)

      IF( istat /= 0 ) GO TO 999

!      CALL pakread(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp,        &
!                   uprt ,vprt ,wprt ,ptprt, pprt,                       &
!                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
!                   ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,        &
!                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
!                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
!                   raing,rainc,prcrate,                                 &
!                   radfrc,radsw,rnflx,                                  &
!                   usflx,vsflx,ptsflx,qvsflx,                           &
!                   ireturn, tem1,tem2)
      CALL pakread

      CLOSE(UNIT=nchanl)
      CALL retunit( nchanl )

    ELSE IF( hinfmt == 6 ) THEN

      CALL getunit( nchanl )
      OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))),                   &
           STATUS='old',FORM='unformatted',IOSTAT=istat)

      IF( istat /= 0 ) GO TO 999

      CALL bn2read(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp,        &
                   uprt ,vprt ,wprt ,ptprt, pprt,                       &
                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
                   ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,        &
                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
                   raing,rainc,prcrate,                                 &
                   radfrc,radsw,rnflx,                                  &
                   usflx,vsflx,ptsflx,qvsflx,                           &
                   ireturn)

      CLOSE(UNIT=nchanl)
      CALL retunit( nchanl )

    ELSE IF (hinfmt == 7) THEN     ! NetCDF format

      packed = 0
!      CALL netread(trim(grdbasfn(1:lengbf)),nx,ny,nz,nstyps,packed,           &
!                   grdbas, x, y, z, zp,                                 &
!                   uprt ,vprt ,wprt ,ptprt, pprt ,                      &
!                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
!                   ubar, vbar, wbar, ptbar, pbar,                       &
!                   rhobar, qvbar,                                       &
!                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
!                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
!                   raing,rainc,prcrate,                                 &
!                   radfrc,radsw,rnflx,                                  &
!                   usflx,vsflx,ptsflx,qvsflx,                           &
!                   tem1, tem2, tem3)
      CALL netread

    ELSE IF (hinfmt == 8) THEN     ! NetCDF packed format

      packed = 1
!      CALL netread(trim(grdbasfn(1:lengbf)),nx,ny,nz,nstyps,packed,           &
!                   grdbas, x, y, z, zp,                                 &
!                   uprt ,vprt ,wprt ,ptprt, pprt ,                      &
!                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
!                   ubar, vbar, wbar, ptbar, pbar,                       &
!                   rhobar, qvbar,                                       &
!                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
!                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
!                   raing,rainc,prcrate,                                 &
!                   radfrc,radsw,rnflx,                                  &
!                   usflx,vsflx,ptsflx,qvsflx,                           &
!                   tem1, tem2, tem3)
      CALL netread

    ELSE IF( hinfmt == 10 ) THEN
!
!-----------------------------------------------------------------------
!
!  Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
      CALL asnctl ('NEWLOCAL', 1, ierr)
      CALL asnfile(trim(grdbasfn(1:lengbf)), '-F f77 -N ieee', ierr)

      CALL getunit( nchanl )

      OPEN(UNIT=nchanl,FILE=trim(trim(grdbasfn(1:lengbf))),STATUS='old',      &
           FORM='unformatted',IOSTAT= istat )

      CALL gribread(nx,ny,nz,nstyps,grdbas,nchanl,btime,x,y,z,zp,       &
                   uprt ,vprt ,wprt ,ptprt, pprt ,                      &
                   qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,              &
                   ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,        &
                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
                   raing,rainc,prcrate,                                 &
                   radfrc,radsw,rnflx,                                  &
                   usflx,vsflx,ptsflx,qvsflx)

      CLOSE(UNIT=nchanl)
      CALL retunit( nchanl )

    ELSE

      WRITE(6,'(a,i3,a)')                                               &
          ' Data format flag had an invalid value ',                    &
            hinfmt ,' program stopped.'
      CALL arpsstop('arpsstop called from dtaread during base state     &
        &read',1)

    END IF

    grdread = 1

    runname = rname

  END IF

  500   CONTINUE

!
!-----------------------------------------------------------------------
!
!  Read data fields.
!
!-----------------------------------------------------------------------
!
  grdbas = 0

  INQUIRE(FILE=trim(datafn(1:lendtf)), EXIST = fexist )
  IF( fexist ) GO TO 100

  INQUIRE(FILE=trim(datafn(1:lendtf))//'.Z', EXIST = fexist )
  IF( fexist ) THEN
    CALL uncmprs( trim(datafn(1:lendtf))//'.Z' )
    GO TO 100
  END IF

  INQUIRE(FILE=trim(datafn(1:lendtf))//'.gz', EXIST = fexist )
  IF( fexist ) THEN
    CALL uncmprs( trim(datafn(1:lendtf))//'.gz' )
    GO TO 100
  END IF

  WRITE(6,'(/1x,a,/1x,a/)')                                             &
       'File '//trim(datafn(1:lendtf))                                  &
       //' or its compressed version not found.',                       &
       'Program stopped in DTAREAD.'
  CALL arpsstop('arpsstop called from dtaread during base read-2',1)

  100   CONTINUE

  IF( hinfmt == 1 ) THEN

!
!-----------------------------------------------------------------------
!
!  Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!

    CALL asnctl ('NEWLOCAL', 1, ierr)
    CALL asnfile(trim(datafn(1:lendtf)), '-F f77 -N ieee', ierr)

    CALL getunit( nchanl )
    OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))),                       &
         STATUS='old',FORM='unformatted',IOSTAT=istat)

    IF( istat /= 0 ) GO TO 998

    CALL binread(nx,ny,nz,nstyps,grdbas, nchanl, time, x,y,z,zp,        &
                 uprt ,vprt ,wprt ,ptprt, pprt,                         &
                 qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                &
                 ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,          &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 ireturn)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE IF( hinfmt == 2 ) THEN

    CALL getunit( nchanl )
    OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))),                       &
         STATUS='old',FORM='formatted',IOSTAT=istat)

    IF( istat /= 0 ) GO TO 998

    CALL ascread(nx,ny,nz,nstyps,grdbas,nchanl,time,x,y,z,zp,           &
                 uprt ,vprt ,wprt ,ptprt, pprt,                         &
                 qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                &
                 ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,          &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 ireturn)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE IF( hinfmt == 3 ) THEN

    CALL hdfread(nx,ny,nz,nstyps,grdbas,trim(datafn(1:lendtf)),               &
                 time,x,y,z,zp,                                         &
                 uprt ,vprt ,wprt ,ptprt, pprt,                         &
                 qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                &
                 ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,          &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 ireturn, tem1)

  ELSE IF( hinfmt == 4 ) THEN

    CALL getunit( nchanl )
    OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))),                       &
         STATUS='old',FORM='unformatted',IOSTAT=istat)

    IF( istat /= 0 ) GO TO 998

!    CALL pakread(nx,ny,nz,nstyps,grdbas,nchanl,time,x,y,z,zp,           &
!                 uprt ,vprt ,wprt ,ptprt, pprt,                         &
!                 qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                &
!                 ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,          &
!                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
!                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
!                 raing,rainc,prcrate,                                   &
!                 radfrc,radsw,rnflx,                                    &
!                 usflx,vsflx,ptsflx,qvsflx,                             &
!                 ireturn, tem1,tem2)
    CALL pakread

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE IF( hinfmt == 6 ) THEN

    CALL getunit( nchanl )
    OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))),                       &
            STATUS='old',FORM='unformatted',IOSTAT=istat)

    IF( istat /= 0 ) GO TO 998

    CALL bn2read(nx,ny,nz,nstyps,grdbas, nchanl, time, x,y,z,zp,        &
                 uprt ,vprt ,wprt ,ptprt, pprt ,                        &
                 qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                &
                 ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,          &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 ireturn)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE IF (hinfmt == 7) THEN     ! NetCDF format

    packed = 0
!    CALL netread(trim(datafn(1:lendtf)),nx,ny,nz,nstyps,packed,               &
!                 grdbas, x, y, z, zp,                                   &
!                 uprt ,vprt ,wprt ,ptprt, pprt ,                        &
!                 qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                &
!                 ubar, vbar, wbar, ptbar, pbar,                         &
!                 rhobar, qvbar,                                         &
!                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
!                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
!                 raing,rainc,prcrate,                                   &
!                 radfrc,radsw,rnflx,                                    &
!                 usflx,vsflx,ptsflx,qvsflx,                             &
!                 tem1, tem2, tem3)
    CALL netread


  ELSE IF (hinfmt == 8) THEN     ! NetCDF packed format

    packed = 1
!    CALL netread(trim(datafn(1:lendtf)),nx,ny,nz,nstyps,packed,               &
!                 grdbas, x, y, z, zp,                                   &
!                 uprt ,vprt ,wprt ,ptprt, pprt ,                        &
!                 qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,                &
!                 ubar, vbar, wbar, ptbar, pbar,                         &
!                 rhobar, qvbar,                                         &
!                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
!                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
!                 raing,rainc,prcrate,                                   &
!                 radfrc,radsw,rnflx,                                    &
!                 usflx,vsflx,ptsflx,qvsflx,                             &
!                 tem1, tem2, tem3)
    CALL netread

  ELSE IF( hinfmt == 9 ) THEN
!
!-----------------------------------------------------------------------
!
!  hinfmt = 9, read GrADS format data file. Since the DTAREAD
!  subroutine doesn't define u, v, w, pt, p, qv, we have to use
!  those base variables such ubar, vbar, wbar, ptbar, pbar, and
!  qvbar as temporary arrays to store u, v, w, pt, p, qv.
!
!-----------------------------------------------------------------------
!
    CALL gradsread(nx,ny,nz,nstyps,                                     &
                   nchanl, trim(datafn(1:lendtf)), time,                      &
                   x,y,z,zp, uprt,vprt,wprt,ptprt,pprt,                 &
                   qvprt,qc,qr,qi,qs,qh, tke,kmh,kmv,                   &
                   ubar,vbar,ptbar,pbar,rhobar,qvbar,                   &
                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
                   raing,rainc,prcrate,                                 &
                   radfrc,radsw,rnflx,                                  &
                   usflx,vsflx,ptsflx,qvsflx,                           &
                   ireturn, tem1)

    DO k = 1, nz
      DO j = 1, ny
        DO i = 1, nx
          ubar (i,j,k) = ubar (i,j,k) - uprt (i,j,k)
          vbar (i,j,k) = vbar (i,j,k) - vprt (i,j,k)
          wbar (i,j,k) = 0.0
          pbar (i,j,k) = pbar (i,j,k) - pprt (i,j,k)
          ptbar(i,j,k) = ptbar(i,j,k) - ptprt(i,j,k)
          qvbar(i,j,k) = qvbar(i,j,k) - qvprt(i,j,k)
        END DO
      END DO
    END DO

    grdin = 1

  ELSE IF( hinfmt == 10 ) THEN
!
!-----------------------------------------------------------------------
!
!  hinfmt = 10, read GRIB format data file.
!
!-----------------------------------------------------------------------
!
    CALL asnctl ('NEWLOCAL', 1, ierr)
    CALL asnfile(trim(datafn(1:lendtf)), '-F f77 -N ieee', ierr)

    CALL getunit( nchanl )

    OPEN(UNIT=nchanl,FILE=trim(trim(datafn(1:lendtf))),STATUS='old',          &
         FORM='unformatted',IOSTAT= istat )

    CALL gribread(nx,ny,nz,nstyps, grdbas, nchanl, time, x,y,z,zp,      &
                  uprt ,vprt ,wprt ,ptprt, pprt ,                       &
                  qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv,               &
                  ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,         &
                  soiltyp,stypfrct,vegtyp,lai,roufns,veg,               &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,             &
                  raing,rainc,prcrate,                                  &
                  radfrc,radsw,rnflx,                                   &
                  usflx,vsflx,ptsflx,qvsflx)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE

    WRITE(6,'(a,i3,a)')                                                 &
          ' Data format flag had an invalid value ',                    &
           hinfmt ,' program stopped.'
    CALL arpsstop('arpsstop called from dtaread during read',1)
  END IF

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        rhobar(i,j,k)= pbar(i,j,k)/                                     &
               ( rd * ptbar(i,j,k)*(pbar(i,j,k)*p0inv)**rddcp )
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Set dx, dy and dz from read-in data
!
!-----------------------------------------------------------------------
!
  dx = x(2) - x(1)
  dy = y(2) - y(1)
  dz = z(2) - z(1)
!
  IF (myproc == 0)  &
     WRITE(6,'(/1x,a/)') 'Min. and max. of the data arrays read in:'

  CALL edgfill(x,1,nx,1,nx,1,1,1,1,1,1,1,1)
  CALL a3dmax0lcl(x,1,nx,1,nx,1,1,1,1,1,1,1,1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(/1x,2(a,e13.6))') 'xmin    = ', amin,',  xmax    =',amax

  CALL edgfill(y,1,ny,1,ny,1,1,1,1,1,1,1,1)
  CALL a3dmax0lcl(y,1,ny,1,ny,1,1,1,1,1,1,1,1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'ymin    = ', amin,',  ymax    =',amax

  CALL edgfill(z,1,nz,1,nz,1,1,1,1,1,1,1,1)
  CALL a3dmax0lcl(z,1,nz,1,nz,1,1,1,1,1,1,1,1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'zmin    = ', amin,',  zmax    =',amax

  CALL edgfill(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz)
  CALL a3dmax0lcl(zp,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'zpmin   = ', amin,',  zpmax   =',amax

  CALL edgfill(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1)
  CALL a3dmax0lcl(ubar,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'ubarmin = ', amin,',  ubarmax =',amax

  CALL edgfill(vbar,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1)
  CALL a3dmax0lcl(vbar,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'vbarmin = ', amin,',  vbarmax =',amax

  CALL edgfill(wbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz)
  CALL a3dmax0lcl(wbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'wbarmin = ', amin,',  wbarmax =',amax

  CALL edgfill(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
  CALL a3dmax0lcl(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'ptbarmin= ', amin,',  ptbarmax=',amax

  CALL edgfill(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
  CALL a3dmax0lcl(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'pbarmin = ', amin,',  pbarmax =',amax

  IF( mstin == 1) THEN

    CALL edgfill(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.6))') 'qvbarmin= ', amin, ',  qvbarmax=',amax
  END IF

  CALL edgfill(uprt,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1)
  CALL a3dmax0lcl(uprt,1,nx,1,nx,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'uprtmin = ', amin,',  uprtmax =',amax

  CALL edgfill(vprt,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1)
  CALL a3dmax0lcl(vprt,1,nx,1,nx-1,1,ny,1,ny,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'vprtmin = ', amin,',  vprtmax =',amax

  CALL edgfill(wprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz)
  CALL a3dmax0lcl(wprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'wprtmin = ', amin,',  wprtmax =',amax

  CALL edgfill(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
  CALL a3dmax0lcl(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'ptprtmin= ', amin,',  ptprtmax=',amax

  CALL edgfill(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
  CALL a3dmax0lcl(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
  IF (myproc == 0)  &
     WRITE(6,'(1x,2(a,e13.6))') 'pprtmin = ', amin,',  pprtmax =',amax

  IF( mstin == 1) THEN

    CALL edgfill(qvprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(qvprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.6))') 'qvprtmin= ', amin,                      &
                              ',  qvprtmax=',amax

    CALL edgfill(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.6))') 'qcmin   = ', amin,                      &
                               ',  qcmax   =',amax

    CALL edgfill(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.6))') 'qrmin   = ', amin,                      &
                               ',  qrmax   =',amax

    IF( icein == 1) THEN

      CALL edgfill(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
      CALL a3dmax0lcl(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))') 'qimin   = ', amin,                    &
                                 ',  qimax   =',amax

      CALL edgfill(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
      CALL a3dmax0lcl(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))') 'qsmin   = ', amin,                    &
                                 ',  qsmax   =',amax

      CALL edgfill(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
      CALL a3dmax0lcl(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))') 'qhmin   = ', amin,                    &
                                 ',  qhmax   =',amax

    END IF

    IF(rainin == 1) THEN

      CALL edgfill(raing,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(raing,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'raingmin= ', amin,',  raingmax=',amax

      CALL edgfill(rainc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(rainc,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'raincmin= ', amin,',  raincmax=',amax

    END IF

    IF( prcin == 1 ) THEN

      CALL edgfill(prcrate(1,1,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(prcrate(1,1,1),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'prcr1min= ', amin,',  prcr1max=',amax

      CALL edgfill(prcrate(1,1,2),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(prcrate(1,1,2),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'prcr2min= ', amin,',  prcr2max=',amax

      CALL edgfill(prcrate(1,1,3),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(prcrate(1,1,3),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'prcr3min= ', amin,',  prcr3max=',amax

      CALL edgfill(prcrate(1,1,4),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(prcrate(1,1,4),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF (myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'prcr4min= ', amin,',  prcr4max=',amax

    END IF

  END IF

  IF (tkein == 1) THEN

    CALL edgfill(tke,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(tke,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'tkemin  = ', amin,',  tkemax  =',amax

  END IF

  IF (trbin == 1) THEN

    CALL edgfill(kmh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(kmh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
      WRITE(6,'(1x,2(a,e13.4))')                                          &
        'kmhmin  = ', amin,',  kmhmax  =',amax

    CALL edgfill(kmv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(kmv,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
      WRITE(6,'(1x,2(a,e13.4))')                                          &
        'kmvmin  = ', amin,',  kmvmax  =',amax

  END IF

  IF( sfcin == 1 ) THEN

    DO is = 0, nstyp

      IF(myproc == 0)  &
         PRINT*,'Max/min of soil model variables for type index ', is

      CALL edgfill(tsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(tsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF(myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'tsfcmin = ', amin,',  tsfcmax =',amax

      CALL edgfill(tsoil(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(tsoil(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF(myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'tsoilmin= ', amin,',  tsoilmax=',amax

      CALL edgfill(wetsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(wetsfc(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF(myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'wetsmin = ', amin,',  wetsmax =',amax

      CALL edgfill(wetdp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(wetdp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF(myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'wetdmin = ', amin,',  wetdmax =',amax

      CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
      CALL a3dmax0lcl(wetcanp(1,1,is),1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
      IF(myproc == 0)  &
         WRITE(6,'(1x,2(a,e13.6))')                                        &
          'wetcmin = ', amin,',  wetcmax =',amax

    END DO

  END IF

  IF ( radin == 1 ) THEN

    CALL edgfill(radfrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
    CALL a3dmax0lcl(radfrc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'radfmin = ', amin,',  radfmax =',amax

    CALL edgfill(radsw,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
    CALL a3dmax0lcl(radsw,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'radswmin= ', amin,',  radswmax=',amax

    CALL edgfill(rnflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
    CALL a3dmax0lcl(rnflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'rnflxmin= ', amin,',  rnflxmax=',amax

  END IF

  IF ( flxin == 1 ) THEN

    CALL edgfill(usflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
    CALL a3dmax0lcl(usflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'usflxmin= ', amin,',  usflxmax=',amax

    CALL edgfill(vsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
    CALL a3dmax0lcl(vsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'vsflxmin= ', amin,',  vsflxmax=',amax

    CALL edgfill(ptsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
    CALL a3dmax0lcl(ptsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'ptflxmin= ', amin,',  ptflxmax=',amax

    CALL edgfill(qvsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1)
    CALL a3dmax0lcl(qvsflx,1,nx,1,nx-1,1,ny,1,ny-1,1,1,1,1,amax,amin)
    IF (myproc == 0)  &
       WRITE(6,'(1x,2(a,e13.4))')                                          &
        'qvflxmin= ', amin,',  qvflxmax=',amax

  END IF

  RETURN


!##################################################################
!##################################################################
!######                                                      ######
!######                ENTRY SETGBRD                         ######
!######                                                      ######
!##################################################################
!##################################################################

  ENTRY setgbrd (iread)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Set the gridread parameter so that the grid/base state file
!  will (iread=0) or will not (iread=1) be read. This is useful
!  if more than one history file is to be accessed during a run
!  (primarily in data assimilation experiments).
!
!  As a default, the grid/base file is read on the first call to
!  dtaread and not on subsequent calls.  dtaread resets grdread
!  to 1 after each time a grid/base file is read.
!
!-----------------------------------------------------------------------
!

  grdread=iread

  RETURN

!  997   CONTINUE
  WRITE(6,'(1x,a,a,/1x,i3,a)')                                          &
      'Error occured when opening GRIB file ',trim(datafn(1:lendtf)),         &
      'with the FILE pointer ',nchanl,' Program stopped in DTAREAD.'
  CALL arpsstop('arpsstop called from dtaread during Grid read',1)

  998   CONTINUE
  WRITE(6,'(1x,a,a,/1x,i3,a)')                                          &
      'Error occured when opening file ',trim(datafn(1:lendtf)),              &
      'using FORTRAN unit ',nchanl,' Program stopped in DTAREAD.'
  CALL arpsstop('arpsstop called from dtaread during file read',1)

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

  CALL arpsstop('arpsstop called from dtaread during file read',1)

END SUBROUTINE dtaread
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE DTADUMP                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE dtadump(nx,ny,nz,nstyps,                                     & 24,35
           houtfmt,nchout,filnam,grdbas,flcmprs,                        &
           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv,              &
           ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,                      &
           x,y,z,zp,hterain, j1,j2,j3,                                  &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,                    &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,                                          &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           tem1,tem2,tem3)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Coordinate the history data dump in various data formats by
!  calling the appropriate history dump routine.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  8/23/1992
!
!  MODIFICATION HISTORY:
!
!  4/4/93  (M. Xue)
!  Modified, so that data on the original staggered grid are written
!  out. Averaging to the volume center is no longer done.
!
!  9/1/94 (Y. Lu)
!  Cleaned up documentation.
!
!  4/22/1995 (M. Xue)
!  Modified external boundary condition file naming format, so
!  that the files can be directed to directory 'dirname'.
!
!  4/24/1995 (M. Xue)
!  Changed boundary file naming convention.
!
!  12/07/95 (Yuhe Liu)
!  Added a new data dumping format, GRIB. The parameter definitions
!  is mostly in consistance with NMC standard, despite those
!  variables which are not in the WMO and NMC's tables (Table 2)
!
!  01/22/1996 (Yuhe Liu)
!  Changed the ARPS history dumpings to dump the total values of
!  ARPS variables, instead of the perturbated values.
!
!  04/30/1997 (Fanyou Kong -- CMRP)
!  Add Vis5D format conversion
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    nstyps   number of soil type
!
!    houtfmt  The format flag of history data dump
!    nchout   FORTRAN I/O channel number for history data output.
!    filnam   The name of history data dump file
!    grdbas   Flag indicating if this is a call for the data dump
!             of grid and base state arrays only. If so, grdbas=1.
!    flcmprs  Option for file compression (0 or 1), not used in this
!             routine. A identical variable, filcmprs, is passed
!             through include file globcst.inc
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        Vertical component of Cartesian velocity at a given
!             time level (m/s)
!    ptprt    Perturbation potential temperature at a given time
!             level (K)
!    pprt     Perturbation pressure at  a given time level (Pascal)
!    qv       Water vapor specific humidity at a given time level (kg/kg)
!    qc       Cloud water mixing ratio at a given time level (kg/kg)
!    qr       Rainwater mixing ratio at a given time level (kg/kg)
!    qi       Cloud ice mixing ratio at a given time level (kg/kg)
!    qs       Snow mixing ratio at a given time level (kg/kg)
!    qh       Hail mixing ratio at a given time level (kg/kg)
!    tke      Turbulent Kinetic Energy ((m/s)**2)
!
!    kmh      Horizontal turb. mixing coef. for momentum ( m**2/s )
!    kmv      Vertical turb. mixing coef. for momentum ( m**2/s )
!
!    ubar     Base state x-velocity component (m/s)
!    vbar     Base state y-velocity component (m/s)
!    wbar     Base state vertical velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhobar   Base state density (kg/m**3)
!    qvbar    Base state water vapor specific humidity (kg/kg)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!    hterain  Terrain height (m)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at ground (K) (in top 1 cm layer)
!    tsoil    Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture in the top 1 cm layer
!    wetdp    Deep soil moisture in the deep 1 m layer
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  OUTPUT:
!
!    None.
!
!  WORK ARRAY:
!
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

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

  INTEGER :: houtfmt           ! Flag of history data dump format
  INTEGER :: nchout            ! FORTRAN I/O channel number for output
  CHARACTER (LEN=80) :: filnam       ! The name of history dump data file
  INTEGER :: grdbas            ! If this is a grid/base state array dump
  INTEGER :: flcmprs           ! Option for file compression (0 or 1).

  REAL :: u     (nx,ny,nz)     ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz)     ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz)     ! Total w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz)     ! Perturbation pressure (Pascal)

  REAL :: qv    (nx,ny,nz)     ! Water vapor specific humidity (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 :: tke   (nx,ny,nz)     ! Turbulent Kinetic Energy ((m/s)**2)

  REAL :: kmh   (nx,ny,nz)     ! Horizontal turb. mixing coef. for
                               ! momentum. ( m**2/s )
  REAL :: kmv   (nx,ny,nz)     ! Vertical turb. mixing coef. for
                               ! momentum. ( m**2/s )

  REAL :: ubar  (nx,ny,nz)     ! Base state x-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state y-velocity (m/s)
  REAL :: wbar  (nx,ny,nz)     ! Base state z-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)

  REAL :: x     (nx)           ! The x-coord. of the physical and
                               ! computational grid. Defined at u-point.
  REAL :: y     (ny)           ! The y-coord. of the physical and
                               ! computational grid. Defined at v-point.
  REAL :: z     (nz)           ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of the staggered grid.

  REAL :: hterain(nx,ny)       ! Terrain height.

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as - d( zp )/d( x )
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as - d( zp )/d( y )
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z )

  INTEGER :: nstyps                  ! Number of soil types
  INTEGER :: soiltyp(nx,ny,nstyps)   ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)  ! Soil type fractions
  INTEGER :: vegtyp (nx,ny)          ! Vegetation type
  REAL :: lai    (nx,ny)          ! Leaf Area Index
  REAL :: roufns (nx,ny)          ! Surface roughness
  REAL :: veg    (nx,ny)          ! Vegetation fraction

  REAL :: tsfc   (nx,ny,0:nstyps)    ! Temperature at surface (K)
  REAL :: tsoil  (nx,ny,0:nstyps)    ! Deep soil temperature (K)
  REAL :: wetsfc (nx,ny,0:nstyps)    ! Surface soil moisture
  REAL :: wetdp  (nx,ny,0:nstyps)    ! Deep soil moisture
  REAL :: wetcanp(nx,ny,0:nstyps)    ! Canopy water amount
  REAL :: snowdpth(nx,ny)            ! Snow depth (m)

  REAL :: raing(nx,ny)         ! Grid supersaturation rain
  REAL :: rainc(nx,ny)         ! Cumulus convective rain
  REAL :: prcrate(nx,ny,4)     ! precipitation rates (kg/(m**2*s))
                               ! prcrate(1,1,1) = total precip. rate
                               ! prcrate(1,1,2) = grid scale precip. rate
                               ! prcrate(1,1,3) = cumulus precip. rate
                               ! prcrate(1,1,4) = microphysics precip. rate

  REAL :: radfrc(nx,ny,nz)     ! Radiation forcing (K/s)
  REAL :: radsw (nx,ny)        ! Solar radiation reaching the surface
  REAL :: rnflx (nx,ny)        ! Net radiation flux absorbed by surface

  REAL :: usflx (nx,ny)        ! Surface flux of u-momentum (kg/(m*s**2))
  REAL :: vsflx (nx,ny)        ! Surface flux of v-momentum (kg/(m*s**2))
  REAL :: ptsflx(nx,ny)        ! Surface heat flux (K*kg/(m*s**2))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array
  REAL :: tem2  (nx,ny,nz)     ! Temporary work array
  REAL :: tem3  (nx,ny,nz)     ! Temporary work array

  CHARACTER (LEN=256) :: filnamr

  INTEGER :: istat, nchout0, nchout1
  INTEGER :: ierr
  INTEGER :: packed

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
  REAL :: zpmax
!
!-----------------------------------------------------------------------
!
!  Declaration for dumping the external boundary variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i, j, k
  INTEGER :: itema, lenstr
  INTEGER :: iyr, imon, idy, ihr, imin, isec
  CHARACTER (LEN=15) :: ctime
  CHARACTER (LEN=80) :: exbcfn
  CHARACTER (LEN=80) :: temchar
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'bndry.inc'
  INCLUDE 'mp.inc'            ! Message passing parameters.
!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
  INCLUDE 'exbc.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  IF (sfcout == 1) snowout = 1   ! snowout only independetly zero
                                 ! for pre-snow cover versions.
!
!  IF( houtfmt.eq.0 ) RETURN   ! No data dump
!
!
  IF( houtfmt == 1 ) THEN     ! Unformatted binary data dump.
!
!-----------------------------------------------------------------------
!
!  History dump of unformatted binary data
!
!-----------------------------------------------------------------------
!

    CALL getunit( nchout )

!
!-----------------------------------------------------------------------
!
!  Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!

    CALL asnctl ('NEWLOCAL', 1, ierr)
    CALL asnfile(filnam, '-F f77 -N ieee', ierr)

    OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new',                    &
         FORM='unformatted',IOSTAT= istat )

    IF( istat /= 0) GO TO 999

    CALL bindump(nx,ny,nz,nstyps,nchout, grdbas,                        &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
                 kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar,            &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 tem1)

    CLOSE(UNIT=nchout)
    CALL retunit( nchout )

  ELSE IF( houtfmt == 2 ) THEN     ! Formatted ASCII data dump
!
!-----------------------------------------------------------------------
!
!  History dump of formatted ASCII data
!
!-----------------------------------------------------------------------
!

    CALL getunit( nchout )
    OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new',                    &
         FORM='formatted',IOSTAT= istat )

    IF( istat /= 0) GO TO 999

    CALL ascdump(nx,ny,nz,nstyps, nchout, grdbas,                       &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
                 kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar,            &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 tem1)

    CLOSE(UNIT=nchout)
    CALL retunit( nchout )

  ELSE IF( houtfmt == 3 ) THEN    ! HDF4 format dump.
!
!-----------------------------------------------------------------------
!
!  History data dump in NCSA HDF4 file format.
!
!-----------------------------------------------------------------------
!

    CALL hdfdump(nx,ny,nz,nstyps,filnam, grdbas,                        &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
                 kmh,kmv,ubar,vbar,ptbar,pbar,rhobar,qvbar,             &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 tem1)

  ELSE IF( houtfmt == 4 ) THEN     ! Packed binary data dump.
!
!-----------------------------------------------------------------------
!
!  History dump of packed unformatted binary data
!
!-----------------------------------------------------------------------
!

    CALL getunit( nchout )
    OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new',                    &
         FORM='unformatted',IOSTAT= istat )

    IF( istat /= 0) GO TO 999

!    CALL pakdump(nx,ny,nz,nstyps,nchout, grdbas,                        &
!                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
!                 kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar,            &
!                 x,y,z,zp,hterain, j1,j2,j3,                            &
!                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
!                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
!                 raing,rainc,prcrate,                                   &
!                 radfrc,radsw,rnflx,                                    &
!                 usflx,vsflx,ptsflx,qvsflx,                             &
!                 tem1, tem2)
    CALL pakdump

    CLOSE(UNIT=nchout)
    CALL retunit( nchout )

  ELSE IF( houtfmt == 5 ) THEN     ! Data dump for Savi3D
!
!-----------------------------------------------------------------------
!
!  History dump for Savi3D.
!  Data at all time levels are dumped into a single file.
!
!-----------------------------------------------------------------------
!
    nchout0 = 79

    CALL svidump(nx,ny,nz,nstyps,nchout0,filnam, grdbas,                &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
                 kmh,kmv,ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,        &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 tem1,tem2)

  ELSE IF( houtfmt == 6 ) THEN     ! Unformatted binary data dump.
!
!-----------------------------------------------------------------------
!
!  History dump of unformatted binary data
!
!-----------------------------------------------------------------------
!

    CALL getunit( nchout )
    OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new',                    &
         FORM='unformatted',IOSTAT= istat )

    IF( istat /= 0) GO TO 999

    CALL bn2dump(nx,ny,nz,nstyps, nchout, grdbas,                       &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
                 kmh,kmv, ubar,vbar,ptbar,pbar,rhobar,qvbar,            &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 tem1)

    CLOSE(UNIT=nchout)
    CALL retunit( nchout )
!
!-----------------------------------------------------------------------
!
!  NetCDF format dump.
!
!-----------------------------------------------------------------------
!
  ELSE IF (houtfmt == 7) THEN    ! NetCDF format

    packed = 0
!    CALL netdump(filnam, nx,ny,nz,nstyps,packed, grdbas,                &
!                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
!                 kmh,kmv, ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,       &
!                 x,y,z,zp,hterain, j1,j2,j3,                            &
!                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
!                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
!                 raing,rainc,prcrate,                                   &
!                 radfrc,radsw,rnflx,                                    &
!                 usflx,vsflx,ptsflx,qvsflx,                             &
!                 tem1,tem2,tem3)
     CALL netdump
!
!-----------------------------------------------------------------------
!
!  Packed NetCDF format dump.
!
!-----------------------------------------------------------------------
!
  ELSE IF (houtfmt == 8) THEN    ! Packed NetCDF format

    packed = 1
!    CALL netdump(filnam, nx,ny,nz,nstyps,packed, grdbas,                &
!                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
!                 kmh,kmv, ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,       &
!                 x,y,z,zp,hterain, j1,j2,j3,                            &
!                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
!                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
!                 raing,rainc,prcrate,                                   &
!                 radfrc,radsw,rnflx,                                    &
!                 usflx,vsflx,ptsflx,qvsflx,                             &
!                 tem1,tem2,tem3)
    CALL netdump

  ELSE IF( houtfmt == 9 ) THEN     ! Unformatted binary data dump.
!
    CALL gradsdump(nx,ny,nz,nstyps,nchout, filnam, istager,             &
                   u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,              &
                   kmh,kmv,ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,      &
                   x,y,z,zp,hterain, j1,j2,j3,                          &
                   soiltyp,stypfrct,vegtyp,lai,roufns,veg,              &
                   tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,            &
                   raing,rainc,prcrate,                                 &
                   radfrc,radsw,rnflx,                                  &
                   usflx,vsflx,ptsflx,qvsflx,                           &
                   tem1,tem2)

  ELSE IF( houtfmt == 10 ) THEN     ! GRIB data dump.
!
!-----------------------------------------------------------------------
!
!  History dump of GRIB format data
!
!-----------------------------------------------------------------------
!
    lenstr = 80
    CALL strlnth( filnam, lenstr )

    CALL asnctl ('NEWLOCAL', 1, ierr)
    CALL asnfile(filnam, '-F f77 -N ieee', ierr)

    CALL getunit( nchout )

    OPEN(UNIT=nchout,FILE=trim(filnam),STATUS='new',                    &
         FORM='unformatted',IOSTAT= istat )

    IF( istat /= 0) GO TO 999
!
    CALL gribdump(nx,ny,nz,nstyps,                                      &
                  nchout, filnam(1:lenstr),grdbas,                      &
                  u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv,       &
                  ubar,vbar,ptbar,pbar,rhobar,qvbar,                    &
                  x,y,z,zp,hterain, j1,j2,j3,                           &
                  soiltyp,stypfrct,vegtyp,lai,roufns,veg,               &
                  tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,             &
                  raing,rainc,prcrate,                                  &
                  radfrc,radsw,rnflx,                                   &
                  usflx,vsflx,ptsflx,qvsflx,                            &
                  tem1,tem2,                                            &
                  tem3(1,1,1),tem3(1,1,2),                              &
                  tem3(1,1,3),tem3(1,1,4))

    CLOSE(UNIT=nchout)

    CALL retunit( nchout )

  ELSE IF( houtfmt == 11 ) THEN     ! Vis5D data dump.

    lenstr = 80
    CALL strlnth( filnam, lenstr )

    CALL v5ddump(nx,ny,nz, nstyps,filnam(1:lenstr), istager,            &
                 u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,                &
                 kmh,kmv,ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,        &
                 x,y,z,zp,hterain, j1,j2,j3,                            &
                 soiltyp,stypfrct,vegtyp,lai,roufns,veg,                &
                 tsfc,tsoil,wetsfc,wetdp,wetcanp,snowdpth,              &
                 raing,rainc,prcrate,                                   &
                 radfrc,radsw,rnflx,                                    &
                 usflx,vsflx,ptsflx,qvsflx,                             &
                 tem1,tem2)


  END IF   ! houtfmt switches
!
!-----------------------------------------------------------------------
!
!  Compress file except for Savi3D dump.
!
!-----------------------------------------------------------------------
!

  IF( houtfmt /= 0 .AND. houtfmt /= 5  .AND.                            &
      houtfmt /= 9 .AND. houtfmt /= 11 .AND. filcmprs == 1 )            &
      CALL cmprs( filnam )

!
!-----------------------------------------------------------------------
!
!  Create ready file, indicating history dump writing is complete
!
!-----------------------------------------------------------------------
!
  IF (readyfl == 1 .AND. houtfmt /= 0) THEN
    WRITE (filnamr,'(a)') trim(filnam) // "_ready"
    CALL getunit( nchout1 )
    OPEN (UNIT=nchout1,FILE=trim(filnamr))
    WRITE (nchout1,'(a)') trim(filnam)
    CLOSE (nchout1)
    CALL retunit ( nchout1 )
  END IF
!
!-----------------------------------------------------------------------
!
!  Write u,v,w,pt,p, and qv to an additional file for
!  external boundary conditions.
!
!-----------------------------------------------------------------------
!
  IF ( exbcdmp /= 0 .AND. grdbas /= 1 ) THEN
!
!-----------------------------------------------------------------------
!
!  Write the ARPS predicted fields into external boundary format
!  files. The purpose for doing so is to produce external boundary
!  conditions from ARPS. They can be used to test the external
!  BC code, as well as to run the experiments with the ARPS generated
!  external external boundary conditions.
!
!-----------------------------------------------------------------------
!
    CALL ctim2abss( year,month,day,hour,minute,second, itema )

    itema = itema + INT(curtim)

    CALL abss2ctim( itema,iyr,imon,idy,ihr,imin,isec )

    WRITE (ctime,'(i4.4,2i2.2,a,3i2.2)')                                &
          iyr,imon,idy,'.',ihr,imin,isec

    DO k = 1, nz-1
      DO j = 1, ny-1
        DO i = 1, nx-1
          tem1(i,j,k) = ptbar(i,j,k) + ptprt(i,j,k)
          tem2(i,j,k) = pbar (i,j,k) + pprt (i,j,k)
        END DO
      END DO
    END DO

    CALL edgfill(u,   1,nx,1,nx,   1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL edgfill(v,   1,nx,1,nx-1, 1,ny,1,ny,   1,nz,1,nz-1)
    CALL edgfill(w,   1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz  )
    CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL edgfill(tem2,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL edgfill(qv,  1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)

    exbcfn = runname(1:lfnkey)//'.'//ctime
    lenstr = lfnkey + 16

    IF (mp_opt > 0) THEN
      WRITE(exbcfn,'(a,a,a,a,2i2.2)')                                   &
          runname(1:lfnkey),'.',ctime,'_',loc_x,loc_y
      lenstr  = lenstr + 5
    END IF

    IF( dirname /= ' ' ) THEN

      temchar = exbcfn
      exbcfn  = dirname(1:ldirnam)//'/'//temchar
      lenstr  = lenstr + ldirnam + 1

    END IF

    CALL fnversn(exbcfn,lenstr)

    WRITE(6,'(1x,a,a)')                                                 &
         'Dumping to the external boundary format file: ',              &
         exbcfn(1:lenstr)

!wdt Copyright (c) 2001 Weather Decision Technologies, Inc.
    ! Calculate rayklow, needed for exbcdmp=4 format, if not computed already
    IF ( exbcdmp == 4 .and. rayklow <= 0) THEN
      rayklow = nz-1
      DO k=nz-1,2,-1
        zpmax = zp(1,1,k)
        DO j=1,ny-1
          DO i=1,nx-1
            zpmax = MAX( zp(i,j,k), zpmax )
          END DO
        END DO
        IF( zpmax < zbrdmp ) THEN
          rayklow = MAX(2, k+1)
          EXIT
        END IF
      END DO
    END IF

    CALL writexbc(nx,ny,nz,exbcfn,lenstr,ctime,                         &
                  1,1,1,1,1,1,                                          &
                  qcexout,qrexout,qiexout,qsexout,qhexout,              &
                  u,v,w,tem1,tem2,qv,qc,qr,qi,qs,qh)

  END IF

  RETURN

  999   CONTINUE
  WRITE(6,'(1x,a,a,a,/1x,i3)')                                          &
      'Error occured when opening file ',filnam,                        &
      'using FORTRAN unit ',nchout
  CALL arpsstop(' Program stopped in DTADUMP.',1)

END SUBROUTINE dtadump
!
!##################################################################
!##################################################################
!######                                                      ######
!######            SUBROUTINE GET_DIMS_FROM_DATA             ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE get_dims_from_data(hinfmt,grdbasfn, nx,ny,nz,nstyps, ireturn) 14,17
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  7/17/2000
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    hinfmt   The format of the history data dump
!    grdbasfn Name of the grid/base state array file
!
!  OUTPUT:
!
!    nx,ny,nz The dimension of data arrays
!    nstyps   The number of soil types
!    ireturn  Return status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: hinfmt            ! The format of the history data dump
  CHARACTER (LEN=*  ) :: grdbasfn  ! Name of the grid/base state array file
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  INTEGER :: nstyps            ! Number of soil types
  INTEGER :: ireturn,lengbf,nchanl,ierr,istat,packed,i
  LOGICAL :: fexist
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ireturn = 0
!
!-----------------------------------------------------------------------
!
!  Open and read grid and base state data file depending on the
!  values of parameters grdin and basin, which are read in from the
!  time dependent data set. If grdin or basin is zero, the grid and
!  base state arrays have to be read in from a separate file.
!
!-----------------------------------------------------------------------
!
!  IF ( hinfmt.eq.9 ) GOTO 500

  INQUIRE(FILE=trim(grdbasfn), EXIST = fexist )
  IF( fexist ) GO TO 200

  INQUIRE(FILE=trim(grdbasfn)//'.Z', EXIST = fexist )
  IF( fexist ) THEN
    CALL uncmprs( trim(grdbasfn)//'.Z' )
    GO TO 200
  END IF

  INQUIRE(FILE=trim(grdbasfn)//'.gz', EXIST = fexist )
  IF( fexist ) THEN
    CALL uncmprs( trim(grdbasfn)//'.gz' )
    GO TO 200
  END IF

  WRITE(6,'(/1x,a,/1x,a/)')                                             &
      'File '//trim(grdbasfn)//                                         &
      ' or its compressed version not found.',                          &
      'Program stopped in get_dims_from_data.'
  CALL arpsstop('arpsstop called from get_dims_from_data',1)

  200     CONTINUE

!
!-----------------------------------------------------------------------
!
!  Read grid and base state fields.
!
!-----------------------------------------------------------------------
!
  IF( hinfmt == 1 .OR. hinfmt == 6 ) THEN

    CALL getunit( nchanl )
!
!-----------------------------------------------------------------------
!
!  Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
    CALL asnctl ('NEWLOCAL', 1, ierr)
    CALL asnfile(trim(grdbasfn), '-F f77 -N ieee', ierr)

    OPEN(UNIT=nchanl,FILE=trim(grdbasfn),                           &
         STATUS='old',FORM='unformatted',IOSTAT=istat)

    IF( istat /= 0 ) GO TO 999

    CALL bin_getdims(nchanl, nx,ny,nz,nstyps, ireturn)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE IF( hinfmt == 2 ) THEN

    CALL getunit( nchanl )
    OPEN(UNIT=nchanl,FILE=trim(grdbasfn),                           &
         STATUS='old',FORM='formatted',IOSTAT=istat)

    IF( istat /= 0 ) GO TO 999

    CALL asc_getdims(nchanl, nx,ny,nz,nstyps, ireturn)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE IF( hinfmt == 3 ) THEN

    CALL hdf_getdims(trim(grdbasfn),nx,ny,nz,nstyps, ireturn)

  ELSE IF( hinfmt == 4 ) THEN

    CALL getunit( nchanl )
    OPEN(UNIT=nchanl,FILE=trim(grdbasfn),                           &
         STATUS='old',FORM='unformatted',IOSTAT=istat)

    IF( istat /= 0 ) GO TO 999

    CALL pak_getdims(nchanl, nx,ny,nz,nstyps, ireturn)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE IF (hinfmt == 7.OR.hinfmt == 8) THEN ! NetCDF format

    PRINT*,'NetCDF format is no longer supported by ARPS.'
    PRINT*,'Please check the specified inupt data format.'
    STOP

    packed = 0
!   CALL net_getdims(nchanl, packed, nx,ny,nz,nstyps, ireturn)

  ELSE IF( hinfmt == 10 ) THEN
!
!-----------------------------------------------------------------------
!
!  Cray routines to force binary data file to be in the IEEE format
!
!-----------------------------------------------------------------------
!
    CALL asnctl ('NEWLOCAL', 1, ierr)
    CALL asnfile(trim(grdbasfn), '-F f77 -N ieee', ierr)

    CALL getunit( nchanl )

    OPEN(UNIT=nchanl,FILE=trim(grdbasfn),STATUS='old',              &
         FORM='unformatted',IOSTAT= istat )

    CALL grib_getdims(nchanl, nx,ny,nz,nstyps, ireturn)

    CLOSE(UNIT=nchanl)
    CALL retunit( nchanl )

  ELSE

    WRITE(6,'(a,i3,a)')                                                 &
        ' Data format flag had an invalid value ',                      &
          hinfmt ,' program stopped.'
    CALL arpsstop('arpsstop called from get_dims_from_data-wrong flag',1)

  END IF

!  500   CONTINUE

  RETURN

  999   CONTINUE
  WRITE(6,'(1x,a,a,/1x,i3,a)')                                          &
      'Error occured when opening file ',trim(grdbasfn),            &
      'using FORTRAN unit ',nchanl,' Program stopped in get_dims_from_data.'

  CALL arpsstop('arpsstop called from get_dims_from_data opening file',1)

END SUBROUTINE get_dims_from_data
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE BIN_GETDIMS                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE bin_getdims(inch, nx,ny,nz,nstyps, ireturn) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  7/17/2000.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    inch     Channel number for binary reading.
!
!  OUTPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    nstyps   Number of soil types
!
!    ireturn  Return status indicator
!             =0, successful read of all data
!             =1, error reading data
!             =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  INTEGER :: nstyps            ! Number of soil types
  INTEGER :: inch              ! Channel number for binary reading
  INTEGER :: ireturn           ! Return status indicator
  CHARACTER (LEN=10) :: tmunit
  CHARACTER (LEN=40) :: fmtverin
  INTEGER :: i,nocmnt
  REAL :: time
  CHARACTER (LEN=80) :: runname

  INTEGER :: idummy,totin
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  READ(inch,ERR=110,END=120) fmtverin
  READ(inch,ERR=110,END=120) runname
  WRITE(6,'(//''  THE NAME OF THIS RUN IS:  '',A//)') runname

  READ(inch,ERR=110,END=120) nocmnt

  IF( nocmnt > 0 ) THEN
    DO i=1,nocmnt
      READ(inch,ERR=110,END=120)
    END DO
  END IF

  READ(inch,ERR=110,END=120) time,tmunit
  READ(inch,ERR=110,END=120) nx, ny, nz

  WRITE(6,'(a,3i5,a)') 'nx,ny,nz read in from data are ',nx,ny,nz,      &
      ' respectively.'

  READ(inch,ERR=110,END=120) idummy,idummy,idummy,idummy,idummy,  &
                             idummy,idummy,idummy,idummy,totin,   &
                             idummy,idummy,idummy,idummy,idummy,  &
                             idummy,idummy,idummy,idummy,idummy 
  IF (totin /= 0) THEN
    READ(inch,ERR=110,END=120) ! block of 20 REALs ...
    READ(inch,ERR=110,END=120) nstyps,idummy,idummy,idummy,idummy,  &
                               idummy,idummy,idummy,idummy,idummy,  &
                               idummy,idummy,idummy,idummy,idummy,  &
                               idummy,idummy,idummy,idummy,idummy
    WRITE (6,*) "nstyps read in as",nstyps
  ELSE
    WRITE (6,*) "nstyps not defined in data, set to 4"
    nstyps = 4
  ENDIF

  ireturn = 0

  RETURN

  110   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in BIN_GETDIMS'
  ireturn=1
  RETURN

  120   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in BIN_GETDIMS'
  ireturn=2
  RETURN

END SUBROUTINE bin_getdims
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ASC_GETDIMS                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE asc_getdims(inch, nx,ny,nz,nstyps, ireturn) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  7/17/2000.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    inch     Channel number for ASCII reading.
!
!  OUTPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    nstyps   Number of soil types
!
!    ireturn  Return status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  INTEGER :: nstyps            ! Number of soil types
  INTEGER :: inch              ! Channel number for binary reading
  INTEGER :: ireturn           ! Return status indicator
  CHARACTER (LEN=40) :: fmtverin
  CHARACTER (LEN=10) :: tmunit
  INTEGER :: i,nocmnt
  REAL :: time
  CHARACTER (LEN=80) :: runname

  INTEGER :: idummy,totin
  REAL::     rdummy
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  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

  READ(inch,'(1x,i4)',ERR=110,END=120) nocmnt
  IF( nocmnt > 0 ) THEN
    DO i=1,nocmnt
      READ(inch,'(1x,a80)',ERR=110,END=120)
    END DO
  END IF

  READ(inch,'(1x,e16.8,1x,a10)',ERR=110,END=120) time,tmunit
  READ(inch,'(1x, 3i12)',ERR=110,END=120) nx,ny,nz

  WRITE(6,'(a,3i5,a)') 'nx,ny,nz read in from data are ',nx,ny,nz,      &
      ' respectively.'

  READ(inch,'(1x,10i8)',ERR=110,END=120)                          &
                             idummy,idummy,idummy,idummy,idummy,  &
                             idummy,idummy,idummy,idummy,totin,   &   !should be modified
                             idummy,idummy,idummy,idummy,idummy,  &
                             idummy,idummy,idummy,idummy,idummy
    

  READ(inch,'(1x,8E16.10)',ERR=110,END=120)                       &
                             rdummy,rdummy,rdummy,rdummy,rdummy,  &
                             rdummy,rdummy,rdummy,rdummy,rdummy,  &
                             rdummy,rdummy,rdummy,rdummy,rdummy,  &
                             rdummy,rdummy,rdummy,rdummy,rdummy
    
  IF (totin /= 0) THEN
    READ(inch,'(1x,10i8)',ERR=110,END=120)                          &
                               nstyps,idummy,idummy,idummy,idummy,  &
                               idummy,idummy,idummy,idummy,idummy,  &
                               idummy,idummy,idummy,idummy,idummy,  &
                               idummy,idummy,idummy,idummy,idummy
    READ(inch,'(1x,8E16.10)',ERR=110,END=120) ! block of 20 REALs ...
    WRITE (6,*) "nstyps read in as",nstyps
  ELSE
    WRITE (6,*) "nstyps not defined in data, set to 4"
    nstyps = 4
  ENDIF

  ireturn = 0

  RETURN

  110   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in ASC_GETDIMS.'
  ireturn=1
  RETURN

  120   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in ASC_GETDIMS.'
  ireturn=2
  RETURN
END SUBROUTINE asc_getdims

!##################################################################
!##################################################################
!######                                                      ######
!######             SUBROUTINE HDF_GETDIMS                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################


SUBROUTINE hdf_getdims(filename, nx,ny,nz,nstyps, ireturn),6
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  7/17/2000.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    filename Channel number for binary reading.
!
!  OUTPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    ireturn  Return status indicator
!             =0, successful read of all data
!             =1, error reading data
!             =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: stat, sd_id, sds_id

  CHARACTER (LEN=*) :: filename

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  INTEGER :: nstyps
  INTEGER :: ireturn           ! Return status indicator
  CHARACTER (LEN=10) :: tmunit
  CHARACTER (LEN=40) :: fmtverin
  INTEGER :: i,nocmnt
  REAL :: time
  CHARACTER (LEN=80) :: runname
  CHARACTER (LEN=80 ) :: cmnt(50)  ! String of comments on this model run
  INTEGER istat

!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  CALL hdfopen(filename,1,sd_id)

  IF (sd_id < 0) THEN
    WRITE (6,*) "HDFREAD: ERROR opening ",                              &
                 trim(filename)," for reading."
    GO TO 110
  END IF

!-----------------------------------------------------------------------
!
!  Get dimensions of data in binary file and check against
!  the dimensions passed to HDFREAD
!
!-----------------------------------------------------------------------

  CALL hdfrdi(sd_id,"nx",nx,istat)
  CALL hdfrdi(sd_id,"ny",ny,istat)
  CALL hdfrdi(sd_id,"nz",nz,istat)

  CALL hdfrdi(sd_id,"nstyp",nstyps,istat)

  WRITE(6,'(a,3i5,a)') 'nx,ny,nz read in from data are ',nx,ny,nz,      &
      ' respectively.'
  WRITE (6,*) "nstyps read in as",nstyps

  ireturn = 0
  GO TO 130

!-----------------------------------------------------------------------
!
!  Error during read
!
!-----------------------------------------------------------------------

  110   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in HDF_GETDIMS.'
  ireturn=1

  GO TO 130

!-----------------------------------------------------------------------
!
!  End-of-file during read
!
!-----------------------------------------------------------------------

!  120   CONTINUE
!  WRITE(6,'(/a/)') ' End of file reached in HDF_GETDIMS.'
!  ireturn=2

  130   CONTINUE

!tmp  stat = sfendacc(sd_id)   ! is this necessary?
  CALL hdfclose(sd_id,stat)

  RETURN
END SUBROUTINE hdf_getdims
!
!##################################################################
!##################################################################
!######                                                      ######
!######             SUBROUTINE PAK_GETDIMS                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE pak_getdims(inch, nx,ny,nz,nstyps, ireturn) 1,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  7/17/2000.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    inch     Channel number for binary reading.
!
!  OUTPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    ireturn  Return status indicator
!             =0, successful read of all data
!             =1, error reading data
!             =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  INTEGER :: nstyps            ! Number of soil types
  INTEGER :: inch              ! Channel number for binary reading
  INTEGER :: ireturn           ! Return status indicator
  CHARACTER (LEN=10) :: tmunit
  CHARACTER (LEN=40) :: fmtverin
  INTEGER :: i,nocmnt
  REAL :: time
  CHARACTER (LEN=80) :: runname
  INTEGER :: bgrdin,bbasin,bvarin,mstin,brainin,bicein,btkein,          &
           btrbin,bsfcin,landin,totin,prcin,radin,flxin,                &
           snowcin,snowin
  INTEGER :: ihdlen,ilen
  PARAMETER (ihdlen=5000)
  INTEGER :: ihead(ihdlen)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  READ(inch,ERR=110,END=120) ihead

!  CALL decdhdr(nx,ny,nz,ihdlen,ihead,ilen,time,                         &
!               bgrdin,bbasin,bvarin,mstin,brainin,bicein,btkein,        &
!               btrbin,bsfcin,landin,totin,prcin,radin,flxin,            &
!               snowcin,snowin,fmtverin,tmunit)
   CALL decdhdr

  WRITE (6,*) "nstyps not defined in data, set to 4"
  nstyps = 4

  ireturn =0

  RETURN

  110   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in PAK_GETDIMS.'
  ireturn=1
  RETURN

  120   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in PAK_GETDIMS.'
  ireturn=2

  RETURN
END SUBROUTINE pak_getdims
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GRIB_GETDIMS               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE grib_getdims(inch, nx,ny,nz,nstyps, ireturn) 1,1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!  Read in grid dimensions from base state/grid history data.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  7/17/2000.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT :
!
!    inch     Channel number for binary reading.
!
!  OUTPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the vertical
!
!    nstyps   NUmber of soil types.
!
!    ireturn  Return status indicator
!             =0, successful read of all data
!             =1, error reading data
!             =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: inch
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  INTEGER :: nstyps            ! Number of soil types
  INTEGER :: ireturn           ! Return status indicator

  CHARACTER (LEN=10) :: tmunit
  CHARACTER (LEN=40) :: fmtverin
  INTEGER :: i,j,nocmnt
  REAL :: time
  CHARACTER (LEN=80) :: runname
  CHARACTER (LEN=80 ) :: cmnt(50)  ! String of comments on this model run
!
!-----------------------------------------------------------------------
!
!  Parameters for GRIB packing
!
!-----------------------------------------------------------------------
!
  INTEGER :: nbufsz
  PARAMETER ( nbufsz = 800000 )  ! Size of GRIB buffer
  CHARACTER (LEN=1) :: mgrib(nbufsz) ! Buffer to carry GRIB messages
  INTEGER :: npos,hhdlen
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  READ (inch,ERR=110,END=120) (mgrib(i),i=1,3)

  npos = 0

  hhdlen = ICHAR(mgrib(npos+1))*65536                                   &
         + ICHAR(mgrib(npos+2))*256                                     &
         + ICHAR(mgrib(npos+3))

  BACKSPACE (inch)

  READ (inch,ERR=110,END=120) (mgrib(i),i=1,hhdlen)

  npos = npos + 3

  DO i=1,40
    fmtverin(i:i) = mgrib(npos+i)
  END DO

  npos = npos + 40

  DO i=1,80
    runname(i:i) = mgrib(npos+i)
  END DO

  npos = npos + 80

  nocmnt = ICHAR(mgrib(npos+1))
  npos = npos + 1

  IF( nocmnt > 0 ) THEN
    DO j=1,nocmnt
      DO  i=1,80
        cmnt(j)(i:i) = mgrib(npos+i)
      END DO
      npos = npos + 80
    END DO
  END IF

  CALL ibm2flt( mgrib(npos+1), time )
  npos = npos + 4

  nx = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
  npos = npos + 2
  ny = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
  npos = npos + 2
  nz = ICHAR(mgrib(npos+1))*256 + ICHAR(mgrib(npos+2))
  npos = npos + 2

  WRITE (6,*) "nstyps not defined in data, set to 4"
  nstyps = 4

  ireturn = 0

  RETURN

  110   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in GRIB_GETDIMS'
  ireturn=1
  RETURN

  120   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in GRIB_GETDIMS'
  ireturn=2
  RETURN

  RETURN
END SUBROUTINE grib_getdims
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE EDGFILL                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE edgfill(a,nx1,nx2,ibgn,iend,ny1,ny2,jbgn,jend,               & 451
           nz1,nz2,kbgn,kend)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Fill in the edges of a data array from the valid interior grid
!  points so that the arrays are completely filled. This is done primarily
!  for the benefit of the compression routine which must have a valid
!  value in each array location.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  7/14/92.
!
!  9/1/94 (Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    a        Array to have edges filled in
!
!    nx1:nx2  Dimensioned range of the  first index of array "a"
!    ny1:ny2  Dimensioned range of the second index of array "a"
!    nz1:nz2  Dimensioned range of the  third index of array "a"
!
!    ibgn,iend  Valued range of the  first index of array "a"
!    jgbn,jend  Valued range of the second index of array "a"
!    kbgn,kend  Valued range of the  third index of array "a"
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx1,nx2,ny1,ny2,nz1,nz2
  REAL :: a(nx1:nx2,ny1:ny2,nz1:nz2)
  INTEGER :: ibgn,iend,jbgn,jend,kbgn,kend
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Fill in top and bottom grid levels, as needed.
!
!-----------------------------------------------------------------------
!
  IF( kbgn > nz1 ) THEN
    DO k=nz1,(kbgn-1)
      DO j=jbgn,jend
        DO i=ibgn,iend
          a(i,j,k)=a(i,j,kbgn)
        END DO
      END DO
    END DO
  END IF

  IF( kend < nz2 ) THEN
    DO k=(kend+1),nz2
      DO j=jbgn,jend
        DO i=ibgn,iend
          a(i,j,k)=a(i,j,kend)
        END DO
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Fill in north and south edges, as needed.
!
!-----------------------------------------------------------------------
!
  IF( jbgn > ny1 ) THEN
    DO k=nz1,nz2
      DO j=ny1,(jbgn-1)
        DO i=ibgn,iend
          a(i,j,k)=a(i,jbgn,k)
        END DO
      END DO
    END DO
  END IF

  IF( jend < ny2 ) THEN
    DO k=nz1,nz2
      DO j=(jend+1),ny2
        DO i=ibgn,iend
          a(i,j,k)=a(i,jend,k)
        END DO
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Fill in east and west edges, as needed.
!
!-----------------------------------------------------------------------
!
  IF( ibgn > nx1 ) THEN
    DO k=nz1,nz2
      DO j=ny1,ny2
        DO i=nx1,(ibgn-1)
          a(i,j,k)=a(ibgn,j,k)
        END DO
      END DO
    END DO
  END IF

  IF( iend < nx2 ) THEN
    DO k=nz1,nz2
      DO j=ny1,ny2
        DO i=(iend+1),nx2
          a(i,j,k)=a(iend,j,k)
        END DO
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE edgfill
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE IEDGFILL                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE iedgfill(a,nx1,nx2,ibgn,iend,ny1,ny2,jbgn,jend,              & 17
           nz1,nz2,kbgn,kend)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Fill in the edges of a data array from the valid interior grid
!  points so that the arrays are completely filled. This is done primarily
!  for the benefit of the compression routine which must have a valid
!  value in each array location.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  7/14/92.
!
!  9/1/94 (Y. Lu)
!  Cleaned up documentation.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    a        Array to have edges filled in
!
!    nx1:nx2  Dimensioned range of the  first index of array "a"
!    ny1:ny2  Dimensioned range of the second index of array "a"
!    nz1:nz2  Dimensioned range of the  third index of array "a"
!
!    ibgn,iend  Valued range of the  first index of array "a"
!    jgbn,jend  Valued range of the second index of array "a"
!    kbgn,kend  Valued range of the  third index of array "a"
!
  IMPLICIT NONE
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx1,nx2,ny1,ny2,nz1,nz2
  INTEGER :: a(nx1:nx2,ny1:ny2,nz1:nz2)
  INTEGER :: ibgn,iend,jbgn,jend,kbgn,kend
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Fill in top and bottom grid levels, as needed.
!
!-----------------------------------------------------------------------
!
  IF( kbgn > nz1 ) THEN
    DO k=nz1,(kbgn-1)
      DO j=jbgn,jend
        DO i=ibgn,iend
          a(i,j,k)=a(i,j,kbgn)
        END DO
      END DO
    END DO
  END IF

  IF( kend < nz2 ) THEN
    DO k=(kend+1),nz2
      DO j=jbgn,jend
        DO i=ibgn,iend
          a(i,j,k)=a(i,j,kend)
        END DO
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Fill in north and south edges, as needed.
!
!-----------------------------------------------------------------------
!
  IF( jbgn > ny1 ) THEN
    DO k=nz1,nz2
      DO j=ny1,(jbgn-1)
        DO i=ibgn,iend
          a(i,j,k)=a(i,jbgn,k)
        END DO
      END DO
    END DO
  END IF

  IF( jend < ny2 ) THEN
    DO k=nz1,nz2
      DO j=(jend+1),ny2
        DO i=ibgn,iend
          a(i,j,k)=a(i,jend,k)
        END DO
      END DO
    END DO
  END IF
!
!-----------------------------------------------------------------------
!
!  Fill in east and west edges, as needed.
!
!-----------------------------------------------------------------------
!
  IF( ibgn > nx1 ) THEN
    DO k=nz1,nz2
      DO j=ny1,ny2
        DO i=nx1,(ibgn-1)
          a(i,j,k)=a(ibgn,j,k)
        END DO
      END DO
    END DO
  END IF

  IF( iend < nx2 ) THEN
    DO k=nz1,nz2
      DO j=ny1,ny2
        DO i=(iend+1),nx2
          a(i,j,k)=a(iend,j,k)
        END DO
      END DO
    END DO
  END IF

  RETURN
END SUBROUTINE iedgfill