!
!##################################################################
!##################################################################
!######                                                      ######
!######                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,9
           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,nzsoilin
  INTEGER :: bgrdin,bbasin,bvarin,bicein,btkein,btrbin
  INTEGER :: idummy

  REAL(4) :: rdummy
  REAL(4) :: umove4,   vmove4
  REAL(4) :: xgrdorg4, ygrdorg4
  REAL(4) :: trulat14, trulat24, trulon4
  REAL(4) :: sclfct4
  REAL(4) :: tstop4,   thisdmp4
  REAL(4) :: latitud4, ctrlat4,  ctrlon4

  REAL(4), ALLOCATABLE :: invar1d(:)
  REAL(4), ALLOCATABLE :: invar2d(:,:)
  REAL(4), ALLOCATABLE :: invar3d(:,:,:)
  REAL(4), ALLOCATABLE :: invar3dsoil(:,:,:)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'indtflg.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'mp.inc'            ! mpi parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  ALLOCATE(invar1d(MAX(nx,ny,nz)), STAT = idummy)
  CALL check_alloc_status(idummy, "BINREAD:invar1d")
  ALLOCATE(invar2d(nx,ny),         STAT = idummy)
  CALL check_alloc_status(idummy, "BINREAD:invar2d")
  ALLOCATE(invar3d(nx,ny,nz),      STAT = idummy)
  CALL check_alloc_status(idummy, "BINREAD:invar3d")
  ALLOCATE(invar3dsoil(nx,ny,nzsoil), STAT = idummy)
  CALL check_alloc_status(idummy, "BINREAD:invar3dsoil")
!
!-----------------------------------------------------------------------
!
!  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) rdummy,tmunit
  time = rdummy
