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


SUBROUTINE binread(nx,ny,nz,nzsoil,nstyps,grdbas,inch,time,x,y,z,zp,    & 2,5
           zpsoil,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,                      &
           tsoil,qsoil,wetcanp,snowdpth,                                &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,radswnet,radlwin,                         &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           ireturn)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in binary data set created by ARPS using history dump format
!  No. 1.
!  All data read in are located at the original staggered grid points
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  2/27/92.
!
!  MODIFICATION HISTORY:
!
!  6/08/92
!  Added full documentation (K. Brewster)
!
!  7/14/92 (K. Brewster)
!  Added runname, comment and version number reading
!
!  8/20/92 (M. Xue)
!  Added data reading of computational z coordinate array z.
!
!  4/23/93 (M. Xue)
!  New data format.
!
!  02/06/95 (Y. Liu)
!  Added map projection parameters into the binary dumping
!
!  03/26/96 (G. Bassett)
!  Backwards compatibility added for ARPS 3.2 and ARPS 4.0 binary
!  history dump formats.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  05/15/2002 (J. Brotzge)
!  Added additional variables to allow for multiple soil schemes  
!
!  1 June 2002 Eric Kemp
!  Bug fixes for new soil variables.
!
!-----------------------------------------------------------------------
!
!  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
!    nzsoil   Number of grid points in the soil  
!
!    grdbas   Data read flag.
!             =1, only grid and base state arrays will be read
!             =0, all arrays will be read based on data
!                 parameter setting.
!    inch     Channel number for binary reading.
!             This channel must be opened for unformatted reading
!             by the calling routine.
!
!  OUTPUT:
!
!    time     Time in seconds of data read from "filename"
!
!    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)
!    zpsoil   z coordinate of grid points in the soil (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
!
!    tsoil    Soil temperature (K)
!    qsoil    Soil moisture (m**3/m**3) 
!    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
!    radswnet Net shortwave radiation
!    radlwin  Incoming longwave radiation
!
!    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
!             =0, successful read of all data
!             =1, error reading data
!             =2, end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

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

  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 :: zpsoil(nx,ny,nzsoil) ! Physical height coordinate defined at
                               ! w-point of the soil (m)  


  INTEGER :: grdbas            ! Data read flag.
  INTEGER :: inch              ! Channel number for binary reading
  REAL :: time                 ! Time in seconds of data read
                               ! from "filename"

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

  REAL :: qvprt (nx,ny,nz)     ! Perturbation water vapor mixing
                               ! ratio (kg/kg)
  REAL :: qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
  REAL :: 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 u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: wbar  (nx,ny,nz)     ! Base state w-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor mixing ratio

  INTEGER :: nstyps                    ! Number of soil type
  INTEGER :: soiltyp(nx,ny,nstyps)     ! Soil type
  REAL :: stypfrct(nx,ny,nstyps)    ! Soil type fraction
  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 :: tsoil  (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
  REAL :: qsoil  (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) 
  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 :: radswnet(nx,ny)      ! Net shortwave radiation
  REAL :: radlwin(nx,ny)       ! Incoming longwave radiation

  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**2*s))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  INTEGER :: ireturn           ! Return status indicator
!
!-----------------------------------------------------------------------
!
!  Parameters describing routine that wrote the gridded data
!
!-----------------------------------------------------------------------
!
! 06/28/2002 Zuwen He
!
! fmtver??: to label each data a version. 
! intver??: an integer to allow faster comparison than fmtver??, 
!           which are strings. 
!
! Verion 5.00: significant change in soil variables since version 4.10. 
! 
  CHARACTER (LEN=40) :: fmtver320,fmtver400,fmtver410,fmtver500
  INTEGER  :: intver,intver320,intver400,intver410,intver500

  PARAMETER (fmtver320='003.20 Binary Data',intver320=320)
  PARAMETER (fmtver400='004.00 Binary Data',intver400=400)
  PARAMETER (fmtver410='004.10 Binary Data',intver410=410)
  PARAMETER (fmtver500='005.00 Binary Data',intver500=500)   

  CHARACTER (LEN=40) :: fmtverin
!
  CHARACTER (LEN=10) :: tmunit
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: lchanl
  PARAMETER (lchanl=6)      ! Channel number for formatted printing.

  INTEGER :: i,j,k,is
  INTEGER :: nstyp1
  CHARACTER (LEN=12) :: label
!  INTEGER :: nxin,nyin,nzin
  INTEGER :: nxin,nyin,nzin,nzsoilin
  INTEGER :: bgrdin,bbasin,bvarin,bicein,btkein,btrbin
  INTEGER :: idummy
  REAL :: rdummy
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'indtflg.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'mp.inc'            ! mpi parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Read header info
!
!-----------------------------------------------------------------------
!
  READ(inch,ERR=110,END=120) fmtverin

  IF (fmtverin == fmtver320) THEN 
    intver=intver320
  ELSE IF (fmtverin == fmtver400) THEN 
    intver=intver400
  ELSE IF (fmtverin == fmtver410) THEN 
    intver=intver410
  ELSE IF (fmtverin == fmtver500) THEN 
    intver=intver500
  ELSE 
    IF (myproc == 0) WRITE(6,'(/1x,a,a,a/)')                        &
        'Incoming data format, fmtverin=',fmtverin,                 & 
        ', not found. The Job stopped.'
    CALL arpsstop('arpstop called from binread. ',1)
  END IF 

  IF (myproc == 0) WRITE(6,'(/1x,a,a/)')                        &
      'Incoming data format, fmtverin=',fmtverin

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

  IF (myproc == 0)  &
     WRITE(6,'(//''  THE NAME OF THIS RUN IS:  '',A//)') runname

  IF( nocmnt > 0 ) THEN
    IF(myproc == 0)THEN
      DO i=1,nocmnt
        WRITE(6,'(1x,a)') cmnt(i)
      END DO
    END IF
  END IF

  READ(inch,ERR=110,END=120) time,tmunit
!
!-----------------------------------------------------------------------
!
!  Get dimensions of data in binary file and check against
!  the dimensions passed to BINREAD
!
!-----------------------------------------------------------------------
!
  IF (intver <= intver410) THEN 

    READ(inch,ERR=110,END=120) nxin, nyin, nzin
    nzsoilin = 2  ! for version prior to 410, it is a two-layer soil model

  ELSE IF (intver >= intver500) THEN 

    READ(inch,ERR=110,END=120) nxin, nyin, nzin,nzsoilin

  END IF 
!
! Data validation: dimensions 
!
  IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz) THEN
    IF(myproc == 0)THEN
      WRITE(6,'(1x,a)')                                                   &
         ' Dimensions in BINREAD inconsistent with data.'
      WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
      WRITE(6,'(1x,a,3I15)') ' Expected:  ', nx, ny, nz
      WRITE(6,'(1x,a)')      ' Program aborted in BINREAD.'
    END IF
    CALL arpsstop('arpstop called from binread nx-ny-nz read ',1)
  END IF
 
  IF(nzsoilin /= nzsoil) THEN

    IF (intver <= intver410) THEN 
      IF(myproc == 0) WRITE(6,'(1x,a,a/,2(1x,a/))')           & 
        ' The incoming data version is ', fmtverin,            & 
        ' In the input file, nzsoil must be set to 2. ',       & 
        ' Program aborted in BINREAD.'

    ELSE IF (intver >= intver500) THEN 

      IF(myproc == 0)THEN
        WRITE(6,'(1x,a)') & 
                ' Dimensions in BINREAD inconsistent with data.'
        WRITE(6,'(1x,a,I15)') ' Read were: ', nzsoilin
        WRITE(6,'(1x,a,I15)') ' Expected:  ', nzsoil
        WRITE(6,'(1x,a)') ' Program aborted in BINREAD.'
      END IF
    END IF
    CALL arpsstop('arpstop called from binread nzsoil read ',1)

  END IF 
!
!-----------------------------------------------------------------------
!
!  Read in flags for different data groups
!
!-----------------------------------------------------------------------
!
  IF( grdbas == 1 ) THEN ! Read grid and base state arrays

    IF (myproc == 0)  &
       WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')                               &
         'To read grid and base state data at time ', time,             &
         ' secs = ',(time/60.),' mins.'

    READ(inch,ERR=110,END=120)                                          &
         bgrdin,bbasin,bvarin,mstin,bicein,                             &
         btrbin,idummy,idummy,landin,totin,                             &
         btkein,idummy,idummy,mapproj,month,                            &
         day, year, hour, minute, second

  ELSE ! Normal data reading

    IF (myproc == 0)  &
       WRITE(lchanl,'(1x,a,f8.1,a,f8.3,a/)')'To read data for time:',      &
         time,' secs = ',(time/60.),' mins.'

    READ(inch,ERR=110,END=120)                                          &
         grdin,basin,varin,mstin,icein,                                 &
         trbin, sfcin,rainin,landin,totin,                              &
         tkein,idummy,idummy,mapproj, month,                            &
         day, year, hour, minute, second

  END IF

  READ(inch,ERR=110,END=120)                                            &
                  umove,vmove,xgrdorg,ygrdorg,trulat1,                  &
                  trulat2,trulon,sclfct,rdummy,rdummy,                  &
                  rdummy,rdummy,rdummy,rdummy,rdummy,                   &
                  tstop,thisdmp,latitud,ctrlat,ctrlon

  IF ( totin /= 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  Read in additional parameters for ARPS history dump 4.0 or later
!  version.
!
!-----------------------------------------------------------------------
!
    READ(inch,ERR=110,END=120)                                          &
         nstyp1,  prcin, radin, flxin,snowcin,                           &
         snowin,idummy,idummy,idummy,idummy,                            &
         idummy,idummy,idummy,idummy,idummy,                            &
         idummy,idummy,idummy,idummy,idummy

    IF ( nstyp1 < 1 ) THEN
      nstyp1 = 1
    END IF

    READ(inch,ERR=110,END=120)                                          &
         rdummy,rdummy,rdummy,rdummy,rdummy,                            &
         rdummy,rdummy,rdummy,rdummy,rdummy,                            &
         rdummy,rdummy,rdummy,rdummy,rdummy,                            &
         rdummy,rdummy,rdummy,rdummy,rdummy
  END IF
!
!-----------------------------------------------------------------------
!
!  Read in x, y, z and zp arrays.
!
!----------------------------------------------------------------------
!
  IF( grdin == 1 .OR. grdbas == 1 ) THEN
    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) x
    IF (myproc == 0) WRITE(lchanl,910) label,' x.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) y
    IF (myproc == 0) WRITE(lchanl,910) label,' y.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) z
    IF (myproc == 0) WRITE(lchanl,910) label,' z.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) zp
    IF (myproc == 0) WRITE(lchanl,910) label,' zp.'

    IF (intver <= intver410) THEN 
