!
!##################################################################
!##################################################################
!######                                                      ######
!######                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,2)=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 BINREADSPLIT               ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE binreadsplit(nx,ny,nz,nzsoil,nstyps,grdbas,inch,time,x,y,z,  & 2,120
           zp,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: Yunheng Wang
!  9/03/02.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  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
!
!-----------------------------------------------------------------------
!
! 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

  INTEGER :: nxlg, nylg, nzlg, nzsoillg, n3d
  INTEGER, ALLOCATABLE :: var3di(:,:,:)
  REAL, ALLOCATABLE :: var2d(:,:), var3d(:,:,:)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'indtflg.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'mp.inc'            ! mpi parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nxlg = (nx-3)*nproc_x+3
  nylg = (ny-3)*nproc_y+3
  nzlg = nz
  nzsoillg = nzsoil
  n3d = MAX(nz, nzsoil, nstyps+1, 4)

  ALLOCATE(var2d(nxlg,nylg))
  ALLOCATE(var3d(nxlg,nylg,n3d))
  ALLOCATE(var3di(nxlg,nylg,n3d))

  IF (myproc == 0) THEN

!-----------------------------------------------------------------------
!
!  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 
      WRITE(6,'(/1x,a,a,a/)')                                     &
          'Incoming data format, fmtverin=',fmtverin,                 & 
          ', not found. The Job stopped.'
      CALL arpsstop('arpstop called from BINREADSPLIT. ',1)
    END IF 
  
    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
  
    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 BINREADSPLIT
!
!-----------------------------------------------------------------------
!
    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 /= nxlg .OR. nyin /= nylg .OR. nzin /= nzlg) THEN
      WRITE(6,'(1x,a)')                                                   &
           ' Dimensions in BINREADSPLIT 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 BINREADSPLIT.'
      CALL arpsstop('arpstop called from BINREADSPLIT nx-ny-nz read ',1)
    END IF
   
    IF(nzsoilin /= nzsoillg) THEN
  
      IF (intver <= intver410) THEN 
        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 BINREADSPLIT.'
      ELSE IF (intver >= intver500) THEN 
        WRITE(6,'(1x,a)') & 
             ' Dimensions in BINREADSPLIT inconsistent with data.'
        WRITE(6,'(1x,a,I15)') ' Read were: ', nzsoilin
        WRITE(6,'(1x,a,I15)') ' Expected:  ', nzsoillg
        WRITE(6,'(1x,a)') ' Program aborted in BINREADSPLIT.'
      END IF
      CALL arpsstop('arpstop called from BINREADSPLIT nzsoil read ',1)
  
    END IF 

    nxin = (nxin-3)/nproc_x+3
    nyin = (nyin-3)/nproc_y+3

  END IF

  CALL mpupdatei(intver, 1)
  CALL mpupdatec(runname, 40)
  CALL mpupdater(time, 1)
  CALL mpupdatec(tmunit, 10)

  CALL mpupdatei(nxin, 1)
  CALL mpupdatei(nyin, 1)
  CALL mpupdatei(nzin, 1)
  CALL mpupdatei(nzsoilin, 1)