!
!-----------------------------------------------------------------------
!
!  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)                                            &
                  umove4,  vmove4,  xgrdorg4,ygrdorg4,trulat14,         &
                  trulat24,trulon4, sclfct4, rdummy,  rdummy,           &
                  rdummy,  rdummy,  rdummy,  rdummy,  rdummy,           &
                  tstop4,  thisdmp4,latitud4,ctrlat4, ctrlon4

  umove = umove4
  vmove = vmove4
  xgrdorg = xgrdorg4
  ygrdorg = ygrdorg4
  trulat1 = trulat14
  trulat2 = trulat24
  trulon  = trulon4
  sclfct  = sclfct4
  tstop   = tstop4
  thisdmp = thisdmp4
  latitud = latitud4
  ctrlat  = ctrlat4
  ctrlon  = ctrlon4

  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) invar1d(1:nx)
    x(:) = invar1d(1:nx)
    IF (myproc == 0) WRITE(lchanl,910) label,' x.'

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

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

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) invar3d
    zp(:,:,:) = invar3d(:,:,:)
    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) invar3dsoil
      zpsoil(:,:,:) = invar3dsoil(:,:,:)
      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) invar3d
    ubar(:,:,:) = invar3d(:,:,:)
    IF (myproc == 0) WRITE(lchanl,910) label,' ubar.'

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

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

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

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) invar3d
    pbar(:,:,:) = invar3d(:,:,:)
    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) invar3d
      qvbar(:,:,:) = invar3d(:,:,:)
      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) invar2d
            stypfrct(:,:,is) = invar2d(:,:)
            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) invar2d
      lai(:,:) = invar2d(:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' lai.'

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

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar2d
      veg(:,:) = invar2d(:,:)
      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) invar3d
      uprt(:,:,:) = invar3d(:,:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' uprt.'

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

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

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar3d
      pprt(:,:,:) = invar3d(:,:,:)
      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) invar3d
      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) = invar3d(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) invar3d
      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) = invar3d(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) invar3d
      wprt(:,:,:) = invar3d(:,:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' w.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar3d
      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) = invar3d(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) invar3d
      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) = invar3d(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

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) invar3d

    IF ( totin == 0 ) THEN
      IF (myproc == 0) WRITE(lchanl,910) label,' qvprt.'
      qvprt(:,:,:) = invar3d(:,:,:)
    ELSE
      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) = invar3d(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) invar3d
    qc(:,:,:) = invar3d(:,:,:)
    IF (myproc == 0) WRITE(lchanl,910) label,' qc.'

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) invar3d
    qr(:,:,:) = invar3d(:,:,:)
    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) invar2d
      raing(:,:) = invar2d(:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' raing.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar2d
      rainc(:,:) = invar2d(:,:)
      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) invar2d
      prcrate(:,:,1) = invar2d(:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' prcrate1.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar2d
      prcrate(:,:,2) = invar2d(:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' prcrate2.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar2d
      prcrate(:,:,3) = invar2d(:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' prcrate3.'

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar2d
      prcrate(:,:,4) = invar2d(:,:)
      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) invar3d
      qi(:,:,:) = invar3d(:,:,:)
      IF (myproc == 0) WRITE(lchanl,910) label,' qi.'

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

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar3d
      qh(:,:,:) = invar3d(:,:,:)
      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) invar3d
    tke(:,:,:) = invar3d(:,:,:)
    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) invar3d
    kmh(:,:,:) = invar3d(:,:,:)
    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) invar3d
      kmv(:,:,:) = invar3d(:,:,:)
      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) invar2d
        tsoil(:,:,1,0) = invar2d(:,:)
        IF (myproc == 0) WRITE(lchanl,910) label,' tsfc.'

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

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) invar2d
        qsoil(:,:,1,0) = invar2d(:,:)
        IF (myproc == 0) WRITE(lchanl,910) label,' wetsfc.'

        READ(inch,ERR=110,END=120) label
        READ(inch,ERR=110,END=120) invar2d
        qsoil(:,:,2,0) = invar2d(:,:)
        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) invar3dsoil
        tsoil(:,:,:,0) = invar3dsoil(:,:,:)
        IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.'

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

      END IF  ! intver

      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar2d
      wetcanp(:,:,0) = invar2d(:,:)
      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) invar2d
            tsoil(:,:,1,is) = invar2d(:,:)
            IF (myproc == 0) WRITE(lchanl,910) label,' tsfc.'
  
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) invar2d
            tsoil(:,:,2,is) = invar2d(:,:)
            IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.'
  
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) invar2d
            qsoil(:,:,1,is) = invar2d(:,:)
            IF (myproc == 0) WRITE(lchanl,910) label,' wetsfc.'

            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) invar2d
            tsoil(:,:,2,is) = invar2d(:,:)
            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) invar3dsoil
            tsoil(:,:,:,is) = invar3dsoil(:,:,:)
            IF (myproc == 0) WRITE(lchanl,910) label,' tsoil.'
  
            READ(inch,ERR=110,END=120) label
            READ(inch,ERR=110,END=120) invar3dsoil
            qsoil(:,:,:,is) = invar3dsoil(:,:,:)
            IF (myproc == 0) WRITE(lchanl,910) label,' qsoil.'

          END IF  ! intver

          READ(inch,ERR=110,END=120) label
          READ(inch,ERR=110,END=120) invar2d
          wetcanp(:,:,is) = invar2d(:,:)
          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) invar2d
      snowdpth(:,:) = invar2d(:,:)
      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) invar3d
    radfrc(:,:,:) = invar3d(:,:,:)
    IF (myproc == 0) WRITE(lchanl,910) label,' radfrc.'

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

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) invar2d
    rnflx(:,:) = invar2d(:,:)
    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) invar2d
      radswnet(:,:) = invar2d(:,:)
      WRITE(lchanl,910) label,' radswnet.'
  
      READ(inch,ERR=110,END=120) label
      READ(inch,ERR=110,END=120) invar2d
      radlwin(:,:) = invar2d(:,:)
      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) invar2d
    usflx(:,:) = invar2d(:,:)
    IF (myproc == 0) WRITE(lchanl,910) label,' usflx.'

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

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

    READ(inch,ERR=110,END=120) label
    READ(inch,ERR=110,END=120) invar2d
    qvsflx(:,:) = invar2d(:,:)
    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.'

  DEALLOCATE(invar1d, STAT = idummy)
  DEALLOCATE(invar2d, STAT = idummy)
  DEALLOCATE(invar3d, STAT = idummy)
  DEALLOCATE(invar3dsoil, STAT = idummy)

  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,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 mpisplit2di(var3di(:,:,1),nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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 mpisplit2d(var2d,nx,ny,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

  DEALLOCATE(var2d, var3d, var3di)

  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(4) :: rdummy
  REAL(4) :: umove4,   vmove4,   xgrdorg4, ygrdorg4 
  REAL(4) :: trulat14, trulat24, trulon4,  sclfct4  
  REAL(4) :: tstop4,   thisdmp4, latitud4
  REAL(4) :: ctrlat4,  ctrlon4  

  REAL(4), ALLOCATABLE :: out1d(:)
  REAL(4), ALLOCATABLE :: out2d(:,:)
  REAL(4), ALLOCATABLE :: out3d(:,:,:)
  REAL(4), ALLOCATABLE :: out3dsoil(:,:,:)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'           ! Grid & map parameters.
  INCLUDE 'mp.inc'             ! mpi parameters.
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!---------------------------------------------------------------
!
! Allocate 4 byte arrays
!
!---------------------------------------------------------------

  ALLOCATE(out1d(MAX(nx,ny,nz)), STAT = idummy)
  ALLOCATE(out2d(nx,ny),         STAT = idummy)
  ALLOCATE(out3d(nx,ny,nz),      STAT = idummy)
  ALLOCATE(out3dsoil(nx,ny,nzsoil), STAT = idummy)


  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

  rdummy = curtim
  WRITE(nchanl) rdummy,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
  umove4 = umove
  vmove4 = vmove
  xgrdorg4 = xgrdorg
  ygrdorg4 = ygrdorg
  trulat14 = trulat1
  trulat24 = trulat2
  trulon4  = trulon
  sclfct4  = sclfct
  tstop4   = tstop
  thisdmp4 = thisdmp
  latitud4 = latitud
  ctrlat4  = ctrlat
  ctrlon4  = ctrlon
  WRITE(nchanl)   umove4,   vmove4, xgrdorg4, ygrdorg4, trulat14,       &
                trulat24,  trulon4,  sclfct4,   rdummy,   rdummy,       &
                  rdummy,   rdummy,   rdummy,   rdummy,   rdummy,       &
                  tstop4, thisdmp4,  latitud4, ctrlat4,  ctrlon4
!
!-----------------------------------------------------------------------
!
!  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

    out1d(1:nx) = x
    WRITE(nchanl) 'x coord r1d1'
    WRITE(nchanl) out1d(1:nx)

    out1d(1:ny) = y
    WRITE(nchanl) 'y coord r1d2'
    WRITE(nchanl) out1d(1:ny)

    out1d(1:nz) = z
    WRITE(nchanl) 'z coord r1d3'
    WRITE(nchanl) out1d(1:nz)

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

    IF (intver == intver500) THEN 
      CALL edgfill(zpsoil,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nzsoil,1,nzsoil)
      out3dsoil(:,:,:) = zpsoil(:,:,:)
      WRITE(nchanl) 'zpsoil rs3d0'
      WRITE(nchanl) out3dsoil
    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)
    out3d(:,:,:) = ubar(:,:,:)
    WRITE(nchanl) 'ubar    r3d1'
    WRITE(nchanl) out3d

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

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

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

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

    IF(mstout == 1) THEN

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

    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)
          out2d(:,:) = stypfrct(:,:,is)
          WRITE(nchanl) 'stypfrc r2d0'
          WRITE(nchanl)  out2d

        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'
      out2d(:,:) = lai(:,:)
      WRITE(nchanl)  out2d

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

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

    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)
      out3d(:,:,:) = tem1(:,:,:)
      WRITE(nchanl) 'uprt    r3d1'
      WRITE(nchanl) out3d

      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)
      out3d(:,:,:) = tem1(:,:,:)
      WRITE(nchanl) 'vprt    r3d2'
      WRITE(nchanl) out3d

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

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

    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)
      out3d(:,:,:) = u(:,:,:)
      WRITE(nchanl) 'u       r3d1'
      WRITE(nchanl) out3d

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

      CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
      out3d(:,:,:) = w(:,:,:)
      WRITE(nchanl) 'w       r3d3'
      WRITE(nchanl) out3d