!
! nzsoil must equal to 2, 06/28/2002, Zuwen 
! assume zpsoil(,,2) is one meter below the surface. 
!
      DO j=1,ny
        DO i=1,nx
          zpsoil(i,j,1)=zp(i,j,2) 
          zpsoil(i,j,1)=zpsoil(i,j,1)-1.
        END DO 
      END DO 

      IF (myproc == 0) THEN 
        WRITE(lchanl,910) ' Assign zpsoil. '
        WRITE(lchanl,910) ' Assume zpsoil(,,1) is zp(,,2). '
        WRITE(lchanl,910) ' Assume zpsoil(,,2) is zp(,,2)-1. '
      END IF 
          
    ELSE IF (intver >= intver500) THEN 

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) zpsoil
      IF (myproc == 0)  WRITE(lchanl,910) label,' zpsoil.'

    END IF  ! intver

  END IF  ! grdin
!
!-----------------------------------------------------------------------
!
!  Read in base state fields
!
!----------------------------------------------------------------------
!
  IF( basin == 1 .OR. grdbas == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ubar
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' ubar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) vbar
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' vbar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) wbar
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' wbar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ptbar
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' ptbar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) pbar
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' pbar.'

    IF( mstin == 1) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) qvbar
      IF (myproc == 0)  &
         WRITE(lchanl,910) label,' qvbar.'
    END IF

    IF (landin == 1) THEN

      IF (nstyp1 <= 1) THEN

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) ((soiltyp(i,j,1),i=1,nx),j=1,ny)
        IF (myproc == 0) WRITE(lchanl,910) label,' soiltyp.'

      ELSE

        DO is=1,nstyp1
          IF (is <= nstyps) THEN
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)                                    &
                  ((soiltyp(i,j,is),i=1,nx),j=1,ny)
            IF (myproc == 0)  &
               WRITE(lchanl,910) label,' soiltyp.'
  
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)                                    &
                  ((stypfrct(i,j,is),i=1,nx),j=1,ny)
            IF (myproc == 0)  &
               WRITE(lchanl,910) label,' stypfrct.'
          ELSE
            READ(inch,ERR=110,END=120) label
            IF (myproc == 0)  &
               WRITE(lchanl,910) label,'skipping soiltyp'
            READ(inch,ERR=110,END=120)
            READ(inch,ERR=110,END=120) label
            IF (myproc == 0)  &
               WRITE(lchanl,910) label,'skipping stypfrct.'
            READ(inch,ERR=110,END=120) 
          ENDIF
        END DO

      END IF

      CALL fix_stypfrct_nstyp(nx,ny,nstyp1,nstyp,stypfrct)

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) vegtyp
      IF (myproc == 0) WRITE(lchanl,910) label,' vegtyp.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) lai
      IF (myproc == 0) WRITE(lchanl,910) label,' lai.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) roufns
      IF (myproc == 0) WRITE(lchanl,910) label,' roufns.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) veg
      IF (myproc == 0) WRITE(lchanl,910) label,' veg.'

    END IF

  END IF

  IF( grdbas == 1 ) GO TO 930

  IF( varin == 1 ) THEN

    IF ( totin == 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  Read in perturbations from history dump
!
!-----------------------------------------------------------------------
!
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) uprt
      IF (myproc == 0) WRITE(lchanl,910) label,' uprt.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) vprt
      IF (myproc == 0) WRITE(lchanl,910) label,' vprt.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) wprt
      IF (myproc == 0) WRITE(lchanl,910) label,' wprt.'
!
!-----------------------------------------------------------------------
!
!  Read in scalars
!
!----------------------------------------------------------------------
!
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ptprt
      IF (myproc == 0) WRITE(lchanl,910) label,' ptprt.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) pprt
      IF (myproc == 0) WRITE(lchanl,910) label,' pprt.'

    ELSE
!
!-----------------------------------------------------------------------
!
!  Read in total values of variables from history dump
!
!----------------------------------------------------------------------
!
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) uprt
      IF (myproc == 0) WRITE(lchanl,910) label,' u.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            uprt(i,j,k) = uprt(i,j,k) - ubar(i,j,k)
          END DO
        END DO
      END DO

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) vprt
      IF (myproc == 0) WRITE(lchanl,910) label,' v.'
      DO k=1,nz-1
        DO j=1,ny
          DO i=1,nx-1
            vprt(i,j,k) = vprt(i,j,k) - vbar(i,j,k)
          END DO
        END DO
      END DO

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) wprt
      IF (myproc == 0) WRITE(lchanl,910) label,' w.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ptprt
      IF (myproc == 0) WRITE(lchanl,910) label,' pt.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            ptprt(i,j,k) = ptprt(i,j,k) - ptbar(i,j,k)
          END DO
        END DO
      END DO

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) pprt
      IF (myproc == 0) WRITE(lchanl,910) label,' p.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            pprt(i,j,k) = pprt(i,j,k) - pbar(i,j,k)
          END DO
        END DO
      END DO

    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Read in moisture variables
!
!-----------------------------------------------------------------------
!
  IF( mstin == 1 ) THEN

    IF ( totin == 0 ) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) qvprt
      IF (myproc == 0) WRITE(lchanl,910) label,' qvprt.'

    ELSE

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) qvprt
      IF (myproc == 0) WRITE(lchanl,910) label,' qv.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            qvprt(i,j,k) = qvprt(i,j,k) - qvbar(i,j,k)
          END DO
        END DO
      END DO

    END IF

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) qc
    IF (myproc == 0) WRITE(lchanl,910) label,' qc.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) qr
    IF (myproc == 0) WRITE(lchanl,910) label,' qr.'

    IF( rainin == 1 ) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) raing
      IF (myproc == 0) WRITE(lchanl,910) label,' raing.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) rainc
      IF (myproc == 0) WRITE(lchanl,910) label,' rainc.'
    END IF

    IF( prcin == 1 ) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((prcrate(i,j,1),i=1,nx),j=1,ny)
      IF (myproc == 0) WRITE(lchanl,910) label,' prcrate1.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((prcrate(i,j,2),i=1,nx),j=1,ny)
      IF (myproc == 0) WRITE(lchanl,910) label,' prcrate2.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((prcrate(i,j,3),i=1,nx),j=1,ny)
      IF (myproc == 0) WRITE(lchanl,910) label,' prcrate3.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((prcrate(i,j,4),i=1,nx),j=1,ny)
      IF (myproc == 0) WRITE(lchanl,910) label,' prcrate4.'
    END IF

    IF( icein == 1 ) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) qi
      IF (myproc == 0) WRITE(lchanl,910) label,' qi.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) qs
      IF (myproc == 0) WRITE(lchanl,910) label,' qs.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) qh
      IF (myproc == 0) WRITE(lchanl,910) label,' qh.'

    END IF
  END IF

  IF( tkein == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) tke
    IF (myproc == 0) WRITE(lchanl,910) label,' tke.'

  END IF

  IF( trbin == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) kmh
    IF (myproc == 0) WRITE(lchanl,910) label,' kmh.'

    IF ( intver >= intver410 ) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) kmv
      IF (myproc == 0) WRITE(lchanl,910) label,' kmv.'

    END IF ! intver

  END IF

  IF( sfcin == 1 ) THEN

    IF (nstyp1 <= 1) THEN

      IF (intver <= intver410) THEN 

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) ((tsoil(i,j,1,0),i=1,nx),j=1,ny)
        IF (myproc == 0) WRITE(lchanl,910) label,' tsfc.'

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) ((tsoil(i,j,2,0),i=1,nx),j=1,ny)
        IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.'

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120)  ((qsoil(i,j,1,0),i=1,nx),j=1,ny)
        IF (myproc == 0) WRITE(lchanl,910) label,' wetsfc.'

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120)  ((qsoil(i,j,2,0),i=1,nx),j=1,ny)
        IF (myproc == 0) WRITE(lchanl,910) label,' wetdp.'

      ELSE IF (intver >= intver500) THEN 

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120)                                  &
            (((tsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) 
        IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.'

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120)                                &
            (((qsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) 
        IF (myproc == 0) WRITE(lchanl,910) label,' qsoil.'

      END IF  ! intver

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
            ((wetcanp(i,j,0),i=1,nx),j=1,ny)
      IF (myproc == 0) WRITE(lchanl,910) label,' wetcanp.'

    ELSE

      DO is=0,nstyp1
        IF (is <= nstyps) THEN

          IF (intver <= intver410) THEN 

            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) ((tsoil(i,j,1,is),i=1,nx),j=1,ny)
            IF (myproc == 0) WRITE(lchanl,910) label,' tsfc.'
  
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) ((tsoil(i,j,2,is),i=1,nx),j=1,ny)
            IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.'
  
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) ((qsoil(i,j,1,is),i=1,nx),j=1,ny)
            IF (myproc == 0) WRITE(lchanl,910) label,' wetsfc.'

            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) ((qsoil(i,j,2,is),i=1,nx),j=1,ny)
            IF (myproc == 0) WRITE(lchanl,910) label,' wetdp.'

          ELSE IF (intver >= intver500) THEN 

            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)                                      &
                  (((tsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) 
            IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.'
  
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)                                      &
                  (((qsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) 
            IF (myproc == 0) WRITE(lchanl,910) label,' qsoil.'

          END IF  ! intver

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)                                      &
                ((wetcanp(i,j,is),i=1,nx),j=1,ny)
          IF (myproc == 0) WRITE(lchanl,910) label,' wetcanp.'

        ELSE

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)
          IF (myproc == 0)  &
             WRITE(lchanl,910) label,'skipping tsoil.'

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)
          IF (myproc == 0)  &
             WRITE(lchanl,910) label,'skipping qsoil.'

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)
          IF (myproc == 0)  &
             WRITE(lchanl,910) label,'skipping wetcanp.'

        ENDIF

      END DO

    END IF

    CALL fix_soil_nstyp(nx,ny,nzsoil,nstyp1,nstyp,tsoil,qsoil,wetcanp)

    IF (snowcin == 1) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)
      IF (myproc == 0) WRITE(lchanl,910) label,' snowcvr -- discarding.'
    END IF

    IF (snowin == 1) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          ((snowdpth(i,j),i=1,nx),j=1,ny)
      IF (myproc == 0) WRITE(lchanl,910) label,' snowdpth.'
    END IF

  END IF

  IF( radin == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) radfrc
    IF (myproc == 0) WRITE(lchanl,910) label,' radfrc.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) radsw
    IF (myproc == 0) WRITE(lchanl,910) label,' radsw.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) rnflx
    IF (myproc == 0) WRITE(lchanl,910) label,' rnflx.'

    IF (intver <= intver410) THEN 