!
!-----------------------------------------------------------------------
!
!  Read in flags for different data groups
!
!-----------------------------------------------------------------------
!
  IF (myproc == 0)  THEN
    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.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

  END IF  ! myproc == 0
  CALL mpupdatei(mstin,1)
  CALL mpupdatei(landin,1)
  CALL mpupdatei(totin,1)
  CALL mpupdatei(mapproj,1)
  CALL mpupdatei(month,1)
  CALL mpupdatei(day,1)
  CALL mpupdatei(year,1)
  CALL mpupdatei(hour,1)
  CALL mpupdatei(minute,1)
  CALL mpupdatei(second,1)
  IF(grdbas == 1) THEN
    CALL mpupdatei(bgrdin,1)
    CALL mpupdatei(bbasin,1)
    CALL mpupdatei(bvarin,1)
    CALL mpupdatei(btrbin,1)
    CALL mpupdatei(btkein,1)
  ELSE
    CALL mpupdatei(grdin,1)
    CALL mpupdatei(basin,1)
    CALL mpupdatei(varin,1)
    CALL mpupdatei(trbin,1)
    CALL mpupdatei(tkein,1)
    CALL mpupdatei(icein,1)
    CALL mpupdatei(sfcin,1)
    CALL mpupdatei(rainin,1)
  END IF
  CALL mpupdater(umove,1)
  CALL mpupdater(vmove,1)
  CALL mpupdater(xgrdorg,1)
  CALL mpupdater(ygrdorg,1)
  CALL mpupdater(trulat1,1)
  CALL mpupdater(trulat2,1)
  CALL mpupdater(trulon,1)
  CALL mpupdater(sclfct,1)
  CALL mpupdater(tstop,1)
  CALL mpupdater(thisdmp,1)
  CALL mpupdater(latitud,1)
  CALL mpupdater(ctrlat,1)
  CALL mpupdater(ctrlon,1)
  IF(totin /= 0) THEN
    CALL mpupdatei(nstyp1,1)
    CALL mpupdatei(prcin,1)
    CALL mpupdatei(radin,1)
    CALL mpupdatei(flxin,1)
    CALL mpupdatei(snowcin,1)
    CALL mpupdatei(snowin,1)
  END IF
!
!-----------------------------------------------------------------------
!
!  Read in x, y, z and zp arrays.
!
!----------------------------------------------------------------------
!
  IF( grdin == 1 .OR. grdbas == 1 ) THEN
    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (var2d(i,1),i=1,nxlg)
      WRITE(lchanl,910) label,' x.'
    END IF
    CALL mpisplit1dx(var2d(:,1),nx,x)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (var2d(1,j),j=1,nylg)
      WRITE(lchanl,910) label,' y.'
    END IF
    CALL mpisplit1dy(var2d(1,:),ny,y)

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

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' zp.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,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,2)=zpsoil(i,j,1)-1.
        END DO 
      END DO 

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

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg)
        WRITE(lchanl,910) label,' zpsoil.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nzsoil,zpsoil)

    END IF  ! intver

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

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' ubar.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,ubar)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' vbar.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,vbar)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' wbar.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,wbar)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' ptbar.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,ptbar)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' pbar.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,pbar)

    IF( mstin == 1) THEN
      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' qvbar.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,qvbar)
    END IF

    IF (landin == 1) THEN

      IF (nstyp1 <= 1) THEN

        IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) ((var3di(i,j,1),i=1,nxlg),j=1,nylg)
          WRITE(lchanl,910) label,' soiltyp.'
        END IF
        CALL mpisplit3di(var3di,nx,ny,1,soiltyp(:,:,1))

      ELSE

        IF(myproc == 0) THEN
          var3di(:,:,:) = 0
          var3d(:,:,:)  = 0.0
          DO is=1,nstyp1
            IF (is <= nstyps) THEN
              READ(inch,ERR=110,END=120) label
              READ(inch,ERR=110,END=120) ((var3di(i,j,is),i=1,nxlg),j=1,nylg)
              WRITE(lchanl,910) label,' soiltyp.'
  
              READ(inch,ERR=110,END=120) label
              READ(inch,ERR=110,END=120) ((var3d(i,j,is),i=1,nxlg),j=1,nylg)
              WRITE(lchanl,910) label,' stypfrct.'
            ELSE
              READ(inch,ERR=110,END=120) label
              WRITE(lchanl,910) label,'skipping soiltyp'
              READ(inch,ERR=110,END=120)
              READ(inch,ERR=110,END=120) label
              WRITE(lchanl,910) label,'skipping stypfrct.'
              READ(inch,ERR=110,END=120) 
            END IF
          END DO
        ENDIF
        CALL mpisplit3di(var3di,nx,ny,nstyps,soiltyp)
        CALL mpisplit3d(var3d,nx,ny,nstyps,stypfrct)
      END IF

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

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) ((var3di(i,j,1),i=1,nxlg),j=1,nylg)
        WRITE(lchanl,910) label,' vegtyp.'
      END IF
      CALL mpisplit3di(var3di,nx,ny,1,vegtyp)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) var2d
        WRITE(lchanl,910) label,' lai.'
      END IF
      CALL mpisplit3d(var2d,nx,ny,1,lai)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) var2d
        WRITE(lchanl,910) label,' roufns.'
      END IF
      CALL mpisplit3d(var2d,nx,ny,1,roufns)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) var2d
        WRITE(lchanl,910) label,' veg.'
      END IF
      CALL mpisplit3d(var2d,nx,ny,1,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
!
!-----------------------------------------------------------------------
!
      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' uprt.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,uprt)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' vprt.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,vprt)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' wprt.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,wprt)