!
!-----------------------------------------------------------------------
!
!  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)
      out3d(:,:,:) = tem1(:,:,:)
      WRITE(nchanl) 'pt      r3d0'
      WRITE(nchanl) out3d

      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)
      out3d(:,:,:) = tem1(:,:,:)
      WRITE(nchanl) 'p       r3d0'
      WRITE(nchanl) out3d

    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)
      out3d(:,:,:) = tem1(:,:,:)
      WRITE(nchanl) 'qvprt   r3d0'
      WRITE(nchanl)  out3d

    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)
      out3d(:,:,:) = qv(:,:,:)
      WRITE(nchanl) 'qv      r3d0'
      WRITE(nchanl) out3d

    END IF

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

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

    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'
      out2d(:,:) = raing(:,:)
      WRITE(nchanl)  out2d

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

    END IF   !rainout

    IF ( prcout == 1 ) THEN
      CALL edgfill(prcrate,1,nx,1,nx-1, 1,ny,1,ny-1, 1,4,1,4)

      out2d(:,:) = prcrate(:,:,1)
      WRITE(nchanl) 'prcrat1 r2d0'
      WRITE(nchanl)  out2d

      out2d(:,:) = prcrate(:,:,2)
      WRITE(nchanl) 'prcrat2 r2d0'
      WRITE(nchanl)  out2d

      out2d(:,:) = prcrate(:,:,3)
      WRITE(nchanl) 'prcrat3 r2d0'
      WRITE(nchanl)  out2d

      out2d(:,:) = prcrate(:,:,4)
      WRITE(nchanl) 'prcrat4 r2d0'
      WRITE(nchanl)  out2d
    END IF

    IF(iceout == 1) THEN

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

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

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

    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)
    out3d(:,:,:) = tke(:,:,:)
    WRITE(nchanl) 'tke     r3d0'
    WRITE(nchanl)  out3d

  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)
    out3d(:,:,:) = kmh(:,:,:)
    WRITE(nchanl) 'kmh     r3d0'
    WRITE(nchanl)  out3d

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

  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. 