!
! 06/28/2002 Zuwen He
!
! Do not know how to initialized radswnet, radlwin, yet, 
! at this moment. 
!
      radswnet = 0. 
      radlwin = 0. 

    ELSE IF (intver >= intver500) THEN 

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) radswnet
      WRITE(lchanl,910) label,' radswnet.'
  
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) radlwin
      WRITE(lchanl,910) label,' radlwin.'

    END IF  ! intver 

  END IF

  IF( flxin == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) usflx
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' usflx.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) vsflx
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' vsflx.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ptsflx
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' ptsflx.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) qvsflx
    IF (myproc == 0)  &
       WRITE(lchanl,910) label,' qvsflx.'

  END IF

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

!
!-----------------------------------------------------------------------
!
!  Friendly exit message
!
!----------------------------------------------------------------------
!
  930   CONTINUE

  IF (myproc == 0)  &
     WRITE(6,'(/a,F8.1,a/)')                                               &
      ' Data at time=', time/60,' (min) were successfully read.'

  ireturn = 0

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

  110   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in BINREAD'
  ireturn=1
  RETURN
!
!-----------------------------------------------------------------------
!
!  End-of-file during read
!
!----------------------------------------------------------------------
!

  120   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in BINREAD'
  ireturn=2
  RETURN
END SUBROUTINE binread
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE BN2READ                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE bn2read(nx,ny,nz,nzsoil,nstyps,grdbas,inch,time,x,y,z,zp,    & 2,4
           zpsoil,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,                      &
           tsoil,qsoil,wetcanp,snowdpth,                                &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,radswnet,radlwin,                         &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           ireturn)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in binary data set created by ARPS using history dump format
!  No.2.
!  All data read in are located at the original staggered grid points
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  2/27/92.
!
!  MODIFICATION HISTORY:
!
!  6/08/92  Added full documentation (K. Brewster)
!
!  7/14/92 (K. Brewster)
!  Added runname, comment and version number reading
!
!  8/20/92 (M. Xue)
!  Added data reading of computational z coordinate array z.
!
!  4/23/93 (M. Xue)
!  New data format.
!
!  02/06/95 (Y. Liu)
!  Added map projection parameters into the second binary dumping
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  05/15/2002 (J. Brotzge)
!  Added variables to allow for multiple soil schemes.  
!
!-----------------------------------------------------------------------
!
!  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
!    nzsoil   Number of grid points in the soil  
!
!    grdbas   Data read flag.
!             =1, only grid and base state arrays will be read
!             =0, all arrays will be read based on data
!                          parameter setting.
!    inch     Channel number for binary reading.
!             This channel must be opened for unformatted reading
!             by the calling routine.
!
!  OUTPUT:
!
!    time     Time in seconds of data read from "filename"
!
!    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)
!    zpsoil   z coordinate of grid points in the soil (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
!
!    tsoil    Soil temperature (K)
!    qsoil    Soil moisture (m**3/m**3) 
!    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
!    radswnet Net shortwave radiation
!    radlwin  Incoming longwave radiation
!
!    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
!             =0   successful read of all data
!             =1   error reading data
!             =2   end-of-file reached during read attempt
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

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

  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 :: zpsoil(nx,ny,nzsoil) ! Physical height coordinate defined at
                               ! w-point of the soil (m) 

  INTEGER :: grdbas            ! Data read flag.
  INTEGER :: inch              ! Channel number for binary reading
  REAL :: time                 ! Time in seconds of data read
                               ! from "filename"

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

  REAL :: qvprt (nx,ny,nz)     ! Perturbation water vapor mixing
                               ! ratio (kg/kg)
  REAL :: qc    (nx,ny,nz)     ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz)     ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz)     ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz)     ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz)     ! Hail mixing ratio (kg/kg)
  REAL :: 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 u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: wbar  (nx,ny,nz)     ! Base state w-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor mixing ratio

  INTEGER :: 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 :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
  REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) 
  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 :: radswnet(nx,ny)      ! Net shortwave radiation
  REAL :: radlwin(nx,ny)       ! Incoming longwave radiation

  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**2*s))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  INTEGER :: ireturn           ! Return status indicator
!
!-----------------------------------------------------------------------
!
!  Parameters describing routine that wrote the gridded data
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=40) :: fmtver0,fmtver1,fmtverin
  PARAMETER (fmtver0='004.10 2nd Binary Data')
  PARAMETER (fmtver1='004.10 2nd Binary Data')
  CHARACTER (LEN=10) :: tmunit
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: lchanl
  PARAMETER (lchanl=6)      ! Channel number for formatted printing.

  INTEGER :: i,j,k,is
  INTEGER :: nstyp1
  CHARACTER (LEN=12) :: label
!  INTEGER :: nxin,nyin,nzin
  INTEGER :: nxin,nyin,nzin,nzsoilin
  INTEGER :: bgrdin,bbasin,bvarin,bicein,btkein,btrbin,idummy
  REAL :: rdummy
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'indtflg.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!  Read header info
!
!-----------------------------------------------------------------------
!
  READ(inch,ERR=110,END=120) fmtverin

  IF( fmtverin /= fmtver0 .AND. fmtverin /= fmtver1 ) THEN
    WRITE(6,'(/1x,a/1x,2a/1x,2a/1x,2a/1x,a)')                           &
        'Data format incompatible with the data reader.',               &
        'Format of data is ',fmtverin,' Format of reader is ',fmtver1,  &
        'compitable to ',fmtver0, '. Job stopped.'
    CALL arpsstop('arpstop called from bn2read header read ',1)
  END IF

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

  WRITE(6,'(//''  THE NAME OF THIS RUN IS:  '',A//)') runname

  IF( nocmnt > 0 ) THEN
    DO i=1,nocmnt
      WRITE(6,'(1x,a)') cmnt(i)
    END DO
  END IF

  READ(inch,ERR=110,END=120) time,tmunit
!
!-----------------------------------------------------------------------
!
!  Get dimensions of data in binary file and check against
!  the dimensions passed to BN2READ
!
!-----------------------------------------------------------------------
!
!  READ(inch,ERR=110,END=120) nxin, nyin, nzin
  READ(inch,ERR=110,END=120) nxin, nyin, nzin, nzsoilin

!  IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
  IF( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz .OR. nzsoilin /= nzsoil) THEN
    WRITE(6,'(1x,a)')                                                   &
         ' Dimensions in BN2READ inconsistent with data.'
!    WRITE(6,'(1x,a,3I15)') ' Read were: ', nxin, nyin, nzin
    WRITE(6,'(1x,a,4I15)') ' Read were: ', nxin, nyin, nzin, nzsoilin
    WRITE(6,'(1x,a)')                                                   &
         ' Program aborted in BN2READ.'
!    CALL arpsstop('arpstop called from bn2read nx-ny-nz read',1)
    CALL arpsstop('arpstop called from bn2read nx-ny-nz-nzsoil read',1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Read in flags for different data groups
!
!-----------------------------------------------------------------------
!
  IF( grdbas == 1 ) THEN ! Read grid and base state arrays

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

    READ(inch,ERR=110,END=120)                                          &
         bgrdin,bbasin,bvarin,mstin,bicein,                             &
         btrbin,idummy,idummy,landin,totin,                             &
         btkein,idummy,idummy,mapproj,month,                            &
         day,year,hour,minute,second

  ELSE ! Normal data reading

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

    READ(inch,ERR=110,END=120)                                          &
         grdin,basin,varin,mstin,icein,                                 &
         trbin,sfcin,rainin,landin,totin,                               &
         tkein,idummy,idummy,mapproj,month,                             &
         day,year,hour,minute,second

  END IF

  READ(inch,ERR=110,END=120)                                            &
                  umove,vmove,xgrdorg,ygrdorg,trulat1,                  &
                  trulat2,trulon,sclfct,rdummy,rdummy,                  &
                  rdummy,rdummy,rdummy,rdummy,rdummy,                   &
                  tstop,thisdmp,latitud,ctrlat,ctrlon

  IF ( totin /= 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  Read in additional parameters for ARPS history dump 4.1 or later
!  version.
!
!-----------------------------------------------------------------------
!
    READ(inch,ERR=110,END=120)                                          &
         nstyp1,  prcin, radin, flxin,snowcin,                           &
         snowin,idummy,idummy,idummy,idummy,                            &
         idummy,idummy,idummy,idummy,idummy,                            &
         idummy,idummy,idummy,idummy,idummy

    IF ( nstyp1 < 1 ) THEN
      nstyp1 = 1
    END IF

    READ(inch,ERR=110,END=120)                                          &
         rdummy,rdummy,rdummy,rdummy,rdummy,                            &
         rdummy,rdummy,rdummy,rdummy,rdummy,                            &
         rdummy,rdummy,rdummy,rdummy,rdummy,                            &
         rdummy,rdummy,rdummy,rdummy,rdummy
  END IF
!
!-----------------------------------------------------------------------
!
!  Read in x,y and z at grid cell centers (scalar points).
!
!----------------------------------------------------------------------
!
  IF( grdin == 1 .OR. grdbas == 1 ) THEN
    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) (x(i),i=1,nx)
    WRITE(lchanl,910) label,' x.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) (y(j),j=1,ny)
    WRITE(lchanl,910) label,' y.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) (z(k),k=1,nz)
    WRITE(lchanl,910) label,' z.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((zp(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' zp.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((zpsoil(i,j,k),i=1,nx),j=1,ny),k=1,nzsoil)
    WRITE(lchanl,910) label,' zpsoil.'

  END IF  ! grdin