!
!-----------------------------------------------------------------------
!
!  Read in scalars
!
!----------------------------------------------------------------------
!
      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' ptprt.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,ptprt)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' pprt.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,pprt)

    ELSE
!
!-----------------------------------------------------------------------
!
!  Read in total values of variables from history dump
!
!----------------------------------------------------------------------
!
      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' u.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,uprt)
      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

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' v.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,vprt)
      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

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' w.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,wprt)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' pt.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,ptprt)
      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

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' p.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,pprt)
      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

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' qvprt.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,qvprt)

    ELSE

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' qv.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,qvprt)
      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

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' qc.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,qc)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' qr.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,qr)

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

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) var2d
        WRITE(lchanl,910) label,' rainc.'
      END IF
      CALL mpisplit3d(var2d,nx,ny,1,rainc)
    END IF

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

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) ((var3d(i,j,2),i=1,nxlg),j=1,nylg)
        WRITE(lchanl,910) label,' prcrate2.'

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) ((var3d(i,j,3),i=1,nxlg),j=1,nylg)
        WRITE(lchanl,910) label,' prcrate3.'

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) ((var3d(i,j,4),i=1,nxlg),j=1,nylg)
        WRITE(lchanl,910) label,' prcrate4.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,4,prcrate)
    END IF

    IF( icein == 1 ) THEN

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' qi.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,qi)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' qs.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,qs)

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' qh.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,qh)

    END IF
  END IF

  IF( tkein == 1 ) THEN

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' tke.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,tke)

  END IF

  IF( trbin == 1 ) THEN

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' kmh.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,kmh)

    IF ( intver >= intver410 ) THEN

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
        WRITE(lchanl,910) label,' kmv.'
      END IF
      CALL mpisplit3d(var3d,nx,ny,nz,kmv)

    END IF ! intver

  END IF

  IF( sfcin == 1 ) THEN

    IF (nstyp1 <= 1) THEN

      IF (intver <= intver410) THEN 

        IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) var2d
          WRITE(lchanl,910) label,' tsfc.'
        END IF
        CALL mpisplit3d(var2d,nx,ny,1,tsoil(:,:,1,0))

        IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) var2d 
          WRITE(lchanl,910) label,' tsoil.'
        END IF
        CALL mpisplit3d(var2d,nx,ny,1,tsoil(:,:,2,0))

        IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) var2d 
          WRITE(lchanl,910) label,' wetsfc.'
        END IF
        CALL mpisplit3d(var2d,nx,ny,1,qsoil(:,:,1,0))

        IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) var2d 
          WRITE(lchanl,910) label,' wetdp.'
        END IF
        CALL mpisplit3d(var2d,nx,ny,1,qsoil(:,:,2,0))

      ELSE IF (intver >= intver500) THEN 

        IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg)
          WRITE(lchanl,910) label,' tsoil.'
        END IF
        CALL mpisplit3d(var3d,nx,ny,nzsoil,tsoil(:,:,:,0))

        IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg)
          WRITE(lchanl,910) label,' qsoil.'
        END IF
        CALL mpisplit3d(var3d,nx,ny,nzsoil,qsoil(:,:,:,0))

      END IF  ! intver

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) var2d
        WRITE(lchanl,910) label,' wetcanp.'
      END IF
      CALL mpisplit3d(var2d,nx,ny,1,wetcanp(:,:,0))

    ELSE

      DO is=0,nstyp1

        IF (is <= nstyps) THEN

          IF (intver <= intver410) THEN 

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

            IF(myproc == 0) THEN
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) var2d
            WRITE(lchanl,910) label,' wetdp.'
            END IF
            CALL mpisplit3d(var2d,nx,ny,1,qsoil(:,:,2,is))

          ELSE IF (intver >= intver500) THEN 

            IF(myproc == 0) THEN
              READ(inch,ERR=110,END=120) label
              READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg)
              WRITE(lchanl,910) label,' tsoil.'
            END IF
            CALL mpisplit3d(var3d,nx,ny,nzsoil,tsoil(:,:,:,is))
  
            IF(myproc == 0) THEN
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzsoillg)
            WRITE(lchanl,910) label,' qsoil.'
            END IF
            CALL mpisplit3d(var3d,nx,ny,nzsoil,qsoil(:,:,:,is))

          END IF  ! intver

          IF(myproc == 0) THEN
          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) var2d
          WRITE(lchanl,910) label,' wetcanp.'
          END IF
          CALL mpisplit3d(var2d,nx,ny,1,wetcanp(:,:,is))

        ELSE

          IF(myproc == 0) THEN
            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.'
          END IF

        ENDIF

      END DO

    END IF

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

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

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

  END IF

  IF( radin == 1 ) THEN

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) (((var3d(i,j,k),i=1,nxlg),j=1,nylg),k=1,nzlg)
      WRITE(lchanl,910) label,' radfrc.'
    END IF
    CALL mpisplit3d(var3d,nx,ny,nz,radfrc)

    IF(myproc == 0) THEN
    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) var2d
    WRITE(lchanl,910) label,' radsw.'
    END IF
    CALL mpisplit3d(var2d,nx,ny,1,radsw)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) var2d
      WRITE(lchanl,910) label,' rnflx.'
    END IF
    CALL mpisplit3d(var2d,nx,ny,1,rnflx)

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

    ELSE IF (intver >= intver500) THEN 

      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) var2d
        WRITE(lchanl,910) label,' radswnet.'
      END IF
      CALL mpisplit3d(var2d,nx,ny,1,radswnet)
  
      IF(myproc == 0) THEN
        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) var2d
        WRITE(lchanl,910) label,' radlwin.'
      END IF
      CALL mpisplit3d(var2d,nx,ny,1,radlwin)

    END IF  ! intver 

  END IF

  IF( flxin == 1 ) THEN

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) var2d
      WRITE(lchanl,910) label,' usflx.'
    END IF
    CALL mpisplit3d(var2d,nx,ny,1,usflx)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) var2d
      WRITE(lchanl,910) label,' vsflx.'
    END IF
    CALL mpisplit3d(var2d,nx,ny,1,vsflx)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) var2d
      WRITE(lchanl,910) label,' ptsflx.'
    END IF
    CALL mpisplit3d(var2d,nx,ny,1,ptsflx)

    IF(myproc == 0) THEN
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) var2d
      WRITE(lchanl,910) label,' qvsflx.'
    END IF
    CALL mpisplit3d(var2d,nx,ny,1,qvsflx)

  END IF

  DEALLOCATE(var2d, var3d, var3di)
  
  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 BINREADSPLIT'
  ireturn=1
  RETURN