!
        out3d(:,:,1)=0.

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

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

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

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

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

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

        out3dsoil(:,:,:) = tsoil(:,:,:,0)
        WRITE(nchanl) 'tsoil  rs3d0'
        WRITE(nchanl)  out3dsoil
        out3dsoil(:,:,:) = qsoil(:,:,:,0)
        WRITE(nchanl) 'qsoil  rs3d0'
        WRITE(nchanl)  out3dsoil
      
      END IF  ! intver

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

    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. 
!
        out3d(:,:,1)=0.

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

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

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

        out3d(:,:,1)=0.
        DO k=2,nzsoil
          DO j=1,ny
            DO i=1,nx
              out3d(i,j,1)=out3d(i,j,1)+qsoil(i,j,k,is)
            END DO 
          END DO 
        END DO 

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

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

        out3dsoil(:,:,:) = tsoil(:,:,:,is)
        WRITE(nchanl) 'tsoil  rs3d0'
        WRITE(nchanl)  out3dsoil
        out3dsoil(:,:,:) = qsoil(:,:,:,is)
        WRITE(nchanl) 'qsoil  rs3d0'
        WRITE(nchanl)  out3dsoil

      END IF  ! intver

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

      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)
      out2d(:,:) = snowdpth(:,:)
      WRITE(nchanl) 'snowdpthr2d0'
      WRITE(nchanl)  out2d
    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)
    out3d(:,:,:) = radfrc(:,:,:)
    WRITE(nchanl) 'radfrc  r3d0'
    WRITE(nchanl) out3d

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

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

    IF (intver >= intver500) THEN 

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

    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)
    out2d(:,:) = usflx(:,:)
    WRITE(nchanl) 'usflx   r2d0'
    WRITE(nchanl) out2d

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

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

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

  END IF   ! flxout
!---------------------------------------------------------------
!
! Deallocate 4 byte working arrays
!
!---------------------------------------------------------------

  DEALLOCATE(out1d)
  DEALLOCATE(out2d)
  DEALLOCATE(out3d)
  DEALLOCATE(out3dsoil)

  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), PARAMETER :: tmunit = 'seconds   '