!
!-----------------------------------------------------------------------
!
!  Read in base state fields
!
!----------------------------------------------------------------------
!
  IF( basin == 1 .OR. grdbas == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((ubar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' ubar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((vbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' vbar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((wbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' wbar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((ptbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' ptbar.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((pbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' pbar.'


    IF( mstin == 1) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((qvbar(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' qvbar.'
    END IF

    IF(landin == 1) THEN

      IF( nstyp1 <= 1 ) THEN

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120)                                      &
              ((soiltyp(i,j,1),i=1,nx),j=1,ny)
        WRITE(lchanl,910) label,' soiltyp.'

      ELSE

        DO is=1,nstyp1
          IF (is <= nstyps) THEN
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)                                  &
                  ((soiltyp(i,j,is),i=1,nx),j=1,ny)
            WRITE(lchanl,910) label,' soiltyp.'

            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)                                  &
                  ((stypfrct(i,j,is),i=1,nx),j=1,ny)
            WRITE(lchanl,910) label,' stypfrct.'
          ELSE
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)
            WRITE(lchanl,910) label,'skipping soiltyp.'

            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120)
            WRITE(lchanl,910) label,'skipping stypfrct.'
          ENDIF
        END DO

      END IF

      CALL fix_stypfrct_nstyp(nx,ny,nstyp1,nstyp,stypfrct)

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((vegtyp (i,j),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' vegtyp.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((lai    (i,j),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' lai.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((roufns (i,j),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' roufns.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((veg    (i,j),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' veg.'

    END IF

  END IF

  IF( grdbas == 1 ) GO TO 930

  IF( varin == 1 ) THEN

    IF ( totin == 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  Read in uprt, vprt, and wprt
!
!----------------------------------------------------------------------
!
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((uprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' uprt.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((vprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' vprt.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((wprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' wprt.'
!
!-----------------------------------------------------------------------
!
!  Read in scalars
!
!----------------------------------------------------------------------
!
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((ptprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' ptprt.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((pprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' pprt.'

    ELSE

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((uprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' u.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            uprt(i,j,k) = uprt(i,j,k) - ubar(i,j,k)
          END DO
        END DO
      END DO

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((vprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' v.'
      DO k=1,nz-1
        DO j=1,ny
          DO i=1,nx-1
            vprt(i,j,k) = vprt(i,j,k) - vbar(i,j,k)
          END DO
        END DO
      END DO

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((wprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' w.'
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((ptprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' pt.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            ptprt(i,j,k) = ptprt(i,j,k) - ptbar(i,j,k)
          END DO
        END DO
      END DO

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((pprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' p.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            pprt(i,j,k) = pprt(i,j,k) - pbar(i,j,k)
          END DO
        END DO
      END DO

    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  Read in moisture variables
!
!----------------------------------------------------------------------
!
  IF( mstin == 1 ) THEN

    IF ( totin == 0 ) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((qvprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' qvprt.'

    ELSE

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((qvprt(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' qv.'
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            qvprt(i,j,k) = qvprt(i,j,k) - qvbar(i,j,k)
          END DO
        END DO
      END DO

    END IF

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((qc(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' qc.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
        (((qr(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' qr.'

    IF( rainin == 1 ) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
            ((raing(i,j),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' raing.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
            ((rainc(i,j),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' rainc.'

    END IF

    IF( prcin == 1 ) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
            ((prcrate(i,j,1),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' prcrate1.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
            ((prcrate(i,j,2),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' prcrate2.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
            ((prcrate(i,j,3),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' prcrate3.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
            ((prcrate(i,j,4),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' prcrate4.'

    END IF

    IF( icein == 1 ) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((qi(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' qi.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((qs(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' qs.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)                                        &
          (((qh(i,j,k),i=1,nx),j=1,ny),k=1,nz)
      WRITE(lchanl,910) label,' qh.'

    END IF
  END IF

  IF( tkein == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
          (((tke(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' tke.'

  END IF

  IF( trbin == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
          (((kmh(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' kmh.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
          (((kmv(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' kmv.'

  END IF

  IF( sfcin == 1 ) THEN

    IF (nstyp1 <= 1) THEN

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((tsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil)
      WRITE(lchanl,910) label,' tsoil.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((qsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil)
      WRITE(lchanl,910) label,' qsoil.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) ((wetcanp(i,j,0),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' wetcanp.'

    ELSE

      DO is=0,nstyp1
        IF (is <= nstyps) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)(((tsoil(i,j,k,is),i=1,nx),j=1,ny), &
               k=1,nzsoil)
          WRITE(lchanl,910) label,' tsoil.'

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)(((qsoil(i,j,k,is),i=1,nx),j=1,ny), &
               k=1,nzsoil)
          WRITE(lchanl,910) label,' qsoil.'

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)((wetcanp(i,j,is),i=1,nx),j=1,ny)
          WRITE(lchanl,910) label,' wetcanp.'
        ELSE
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)
          WRITE(lchanl,910) label,'skipping tsoil.'

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)
          WRITE(lchanl,910) label,'skipping qsoil.'

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120)
          WRITE(lchanl,910) label,'skipping wetcanp.'
        ENDIF
      END DO

    END IF

    CALL fix_soil_nstyp(nx,ny,nzsoil,nstyp1,nstyp,tsoil,qsoil,wetcanp)

    IF(snowcin == 1) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)
      WRITE(lchanl,910) label,' snowcvr -- discarding.'
    END IF

    IF(snowin == 1) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120)((snowdpth(i,j),i=1,nx),j=1,ny)
      WRITE(lchanl,910) label,' snowdpth.'
    END IF

  END IF

  IF( radin == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120)                                          &
          (((radfrc(i,j,k),i=1,nx),j=1,ny),k=1,nz)
    WRITE(lchanl,910) label,' radfrc.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((radsw(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' radsw.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((rnflx(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' rnflx.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((radswnet(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' radswnet.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((radlwin(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' radlwin.'


  END IF

  IF( flxin == 1 ) THEN

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((usflx(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' usflx.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((vsflx(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' vsflx.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((ptsflx(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' ptsflx.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) ((qvsflx(i,j),i=1,nx),j=1,ny)
    WRITE(lchanl,910) label,' qvsflx.'

  END IF

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

!
!-----------------------------------------------------------------------
!
!  Friendly exit message
!
!----------------------------------------------------------------------
!
  930   CONTINUE

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

  ireturn = 0

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

  110   CONTINUE
  WRITE(6,'(/a/)') ' Error reading data in BN2READ'
  ireturn=1
  RETURN
!
!-----------------------------------------------------------------------
!
!  End-of-file during read
!
!----------------------------------------------------------------------
!

  120   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in BN2READ'
  ireturn=2
  RETURN
END SUBROUTINE bn2read
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE BINDUMP                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE bindump(nx,ny,nz,nzsoil,nstyps, nchanl, grdbas,              & 1,54
           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv,              &
           ubar,vbar,ptbar,pbar,rhobar,qvbar,                           &
           x,y,z,zp,zpsoil,                                             &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsoil,qsoil,wetcanp,snowdpth,                                &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,radswnet,radlwin,                         &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write history data into channel nchanl as binary data.
!
!  All data read in are located at the original staggered grid points.
!
!  Note: coordinate fields are dumped as 3 dimensional fields which
!  have been converted from meters to kilometers.  This is for the
!  convenience of the plotting applications.
!
!  The last 4 characters of the 12 character label written out
!  with each 1-,2-, or 3-d array is used by the splitdump and
!  joinfiles subroutines used by the message passing version of the
!  ARPS (and also by some auxiliary ARPS I/O routines)
!  to determine the data type of the array.
!  Key to the labels:
!
!    'nnnnnnn tdds'
!
!     n  - characters containing the name of the variable.
!     t  - type of variable:  "r" for real and "i" for integer.
!     dd - number of dimensions:  "1d" "2d" or "3d".
!     s  - staggered dimension: "0" for centered,
!                               "1" for staggered in x,
!                               "2" for staggered in y,
!                               "3" for staggered in z.
!
!-----------------------------------------------------------------------
!
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/10/92.
!
!  MODIFICATION HISTORY:
!
!  6/06/92 (M. Xue)
!  Added full documentation.
!
!  7/13/92 (K. Brewster)
!  Added runname, comment and version number writing
!
!  8/23/92 (M. Xue)
!  Modify to perform the dumping of both base and t-dependent arrays
!  and added control on grid staggering.
!
!  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.
!
!  02/06/95 (Y. Liu)
!  Added map projection parameters into the binary dumping
!
!  03/26/96 (G. Bassett)
!  Labels were modified to include information about array type.
!  This information is used by splitdump and joinfiles subroutines.
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  05/15/2002 (J. Brotzge)
!  Added to allow for multiple soil schemes.  
!
!-----------------------------------------------------------------------
!
!  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
!    nzsoil   Number of grid points in the soil  
!
!    nchanl   FORTRAN I/O channel number for history data output.
!    grdbas   Flag indicating if this is a call for the data dump
!             of grid and base state arrays only. If so, grdbas=1.
!
!    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 zonal velocity component (m/s)
!    vbar     Base state meridional 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)
!    zpsoil   Vertical coordinate of grid points in the soil (m)
!
!    soiltyp  Soil type
!    stypfrct  Soil type fraction
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsoil    Soil temperature (K) 
!    qsoil    Soil moisture (m**3/m**3) 
!    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
!    radswnet Net shortwave radiation
!    radlwin  Incoming longwave radiation
!
!    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 (kg/(m**2*s))
!
!  OUTPUT:
!
!    None.
!
!  WORK ARRAY:
!
!    tem1     Temporary work array.
!
!
!-----------------------------------------------------------------------
!
!  The following parameters are passed into this subroutine through
!  a common block in globcst.inc, and they determine which
!  variables are output.
!
!  grdout =0 or 1. If grdout=0, grid variables are not dumped.
!  basout =0 or 1. If basout=0, base state variables are not dumped.
!  varout =0 or 1. If varout=0, model perturbation variables are not dumped.
!  mstout =0 or 1. If mstout=0, water variables are not dumped.
!  rainout=0 or 1. If rainout=0, rain variables are not dumped.
!  prcout =0 or 1. If prcout=0, precipitation rates are not dumped.
!  iceout =0 or 1. If iceout=0, qi, qs and qh are not dumped.
!  tkeout =0 or 1. If tkeout=0, tke is not dumped.
!  trbout =0 or 1. If trbout=0, turbulence parameter km is not dumped.
!  sfcout =0 or 1. If sfcout=0, surface variables are not dumped.
!  landout=0 or 1. If landout=0, surface propertty arrays are not dumped.
!  radout =0 or 1. If radout =0, radiation arrays are not dumped.
!  flxout =0 or 1. If flxout =0, surface flux arrays are not dumped.
!
!  These following parameters are also passed in through common
!  blocks in globcst.inc.
!
!  runname,curtim,umove,vmove,xgrdorg,ygrdorg
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

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

  INTEGER :: nchanl            ! FORTRAN I/O channel number for output
  INTEGER :: grdbas            ! If this is a grid/base state array dump

  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 u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)

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

  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 :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
  REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) 
  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 :: radswnet(nx,ny)      ! Net shortwave radiation
  REAL :: radlwin(nx,ny)       ! Incoming longwave radiation

  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**2*s))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Parameters describing this routine
!
!-----------------------------------------------------------------------
!
! 06/28/2002 Zuwen He
!
! fmtver??: to label each data a version. 
! intver??: an integer to allow faster comparison than fmtver??, 
!           which are strings. 
!
! Verion 5.00: significant change in soil variables since version 4.10. 
! 
  CHARACTER (LEN=40) :: fmtver,fmtver410,fmtver500
  INTEGER  :: intver,intver410,intver500
  PARAMETER (fmtver410='004.10 Binary Data',intver410=410)
  PARAMETER (fmtver500='005.00 Binary Data',intver500=500)   
!
  CHARACTER (LEN=10) :: tmunit
  PARAMETER (tmunit='seconds   ')
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,l,is
  INTEGER :: idummy
  REAL :: rdummy
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'           ! Grid & map parameters.
  INCLUDE 'mp.inc'             ! mpi parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  intver = intver500  !  for the time being, in the future, we will
                      !  allow to dump data in the different version
                      !  intver will be assigned from input file
  IF (intver == intver410) THEN 
    fmtver=fmtver410
  ELSE IF (intver == intver500) THEN 
    fmtver=fmtver500
  ELSE 
    IF (myproc == 0) WRITE(6,'(/1x,a,i10,a/)')                        &
        'Data format, intver=',intver,', not found. The Job stopped.'
    CALL arpsstop('arpstop called from bindump.',1)
  END IF 

  IF (myproc == 0) WRITE(6,'(/1x,a,a,a/)')                        &
      'Data dump format, fmtver=',fmtver,'. ' 

  IF (myproc == 0)  &
     WRITE(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim
!
!-----------------------------------------------------------------------
!
!  Write header info
!
!-----------------------------------------------------------------------
!
  WRITE(nchanl) fmtver
  WRITE(nchanl) runname
  WRITE(nchanl) nocmnt
  IF( nocmnt > 0 ) THEN
    DO l=1,nocmnt
      WRITE(nchanl) cmnt(l)
    END DO
  END IF

  WRITE(nchanl) curtim,tmunit

  IF (intver == intver410) THEN 
    WRITE(nchanl) nx,ny,nz
  ELSE IF (intver == intver500) THEN 
    WRITE(nchanl) nx,ny,nz,nzsoil
  END IF 
!
!-----------------------------------------------------------------------
!
!  Write the flags for different data groups.
!
!-----------------------------------------------------------------------
!
  idummy = 0

  IF( grdbas == 1 ) THEN

    WRITE(nchanl)      1,      1,      0, mstout,      0,               &
                       0,      0,      0, landout,totout,               &
                       0, idummy, idummy, mapproj, month,               &
                     day,   year,   hour, minute, second

  ELSE

    WRITE(nchanl) grdout, basout, varout, mstout, iceout,               &
                  trbout, sfcout, rainout,landout,totout,               &
                  tkeout, idummy, idummy, mapproj, month,               &
                     day,   year,   hour, minute, second

  END IF

  rdummy = 0.0
  WRITE(nchanl)   umove,   vmove, xgrdorg, ygrdorg, trulat1,            &
                trulat2,  trulon,  sclfct,  rdummy,  rdummy,            &
                 rdummy,  rdummy,  rdummy,  rdummy,  rdummy,            &
                  tstop, thisdmp,  latitud, ctrlat,  ctrlon
!
!-----------------------------------------------------------------------
!
!  If totout=1, write additional parameters to history dump files.
!  This is for ARPS version 4.1.2 or later.
!
!-----------------------------------------------------------------------
!
  IF ( totout == 1 ) THEN
    WRITE(nchanl) nstyp,  prcout, radout, flxout,      0,  & ! 0 for snowcvr
              snowout,idummy, idummy, idummy, idummy,                   &
                  idummy, idummy, idummy, idummy, idummy,               &
                  idummy, idummy, idummy, idummy, idummy

    WRITE(nchanl) rdummy, rdummy, rdummy, rdummy, rdummy,               &
                  rdummy, rdummy, rdummy, rdummy, rdummy,               &
                  rdummy, rdummy, rdummy, rdummy, rdummy,               &
                  rdummy, rdummy, rdummy, rdummy, rdummy
  END IF
!
!-----------------------------------------------------------------------
!
!  If grdout=1 or grdbas=1, write out grid variables
!
!-----------------------------------------------------------------------
!
  IF(grdout == 1 .OR. grdbas == 1 ) THEN

    WRITE(nchanl) 'x coord r1d1'
    WRITE(nchanl) x

    WRITE(nchanl) 'y coord r1d2'
    WRITE(nchanl) y

    WRITE(nchanl) 'z coord r1d3'
    WRITE(nchanl) z

    CALL edgfill(zp,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
    WRITE(nchanl) 'zp coor r3d0'
    WRITE(nchanl) zp

    IF (intver == intver500) THEN 
      CALL edgfill(zpsoil,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nzsoil,1,nzsoil)
      WRITE(nchanl) 'zpsoil rs3d0'
      WRITE(nchanl) zpsoil
    END IF 

  END IF    ! grdout
!
!-----------------------------------------------------------------------
!
!  If basout=1, write out base state variables.
!
!-----------------------------------------------------------------------
!
  IF(basout == 1 .OR. grdbas == 1 ) THEN

    CALL edgfill(ubar,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'ubar    r3d1'
    WRITE(nchanl) ubar

    CALL edgfill(vbar,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
    WRITE(nchanl) 'vbar    r3d2'
    WRITE(nchanl) vbar

    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem1(i,j,k) = 0.0
        END DO
      END DO
    END DO
    WRITE(nchanl) 'wbar    r3d3'
    WRITE(nchanl) tem1

    CALL edgfill(ptbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'ptbar   r3d0'
    WRITE(nchanl)  ptbar

    CALL edgfill(pbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'pbar    r3d0'
    WRITE(nchanl)  pbar

    IF(mstout == 1) THEN

      CALL edgfill(qvbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'qvbar   r3d0'
      WRITE(nchanl)  qvbar

    END IF

    IF(landout == 1) THEN

      IF( nstyp <= 1 ) THEN

        CALL iedgfill(soiltyp(1,1,1),1,nx,1,nx-1, 1,ny,1,ny-1,          &
                      1,1,1,1)
        WRITE(nchanl) 'soiltyp i2d0'
        WRITE(nchanl) ((soiltyp(i,j,1),i=1,nx),j=1,ny)

      ELSE
        DO is=1,nstyp

          CALL iedgfill(soiltyp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,       &
                        1,1,1,1)
          WRITE(nchanl) 'soiltyp i2d0'
          WRITE(nchanl) ((soiltyp(i,j,is),i=1,nx),j=1,ny)

          CALL edgfill(stypfrct(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,       &
                       1,1,1,1)
          WRITE(nchanl) 'stypfrc r2d0'
          WRITE(nchanl) ((stypfrct(i,j,is),i=1,nx),j=1,ny)

        END DO

      END IF

      CALL iedgfill(vegtyp ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'vegtyp  i2d0'
      WRITE(nchanl)  vegtyp

      CALL edgfill(lai    ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'lai     r2d0'
      WRITE(nchanl)  lai

      CALL edgfill(roufns ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'roufns  r2d0'
      WRITE(nchanl)  roufns

      CALL edgfill(veg    ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'veg     r2d0'
      WRITE(nchanl)  veg

    END IF

  END IF

  IF ( grdbas == 1 ) RETURN
!
!-----------------------------------------------------------------------
!
!  If varout = 1, Write out uprt, vprt, wprt, ptprt, pprt.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Write out u,v and w.
!
!-----------------------------------------------------------------------
!
  IF(varout == 1) THEN

    IF ( totout == 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  Write out perturbations to history dump
!
!-----------------------------------------------------------------------
!
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx
            tem1(i,j,k)=u(i,j,k)-ubar(i,j,k)
          END DO
        END DO
      END DO
      CALL edgfill(tem1,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'uprt    r3d1'
      WRITE(nchanl) tem1

      DO k=1,nz-1
        DO i=1,nx-1
          DO j=1,ny
            tem1(i,j,k)=v(i,j,k)-vbar(i,j,k)
          END DO
        END DO
      END DO
      CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
      WRITE(nchanl) 'vprt    r3d2'
      WRITE(nchanl) tem1

      CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
      WRITE(nchanl) 'wprt    r3d3'
      WRITE(nchanl) w
!
!-----------------------------------------------------------------------
!
!  Write out scalars
!
!-----------------------------------------------------------------------
!
      CALL edgfill(ptprt,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'ptprt   r3d0'
      WRITE(nchanl) ptprt

      CALL edgfill(pprt,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'pprt    r3d0'
      WRITE(nchanl) pprt

    ELSE
!
!-----------------------------------------------------------------------
!
!  Write out total values to history dump
!
!-----------------------------------------------------------------------
!
      CALL edgfill(u,1,nx,1,nx, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'u       r3d1'
      WRITE(nchanl) u

      CALL edgfill(v,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
      WRITE(nchanl) 'v       r3d2'
      WRITE(nchanl) v

      CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
      WRITE(nchanl) 'w       r3d3'
      WRITE(nchanl) w
!
!-----------------------------------------------------------------------
!
!  Write out scalars
!
!-----------------------------------------------------------------------
!
      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)
          END DO
        END DO
      END DO
      CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'pt      r3d0'
      WRITE(nchanl) tem1

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k) = pbar(i,j,k) + pprt(i,j,k)
          END DO
        END DO
      END DO
      CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'p       r3d0'
      WRITE(nchanl) tem1

    END IF

  END IF     ! varout
!
!-----------------------------------------------------------------------
!
!  If mstout = 1, write out moisture scalars.
!
!-----------------------------------------------------------------------
!
  IF(mstout == 1) THEN

    IF( totout == 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  Write out perturbation to history dump
!
!-----------------------------------------------------------------------
!
      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k)=qv(i,j,k)-qvbar(i,j,k)
          END DO
        END DO
      END DO

      CALL edgfill(tem1,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'qvprt   r3d0'
      WRITE(nchanl)  tem1

    ELSE
!
!-----------------------------------------------------------------------
!
!  Write out total values to history dump
!
!-----------------------------------------------------------------------
!
      CALL edgfill(qv,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'qv      r3d0'
      WRITE(nchanl)  qv

    END IF

    CALL edgfill(qc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'qc      r3d0'
    WRITE(nchanl)  qc

    CALL edgfill(qr,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'qr      r3d0'
    WRITE(nchanl)  qr

    IF(rainout == 1) THEN

      CALL edgfill(raing,   1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'raing   r2d0'
      WRITE(nchanl)  raing

      CALL edgfill(rainc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'rainc   r2d0'
      WRITE(nchanl)  rainc

    END IF   !rainout

    IF ( prcout == 1 ) THEN
      CALL edgfill(prcrate,1,nx,1,nx-1, 1,ny,1,ny-1, 1,4,1,4)
      WRITE(nchanl) 'prcrat1 r2d0'
      WRITE(nchanl) ((prcrate(i,j,1),i=1,nx),j=1,ny)
      WRITE(nchanl) 'prcrat2 r2d0'
      WRITE(nchanl) ((prcrate(i,j,2),i=1,nx),j=1,ny)
      WRITE(nchanl) 'prcrat3 r2d0'
      WRITE(nchanl) ((prcrate(i,j,3),i=1,nx),j=1,ny)
      WRITE(nchanl) 'prcrat4 r2d0'
      WRITE(nchanl) ((prcrate(i,j,4),i=1,nx),j=1,ny)
    END IF

    IF(iceout == 1) THEN

      CALL edgfill(qi,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'qi      r3d0'
      WRITE(nchanl)  qi

      CALL edgfill(qs,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'qs      r3d0'
      WRITE(nchanl)  qs

      CALL edgfill(qh,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      WRITE(nchanl) 'qh      r3d0'
      WRITE(nchanl)  qh

    END IF   !iceout

  END IF   !mstout
!
!-----------------------------------------------------------------------
!
!  If tkeout = 1, write out tke.
!
!-----------------------------------------------------------------------
!
  IF( tkeout == 1 ) THEN

    CALL edgfill(tke,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'tke     r3d0'
    WRITE(nchanl) tke

  END IF
!
!-----------------------------------------------------------------------
!
!  If trbout = 1, write out the turbulence parameter, km.
!
!-----------------------------------------------------------------------
!
  IF( trbout == 1 ) THEN

    CALL edgfill(kmh,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'kmh     r3d0'
    WRITE(nchanl) kmh

    CALL edgfill(kmv,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'kmv     r3d0'
    WRITE(nchanl) kmv

  END IF   ! trbout

!
!-----------------------------------------------------------------------
!
!  If sfcout = 1, write out the surface variables,
!  tsoil, qsoil, and wetcanp.
!
!-----------------------------------------------------------------------
!
  IF( sfcout == 1) THEN

    IF( nstyp <= 1 ) THEN

      CALL edgfill(tsoil(1,1,1,0),  1,nx,1,nx-1, 1,ny,1,ny-1,             &
                   1,nzsoil,1,nzsoil)
      CALL edgfill(qsoil(1,1,1,0), 1,nx,1,nx-1, 1,ny,1,ny-1,             &
                   1,nzsoil,1,nzsoil)

      IF (intver == intver410) THEN 
!
! 06/28/2002 Zuwen He
!
! In version 410, tsoil is the average of tsoil from 2 to nzsoil in later 
! version, and wetdp is similar. 
!
        tem1(:,:,1)=0.

        DO k=2,nzsoil
          DO j=1,ny
            DO i=1,nx
              tem1(i,j,1)=tem1(i,j,1)+tsoil(i,j,k,0)
            END DO 
          END DO 
        END DO 

        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) 
          END DO 
        END DO 

        WRITE(nchanl) 'tsfc   r2d0'
        WRITE(nchanl) ((tsoil(i,j,1,0),i=1,nx),j=1,ny)
        WRITE(nchanl) 'tsoil  r2d0'
        WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny)

        tem1(:,:,1)=0.

        DO k=2,nzsoil
          DO j=1,ny
            DO i=1,nx
              tem1(i,j,1)=tem1(i,j,1)+qsoil(i,j,k,0)
            END DO 
          END DO 
        END DO 

        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) 
          END DO 
        END DO 

        WRITE(nchanl) 'wetsfc r2d0'
        WRITE(nchanl) ((qsoil(i,j,1,0),i=1,nx),j=1,ny)
        WRITE(nchanl) 'wetdp  r2d0'
        WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny)
      
      ELSE IF (intver == intver500) THEN 

        WRITE(nchanl) 'tsoil  rs3d0'
        WRITE(nchanl) (((tsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) 
        WRITE(nchanl) 'qsoil  rs3d0'
        WRITE(nchanl) (((qsoil(i,j,k,0),i=1,nx),j=1,ny),k=1,nzsoil) 
      
      END IF  ! intver

      CALL edgfill(wetcanp(1,1,0),1,nx,1,nx-1, 1,ny,1,ny-1,             &
                   1,1,1,1)
      WRITE(nchanl) 'wetcanp r2d0'
      WRITE(nchanl) ((wetcanp(i,j,0),i=1,nx),j=1,ny)

    ELSE

      DO is=0,nstyp

        CALL edgfill(tsoil(1,1,1,is),  1,nx,1,nx-1, 1,ny,1,ny-1,          &
                     1,nzsoil,1,nzsoil)
        CALL edgfill(qsoil(1,1,1,is), 1,nx,1,nx-1, 1,ny,1,ny-1,          &
                     1,nzsoil,1,nzsoil)

      IF (intver == intver410) THEN 
!
! 06/28/2002 Zuwen He
!
! In version 410, tsoil is the average of tsoil from 2 to nzsoil in later 
! version, and wetdp is similar. 
!
        tem1(:,:,1)=0.

        DO k=2,nzsoil
          DO j=1,ny
            DO i=1,nx
              tem1(i,j,1)=tem1(i,j,1)+tsoil(i,j,k,is)
            END DO 
          END DO 
        END DO 

        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) 
          END DO 
        END DO 

        WRITE(nchanl) 'tsfc  r2d0'
        WRITE(nchanl) ((tsoil(i,j,1,is),i=1,nx),j=1,ny)
        WRITE(nchanl) 'tsoil  r2d0'
        WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny)

        tem1(:,:,1)=0.

        DO k=2,nzsoil
          DO j=1,ny
            DO i=1,nx
              tem1(i,j,1)=tem1(i,j,1)+qsoil(i,j,k,is)
            END DO 
          END DO 
        END DO 

        DO j=1,ny
          DO i=1,nx
            tem1(i,j,1)=tem1(i,j,1)/float(nzsoil-1) 
          END DO 
        END DO 

        WRITE(nchanl) 'wetsfc r2d0'
        WRITE(nchanl) ((qsoil(i,j,1,is),i=1,nx),j=1,ny)
        WRITE(nchanl) 'wetdp  r2d0'
        WRITE(nchanl) ((tem1(i,j,1),i=1,nx),j=1,ny)
      
        ELSE IF (intver == intver500) THEN 

          WRITE(nchanl) 'tsoil  rs3d0'
          WRITE(nchanl) (((tsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) 
          WRITE(nchanl) 'qsoil  rs3d0'
          WRITE(nchanl) (((qsoil(i,j,k,is),i=1,nx),j=1,ny),k=1,nzsoil) 

        END IF  ! intver

        CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,          &
                     1,1,1,1)
        WRITE(nchanl) 'wetcanp r2d0'
        WRITE(nchanl) ((wetcanp(i,j,is),i=1,nx),j=1,ny)

      END DO
    END IF

    IF (snowout == 1) THEN

      CALL edgfill(snowdpth,1,nx,1,nx-1, 1,ny,1,ny-1,                   &
                    1,1,1,1)
      WRITE(nchanl) 'snowdpthr2d0'
      WRITE(nchanl) ((snowdpth(i,j),i=1,nx),j=1,ny)
    END IF

  END IF
!
!-----------------------------------------------------------------------
!
!  If radout = 1, write out the radiation arrays
!
!-----------------------------------------------------------------------
!
  IF( radout == 1 ) THEN

    CALL edgfill(radfrc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    WRITE(nchanl) 'radfrc  r3d0'
    WRITE(nchanl) radfrc

    CALL edgfill(radsw,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    WRITE(nchanl) 'radsw   r2d0'
    WRITE(nchanl) radsw

    CALL edgfill(rnflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    WRITE(nchanl) 'rnflx   r2d0'
    WRITE(nchanl) rnflx

    IF (intver >= intver500) THEN 

      CALL edgfill(radswnet,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'radswnetr2d0'
      WRITE(nchanl) radswnet 
  
      CALL edgfill(radlwin,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      WRITE(nchanl) 'radlwin r2d0'
      WRITE(nchanl) radlwin

    END IF  ! intver 

  END IF   ! radout
!
!-----------------------------------------------------------------------
!
!  If flxout = 1, write out the surface fluxes
!
!-----------------------------------------------------------------------
!
  IF( flxout == 1 ) THEN

    CALL edgfill(usflx,1,nx,1,nx, 1,ny,1,ny-1, 1,1,1,1)
    WRITE(nchanl) 'usflx   r2d0'
    WRITE(nchanl) usflx

    CALL edgfill(vsflx,1,nx,1,nx-1, 1,ny,1,ny, 1,1,1,1)
    WRITE(nchanl) 'vsflx   r2d0'
    WRITE(nchanl) vsflx

    CALL edgfill(ptsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    WRITE(nchanl) 'ptsflx  r2d0'
    WRITE(nchanl) ptsflx

    CALL edgfill(qvsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    WRITE(nchanl) 'qvsflx  r2d0'
    WRITE(nchanl) qvsflx

  END IF   ! flxout

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


SUBROUTINE bn2dump(nx,ny,nz,nzsoil,nstyps, nchanl, grdbas,              & 1
           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,tke,kmh,kmv,              &
           ubar,vbar,ptbar,pbar,rhobar,qvbar,                           &
           x,y,z,zp,zpsoil,                                             &
           soiltyp,stypfrct,vegtyp,lai,roufns,veg,                      &
           tsoil,qsoil,wetcanp,snowdpth,                                &
           raing,rainc,prcrate,                                         &
           radfrc,radsw,rnflx,radswnet,radlwin,                         &
           usflx,vsflx,ptsflx,qvsflx,                                   &
           tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write history data into channel nchanl as binary data.
!
!  This routine can dump the data arrays in a model subdomain and
!  at selected data points.
!
!  All data read in are located at the original staggered grid points.
!
!  Note: coordinate fields are dumped as 3 dimensional fields which
!  have been converted from meters to kilometers.  This is for the
!  convenience of the plotting applications.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/10/92.
!
!  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.
!
!  02/06/95 (Y. Liu)
!  Added map projection parameters into the second binary dumping
!
!  12/09/1998 (Donghai Wang)
!  Added the snow cover.
!
!  05/15/2002 (J. Brotzge)
!  Added to allow for multiple soil schemes.  
!
!-----------------------------------------------------------------------
!
!  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
!    nzsoil   Number of grid points in the soil  
!
!    nchanl   FORTRAN I/O channel number for history data output.
!    grdbas   Flag indicating if this is a call for the data dump
!             of grid and base state arrays only. If so, grdbas=1.
!
!    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 zonal velocity component (m/s)
!    vbar     Base state meridional 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)
!    zpsoil   Vertical coordinate of grid points in soil (m)  
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsoil    Soil temperature (K) 
!    qsoil    Soil moisture (m**3/m**3)  
!    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
!    radswnet Net shortwave radiation
!    radlwin  Incoming longwave radiation
!
!    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.
!
!
!-----------------------------------------------------------------------
!
!  The following parameters are passed into this subroutine through
!  a common block in globcst.inc.  These parameters determine which
!  variables are output.
!
!  grdout =0 or 1. If grdout=0, grid variables are not dumped.
!  basout =0 or 1. If basout=0, base state variables are not dumped.
!  varout =0 or 1. If varout=0, model perturbation variables are not dumped.
!  mstout =0 or 1. If mstout=0, water variables are not dumped.
!  rainout=0 or 1. If rainout=0, rain variables are not dumped.
!  prcout =0 or 1. If prcout=0, precipitation rates are not dumped.
!  iceout =0 or 1. If iceout=0, qi, qs and qh are not dumped.
!  trbout =0 or 1. If trbout=0, turbulence parameter km is not dumped.
!  tkeout =0 or 1. If tkeout=0, tke is not dumped.
!  radout =0 or 1. If radout=0, radiation arrays are not dumped.
!  flxout =0 or 1. If flxout=0, surface fluxes are not dumped.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

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

  INTEGER :: nchanl            ! FORTRAN I/O channel number for output
  INTEGER :: grdbas            ! If this is a grid/base state array dump

  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 u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)

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

  INTEGER :: nstyps                  ! Number of soil types
  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 :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K)
  REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) 
  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 :: radswnet(nx,ny)      ! Net shortwave radiation
  REAL :: radlwin(nx,ny)       ! Incoming longwave radiation

  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**2*s))
  REAL :: qvsflx(nx,ny)        ! Surface moisture flux (kg/(m**2*s))

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Parameters describing this routine
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=40) :: fmtver
  PARAMETER (fmtver='004.10 2nd Binary Data')
  CHARACTER (LEN=10) :: tmunit
  PARAMETER (tmunit='seconds   ')
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,l,is, idummy
  REAL :: rdummy

  INTEGER :: nxout,nyout,nzout ! The size of array to be written out.
  INTEGER :: nzsoilout         ! The size of array to be written out.


  INTEGER :: ist ,ind ,isk ,jst ,jnd ,jsk ,kst ,knd ,ksk
  INTEGER :: ist1,ind1,isk1,jst1,jnd1,jsk1,kst1,knd1,ksk1

  INTEGER :: lst, lnd, lsk 

  INTEGER :: setdomn,setskip
  SAVE setdomn, setskip
  SAVE ist,ind,isk,jst,jnd,jsk,kst,knd,ksk,lst,lnd,lsk 

  DATA setdomn, setskip /0,0/
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  IF( setdomn == 0) THEN  ! If these parameters are nevers set ...

    ist = 1
    jst = 1
    kst = 1
    lst = 1 
    ind = nx
    jnd = ny
    knd = nz
    lnd = nzsoil 

  END IF

  IF( setskip == 0) THEN  ! If these parameters are nevers set ...

    isk = 1
    jsk = 1
    ksk = 1
    lsk = 1 

  END IF

  WRITE(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim
!
!-----------------------------------------------------------------------
!
!  Write header info
!
!-----------------------------------------------------------------------
!
  WRITE(nchanl) fmtver
  WRITE(nchanl) runname
  WRITE(nchanl) nocmnt
  IF( nocmnt > 0 ) THEN
    DO l=1,nocmnt
      WRITE(nchanl) cmnt(l)
    END DO
  END IF
  WRITE(nchanl) curtim,tmunit

  nxout = ist+(ind-ist)/isk
  nyout = jst+(jnd-jst)/jsk
  nzout = kst+(knd-kst)/ksk
  nzsoilout = lst+(lnd-lst)/lsk 

  WRITE(nchanl) nxout, nyout, nzout, nzsoilout 
  PRINT*,'nxout= ',nxout,'nyout= ',nyout,'nzout= ',nzout,            &
         'nzsoilout= ',nzsoilout 
!
!-----------------------------------------------------------------------
!
!  Write the flags for different data groups.
!
!-----------------------------------------------------------------------
!
  idummy = 0

  IF( grdbas == 1 ) THEN

    WRITE(nchanl)      1,      1,      0, mstout,      0,               &
                       0, idummy, idummy, landout, totout,              &
                  idummy, idummy, idummy, mapproj, month,               &
                     day,   year,   hour, minute, second

  ELSE

    WRITE(nchanl) grdout, basout, varout, mstout, iceout,               &
                  trbout, sfcout, rainout,landout,totout,               &
                  tkeout, idummy, idummy, mapproj, month,               &
                     day,   year,   hour, minute, second

  END IF

  rdummy = 0.0
  WRITE(nchanl)   umove,   vmove, xgrdorg, ygrdorg, trulat1,            &
                trulat2,  trulon,  sclfct,  rdummy,  rdummy,            &
                 rdummy,  rdummy,  rdummy,  rdummy,  rdummy,            &
                  tstop, thisdmp,  latitud, ctrlat,  ctrlon

  IF ( totout /= 0 ) THEN
!
!-----------------------------------------------------------------------
!
!  Add new parameters for new version of history data dump
!
!-----------------------------------------------------------------------
!
    WRITE(nchanl) nstyp,  prcout, radout, flxout,      0,  & ! 0 for snowcvr
              snowout,idummy, idummy, idummy, idummy,                   &
                  idummy, idummy, idummy, idummy, idummy,               &
                  idummy, idummy, idummy, idummy, idummy

    WRITE(nchanl) rdummy, rdummy, rdummy, rdummy, rdummy,               &
                  rdummy, rdummy, rdummy, rdummy, rdummy,               &
                  rdummy, rdummy, rdummy, rdummy, rdummy,               &
                  rdummy, rdummy, rdummy, rdummy, rdummy

  END IF
!
!-----------------------------------------------------------------------
!
!  If grdout=1 or grdbas=1, write out grid variables
!
!-----------------------------------------------------------------------
!
  IF(grdout == 1 .OR. grdbas == 1 ) THEN

    WRITE(nchanl) 'x coordinate'
    WRITE(nchanl) ( x(i), i=ist,ind,isk)

    WRITE(nchanl) 'y coordinate'
    WRITE(nchanl) ( y(j), j=jst,jnd,jsk)

    WRITE(nchanl) 'z coordinate'
    WRITE(nchanl) ( z(k), k=kst,knd,ksk)

    WRITE(nchanl) 'zp coord    '
    WRITE(nchanl) ((( zp(i,j,k),                                        &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    WRITE(nchanl) 'zpsoil coord    '
    WRITE(nchanl) (((zpsoil(i,j,l),i=ist,ind,isk),j=jst,jnd,jsk),       &
                  l=lst,lnd,lsk)

  END IF    ! grdout
!
!-----------------------------------------------------------------------
!
!  If basout=1, write out base state variables.
!
!-----------------------------------------------------------------------
!
  IF(basout == 1 .OR. grdbas == 1 ) THEN

    WRITE(nchanl) 'ubar        '
    WRITE(nchanl) ((( ubar(i,j,k),                                      &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    WRITE(nchanl) 'vbar        '
    WRITE(nchanl) ((( vbar(i,j,k),                                      &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    DO k=kst,knd,ksk
      DO j=jst,jnd,jsk
        DO i=ist,ind,isk
          tem1(i,j,k) = 0.0
        END DO
      END DO
    END DO
    WRITE(nchanl) 'wbar        '
    WRITE(nchanl) (((tem1(i,j,k),                                       &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    WRITE(nchanl) 'ptbar       '
    WRITE(nchanl) (((ptbar(i,j,k),                                      &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    WRITE(nchanl) 'pbar        '
    WRITE(nchanl) (((pbar(i,j,k),                                       &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    IF(mstout == 1) THEN

      WRITE(nchanl) 'qvbar       '
      WRITE(nchanl) (((qvbar(i,j,k),                                    &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    END IF

    IF(landout == 1) THEN

      IF( nstyp <= 1 ) THEN
        WRITE(nchanl) 'soiltyp     '
        WRITE(nchanl) ((soiltyp(i,j,1),i=ist,ind,isk),                  &
                        j=jst,jnd,jsk)
      ELSE
        DO is=1,nstyp
          WRITE(nchanl) 'soiltyp     '
          WRITE(nchanl) ((soiltyp(i,j,is),i=ist,ind,isk),               &
                        j=jst,jnd,jsk)

          WRITE(nchanl) 'stypfrct    '
          WRITE(nchanl) ((stypfrct(i,j,is),i=ist,ind,isk),              &
                        j=jst,jnd,jsk)
        END DO
      END IF

      WRITE(nchanl) 'vegtyp      '
      WRITE(nchanl) ((vegtyp (i,j),i=ist,ind,isk),j=jst,jnd,jsk)

      WRITE(nchanl) 'lai         '
      WRITE(nchanl) ((lai    (i,j),i=ist,ind,isk),j=jst,jnd,jsk)

      WRITE(nchanl) 'roufns      '
      WRITE(nchanl) ((roufns (i,j),i=ist,ind,isk),j=jst,jnd,jsk)

      WRITE(nchanl) 'veg         '
      WRITE(nchanl) ((veg    (i,j),i=ist,ind,isk),j=jst,jnd,jsk)

    END IF

  END IF


  IF ( grdbas == 1 ) RETURN

!
!-----------------------------------------------------------------------
!
!  If varout = 1, Write out uprt, vprt, wprt, ptprt, pprt.
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Write out u, v and w
!
!-----------------------------------------------------------------------
!
  IF(varout == 1) THEN

    IF ( totout == 0 ) THEN
!
!-----------------------------------------------------------------------
!
!    Write out perturbatios to history dump
!
!-----------------------------------------------------------------------
!
      WRITE(nchanl) 'uprt        '
      WRITE(nchanl) ((( u(i,j,k)-ubar(i,j,k),                           &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'vprt        '
      WRITE(nchanl) ((( v(i,j,k)-vbar(i,j,k),                           &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'wprt        '
      WRITE(nchanl) ((( w(i,j,k),                                       &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)
!
!-----------------------------------------------------------------------
!
!  Write out scalars
!
!-----------------------------------------------------------------------
!
      WRITE(nchanl) 'ptprt       '
      WRITE(nchanl) ((( ptprt(i,j,k),                                   &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'pprt        '
      WRITE(nchanl) ((( pprt(i,j,k),                                    &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    ELSE
!
!-----------------------------------------------------------------------
!
!    Write out total values to history dump
!
!-----------------------------------------------------------------------
!
      WRITE(nchanl) 'u           '
      WRITE(nchanl) ((( u(i,j,k),                                       &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'v           '
      WRITE(nchanl) ((( v(i,j,k),                                       &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'w           '
      WRITE(nchanl) ((( w(i,j,k),                                       &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)
!
!-----------------------------------------------------------------------
!
!  Write out scalars
!
!-----------------------------------------------------------------------
!
      WRITE(nchanl) 'pt          '
      WRITE(nchanl) ((( ptprt(i,j,k)+ptbar(i,j,k),                      &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'p           '
      WRITE(nchanl) ((( pprt(i,j,k)+pbar(i,j,k),                        &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    END IF

  END IF     ! varout
!
!-----------------------------------------------------------------------
!
!  If mstout = 1, write out moisture scalars.
!
!-----------------------------------------------------------------------
!
  IF(mstout == 1) THEN

    IF ( totout == 0 ) THEN
!
!-----------------------------------------------------------------------
!
!    Write out perturbations to history dump
!
!-----------------------------------------------------------------------
!
      WRITE(nchanl) 'qvprt       '
      WRITE(nchanl) ((( qv(i,j,k)-qvbar(i,j,k),                         &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    ELSE
!
!-----------------------------------------------------------------------
!
!    Write out total values to history dump
!
!-----------------------------------------------------------------------
!
      WRITE(nchanl) 'qv          '
      WRITE(nchanl) ((( qv(i,j,k),                                      &
                    i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    END IF

    WRITE(nchanl) 'qc          '
    WRITE(nchanl) ((( qc(i,j,k),                                        &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    WRITE(nchanl) 'qr          '
    WRITE(nchanl) ((( qr(i,j,k),                                        &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    IF(rainout == 1) THEN

      WRITE(nchanl) 'raing       '
      WRITE(nchanl) ((raing(i,j),i=ist,ind,isk),j=jst,jnd,jsk)

      WRITE(nchanl) 'rainc       '
      WRITE(nchanl) ((rainc(i,j),i=ist,ind,isk),j=jst,jnd,jsk)

    END IF   !rainout

    IF ( prcout == 1 ) THEN

      WRITE(nchanl) 'prcrate1    '
      WRITE(nchanl) ((prcrate(i,j,1),i=ist,ind,isk),j=jst,jnd,jsk)
      WRITE(nchanl) 'prcrate2    '
      WRITE(nchanl) ((prcrate(i,j,2),i=ist,ind,isk),j=jst,jnd,jsk)
      WRITE(nchanl) 'prcrate3    '
      WRITE(nchanl) ((prcrate(i,j,3),i=ist,ind,isk),j=jst,jnd,jsk)
      WRITE(nchanl) 'prcrate4    '
      WRITE(nchanl) ((prcrate(i,j,4),i=ist,ind,isk),j=jst,jnd,jsk)

    END IF   ! prcout

    IF(iceout == 1) THEN

      WRITE(nchanl) 'qi          '
      WRITE(nchanl) ((( qi(i,j,k),                                      &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'qs          '
      WRITE(nchanl) ((( qs(i,j,k),                                      &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

      WRITE(nchanl) 'qh          '
      WRITE(nchanl) ((( qh(i,j,k),                                      &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    END IF   !iceout

  END IF   !mstout
!
!-----------------------------------------------------------------------
!
!  If tkeout = 1, write out turbulence parameter, km.
!
!-----------------------------------------------------------------------
!
  IF( tkeout == 1 ) THEN

    WRITE(nchanl) 'tke          '
    WRITE(nchanl) ((( tke(i,j,k),                                       &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

  END IF   ! tkeout

!
!-----------------------------------------------------------------------
!
!  If trbout = 1, write out turbulence parameter, km.
!
!-----------------------------------------------------------------------
!
  IF( trbout == 1 ) THEN

    WRITE(nchanl) 'kmh          '
    WRITE(nchanl) ((( kmh(i,j,k),                                       &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    WRITE(nchanl) 'kmv          '
    WRITE(nchanl) ((( kmv(i,j,k),                                       &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

  END IF   ! trbout
!
!-----------------------------------------------------------------------
!
!  If sfcout = 1, write out the surface variables, tsoil,
!  qsoil, and wetcanp.
!
!-----------------------------------------------------------------------
!
  IF( sfcout == 1) THEN

    IF( nstyp <= 1 ) THEN

      WRITE(nchanl) 'tsoil         '
      WRITE(nchanl) (((tsoil(i,j,l,0),i=ist,ind,isk),j=jst,jnd,jsk),  &
                     l=lst,lnd,lsk) 

      WRITE(nchanl) 'qsoil       '
      WRITE(nchanl) (((qsoil(i,j,l,0),i=ist,ind,isk),j=jst,jnd,jsk),  &
                     l=lst,lnd,lsk) 

      WRITE(nchanl) 'wetcanp      '
      WRITE(nchanl) ((wetcanp(i,j,0),i=ist,ind,isk),j=jst,jnd,jsk)

    ELSE

      DO is=0,nstyp
        WRITE(nchanl) 'tsoil         '
        WRITE(nchanl) (((tsoil(i,j,l,is),i=ist,ind,isk),                   &
                      j=jst,jnd,jsk),l=lst,lnd,lsk) 

        WRITE(nchanl) 'qsoil       '
        WRITE(nchanl) (((qsoil(i,j,l,is),i=ist,ind,isk),                  &
                      j=jst,jnd,jsk),l=lst,lnd,lsk) 

        WRITE(nchanl) 'wetcanp      '
        WRITE(nchanl) ((wetcanp(i,j,is),i=ist,ind,isk),                 &
                      j=jst,jnd,jsk)
      END DO

    END IF

    IF(snowout == 1) THEN

      WRITE(nchanl) 'snowdpth     '
      WRITE(nchanl) ((snowdpth(i,j),i=ist,ind,isk),                     &
                    j=jst,jnd,jsk)
    END IF

  END IF   ! sfcout done
!
!-----------------------------------------------------------------------
!
!  If radout = 1, write out radiation arrays, radfrc, radsw and
!  rnflx.
!
!-----------------------------------------------------------------------
!
  IF( radout == 1 ) THEN

    WRITE(nchanl) 'radfrc       '
    WRITE(nchanl) ((( radfrc(i,j,k),                                    &
                  i=ist,ind,isk),j=jst,jnd,jsk),k=kst,knd,ksk)

    WRITE(nchanl) 'radsw        '
    WRITE(nchanl) (( radsw(i,j),                                        &
                  i=ist,ind,isk),j=jst,jnd,jsk)

    WRITE(nchanl) 'rnflx'
    WRITE(nchanl) (( rnflx(i,j),                                        &
                  i=ist,ind,isk),j=jst,jnd,jsk)

    WRITE(nchanl) 'radswnet        '
    WRITE(nchanl) (( radswnet(i,j),                                     &
                  i=ist,ind,isk),j=jst,jnd,jsk)

    WRITE(nchanl) 'radlwin        '
    WRITE(nchanl) (( radlwin(i,j),                                      &
                  i=ist,ind,isk),j=jst,jnd,jsk)


  END IF   ! radout
!
!-----------------------------------------------------------------------
!
!  If flxout = 1, write out surface fluxes
!
!-----------------------------------------------------------------------
!
  IF( flxout == 1 ) THEN

    WRITE(nchanl) 'usflx        '
    WRITE(nchanl) (( usflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk)

    WRITE(nchanl) 'vsflx        '
    WRITE(nchanl) (( vsflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk)

    WRITE(nchanl) 'ptsflx       '
    WRITE(nchanl) (( ptsflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk)

    WRITE(nchanl) 'qvsflx       '
    WRITE(nchanl) (( qvsflx(i,j),i=ist,ind,isk),j=jst,jnd,jsk)

  END IF   ! flxout

  RETURN

  ENTRY bdmpdomn(ist1,ind1,jst1,jnd1,kst1,knd1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  To set the start and end indicies of the model subdomain
!  in which the data is dumped out.
!
!-----------------------------------------------------------------------
!
  ist = ist1
  jst = jst1
  kst = kst1
  ind = ind1
  jnd = jnd1
  knd = knd1

  setdomn = 1

  RETURN

  ENTRY bdmpskip(isk1, jsk1, ksk1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  To set data skip parameters for data dump.
!
!-----------------------------------------------------------------------
!

  isk = isk1
  jsk = jsk1
  ksk = ksk1

  setskip = 1

  RETURN
END SUBROUTINE bn2dump