!
!-----------------------------------------------------------------------
!
!  End-of-file during read
!
!----------------------------------------------------------------------
!

  120   CONTINUE
  WRITE(6,'(/a/)') ' End of file reached in BINREADSPLIT'
  ireturn=2
  RETURN
END SUBROUTINE binreadsplit
!
!##################################################################
!##################################################################
!######                                                      ######
!######                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 BINJOINDUMP                ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE binjoindump(nx,ny,nz,nzsoil,nstyps, nchanl, grdbas,         & 1,114
           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: Yunheng Wang
!  8/16/02.
!  Based on subroutine bindump.
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  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

  INTEGER :: nxlg, nylg, nzlg, nzsoillg
  INTEGER :: n3d, istat

  REAL, ALLOCATABLE :: out1d(:)
  REAL, ALLOCATABLE :: out3d(:,:,:)
  INTEGER, ALLOCATABLE :: out3di(:,:,:)
  REAL, ALLOCATABLE :: outtsoil(:,:,:,:), outqsoil(:,:,:,:)

!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'           ! Grid & map parameters.
  INCLUDE 'mp.inc'             ! mpi parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  nxlg = nproc_x*(nx-3)+3
  nylg = nproc_y*(ny-3)+3 
  nzlg = nz
  nzsoillg = nzsoil
  n3d = MAX(nzlg, nzsoillg, nstyps+1, 4)   ! 3rd dimenson for out3d

  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 binjoindumpp.',1)
  END IF 

  IF (myproc == 0) THEN

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

     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) nxlg,nylg,nzlg
    ELSE IF (intver == intver500) THEN 
      WRITE(nchanl) nxlg,nylg,nzlg,nzsoillg
    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

  END IF  ! myproc == 0

  ALLOCATE (out1d( MAX(nxlg,nylg,nzlg) ),stat=istat)
  IF (istat /= 0)   &
    CALL arpsstop("BINJOINDUMP: ERROR allocating out1d, returning",1)

  ALLOCATE (out3d( nxlg,nylg, n3d ),stat=istat)
  IF (istat /= 0) &
    CALL arpsstop("BINJOINDUMP: ERROR allocating out3d, returning",1)
  
  ALLOCATE (out3di( nxlg,nylg, nstyps ),stat=istat)
  IF (istat /= 0) &
    CALL arpsstop("BINJOINDUMP: ERROR allocating out3di, returning",1)
  
  ALLOCATE (outtsoil( nxlg,nylg, nzsoillg, 0:nstyps ),stat=istat)
  IF (istat /= 0) &
    CALL arpsstop("BINJOINDUMP: ERROR allocating outtsoil, returning",1)
  
  ALLOCATE (outqsoil( nxlg,nylg, nzsoillg, 0:nstyps ),stat=istat)
  IF (istat /= 0) &
    CALL arpsstop("BINJOINDUMP: ERROR allocating outqsoil, returning",1)

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

    CALL mpimerge1dx(x,nx,out1d)
    IF (myproc == 0) THEN
      WRITE(nchanl) 'x coord r1d1'
      WRITE(nchanl) out1d(1:nxlg)
    END IF 

    CALL mpimerge1dy(y,ny,out1d)
    IF (myproc == 0) THEN
      WRITE(nchanl) 'y coord r1d2'
      WRITE(nchanl) out1d(1:nylg)
    END IF 

    IF (myproc == 0) THEN
      WRITE(nchanl) 'z coord r1d3'
      WRITE(nchanl) z
    END IF 

    CALL mpimerge3d(zp,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg)
      WRITE(nchanl) 'zp coor r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    IF (intver == intver500) THEN 
      CALL mpimerge3d(zpsoil,nx,ny,nzsoil,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzsoillg,1,nzsoillg)
        WRITE(nchanl) 'zpsoil rs3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzsoillg)
      END IF 
    END IF 

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

    CALL mpimerge3d(ubar,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'ubar    r3d1'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    CALL mpimerge3d(vbar,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'vbar    r3d2'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    IF (myproc == 0) THEN
      out3d(:,:,:) = 0.0
      WRITE(nchanl) 'wbar    r3d3'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    CALL mpimerge3d(ptbar,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'ptbar   r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    CALL mpimerge3d(pbar,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'pbar    r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    IF(mstout == 1) THEN

      CALL mpimerge3d(qvbar,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'qvbar   r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

    END IF

    IF(landout == 1) THEN

      IF( nstyp <= 1 ) THEN

        CALL mpimerge3di(soiltyp(:,:,1),nx,ny,1,out3di)
        IF (myproc == 0) THEN
          CALL iedgfill(out3di(1,1,1),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,          &
                      1,1,1,1)
          WRITE(nchanl) 'soiltyp i2d0'
          WRITE(nchanl) out3di(1:nxlg,1:nylg,1)
        END IF 

      ELSE
        CALL mpimerge3di(soiltyp,nx,ny,nstyp,out3di)
        CALL mpimerge3d(stypfrct,nx,ny,nstyp,out3d)

        IF (myproc == 0) THEN
        DO is=1,nstyp
          CALL iedgfill(out3di(1,1,is),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,       &
                        1,1,1,1)
          WRITE(nchanl) 'soiltyp i2d0'
          WRITE(nchanl) ((out3di(i,j,is),i=1,nxlg),j=1,nylg)

          CALL edgfill(out3d(1,1,is),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,       &
                       1,1,1,1)
          WRITE(nchanl) 'stypfrc r2d0'
          WRITE(nchanl) ((out3d(i,j,is),i=1,nxlg),j=1,nylg)

        END DO
        END IF 

      END IF

      CALL mpimerge3di(vegtyp,nx,ny,1,out3di)
      IF (myproc == 0) THEN
        CALL iedgfill(out3di ,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'vegtyp  i2d0'
        WRITE(nchanl) out3di(1:nxlg,1:nylg,1)
      END IF 

      CALL mpimerge3d(lai,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'lai     r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

      CALL mpimerge3d(roufns,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d ,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'roufns  r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

      CALL mpimerge3d(veg,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d    ,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'veg     r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

    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 mpimerge3d(tem1,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'uprt    r3d1'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      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 mpimerge3d(tem1,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'vprt    r3d2'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      CALL mpimerge3d(w,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg)
        WRITE(nchanl) 'wprt    r3d3'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 
!
!-----------------------------------------------------------------------
!
!  Write out scalars
!
!-----------------------------------------------------------------------
!
      CALL mpimerge3d(ptprt,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'ptprt   r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      CALL mpimerge3d(pprt,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'pprt    r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

    ELSE
!
!-----------------------------------------------------------------------
!
!  Write out total values to history dump
!
!-----------------------------------------------------------------------
!
      CALL mpimerge3d(u,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'u       r3d1'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      CALL mpimerge3d(v,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'v       r3d2'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      CALL mpimerge3d(w,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg)
        WRITE(nchanl) 'w       r3d3'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 
!
!-----------------------------------------------------------------------
!
!  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 mpimerge3d(tem1,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'pt      r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      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 mpimerge3d(tem1,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'p       r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

    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 mpimerge3d(tem1,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'qvprt   r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

    ELSE
!
!-----------------------------------------------------------------------
!
!  Write out total values to history dump
!
!-----------------------------------------------------------------------
!
      CALL mpimerge3d(qv,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'qv      r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

    END IF

    CALL mpimerge3d(qc,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'qc      r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    CALL mpimerge3d(qr,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'qr      r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    IF(rainout == 1) THEN

      CALL mpimerge3d(raing,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,   1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'raing   r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

      CALL mpimerge3d(rainc,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'rainc   r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

    END IF   !rainout

    IF ( prcout == 1 ) THEN
      CALL mpimerge3d(prcrate,nx,ny,4,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,4,1,4)
        WRITE(nchanl) 'prcrat1 r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
        WRITE(nchanl) 'prcrat2 r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,2)
        WRITE(nchanl) 'prcrat3 r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,3)
        WRITE(nchanl) 'prcrat4 r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,4)
      END IF 
    END IF

    IF(iceout == 1) THEN

      CALL mpimerge3d(qi,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'qi      r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      CALL mpimerge3d(qs,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'qs      r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

      CALL mpimerge3d(qh,nx,ny,nz,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
        WRITE(nchanl) 'qh      r3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
      END IF 

    END IF   !iceout

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

    CALL mpimerge3d(tke,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'tke     r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

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

    CALL mpimerge3d(kmh,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'kmh     r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    CALL mpimerge3d(kmv,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'kmv     r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

  END IF   ! trbout

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

    IF( nstyp <= 1 ) THEN

      CALL mpimerge3d(tsoil(:,:,:,0),nx,ny,nzsoil,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,  1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,             &
                   1,nzsoillg,1,nzsoillg)
        WRITE(nchanl) 'tsoil  rs3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzsoillg)
      END IF 

      CALL mpimerge3d(qsoil(:,:,:,0),nx,ny,nzsoil,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d, 1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,             &
                   1,nzsoillg,1,nzsoillg)
        WRITE(nchanl) 'qsoil  rs3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzsoillg)
      END IF 
      

      CALL mpimerge3d(wetcanp(:,:,0),nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,             &
                   1,1,1,1)
        WRITE(nchanl) 'wetcanp r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

    ELSE

      CALL mpimerge4d(tsoil,nx,ny,nzsoil,nstyp+1,outtsoil)
      CALL mpimerge4d(qsoil,nx,ny,nzsoil,nstyp+1,outqsoil)
      CALL mpimerge3d(wetcanp,nx,ny,nstyp+1,out3d)

      IF (myproc == 0) THEN

        DO is=0,nstyp
          CALL edgfill(outtsoil(1,1,1,is),  1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,          &
                     1,nzsoillg,1,nzsoillg)
          CALL edgfill(outqsoil(1,1,1,is), 1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,          &
                     1,nzsoillg,1,nzsoillg)

          WRITE(nchanl) 'tsoil  rs3d0'
          WRITE(nchanl) outtsoil(1:nxlg,1:nylg,1:nzsoillg,is)
          WRITE(nchanl) 'qsoil  rs3d0'
          WRITE(nchanl) outqsoil(1:nxlg,1:nylg,1:nzsoillg,is)

          CALL edgfill(out3d(1,1,is+1),1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,          &
                     1,1,1,1)
          WRITE(nchanl) 'wetcanp r2d0'
          WRITE(nchanl) out3d(1:nxlg,1:nylg,is+1)
        END DO

      END IF

    END IF 

    IF (snowout == 1) THEN

      CALL mpimerge3d(snowdpth,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1,                   &
                    1,1,1,1)
        WRITE(nchanl) 'snowdpthr2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

    END IF

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

    CALL mpimerge3d(radfrc,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,nzlg,1,nzlg-1)
      WRITE(nchanl) 'radfrc  r3d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzlg)
    END IF 

    CALL mpimerge3d(radsw,nx,ny,1,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
      WRITE(nchanl) 'radsw   r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 

    CALL mpimerge3d(rnflx,nx,ny,1,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
      WRITE(nchanl) 'rnflx   r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 

    IF (intver >= intver500) THEN 

      CALL mpimerge3d(radswnet,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'radswnetr2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 
  
      CALL mpimerge3d(radlwin,nx,ny,1,out3d)
      IF (myproc == 0) THEN
        CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
        WRITE(nchanl) 'radlwin r2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 

    END IF  ! intver 

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

    CALL mpimerge3d(usflx,nx,ny,1,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg, 1,nylg,1,nylg-1, 1,1,1,1)
      WRITE(nchanl) 'usflx   r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 

    CALL mpimerge3d(vsflx,nx,ny,1,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg, 1,1,1,1)
      WRITE(nchanl) 'vsflx   r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 
  
    CALL mpimerge3d(ptsflx,nx,ny,1,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
      WRITE(nchanl) 'ptsflx  r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 
  
    CALL mpimerge3d(qvsflx,nx,ny,1,out3d)
    IF (myproc == 0) THEN
      CALL edgfill(out3d,1,nxlg,1,nxlg-1, 1,nylg,1,nylg-1, 1,1,1,1)
      WRITE(nchanl) 'qvsflx  r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 

  END IF   ! flxout

  DEALLOCATE(out1d)
  DEALLOCATE(out3d, out3di)
  DEALLOCATE(outtsoil, outqsoil)

  RETURN
END SUBROUTINE binjoindump
!
!##################################################################
!##################################################################
!######                                                      ######
!######                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