!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,l,is
  INTEGER :: idummy
  REAL    :: rdummy

  INTEGER :: nxlg, nylg
  INTEGER :: n3d, istat

  REAL,    ALLOCATABLE :: out1d(:)
  REAL,    ALLOCATABLE :: out3d(:,:,:)
  INTEGER, ALLOCATABLE :: out3di(:,:,:)
  REAL,    ALLOCATABLE :: outtsoil(:,:,:,:)
  REAL,    ALLOCATABLE :: 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 

  n3d = MAX(nz, nzsoil, 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,nz
    ELSE IF (intver == intver500) THEN 
      WRITE(nchanl) nxlg,nylg,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

  END IF  ! myproc == 0

  ALLOCATE (out1d( MAX(nxlg,nylg) ),stat=istat)
  CALL check_alloc_status(istat, "binjoindump:out1d")

  ALLOCATE (out3d( nxlg,nylg, n3d ),stat=istat)
  CALL check_alloc_status(istat, "binjoindump:out3d")
  
  ALLOCATE (out3di( nxlg,nylg, nstyps ),stat=istat)
  CALL check_alloc_status(istat, "binjoindump:out3di")
  
  ALLOCATE (outtsoil( nxlg,nylg, nzsoil, 0:nstyps ),stat=istat)
  CALL check_alloc_status(istat, "binjoindump:outtsoil")
  
  ALLOCATE (outqsoil( nxlg,nylg, nzsoil, 0:nstyps ),stat=istat)
  CALL check_alloc_status(istat, "binjoindump:outqsoil")

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

    IF (intver == intver500) THEN 
      CALL edgfill(zpsoil,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nzsoil,1,nzsoil)
      CALL mpimerge3d(zpsoil,nx,ny,nzsoil,out3d)
      IF (myproc == 0) THEN
        WRITE(nchanl) 'zpsoil rs3d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nzsoil)
      END IF 
    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)
    CALL mpimerge3d(ubar,nx,ny,nz,out3d)
    IF (myproc == 0) THEN
      WRITE(nchanl) 'ubar    r3d1'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1:nz)
    END IF 

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

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

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

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

    IF(mstout == 1) THEN

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

    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)

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

      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)

          CALL edgfill(stypfrct(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,       &
                       1,1,1,1)
        END DO

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

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

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

          END DO
        END IF 

      END IF

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

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

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

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

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

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

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

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

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

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

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

    END IF

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

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

    IF(rainout == 1) THEN

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

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

    END IF   !rainout

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

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

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

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

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

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

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

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

    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)
         CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,          &
                      1,1,1,1)

      END DO

      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

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

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

      END IF

    END IF 

    IF (snowout == 1) THEN

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

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

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

    IF (intver >= intver500) THEN 

      CALL edgfill(radswnet,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL mpimerge2d(radswnet,nx,ny,out3d)
      IF (myproc == 0) THEN
        WRITE(nchanl) 'radswnetr2d0'
        WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
      END IF 
  
      CALL edgfill(radlwin,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL mpimerge2d(radlwin,nx,ny,out3d)
      IF (myproc == 0) THEN
        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 edgfill(usflx,1,nx,1,nx, 1,ny,1,ny-1, 1,1,1,1)
    CALL mpimerge2d(usflx,nx,ny,out3d)
    IF (myproc == 0) THEN
      WRITE(nchanl) 'usflx   r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 

    CALL edgfill(vsflx,1,nx,1,nx-1, 1,ny,1,ny, 1,1,1,1)
    CALL mpimerge2d(vsflx,nx,ny,out3d)
    IF (myproc == 0) THEN
      WRITE(nchanl) 'vsflx   r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 
  
    CALL edgfill(ptsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL mpimerge2d(ptsflx,nx,ny,out3d)
    IF (myproc == 0) THEN
      WRITE(nchanl) 'ptsflx  r2d0'
      WRITE(nchanl) out3d(1:nxlg,1:nylg,1)
    END IF 
  
    CALL edgfill(qvsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL mpimerge2d(qvsflx,nx,ny,out3d)
    IF (myproc == 0) THEN
      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