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


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

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in history data in the NCSA HDF4 format.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/15
!
!  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
!
!    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.
!    filename  Character variable nhming the input HDF file
!
!  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)
!
!    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
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at ground (K) (in top 1 cm layer)
!    tsoil     Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture in the top 1 cm layer
!    wetdp    Deep soil moisture in the deep 1 m layer
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!    ireturn  Return status indicator
!
!    WORK ARRAY
!
!    tem1
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

  INTEGER :: nx,ny,nz

  INTEGER :: grdbas
  CHARACTER (LEN=*) :: filename

  REAL :: x     (nx)           ! x coord.
  REAL :: y     (ny)           ! y coord.
  REAL :: z     (nz)           ! z coord.
  REAL :: zp    (nx,ny,nz)     ! physical x coord.

  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)    ! Fraction of soil types
  INTEGER :: vegtyp(nx,ny)     ! Vegetation type
  REAL :: lai    (nx,ny)    ! Leaf Area Index
  REAL :: roufns (nx,ny)    ! Surface roughness
  REAL :: veg    (nx,ny)    ! Vegetation fraction

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

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

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

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

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

  INTEGER (KIND=selected_int_kind(4)), ALLOCATABLE :: itmp(:,:,:) ! Temporary array
  REAL, ALLOCATABLE :: hmax(:), hmin(:) ! Temporary array

  INTEGER :: ireturn
  REAL :: time

!-----------------------------------------------------------------------
!
!  Parameters describing routine that wrote the gridded data
!
!-----------------------------------------------------------------------

  CHARACTER (LEN=40) :: fmtver,fmtverin
  PARAMETER (fmtver='004.10 HDF4 Coded Data')
  CHARACTER (LEN=10) :: tmunit

!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------

  INTEGER :: lchanl
  PARAMETER (lchanl=6)      ! Channel number for formatted printing.

  INTEGER :: isizes(3), isize2
  INTEGER :: i,j,k,iret,ndim,itime,is
  INTEGER :: nxin,nyin,nzin

  INTEGER :: bgrdin,bbasin,bvarin,bicein,btrbin,btkein

  INTEGER :: l

  INTEGER :: istat, sd_id, sds_id

  INTEGER :: nstyp1,nstypin

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'        ! Grid parameters
  INCLUDE 'indtflg.inc'
  INCLUDE 'alloc.inc'       ! allocation parameters & declarations

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

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

  ALLOCATE (itmp(nx,ny,nz),stat=istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFDUMP: ERROR allocating itmp, returning"
    RETURN
  END IF
  ALLOCATE (hmax(nz),stat=istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFDUMP: ERROR allocating hmax, returning"
    RETURN
  END IF
  ALLOCATE (hmin(nz),stat=istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFDUMP: ERROR allocating hmin, returning"
    RETURN
  END IF

!-----------------------------------------------------------------------
!
!  Read header info
!
!-----------------------------------------------------------------------

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

  CALL hdfrdc(sd_id,40,"fmtver",fmtverin,istat)

  IF ( fmtverin /= fmtver ) THEN
    WRITE(6,'(/1x,a/1x,2a/1x,2a)')                                      &
        'Data format incompatible with the data reader.',               &
        'Format of data is ',fmtverin,' Format of reader is ',fmtver
    CALL arpsstop('Job stopped in HDFREAD',1)
  END IF

  CALL hdfrdc(sd_id,40,"runname",runname,istat)
  CALL hdfrdi(sd_id,"nocmnt",nocmnt,istat)
  IF( nocmnt > 0 ) THEN
    CALL hdfrdc(sd_id,80*nocmnt,"cmnt",cmnt,istat)
  END IF

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

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

  CALL hdfrdc(sd_id,10,"tmunit",tmunit,istat)
  CALL hdfrdr(sd_id,"time",time,istat)

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

  CALL hdfrdi(sd_id,"nx",nxin,istat)
  CALL hdfrdi(sd_id,"ny",nyin,istat)
  CALL hdfrdi(sd_id,"nz",nzin,istat)

  IF ( nxin /= nx .OR. nyin /= ny .OR. nzin /= nz ) THEN
    WRITE(6,'(1x,a)') ' Dimensions in HDFREAD 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 HDFREAD.'
    CALL arpsstop('Job stopped in HDFREAD',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.'

    CALL hdfrdi(sd_id,"grdflg",bgrdin,istat)
    CALL hdfrdi(sd_id,"basflg",bbasin,istat)
    CALL hdfrdi(sd_id,"varflg",bvarin,istat)
    CALL hdfrdi(sd_id,"mstflg",mstin,istat)
    CALL hdfrdi(sd_id,"iceflg",bicein,istat)
    CALL hdfrdi(sd_id,"trbflg",btrbin,istat)
    CALL hdfrdi(sd_id,"landflg",landin,istat)
    CALL hdfrdi(sd_id,"totflg",totin,istat)
    CALL hdfrdi(sd_id,"tkeflg",btkein,istat)

  ELSE ! Normal data reading

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

    CALL hdfrdi(sd_id,"grdflg",grdin,istat)
    CALL hdfrdi(sd_id,"basflg",basin,istat)
    CALL hdfrdi(sd_id,"varflg",varin,istat)
    CALL hdfrdi(sd_id,"mstflg",mstin,istat)
    CALL hdfrdi(sd_id,"iceflg",icein,istat)
    CALL hdfrdi(sd_id,"trbflg",trbin,istat)
    CALL hdfrdi(sd_id,"sfcflg",sfcin,istat)
    CALL hdfrdi(sd_id,"rainflg",rainin,istat)
    !Don't read landin (land data only written to base file,
    !landin flag set there).
    !CALL hdfrdi(sd_id,"landflg",landin,istat) 
    CALL hdfrdi(sd_id,"totflg",totin,istat)
    CALL hdfrdi(sd_id,"tkeflg",tkein,istat)

  END IF

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

  IF ( nstyp1 < 1 ) THEN
    nstyp1 = 1
  END IF
  IF (nstyp1 > nstyp) THEN
    WRITE (6,*) "HDFREAD: WARNING, nstyp in file (",nstyp1,  &
       ") greater than that specified in input file (",nstyp,  &
       "), using only",nstyp 
    nstypin = nstyp
  ELSE
    nstypin = nstyp1
  ENDIF

  CALL hdfrdi(sd_id,"prcflg",prcin,istat)
  CALL hdfrdi(sd_id,"radflg",radin,istat)
  CALL hdfrdi(sd_id,"flxflg",flxin,istat)
  CALL hdfrdi(sd_id,"snowflg",snowin,istat)

  CALL hdfrdi(sd_id,"month",month,istat)
  CALL hdfrdi(sd_id,"day",day,istat)
  CALL hdfrdi(sd_id,"year",year,istat)
  CALL hdfrdi(sd_id,"hour",hour,istat)
  CALL hdfrdi(sd_id,"minute",minute,istat)
  CALL hdfrdi(sd_id,"second",second,istat)

  CALL hdfrdr(sd_id,"umove",umove,istat)
  CALL hdfrdr(sd_id,"vmove",vmove,istat)
  CALL hdfrdr(sd_id,"xgrdorg",xgrdorg,istat)
  CALL hdfrdr(sd_id,"ygrdorg",ygrdorg,istat)

  CALL hdfrdi(sd_id,"mapproj",mapproj,istat)
  CALL hdfrdr(sd_id,"trulat1",trulat1,istat)
  CALL hdfrdr(sd_id,"trulat2",trulat2,istat)
  CALL hdfrdr(sd_id,"trulon",trulon,istat)
  CALL hdfrdr(sd_id,"sclfct",sclfct,istat)
  CALL hdfrdr(sd_id,"tstop",tstop,istat)
  CALL hdfrdr(sd_id,"thisdmp",thisdmp,istat)
  CALL hdfrdr(sd_id,"latitud",latitud,istat)
  CALL hdfrdr(sd_id,"ctrlat",ctrlat,istat)
  CALL hdfrdr(sd_id,"ctrlon",ctrlon,istat)

!-----------------------------------------------------------------------
!
!  Read in x,y and z at grid cell centers (scalar points).
!
!-----------------------------------------------------------------------

  IF( grdin == 1 .OR. grdbas == 1 ) THEN

    CALL hdfrd1d(sd_id,"x",nx,x,istat)
    IF (istat /= 0) GO TO 110

    CALL hdfrd1d(sd_id,"y",ny,y,istat)
    IF (istat /= 0) GO TO 110

    CALL hdfrd1d(sd_id,"z",nz,z,istat)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"zp",nx,ny,nz,zp,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

  END IF  ! grdin

!-----------------------------------------------------------------------
!
!  Read in base state fields
!
!-----------------------------------------------------------------------

  IF( basin == 1 .OR. grdbas == 1 ) THEN

    CALL hdfrd3d(sd_id,"ubar",nx,ny,nz,ubar,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"vbar",nx,ny,nz,vbar,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"wbar",nx,ny,nz,wbar,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"ptbar",nx,ny,nz,ptbar,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"pbar",nx,ny,nz,pbar,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    IF( mstin == 1) THEN
      CALL hdfrd3d(sd_id,"qvbar",nx,ny,nz,qvbar,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110
    END IF

    IF (landin == 1) THEN

      CALL hdfrd3di(sd_id,"soiltyp",nx,ny,nstypin,soiltyp(1,1,1),istat)
      IF (istat /= 0) GO TO 110

      CALL hdfrd3d(sd_id,"stypfrct",nx,ny,nstypin,                      &
          stypfrct(1,1,1),istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

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

      CALL hdfrd2di(sd_id,"vegtyp",nx,ny,vegtyp,istat)
      IF (istat /= 0) GO TO 110

      CALL hdfrd2d(sd_id,"lai",nx,ny,lai,istat,itmp)
      IF (istat /= 0) GO TO 110

      CALL hdfrd2d(sd_id,"roufns",nx,ny,roufns,istat,itmp)
      IF (istat /= 0) GO TO 110

      CALL hdfrd2d(sd_id,"veg",nx,ny,veg,istat,itmp)
      IF (istat /= 0) GO TO 110

    END IF

  END IF

  IF( grdbas == 1 ) GO TO 930

  IF( varin == 1 ) THEN

    IF ( totin == 0 ) THEN

!-----------------------------------------------------------------------
!
!  Read in perturbations from history dump
!
!-----------------------------------------------------------------------

      CALL hdfrd3d(sd_id,"uprt",nx,ny,nz,uprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

      CALL hdfrd3d(sd_id,"vprt",nx,ny,nz,vprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

      CALL hdfrd3d(sd_id,"wprt",nx,ny,nz,wprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

!-----------------------------------------------------------------------
!
!  Read in scalars
!
!-----------------------------------------------------------------------

      CALL hdfrd3d(sd_id,"ptprt",nx,ny,nz,ptprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

      CALL hdfrd3d(sd_id,"pprt",nx,ny,nz,pprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

    ELSE

!-----------------------------------------------------------------------
!
!  Read in total values of variables from history dump
!
!-----------------------------------------------------------------------

      CALL hdfrd3d(sd_id,"u",nx,ny,nz,uprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110
      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

      CALL hdfrd3d(sd_id,"v",nx,ny,nz,vprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110
      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

      CALL hdfrd3d(sd_id,"w",nx,ny,nz,wprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

      CALL hdfrd3d(sd_id,"pt",nx,ny,nz,ptprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110
      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

      CALL hdfrd3d(sd_id,"p",nx,ny,nz,pprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110
      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

      CALL hdfrd3d(sd_id,"qvprt",nx,ny,nz,qvprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

    ELSE

      CALL hdfrd3d(sd_id,"qv",nx,ny,nz,qvprt,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110
      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

    CALL hdfrd3d(sd_id,"qc",nx,ny,nz,qc,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"qr",nx,ny,nz,qr,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    IF( rainin == 1 ) THEN
      CALL hdfrd2d(sd_id,"raing",nx,ny,raing,istat,itmp)
      IF (istat /= 0) GO TO 110

      CALL hdfrd2d(sd_id,"rainc",nx,ny,rainc,istat,itmp)
      IF (istat /= 0) GO TO 110
    END IF

    IF( prcin == 1 ) THEN
      CALL hdfrd2d(sd_id,"prcrate1",nx,ny,prcrate(1,1,1),istat,itmp)
      IF (istat /= 0) GO TO 110

      CALL hdfrd2d(sd_id,"prcrate2",nx,ny,prcrate(1,1,2),istat,itmp)
      IF (istat /= 0) GO TO 110

      CALL hdfrd2d(sd_id,"prcrate3",nx,ny,prcrate(1,1,3),istat,itmp)
      IF (istat /= 0) GO TO 110

      CALL hdfrd2d(sd_id,"prcrate4",nx,ny,prcrate(1,1,4),istat,itmp)
      IF (istat /= 0) GO TO 110
    END IF

    IF( icein == 1 ) THEN

      CALL hdfrd3d(sd_id,"qi",nx,ny,nz,qi,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

      CALL hdfrd3d(sd_id,"qs",nx,ny,nz,qs,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

      CALL hdfrd3d(sd_id,"qh",nx,ny,nz,qh,istat,itmp,hmax,hmin)
      IF (istat /= 0) GO TO 110

    END IF
  END IF

  IF( tkein == 1 ) THEN

    CALL hdfrd3d(sd_id,"tke",nx,ny,nz,tke,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

  END IF

  IF( trbin == 1 ) THEN

    CALL hdfrd3d(sd_id,"kmh",nx,ny,nz,kmh,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"kmv",nx,ny,nz,kmv,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

  END IF

  IF( sfcin == 1 ) THEN

    CALL hdfrd3d(sd_id,"tsfc",nx,ny,nstypin+1,tsfc,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"tsoil",nx,ny,nstypin+1,tsoil,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"wetsfc",nx,ny,nstypin+1,wetsfc,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"wetdp",nx,ny,nstypin+1,wetdp,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd3d(sd_id,"wetcanp",nx,ny,nstypin+1,wetcanp,istat,itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL fix_soil_nstyp(nx,ny,nstyp1,nstyp,tsfc,tsoil,wetsfc,wetdp,wetcanp)

    IF (snowin == 1) THEN
      CALL hdfrd2d(sd_id,"snowdpth",nx,ny,snowdpth,istat,itmp)
      IF (istat /= 0) GO TO 110
    END IF

  END IF

  IF( radin == 1 ) THEN

    CALL hdfrd3d(sd_id,"radfrc",nx,ny,nz,radfrc,istat,                   &
                   itmp,hmax,hmin)
    IF (istat /= 0) GO TO 110

    CALL hdfrd2d(sd_id,"radsw",nx,ny,radsw,istat,itmp)
    IF (istat /= 0) GO TO 110

    CALL hdfrd2d(sd_id,"rnflx",nx,ny,rnflx,istat,itmp)
    IF (istat /= 0) GO TO 110

  END IF

  IF( flxin == 1 ) THEN

    CALL hdfrd2d(sd_id,"usflx",nx,ny,usflx,istat,itmp)
    IF (istat /= 0) GO TO 110

    CALL hdfrd2d(sd_id,"vsflx",nx,ny,vsflx,istat,itmp)
    IF (istat /= 0) GO TO 110

    CALL hdfrd2d(sd_id,"ptsflx",nx,ny,ptsflx,istat,itmp)
    IF (istat /= 0) GO TO 110

    CALL hdfrd2d(sd_id,"qvsflx",nx,ny,qvsflx,istat,itmp)
    IF (istat /= 0) GO TO 110

  END IF

!-----------------------------------------------------------------------
!
!  Friendly exit message
!
!-----------------------------------------------------------------------

  930   CONTINUE

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

  ireturn = 0

  GO TO 130

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

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

  GO TO 130

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

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

  130   CONTINUE

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

  DEALLOCATE (itmp,stat=istat)
  DEALLOCATE (hmax,stat=istat)
  DEALLOCATE (hmin,stat=istat)

  RETURN
END SUBROUTINE hdfread

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


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

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Produces a history data file "filename" in the NCSA HDF4 format by
!  calling HDF library subroutines.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  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
!
!    filename File name of history dump data.
!    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)
!    hterain  Terrain height (m)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!    soiltyp  Soil type
!    vegtyp   Vegetation type
!    lai      Leaf Area Index
!    roufns   Surface roughness
!    veg      Vegetation fraction
!
!    tsfc     Temperature at ground (K) (in top 1 cm layer)
!    tsoil     Deep soil temperature (K) (in deep 1 m layer)
!    wetsfc   Surface soil moisture in the top 1 cm layer
!    wetdp    Deep soil moisture in the deep 1 m layer
!    wetcanp  Canopy water amount
!
!    raing    Grid supersaturation rain
!    rainc    Cumulus convective rain
!    prcrate  Precipitation rates
!
!    radfrc   Radiation forcing (K/s)
!    radsw    Solar radiation reaching the surface
!    rnflx    Net radiation flux absorbed by surface
!
!    usflx    Surface flux of u-momentum (kg/(m*s**2))
!    vsflx    Surface flux of v-momentum (kg/(m*s**2))
!    ptsflx   Surface heat flux (K*kg/(m**2 * s ))
!    qvsflx   Surface moisture flux of (kg/(m**2 * s))
!
!  OUTPUT:
!
!    None.
!
!  WORK ARRAY:
!
!    tem1     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, 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, the eddy viscosity km is not dumped.
!  radout =0 or 1. If radout=0, the radiation arrays are not dumped.
!  flxout =0 or 1. If flxout=0, the surface fluxes are not dumped.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE

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

  CHARACTER (LEN=*) :: filename
  INTEGER :: grdbas            ! If this is a grid/base state 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 :: hterain(nx,ny)       ! Terrain height.

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

  INTEGER :: nstyps             ! Number of soil type

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

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

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

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

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

  REAL :: tem1  (nx,ny,nz)     ! Temporary work array
  INTEGER(2), allocatable :: itmp(:,:,:) ! Temporary array
  REAL, allocatable :: hmax(:), hmin(:) ! Temporary array

  REAL :: dx_out,dy_out

!-----------------------------------------------------------------------
!
!  Parameters describing this routine
!
!-----------------------------------------------------------------------

  CHARACTER (LEN=40) :: fmtver
  PARAMETER (fmtver='004.10 HDF4 Coded Data')
  CHARACTER (LEN=10) :: tmunit
  PARAMETER (tmunit='seconds   ')

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: i,j,k,l,is
  INTEGER :: nstypout
  INTEGER :: istat, sd_id
 
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'        ! Grid parameters

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

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

  IF (hdfcompr > 3) THEN
    ALLOCATE (itmp(nx,ny,nz),stat=istat)
    IF (istat /= 0) THEN
      WRITE (6,*) "HDFDUMP: ERROR allocating itmp, returning"
      RETURN
    END IF
    ALLOCATE (hmax(nz),stat=istat)
    IF (istat /= 0) THEN
      WRITE (6,*) "HDFDUMP: ERROR allocating hmax, returning"
      RETURN
    END IF
    ALLOCATE (hmin(nz),stat=istat)
    IF (istat /= 0) THEN
      WRITE (6,*) "HDFDUMP: ERROR allocating hmin, returning"
      RETURN
    END IF
  END IF

  WRITE(6,'(1x,a,f13.3,a,a/)')                                          &
       'Writing HDF4 data at time=', curtim,' into file ',filename

!-----------------------------------------------------------------------
!
!  Create the HDF4 file.
!
!-----------------------------------------------------------------------

  CALL hdfopen(filename,2,sd_id)
  IF (sd_id < 0) THEN
    WRITE (6,*) "HDFDUMP: ERROR creating HDF4 file: ",                  &
                trim(filename)
    GO TO 600
  END IF

  CALL hdfwrtc(sd_id, 40, 'fmtver', fmtver, istat)
  CALL hdfwrtc(sd_id, 40, 'runname', runname, istat)
  CALL hdfwrti(sd_id, 'nocmnt', nocmnt, istat)

  IF( nocmnt > 0 ) THEN
    CALL hdfwrtc(sd_id, 80*nocmnt, 'cmnt', cmnt, istat)
  END IF

  CALL hdfwrtc(sd_id, 7, 'tmunit', 'seconds', istat)
  CALL hdfwrtr(sd_id, 'time', curtim, istat)

  CALL hdfwrti(sd_id, 'nx', nx, istat)
  CALL hdfwrti(sd_id, 'ny', ny, istat)
  CALL hdfwrti(sd_id, 'nz', nz, istat)

  IF( grdbas == 1 ) THEN

    CALL hdfwrti(sd_id, 'grdflg', 1, istat)
    CALL hdfwrti(sd_id, 'basflg', 1, istat)
    CALL hdfwrti(sd_id, 'varflg', 0, istat)
    CALL hdfwrti(sd_id, 'mstflg', 1, istat)
    CALL hdfwrti(sd_id, 'iceflg', 0, istat)
    CALL hdfwrti(sd_id, 'trbflg', 0, istat)
    CALL hdfwrti(sd_id, 'sfcflg', 0, istat)
    CALL hdfwrti(sd_id, 'rainflg', 0, istat)
    CALL hdfwrti(sd_id, 'landflg', 1, istat)
    CALL hdfwrti(sd_id, 'totflg', 1, istat)
    CALL hdfwrti(sd_id, 'tkeflg', 0, istat)

  ELSE

    CALL hdfwrti(sd_id, 'grdflg', grdout, istat)
    CALL hdfwrti(sd_id, 'basflg', basout, istat)
    CALL hdfwrti(sd_id, 'varflg', varout, istat)
    CALL hdfwrti(sd_id, 'mstflg', mstout, istat)
    CALL hdfwrti(sd_id, 'iceflg', iceout, istat)
    CALL hdfwrti(sd_id, 'trbflg', trbout, istat)
    CALL hdfwrti(sd_id, 'sfcflg', sfcout, istat)
    CALL hdfwrti(sd_id, 'rainflg', rainout, istat)
    CALL hdfwrti(sd_id, 'landflg', landout*basout, istat)
    CALL hdfwrti(sd_id, 'totflg', totout, istat)
    CALL hdfwrti(sd_id, 'tkeflg', tkeout, istat)

  END IF

  nstypout = max(1,nstyp)
  CALL hdfwrti(sd_id, 'nstyp', nstypout, istat)
  CALL hdfwrti(sd_id, 'prcflg', prcout, istat)
  CALL hdfwrti(sd_id, 'radflg', radout, istat)
  CALL hdfwrti(sd_id, 'flxflg', flxout, istat)
  CALL hdfwrti(sd_id, 'snowflg', snowout, istat)

  CALL hdfwrti(sd_id, 'day', day, istat)
  CALL hdfwrti(sd_id, 'year', year, istat)
  CALL hdfwrti(sd_id, 'month', month, istat)
  CALL hdfwrti(sd_id, 'hour', hour, istat)
  CALL hdfwrti(sd_id, 'minute', minute, istat)
  CALL hdfwrti(sd_id, 'second', second, istat)

  CALL hdfwrtr(sd_id, 'umove', umove, istat)
  CALL hdfwrtr(sd_id, 'vmove', vmove, istat)
  CALL hdfwrtr(sd_id, 'xgrdorg', xgrdorg, istat)
  CALL hdfwrtr(sd_id, 'ygrdorg', ygrdorg, istat)

  CALL hdfwrti(sd_id, 'mapproj', mapproj, istat)
  CALL hdfwrtr(sd_id, 'trulat1', trulat1, istat)
  CALL hdfwrtr(sd_id, 'trulat2', trulat2, istat)
  CALL hdfwrtr(sd_id, 'trulon', trulon, istat)
  CALL hdfwrtr(sd_id, 'sclfct', sclfct, istat)
  CALL hdfwrtr(sd_id, 'tstop', tstop, istat)
  CALL hdfwrtr(sd_id, 'thisdmp', thisdmp, istat)
  CALL hdfwrtr(sd_id, 'latitud', latitud, istat)
  CALL hdfwrtr(sd_id, 'ctrlat', ctrlat, istat)
  CALL hdfwrtr(sd_id, 'ctrlon', ctrlon, istat)

  dx_out = x(2) - x(1)
  dy_out = y(2) - y(1)
  CALL hdfwrtr(sd_id, 'dx', dx_out, istat)
  CALL hdfwrtr(sd_id, 'dy', dy_out, istat)

!-----------------------------------------------------------------------
!
!  If grdout=1 or grdbas=1, write out grid variables
!
!-----------------------------------------------------------------------

  IF(grdout == 1 .OR. grdbas == 1 ) THEN

    CALL hdfwrt1d(x,nx,sd_id,'x','x coordinate','m')
    CALL hdfwrt1d(y,ny,sd_id,'y','y coordinate','m')
    CALL hdfwrt1d(z,nz,sd_id,'z','z coordinate','m')
    CALL hdfwrt3d(zp,nx,ny,nz,sd_id,0,hdfcompr,                         &
                  'zp','Physical height coordinate','m',                &
                   itmp,hmax,hmin)

  END IF

!-----------------------------------------------------------------------
!
!  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 hdfwrt3d(ubar,nx,ny,nz,sd_id,1,hdfcompr,                       &
                  'ubar','Base state u-velocity','m/s',                 &
                   itmp,hmax,hmin)

    CALL edgfill(vbar,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
    CALL hdfwrt3d(vbar,nx,ny,nz,sd_id,2,hdfcompr,                       &
                  'vbar','Base state v-velocity','m/s',                 &
                   itmp,hmax,hmin)

    DO k=1,nz
      DO j=1,ny
        DO i=1,nx
          tem1(i,j,k) = 0.0
        END DO
      END DO
    END DO
    CALL hdfwrt3d(tem1,nx,ny,nz,sd_id,3,hdfcompr,                       &
                  'wbar','Base state w-velocity','m/s',                 &
                   itmp,hmax,hmin)

    CALL edgfill(ptbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL hdfwrt3d(ptbar,nx,ny,nz,sd_id,0,hdfcompr,                      &
                  'ptbar','Base state potential temperature','K',       &
                   itmp,hmax,hmin)

    CALL edgfill(pbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL hdfwrt3d(pbar,nx,ny,nz,sd_id,0,hdfcompr,                       &
                  'pbar','Base state pressure','Pascal',                &
                   itmp,hmax,hmin)

    IF(mstout == 1) THEN

      CALL edgfill(qvbar,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL hdfwrt3d(qvbar,nx,ny,nz,sd_id,0,hdfcompr,                    &
          'qvbar','Base state water vapor specific humidity','kg/kg',   &
                   itmp,hmax,hmin)

    END IF

    IF(landout == 1) THEN

      CALL iedgfill(soiltyp(1,1,1),1,nx,1,nx-1, 1,ny,1,ny-1,          &
                      1,nstypout,1,nstypout)
      CALL hdfwrt3di(soiltyp,nx,ny,nstypout,sd_id,0,0,                &
                'soiltyp','Soil type','index')
      CALL edgfill(stypfrct(1,1,1),1,nx,1,nx-1, 1,ny,1,ny-1,          &
                     1,nstypout,1,nstypout)
      CALL hdfwrt3d(stypfrct(1,1,1),nx,ny,nstypout,sd_id,0,hdfcompr,  &
          'stypfrct','Soil type fractional coverage','fraction',      &
                 itmp,hmax,hmin)

      CALL iedgfill(vegtyp ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL hdfwrt2di(vegtyp,nx,ny,sd_id,0,0,                            &
                  'vegtyp','Vegetation type','index')

      CALL edgfill(lai    ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL hdfwrt2d(lai,nx,ny,sd_id,0,hdfcompr,                         &
                  'lai','Leaf Area Index','index',itmp)

      CALL edgfill(roufns ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL hdfwrt2d(roufns,nx,ny,sd_id,0,hdfcompr,                      &
                  'roufns','Surface roughness','0-1',itmp)

      CALL edgfill(veg    ,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL hdfwrt2d(veg,nx,ny,sd_id,0,hdfcompr,                         &
                  'veg','Vegetation fraction','fraction',itmp)

    END IF

  END IF

  IF ( grdbas == 1 ) GO TO 600

!-----------------------------------------------------------------------
!
!  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 hdfwrt3d(tem1,nx,ny,nz,sd_id,1,hdfcompr,                     &
                    'uprt','Perturbation u-velocity','m/s',             &
                     itmp,hmax,hmin)

      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 hdfwrt3d(tem1,nx,ny,nz,sd_id,2,hdfcompr,                     &
                    'vprt','Perturbation v-velocity','m/s',             &
                     itmp,hmax,hmin)

      CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
      CALL hdfwrt3d(w,nx,ny,nz,sd_id,3,hdfcompr,                        &
                    'wprt','Perturbation w-velocity','m/s',             &
                     itmp,hmax,hmin)

!-----------------------------------------------------------------------
!
!  Write out scalars
!
!-----------------------------------------------------------------------

      CALL edgfill(ptprt,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL hdfwrt3d(ptprt,nx,ny,nz,sd_id,0,hdfcompr,                    &
                    'ptprt','Perturbation potential temperature','K',   &
                     itmp,hmax,hmin)

      CALL edgfill(pprt,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL hdfwrt3d(pprt,nx,ny,nz,sd_id,0,hdfcompr,                     &
                    'pprt','Perturbation pressure','Pascal',            &
                     itmp,hmax,hmin)

    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 hdfwrt3d(u,nx,ny,nz,sd_id,1,hdfcompr,                        &
                    'u','u-velocity','m/s',                             &
                     itmp,hmax,hmin)

      CALL edgfill(v,1,nx,1,nx-1, 1,ny,1,ny, 1,nz,1,nz-1)
      CALL hdfwrt3d(v,nx,ny,nz,sd_id,2,hdfcompr,                        &
                    'v','v-velocity','m/s',                             &
                     itmp,hmax,hmin)

      CALL edgfill(w,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz)
      CALL hdfwrt3d(w,nx,ny,nz,sd_id,3,hdfcompr,                        &
                    'w','w-velocity','m/s',                             &
                     itmp,hmax,hmin)

!-----------------------------------------------------------------------
!
!  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 hdfwrt3d(tem1,nx,ny,nz,sd_id,0,hdfcompr,                     &
                    'pt','Potential temperature','K',                   &
                     itmp,hmax,hmin)

      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 hdfwrt3d(tem1,nx,ny,nz,sd_id,0,hdfcompr,                     &
                    'p','Pressure','Pascal',                            &
                     itmp,hmax,hmin)

    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 hdfwrt3d(tem1,nx,ny,nz,sd_id,0,hdfcompr,                     &
          'qvprt','Pert. water vapor specific humidity','kg/kg',        &
                     itmp,hmax,hmin)

    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 hdfwrt3d(qv,nx,ny,nz,sd_id,0,hdfcompr,                       &
                    'qv','Water vapor specific humidity','kg/kg',       &
                     itmp,hmax,hmin)

    END IF

    CALL edgfill(qc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL hdfwrt3d(qc,nx,ny,nz,sd_id,0,hdfcompr,                         &
                  'qc','Cloud water mixing ratio','kg/kg',              &
                   itmp,hmax,hmin)

    CALL edgfill(qr,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL hdfwrt3d(qr,nx,ny,nz,sd_id,0,hdfcompr,                         &
                  'qr','Rain water mixing ratio','kg/kg',               &
                   itmp,hmax,hmin)

    IF(rainout == 1) THEN

      CALL edgfill(raing,   1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL hdfwrt2d(raing,nx,ny,sd_id,0,hdfcompr,                       &
                  'raing','Grid supersaturation rain','mm',itmp)

      CALL edgfill(rainc,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
      CALL hdfwrt2d(rainc,nx,ny,sd_id,0,hdfcompr,                       &
                  'rainc','Cumulus convective rain','mm',itmp)

    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 hdfwrt2d(prcrate(1,1,1),nx,ny,sd_id,0,hdfcompr,              &
          'prcrate1','Total precip. rate','kg/(m**2*s)',itmp)
      CALL hdfwrt2d(prcrate(1,1,2),nx,ny,sd_id,0,hdfcompr,              &
          'prcrate2','Grid scale precip. rate','kg/(m**2*s)',itmp)
      CALL hdfwrt2d(prcrate(1,1,3),nx,ny,sd_id,0,hdfcompr,              &
          'prcrate3','Cumulative precip. rate','kg/(m**2*s)',itmp)
      CALL hdfwrt2d(prcrate(1,1,4),nx,ny,sd_id,0,hdfcompr,              &
          'prcrate4','Microphysics precip. rate','kg/(m**2*s)',itmp)
    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 hdfwrt3d(qi,nx,ny,nz,sd_id,0,hdfcompr,                       &
                  'qi','Cloud ice mixing ratio','kg/kg',                &
                   itmp,hmax,hmin)

      CALL edgfill(qs,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL hdfwrt3d(qs,nx,ny,nz,sd_id,0,hdfcompr,                       &
                  'qs','Snow mixing ratio','kg/kg',                     &
                   itmp,hmax,hmin)

      CALL edgfill(qh,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
      CALL hdfwrt3d(qh,nx,ny,nz,sd_id,0,hdfcompr,                       &
                  'qh','Hail mixing ratio','kg/kg',                     &
                   itmp,hmax,hmin)

    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 hdfwrt3d(tke,nx,ny,nz,sd_id,0,hdfcompr,                        &
                  'tke','Turbulent Kinetic Energy','(m/s)**2',          &
                   itmp,hmax,hmin)

  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 hdfwrt3d(kmh,nx,ny,nz,sd_id,0,hdfcompr,                        &
        'kmh','Hori. turb. mixing coef. for momentum','m**2/s',         &
                   itmp,hmax,hmin)

    CALL edgfill(kmv,1,nx,1,nx-1, 1,ny,1,ny-1, 1,nz,1,nz-1)
    CALL hdfwrt3d(kmv,nx,ny,nz,sd_id,0,hdfcompr,                        &
        'kmv','Vert. turb. mixing coef. for momentum','m**2/s',         &
                   itmp,hmax,hmin)

  END IF   ! trbout


!-----------------------------------------------------------------------
!
!  If sfcout = 1, write out the surface variables,
!  tsfc, tsoil, wetsfc, wetdp, and wetcanp.
!
!-----------------------------------------------------------------------

  IF( sfcout == 1) THEN

    DO is=0,nstypout

      CALL edgfill(tsfc(1,1,is),   1,nx,1,nx-1, 1,ny,1,ny-1,          &
                   1,1,1,1)
      CALL edgfill(tsoil(1,1,is),  1,nx,1,nx-1, 1,ny,1,ny-1,          &
                   1,1,1,1)
      CALL edgfill(wetsfc(1,1,is), 1,nx,1,nx-1, 1,ny,1,ny-1,          &
                   1,1,1,1)
      CALL edgfill(wetdp(1,1,is),  1,nx,1,nx-1, 1,ny,1,ny-1,          &
                   1,1,1,1)
      CALL edgfill(wetcanp(1,1,is),1,nx,1,nx-1, 1,ny,1,ny-1,          &
                   1,1,1,1)

    END DO

    CALL hdfwrt3d(tsfc,nx,ny,nstypout+1,sd_id,0,hdfcompr,                &
                'tsfc','Surface ground temperature','K',              &
                 itmp,hmax,hmin)
    CALL hdfwrt3d(tsoil,nx,ny,nstypout+1,sd_id,0,hdfcompr,               &
                'tsoil','Deep soil temperature','K',                  &
                 itmp,hmax,hmin)
    CALL hdfwrt3d(wetsfc,nx,ny,nstypout+1,sd_id,0,hdfcompr,              &
                'wetsfc','Surface soil moisture','fraction',          &
                 itmp,hmax,hmin)
    CALL hdfwrt3d(wetdp,nx,ny,nstypout+1,sd_id,0,hdfcompr,               &
                'wetdp','Deep soil moisture','fraction',              &
                 itmp,hmax,hmin)
    CALL hdfwrt3d(wetcanp,nx,ny,nstypout+1,sd_id,0,hdfcompr,             &
                'wetcanp','Canopy water amount','fraction',           &
                 itmp,hmax,hmin)

    IF (snowout == 1) THEN

      CALL edgfill(snowdpth,1,nx,1,nx-1, 1,ny,1,ny-1,1,1,1,1)
      CALL hdfwrt2d(snowdpth,nx,ny,sd_id,0,hdfcompr,                    &
                  'snowdpth','Snow depth','m',itmp)
    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 hdfwrt3d(radfrc,nx,ny,nz,sd_id,0,hdfcompr,                     &
                  'radfrc','Radiation forcing','K/s',                   &
                   itmp,hmax,hmin)

    CALL edgfill(radsw,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL hdfwrt2d(radsw,nx,ny,sd_id,0,hdfcompr,                         &
        'radsw','Solar radiation reaching the surface','W/m**2',itmp)

    CALL edgfill(rnflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL hdfwrt2d(rnflx,nx,ny,sd_id,0,hdfcompr,                         &
        'rnflx','Net radiation flux absorbed by surface','W/m**2',      &
        itmp)

  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 hdfwrt2d(usflx,nx,ny,sd_id,0,hdfcompr,                         &
        'usflx','Surface flux of u-momentum','kg/(m*s**2)',itmp)

    CALL edgfill(vsflx,1,nx,1,nx-1, 1,ny,1,ny, 1,1,1,1)
    CALL hdfwrt2d(vsflx,nx,ny,sd_id,0,hdfcompr,                         &
        'vsflx','Surface flux of v-momentum','kg/(m*s**2)',itmp)

    CALL edgfill(ptsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL hdfwrt2d(ptsflx,nx,ny,sd_id,0,hdfcompr,                        &
        'ptsflx','Surface heat flux','K*kg/(m*s**2)',itmp)

    CALL edgfill(qvsflx,1,nx,1,nx-1, 1,ny,1,ny-1, 1,1,1,1)
    CALL hdfwrt2d(qvsflx,nx,ny,sd_id,0,hdfcompr,                        &
        'qvsflx','Surface moisture flux','kg/(m**2*s)',itmp)

  END IF   ! flxout

  600   CONTINUE

  CALL hdfclose(sd_id,istat)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFDUMP: ERROR on closing file ",trim(filename),       &
                " (status",istat,")"
  END IF

  IF (hdfcompr > 3) THEN
    DEALLOCATE (itmp,stat=istat)
    DEALLOCATE (hmax,stat=istat)
    DEALLOCATE (hmin,stat=istat)
  END IF

  RETURN
END SUBROUTINE hdfdump

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


SUBROUTINE hdfwrt3d(var,nx,ny,nz,sd_id,stag_dim,hdfcompr,               & 45,1
           name,comment,units,itmp,hmax,hmin)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a 3-D real array to an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    var      Array to be written to the file
!
!    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
!
!    sd_id    HDF id of the output file
!
!    stag_dim Dimension of grid staggering (0-none, 1-x, 2-y, 3-z)
!
!    hdfcompr Compression flag (0-none)
!
!    name     Variable name
!    comment  Destriptive string
!    units    String destribing units of var
!
!    itmp     Scratch array for mapping reals to integers (used for
!             some values of hdfcompr)
!    hmax     Used to store maximum values as a function of z
!    hmin     Used to store minimum values as a function of z
!
!  OUTPUT:
!
!    None.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: var(nx,ny,nz)
  INTEGER :: sd_id, stag_dim, hdfcompr
  CHARACTER (LEN=*) :: name, comment, units
  INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz)
  REAL :: hmax(nz), hmin(nz)

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: istat
  INTEGER :: dims(3),start(3),stride(3)
  INTEGER :: comp_prm(1)
  INTEGER :: sds_id

  INTEGER :: i,j,k

  REAL :: scale
  INTEGER :: itmp1

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfcreate, sfscompress, sfscatt, sfsnatt,                   &
          sfwdata, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  dims(3) = nz
  start(1) = 0
  start(2) = 0
  start(3) = 0
  stride(1) = 1
  stride(2) = 1
  stride(3) = 1

!-----------------------------------------------------------------------
!
!  Create an entry for the variable, set compression and add attributes.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFWRT3D: Writing variable ",trim(name)
  IF (hdfcompr <= 3) THEN
    sds_id = sfcreate(sd_id, trim(name), dfnt_float32, 3, dims)
  ELSE
    sds_id = sfcreate(sd_id, trim(name), dfnt_int16, 3, dims)
  END IF

  comp_prm(1) = 0
  IF (hdfcompr == 1 .OR. hdfcompr == 5) THEN  ! quick gzip
    comp_prm(1) = 1
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 2 .OR. hdfcompr == 6) THEN  ! high gzip
    comp_prm(1) = 6
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 3 .OR. hdfcompr == 7) THEN  ! huffman
    comp_prm(1) = 4   ! this may be a problem on a Cray
    istat = sfscompress(sds_id, 3, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 3)
  END IF
  istat = sfscatt(sds_id, 'hdf_comp_prm', dfnt_int32, 1, comp_prm)

  IF (len_trim(comment) > 0) THEN
    istat = sfscatt(sds_id, 'comment', dfnt_char8,                       &
              len_trim(comment), comment)
  ELSE
    istat = sfscatt(sds_id, 'comment', dfnt_char8, 1, ' ')
  END IF
  IF (len_trim(units) > 0) THEN
    istat = sfscatt(sds_id, 'units', dfnt_char8,                         &
            len_trim(units), units)
  ELSE
    istat = sfscatt(sds_id, 'units', dfnt_char8, 1, ' ')
  END IF
  istat = sfsnatt(sds_id, 'stag_dim', dfnt_int32, 1, stag_dim)

!-----------------------------------------------------------------------
!
!  If called for, map reals to 16 bit integers
!
!-----------------------------------------------------------------------

  IF (hdfcompr > 3) THEN
    DO k=1,nz
      CALL a3dmax0lcl(var(1,1,k),1,nx,1,nx,1,ny,1,ny,1,1,1,1,           &
                   hmax(k),hmin(k))
      IF (hmax(k) == hmin(k)) hmax(k) = hmin(k) + 1.
      scale = 65534. / (hmax(k) - hmin(k))
      DO j=1,ny
        DO i=1,nx
          itmp1 = nint(scale * (var(i,j,k) - hmin(k))) - 32767
          itmp(i,j,k) = itmp1
        END DO
      END DO
    END DO
    istat = sfsnatt(sds_id, 'packed16', dfnt_int32, 1, 1)
    istat = sfsnatt(sds_id, 'max', dfnt_float32, nz, hmax)
    istat = sfsnatt(sds_id, 'min', dfnt_float32, nz, hmin)
  END IF

!-----------------------------------------------------------------------
!
!  Write the data in var out to the file
!
!-----------------------------------------------------------------------

  IF (hdfcompr <= 3) THEN
    istat = sfwdata(sds_id, start, stride, dims, var)
  ELSE
    istat = sfwdata(sds_id, start, stride, dims, itmp)
  END IF

  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT3D: ERROR writing variable ",trim(name)
  END IF

  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT3D: ERROR writing variable ",trim(name)
  END IF

  RETURN
END SUBROUTINE hdfwrt3d

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


SUBROUTINE hdfwrt3di(var,nx,ny,nz,sd_id,stag_dim,hdfcompr,              & 1
           name,comment,units)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a 3-D integer array to an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    var      Array to be written to the file
!
!    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
!
!    sd_id    HDF id of the output file
!
!    stag_dim Dimension of grid staggering (0-none, 1-x, 2-y, 3-z)
!
!    hdfcompr Compression flag (0-none)
!
!    name     Variable name
!    comment  Destriptive string
!    units    String destribing units of var
!
!  OUTPUT:
!
!    None.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: var(nx,ny,nz)
  INTEGER :: sd_id, stag_dim, hdfcompr
  CHARACTER (LEN=*) :: name, comment, units

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: istat
  INTEGER :: dims(3),start(3),stride(3)
  INTEGER :: sds_id
  INTEGER :: comp_prm(1)

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfcreate, sfscompress, sfscatt, sfsnatt,                   &
          sfwdata, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  dims(3) = nz
  start(1) = 0
  start(2) = 0
  start(3) = 0
  stride(1) = 1
  stride(2) = 1
  stride(3) = 1

!-----------------------------------------------------------------------
!
!  Create an entry for the variable, set compression and add attributes.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFWRT3DI: Writing variable ",trim(name)
  sds_id = sfcreate(sd_id, name, dfnt_int32, 3, dims)

  comp_prm(1) = 0
  IF (hdfcompr == 1 .OR. hdfcompr == 5) THEN  ! quick gzip
    comp_prm(1) = 1
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 2 .OR. hdfcompr == 6) THEN  ! high gzip
    comp_prm(1) = 6
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 3 .OR. hdfcompr == 7) THEN  ! huffman
    comp_prm(1) = 4   ! this may be a problem on a Cray
    istat = sfscompress(sds_id, 3, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 3)
  END IF
  istat = sfscatt(sds_id, 'hdf_comp_prm', dfnt_int32, 1, comp_prm)

  IF (len_trim(comment) > 0) THEN
    istat = sfscatt(sds_id, 'comment', dfnt_char8,                       &
              len_trim(comment), comment)
  ELSE
    istat = sfscatt(sds_id, 'comment', dfnt_char8, 1, ' ')
  END IF
  IF (len_trim(units) > 0) THEN
    istat = sfscatt(sds_id, 'units', dfnt_char8,                         &
            len_trim(units), units)
  ELSE
    istat = sfscatt(sds_id, 'units', dfnt_char8, 1, ' ')
  END IF
  istat = sfsnatt(sds_id, 'stag_dim', dfnt_int32, 1, stag_dim)

!-----------------------------------------------------------------------
!
!  Write the data in var out to the file
!
!-----------------------------------------------------------------------

  istat = sfwdata(sds_id, start, stride, dims, var)

  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT3DI: ERROR writing variable ",trim(name)
  END IF
  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT3DI: ERROR writing variable ",trim(name)
  END IF

  RETURN
END SUBROUTINE hdfwrt3di

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


SUBROUTINE hdfwrt2d(var,nx,ny,sd_id,stag_dim,hdfcompr,                  & 16,1
           name,comment,units,itmp)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a 2-D real array to an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    var      Array to be written to the file
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!
!    sd_id    HDF id of the output file
!
!    stag_dim Dimension of grid staggering (0-none, 1-x, 2-y, 3-z)
!
!    hdfcompr Compression flag (0-none)
!
!    name     Variable name
!    comment  Destriptive string
!    units    String destribing units of var
!
!    itmp     Scratch array for mapping reals to integers (used for
!             some values of hdfcompr)
!
!  OUTPUT:
!
!    None.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: nx,ny
  REAL :: var(nx,ny)
  INTEGER :: sd_id, stag_dim, hdfcompr
  CHARACTER (LEN=*) :: name, comment, units
  INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny)

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: istat
  INTEGER :: dims(2),start(2),stride(2)
  INTEGER :: sds_id
  INTEGER :: comp_prm(1)
  REAL :: amax,amin,scale
  INTEGER :: i,j,k
  INTEGER :: itmp1

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfcreate, sfscompress, sfscatt, sfsnatt,                   &
          sfwdata, sfendacc

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

  IF (myproc == 0) WRITE (6,*) "HDFWRT2D: Writing variable ",trim(name)

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  start(1) = 0
  start(2) = 0
  stride(1) = 1
  stride(2) = 1

!-----------------------------------------------------------------------
!
!  Create an entry for the variable, set compression and add attributes.
!
!-----------------------------------------------------------------------

  IF (hdfcompr <= 3) THEN
    sds_id = sfcreate(sd_id, trim(name), dfnt_float32, 2, dims)
  ELSE
    sds_id = sfcreate(sd_id, trim(name), dfnt_int16, 2, dims)
  END IF

  comp_prm(1) = 0
  IF (hdfcompr == 1 .OR. hdfcompr == 5) THEN  ! quick gzip
    comp_prm(1) = 1
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 2 .OR. hdfcompr == 6) THEN  ! high gzip
    comp_prm(1) = 6
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 3 .OR. hdfcompr == 7) THEN  ! huffman
    comp_prm(1) = 4   ! this may be a problem on a Cray
    istat = sfscompress(sds_id, 3, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 3)
  END IF
  istat = sfscatt(sds_id, 'hdf_comp_prm', dfnt_int32, 1, comp_prm)

  IF (len_trim(comment) > 0) THEN
    istat = sfscatt(sds_id, 'comment', dfnt_char8,                       &
              len_trim(comment), comment)
  ELSE
    istat = sfscatt(sds_id, 'comment', dfnt_char8, 1, ' ')
  END IF
  IF (len_trim(units) > 0) THEN
    istat = sfscatt(sds_id, 'units', dfnt_char8,                         &
            len_trim(units), units)
  ELSE
    istat = sfscatt(sds_id, 'units', dfnt_char8, 1, ' ')
  END IF
  istat = sfsnatt(sds_id, 'stag_dim', dfnt_int32, 1, stag_dim)

!-----------------------------------------------------------------------
!
!  If called for, map reals to 16 bit integers
!
!-----------------------------------------------------------------------

  IF (hdfcompr > 3) THEN
    CALL a3dmax0lcl(var,1,nx,1,nx,1,ny,1,ny,1,1,1,1,amax,amin)
    IF (amax == amin) amax = amin + 1.
    scale = 65534. / (amax - amin)
    DO j=1,ny
      DO i=1,nx
        itmp1 = nint(scale * (var(i,j) - amin)) - 32767
        itmp(i,j) = itmp1
      END DO
    END DO
    istat = sfsnatt(sds_id, 'packed16', dfnt_int32, 1, 1)
    istat = sfsnatt(sds_id, 'max', dfnt_float32, 1, amax)
    istat = sfsnatt(sds_id, 'min', dfnt_float32, 1, amin)
  END IF

!-----------------------------------------------------------------------
!
!  Write the data in var out to the file
!
!-----------------------------------------------------------------------

  IF (hdfcompr <= 3) THEN
    istat = sfwdata(sds_id, start, stride, dims, var)
  ELSE
    istat = sfwdata(sds_id, start, stride, dims, itmp)
  END IF

  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT2D: ERROR writing variable ",trim(name)
  END IF
  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT2D: ERROR writing variable ",trim(name)
  END IF

  RETURN
END SUBROUTINE hdfwrt2d

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


SUBROUTINE hdfwrt2di(var,nx,ny,sd_id,stag_dim,hdfcompr,                 & 1
           name,comment,units)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a 2-D integer array to an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    var      Array to be written to the file
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!
!    sd_id    HDF id of the output file
!
!    stag_dim Dimension of grid staggering (0-none, 1-x, 2-y, 3-z)
!
!    hdfcompr Compression flag (0-none)
!
!    name     Variable name
!    comment  Destriptive string
!    units    String destribing units of var
!
!  OUTPUT:
!
!    None.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
  INTEGER :: nx,ny
  INTEGER :: var(nx,ny)
  INTEGER :: sd_id, stag_dim, hdfcompr
  CHARACTER (LEN=*) :: name, comment, units

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: istat
  INTEGER :: dims(2),start(2),stride(2)
  INTEGER :: sds_id
  INTEGER :: comp_prm(1)

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfcreate, sfscompress, sfscatt, sfsnatt,                   &
          sfwdata, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  start(1) = 0
  start(2) = 0
  stride(1) = 1
  stride(2) = 1

!-----------------------------------------------------------------------
!
!  Create an entry for the variable, set compression and add attributes.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFWRT2DI: Writing variable ",trim(name)

  sds_id = sfcreate(sd_id, trim(name), dfnt_int32, 2, dims)

  comp_prm(1) = 0
  IF (hdfcompr == 1 .OR. hdfcompr == 5) THEN  ! quick gzip
    comp_prm(1) = 1
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 2 .OR. hdfcompr == 6) THEN  ! high gzip
    comp_prm(1) = 6
    istat = sfscompress(sds_id, 4, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 4)
  ELSE IF (hdfcompr == 3 .OR. hdfcompr == 7) THEN  ! huffman
    comp_prm(1) = 4   ! this may be a problem on a Cray
    istat = sfscompress(sds_id, 3, comp_prm)
    istat = sfscatt(sds_id, 'hdf_comp_code', dfnt_int32, 1, 3)
  END IF
  istat = sfscatt(sds_id, 'hdf_comp_prm', dfnt_int32, 1, comp_prm)

  IF (len_trim(comment) > 0) THEN
    istat = sfscatt(sds_id, 'comment', dfnt_char8,                       &
              len_trim(comment), comment)
  ELSE
    istat = sfscatt(sds_id, 'comment', dfnt_char8, 1, ' ')
  END IF
  IF (len_trim(units) > 0) THEN
    istat = sfscatt(sds_id, 'units', dfnt_char8,                         &
            len_trim(units), units)
  ELSE
    istat = sfscatt(sds_id, 'units', dfnt_char8, 1, ' ')
  END IF
  istat = sfsnatt(sds_id, 'stag_dim', dfnt_int32, 1, stag_dim)

!-----------------------------------------------------------------------
!
!  Write the data in var out to the file
!
!-----------------------------------------------------------------------

  istat = sfwdata(sds_id, start, stride, dims, var)

  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT2DI: ERROR writing variable ",trim(name)
  END IF
  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT2DI: ERROR writing variable ",trim(name)
  END IF

  RETURN
END SUBROUTINE hdfwrt2di

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


SUBROUTINE hdfwrt1d(var,num,sd_id,name,comment,units) 3

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a 1-D real array to an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    var      Array to be written to the file
!
!    num      Number of grid points
!
!    sd_id    HDF id of the output file
!
!    name     Variable name
!    comment  Destriptive string
!    units    String destribing units of var
!
!  OUTPUT:
!
!    None.
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: num
  REAL :: var(num)
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name, comment, units

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: istat
  INTEGER :: dims(1),start(1),stride(1)
  INTEGER :: sds_id

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfcreate, sfscompress, sfscatt, sfsnatt,                   &
          sfwdata, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = num
  start(1) = 0
  stride(1) = 1

!-----------------------------------------------------------------------
!
!  Create an entry for the variable, set compression and add attributes.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFWRT1D: Writing variable ",trim(name)

  sds_id = sfcreate(sd_id, trim(name), dfnt_float32, 1, dims)

  IF (len_trim(comment) > 0) THEN
    istat = sfscatt(sds_id, 'comment', dfnt_char8,                       &
              len_trim(comment), comment)
  ELSE
    istat = sfscatt(sds_id, 'comment', dfnt_char8, 1, ' ')
  END IF
  IF (len_trim(units) > 0) THEN
    istat = sfscatt(sds_id, 'units', dfnt_char8,                         &
            len_trim(units), units)
  ELSE
    istat = sfscatt(sds_id, 'units', dfnt_char8, 1, ' ')
  END IF

!-----------------------------------------------------------------------
!
!  Write the data in var out to the file
!
!-----------------------------------------------------------------------

  istat = sfwdata(sds_id, start, stride, dims, var)

  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT1D: ERROR writing variable ",trim(name)
  END IF
  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFWRT1D: ERROR writing variable ",trim(name)
  END IF

  RETURN
END SUBROUTINE hdfwrt1d

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


SUBROUTINE hdfwrtr(sd_id,name,val,istat) 31

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a real attribute
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/05
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the file or variable containing the
!             named attribute
!
!  OUTPUT:
!
!    val      The value of the attribute
!    istat     Status of the read (0-okay, 1-write error)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  REAL :: val
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  INTEGER :: istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfsnatt

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

  istat = sfsnatt(sd_id, trim(name), dfnt_float32, 1, val)

  IF (istat == -1) THEN
    WRITE (6,*) "HDFWRTR: ERROR writing variable ",trim(name),"."
    istat = 1
  ELSE
    WRITE (6,*) "HDFWRTR: Wrote variable ",trim(name)," value:",val
  END IF

  RETURN
END SUBROUTINE hdfwrtr

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


SUBROUTINE hdfwrti(sd_id,name,val,istat) 50

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out an integer attribute
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/05
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the file or variable containing the
!             named attribute
!
!  OUTPUT:
!
!    val      The value of the attribute
!    istat     Status of the read (0-okay, 1-write error)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: val
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  INTEGER :: istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfsnatt

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

  istat = sfsnatt(sd_id, trim(name), dfnt_int32, 1, val)

  IF (istat == -1) THEN
    WRITE (6,*) "HDFWRTI: ERROR writing variable ",trim(name),"."
    istat = 1
  ELSE
    WRITE (6,*) "HDFWRTI: Wrote variable ",trim(name)," value:",val
  END IF

  RETURN
END SUBROUTINE hdfwrti

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


SUBROUTINE hdfwrtc(sd_id,strlen,name,string,istat) 4

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Write out a string attribute
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/11
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    strlen   Length of string to be written out (set to len(string)
!             if strlen passed in as 0)
!
!    sd_id    HDF id of the file or variable containing the
!             named attribute
!
!  OUTPUT:
!
!    string   The string to be written out
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  CHARACTER (LEN=*) :: string
  INTEGER :: strlen, istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfscatt

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


  IF (strlen == 0) strlen = LEN(string)
  istat = sfscatt(sd_id, trim(name), dfnt_char8, strlen, string)

  IF (istat == -1) THEN
    WRITE (6,*) "HDFWRTC: ERROR writing variable ",trim(name),"."
    istat = 1
  ELSE
    WRITE (6,*) "HDFWRTC: Wrote variable ",trim(name)," value:",        &
                trim(string)
  END IF

  RETURN
END SUBROUTINE hdfwrtc

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


SUBROUTINE hdfrd3d(sd_id,name,nx,ny,nz,var,istat,itmp,hmax,hmin) 45

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in a 3-D real array from an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the output file
!
!    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
!
!    itmp     Scratch array for mapping reals to integers (used for
!             some values of hdfcompr)
!    hmax     Used to store maximum values as a function of z
!    hmin     Used to store minimum values as a function of z
!
!  OUTPUT:
!
!    var      Array to be read in
!
!    istat     Status of read (0-okay, 1-error when reading)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  REAL :: var(nx,ny,nz)
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  INTEGER :: istat
  INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny,nz)
  REAL :: hmax(nz),hmin(nz)

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: dims(3),start(3),stride(3)
  INTEGER :: sds_index,sds_id,attr_index
  INTEGER :: packed16
  INTEGER :: istat1,istat2,istat3
  REAL :: scale
  INTEGER :: i,j,k

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfn2index, sfselect, sfrdata, sffattr, sfrnatt, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  dims(3) = nz
  start(1) = 0
  start(2) = 0
  start(3) = 0
  stride(1) = 1
  stride(2) = 1
  stride(3) = 1

!-----------------------------------------------------------------------
!
!  Get the SDS ID for the variable.
!
!-----------------------------------------------------------------------

  sds_index = sfn2index(sd_id, trim(name))

  IF (sds_index == -1) THEN
    WRITE (6,*) "HDFRD3D: WARNING, variable ",                            &
                 trim(name)," not found in file."
    istat = 1
    RETURN
  END IF

  sds_id = sfselect(sd_id, sds_index)

  attr_index = sffattr(sds_id, "packed16")
  IF (attr_index >= 0) THEN
    istat1 = sfrnatt(sds_id, attr_index, packed16)
    attr_index = sffattr(sds_id, "max")
    istat2 = sfrnatt(sds_id, attr_index, hmax)
    attr_index = sffattr(sds_id, "min")
    istat3 = sfrnatt(sds_id, attr_index, hmin)
    IF (istat1 == -1 .OR. istat2 == -1 .OR. istat3 == -1) THEN
      WRITE (6,*) "HDFRD3D: ERROR reading max/min for ",trim(name)
      istat = 2
      RETURN
    END IF
  ELSE
    packed16 = 0
  END IF

!-----------------------------------------------------------------------
!
!  Read data into var.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFRD3D: Reading variable ",trim(name)
  IF (packed16 == 0) THEN
    istat = sfrdata(sds_id, start, stride, dims, var)
  ELSE
    istat = sfrdata(sds_id, start, stride, dims, itmp)
  END IF
  IF (istat == -1) THEN
    WRITE (6,*) "HDFRD3D: ERROR reading variable ",trim(name),"."
    istat = 2
    RETURN
  END IF

!-----------------------------------------------------------------------
!
!  If called for, map 16 bit integers to reals
!
!-----------------------------------------------------------------------

  IF (packed16 /= 0) THEN
    DO k=1,nz
      scale = (hmax(k) - hmin(k)) / 65534.
      DO j=1,ny
        DO i=1,nx
          var(i,j,k) = scale * (itmp(i,j,k) + 32767) + hmin(k)
        END DO
      END DO
    END DO
  END IF

  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFRD3D: ERROR reading variable ",trim(name)
    istat = 2
  END IF

  RETURN
END SUBROUTINE hdfrd3d

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


SUBROUTINE hdfrd3di(sd_id,name,nx,ny,nz,var,istat) 1

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in a 3-D integer array from an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the output file
!
!    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
!
!  OUTPUT:
!
!    var      Array to be read in
!
!    istat     Status of read (0-okay, 1-error when reading)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: nx,ny,nz
  INTEGER :: var(nx,ny,nz)
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  INTEGER :: istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: dims(3),start(3),stride(3)
  INTEGER :: sds_index,sds_id

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfn2index, sfselect, sfrdata, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  dims(3) = nz
  start(1) = 0
  start(2) = 0
  start(3) = 0
  stride(1) = 1
  stride(2) = 1
  stride(3) = 1

!-----------------------------------------------------------------------
!
!  Get the SDS ID for the variable.
!
!-----------------------------------------------------------------------

  sds_index = sfn2index(sd_id, trim(name))

  IF (sds_index == -1) THEN
    WRITE (6,*) "HDFRD3DI: WARNING, variable ",                           &
                 trim(name)," not found in file."
    istat = 1
    RETURN
  END IF

  sds_id = sfselect(sd_id, sds_index)

!-----------------------------------------------------------------------
!
!  Read data into var.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFRD3DI: Reading variable ",trim(name)
  istat = sfrdata(sds_id, start, stride, dims, var)
  IF (istat == -1) THEN
    WRITE (6,*) "HDFRD3DI: ERROR reading variable ",trim(name),"."
    istat = 2
  END IF

  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFRD3DI: ERROR reading variable ",trim(name)
    istat = 2
  END IF

  RETURN
END SUBROUTINE hdfrd3di

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


SUBROUTINE hdfrd2d(sd_id,name,nx,ny,var,istat,itmp) 16

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in a 2-D real array from an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the output file
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!
!    itmp     Scratch array for mapping reals to integers (used for
!             some values of hdfcompr)
!
!  OUTPUT:
!
!    var      Array to be read in
!
!    istat     Status of read (0-okay, 1-error when reading)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: nx,ny
  REAL :: var(nx,ny)
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  INTEGER :: istat
  INTEGER (KIND=selected_int_kind(4)) :: itmp(nx,ny)

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: dims(2),start(2),stride(2)
  INTEGER :: sds_index,sds_id,attr_index
  REAL :: amax,amin, scale
  INTEGER :: istat1,istat2,istat3
  INTEGER :: i,j
  INTEGER :: packed16

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfn2index, sfselect, sfrdata, sffattr, sfrnatt, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  start(1) = 0
  start(2) = 0
  stride(1) = 1
  stride(2) = 1

!-----------------------------------------------------------------------
!
!  Get the SDS ID for the variable.
!
!-----------------------------------------------------------------------

  sds_index = sfn2index(sd_id, trim(name))

  IF (sds_index == -1) THEN
    WRITE (6,*) "HDFRD2D: WARNING, variable ",                            &
                trim(name)," not found in file."
    istat = 1
    RETURN
  END IF

  sds_id = sfselect(sd_id, sds_index)

  attr_index = sffattr(sds_id, "packed16")
  IF (attr_index >= 0) THEN
    istat1 = sfrnatt(sds_id, attr_index, packed16)
    attr_index = sffattr(sds_id, "max")
    istat2 = sfrnatt(sds_id, attr_index, amax)
    attr_index = sffattr(sds_id, "min")
    istat3 = sfrnatt(sds_id, attr_index, amin)
    IF (istat1 == -1 .OR. istat2 == -1 .OR. istat3 == -1) THEN
      WRITE (6,*) "HDFRD2D: ERROR reading max/min for ",trim(name)
      istat = sfendacc(sds_id)
      istat = 2
      RETURN
    END IF
  ELSE
    packed16 = 0
  END IF

!-----------------------------------------------------------------------
!
!  Read data into var.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFRD2D: Reading variable ",trim(name)
  IF (packed16 == 0) THEN
    istat = sfrdata(sds_id, start, stride, dims, var)
  ELSE
    istat = sfrdata(sds_id, start, stride, dims, itmp)
  END IF
  IF (istat == -1) THEN
    WRITE (6,*) "HDFRD2D: ERROR reading variable ",trim(name),"."
    istat = sfendacc(sds_id)
    istat = 2
    RETURN
  END IF

!-----------------------------------------------------------------------
!
!  If called for, map 16 bit integers to reals
!
!-----------------------------------------------------------------------

  IF (packed16 /= 0) THEN
    scale = (amax - amin) / 65534.
    DO j=1,ny
      DO i=1,nx
        var(i,j) = scale * (itmp(i,j) + 32767) + amin
      END DO
    END DO
  END IF

  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFRD2D: ERROR reading variable ",trim(name)
    istat = 2
  END IF

  RETURN
END SUBROUTINE hdfrd2d

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


SUBROUTINE hdfrd2di(sd_id,name,nx,ny,var,istat) 1

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in a 2-D integer array from an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the output file
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!
!  OUTPUT:
!
!    var      Array to be read in
!
!    istat     Status of read (0-okay, 1-error when reading)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: nx,ny
  INTEGER :: var(nx,ny)
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  INTEGER :: istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: dims(2),start(2),stride(2)
  INTEGER :: sds_index,sds_id

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfn2index, sfselect, sfrdata, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = nx
  dims(2) = ny
  start(1) = 0
  start(2) = 0
  stride(1) = 1
  stride(2) = 1

!-----------------------------------------------------------------------
!
!  Get the SDS ID for the variable.
!
!-----------------------------------------------------------------------

  sds_index = sfn2index(sd_id, trim(name))

  IF (sds_index == -1) THEN
    WRITE (6,*) "HDFRD2DI: WARNING, variable ",                           &
                trim(name)," not found in file."
    istat = 1
    RETURN
  END IF

  sds_id = sfselect(sd_id, sds_index)

!-----------------------------------------------------------------------
!
!  Read data into var.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFRD2DI: Reading variable ",trim(name)
  istat = sfrdata(sds_id, start, stride, dims, var)
  IF (istat == -1) THEN
    WRITE (6,*) "HDFRD2DI: ERROR reading variable ",trim(name),"."
    istat = sfendacc(sds_id)
    istat = 2
    RETURN
  END IF

  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFRD2DI: ERROR reading variable ",trim(name)
    istat = 2
  END IF

  RETURN
END SUBROUTINE hdfrd2di

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


SUBROUTINE hdfrd1d(sd_id,name,num,var,istat) 3

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in a 2-D real array from an HDF4 file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the output file
!
!    num      Number of grid points
!
!  OUTPUT:
!
!    var      Array to be read in
!
!    istat     Status of read (0-okay, 1-error when reading)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: num
  REAL :: var(num)
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  INTEGER :: istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: dims(1),start(1),stride(1)
  INTEGER :: sds_index,sds_id

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfn2index, sfselect, sfrdata, sfendacc

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

!-----------------------------------------------------------------------
!
!  Initialize dimension parameters
!
!-----------------------------------------------------------------------

  dims(1) = num
  start(1) = 0
  stride(1) = 1

!-----------------------------------------------------------------------
!
!  Get the SDS ID for the variable.
!
!-----------------------------------------------------------------------

  sds_index = sfn2index(sd_id, trim(name))

  IF (sds_index == -1) THEN
    WRITE (6,*) "HDFRD1D: WARNING, variable ",                            &
                trim(name)," not found in file."
    istat = 1
    RETURN
  END IF

  sds_id = sfselect(sd_id, sds_index)

!-----------------------------------------------------------------------
!
!  Read data into var.
!
!-----------------------------------------------------------------------

  IF (myproc == 0) WRITE (6,*) "HDFRD1D: Reading variable ",trim(name)
  istat = sfrdata(sds_id, start, stride, dims, var)
  IF (istat == -1) THEN
    WRITE (6,*) "HDFRD1D: ERROR reading variable ",trim(name),"."
    istat = sfendacc(sds_id)
    istat = 2
    RETURN
  END IF

  istat = sfendacc(sds_id)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFRD1D: ERROR reading variable ",trim(name)
    istat = 2
  END IF

  RETURN
END SUBROUTINE hdfrd1d

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


SUBROUTINE hdfrdr(sd_id,name,val,istat) 29

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in a real attribute
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the file or variable containing the
!             named attribute
!
!  OUTPUT:
!
!    val      The value of the attribute
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  REAL :: val
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: attr_index, istat

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sffattr, sfrnatt

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

  attr_index = sffattr(sd_id, trim(name))
  IF (attr_index < 0) THEN
    WRITE (6,*) "HDFRDR: WARNING, variable ",                             &
                trim(name)," not found in file."
    istat = 1
    RETURN
  END IF
  istat = sfrnatt(sd_id, attr_index, val)

  IF (istat == -1) THEN
    WRITE (6,*) "HDFRDR: ERROR reading variable ",trim(name),"."
    istat = 2
  ELSE
    IF (myproc == 0)  &
       WRITE (6,*) "HDFRDR: Read variable ",trim(name)," value:",val
  END IF

  RETURN
END SUBROUTINE hdfrdr

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


SUBROUTINE hdfrdi(sd_id,name,val,istat) 62

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in an integer attribute
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    sd_id    HDF id of the file or variable containing the
!             named attribute
!
!  OUTPUT:
!
!    val      The value of the attribute
!    istat    Status indicator (0-okay, 1-variable not found, 2-read error)
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: val
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: attr_index, istat

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sffattr, sfrnatt

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

  attr_index = sffattr(sd_id, trim(name))
  IF (attr_index < 0) THEN
    WRITE (6,*) "HDFRDI: WARNING, variable ",                             &
                trim(name)," not found in file."
    istat = 1
    RETURN
  END IF
  istat = sfrnatt(sd_id, attr_index, val)

  IF (istat == -1) THEN
    WRITE (6,*) "HDFRDI: ERROR reading variable ",trim(name),"."
    istat = 2
    val = -1
  ELSE
    IF (myproc == 0)  &
       WRITE (6,*) "HDFRDI: Read variable ",trim(name)," value:",val
  END IF

  RETURN
END SUBROUTINE hdfrdi

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


SUBROUTINE hdfrdc(sd_id,max_len,name,string,istat) 4

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Read in a string attribute
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/03/15
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    name     Variable name
!
!    max_len  Maximum allowable length of string to be read in
!
!    sd_id    HDF id of the file or variable containing the
!             named attribute
!
!  OUTPUT:
!
!    string   The value of the attribute
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: name
  CHARACTER (LEN=*) :: string
  INTEGER :: max_len

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: attr_index, istat
  INTEGER :: data_type, n_values
  CHARACTER (LEN=128) :: attr_name
  INTEGER :: i

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file
  INCLUDE 'mp.inc'

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sffattr, sfgainfo, sfrnatt

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


  attr_index = sffattr(sd_id, trim(name))

  IF (attr_index < 0) THEN
    WRITE (6,*) "HDFRDC: WARNING, variable ",trim(name),                  &
                " not located."
    istat = 1
    RETURN
  END IF

  istat = sfgainfo(sd_id, attr_index, attr_name, data_type,              &
                    n_values)
  IF (istat /= 0) THEN
    WRITE (6,*) "HDFRDC: ERROR in looking up attributes for ",          &
                trim(name)
    istat = 2
    RETURN
  END IF

  IF (n_values <= max_len) THEN
    istat = sfrnatt(sd_id, attr_index, string)
  ELSE
    WRITE (6,*) "HDFRDC: ERROR: string length for variable ",           &
                trim(name),", ",max_len,                                &
                ", is less than string in file:",n_values,","
    WRITE (6,*) "   value not read in."
  END IF

  DO i=n_values+1,len(string)
    string(i:i) = " "
  END DO

  IF (istat == -1) THEN
    WRITE (6,*) "HDFRDC: ERROR reading variable ",trim(name),"."
    istat = 1
  ELSE
    IF (myproc == 0)  &
       WRITE (6,*) "HDFRDC: Read variable ",trim(name)," value:",  &
       trim(string(1:n_values))
  END IF

  RETURN
END SUBROUTINE hdfrdc

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


SUBROUTINE hdfopen(filename,rdwrtflg,sd_id) 6

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Open a HDF file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/11
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    filename     File name
!
!    rdwrtflg Read/write flag (1-read, 2-write)
!
!  OUTPUT:
!
!    sd_id    HDF id of the file or variable containing the
!             named attribute
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: sd_id
  CHARACTER (LEN=*) :: filename
  INTEGER :: rdwrtflg

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfstart

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

  IF (rdwrtflg == 1) THEN
    sd_id = sfstart(trim(filename), dfacc_read)
  ELSE IF (rdwrtflg == 2) THEN
    sd_id = sfstart(trim(filename), dfacc_create)
  ELSE
    sd_id = -1
    WRITE (6,*) "HDFOPEN: ERROR, unsupported rdwrtflg value of",        &
          rdwrtflg
  END IF

  IF (sd_id <= 0) THEN
    WRITE (6,*) "HDFOPEN: ERROR opening ",trim(filename)
  ENDIF

  RETURN
END SUBROUTINE hdfopen

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


SUBROUTINE hdfclose(sd_id,istat) 6

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Close a HDF file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/04/11
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    sd_id    HDF id of file to close
!
!  OUTPUT:
!
!    istat     Status returned by close
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: sd_id
  INTEGER :: istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfend

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

  istat = sfend(sd_id)

  RETURN
END SUBROUTINE hdfclose

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


SUBROUTINE hdfinfo(sd_id,ndata,nattr,istat)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Get the number of data sets and attributes contained in a HDF data file.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/10/18
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    sd_id     HDF id of the data file
!
!  OUTPUT:
!
!    ndata     Number of data sets in the file
!    nattr     Number of attributes in the file
!    istat     Status returned by sffinfo
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: sd_id
  INTEGER :: ndata,nattr,istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sffinfo

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

  istat = sffinfo(sd_id,ndata,nattr)

  RETURN
END SUBROUTINE hdfinfo

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


SUBROUTINE hdfdinfo(sds_id,name,rank,dims,dtype,nattr,istat)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Get the information about a data set
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/10/18
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    sds_id    HDF id of the data set
!
!  OUTPUT:
!
!    name      Name of the data set
!    rank      Number of dimensions in the data set
!    dims      Number of points for each dimension
!    dtype     Data type
!    nattr     Number of attributes in the data set
!    istat     Status returned by sfginfo
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: sds_id
  CHARACTER (LEN=*) :: name
  INTEGER :: rank,dims(6),dtype,nattr,istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfginfo

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

  istat = sfginfo(sds_id,name,rank,dims,dtype,nattr)

  RETURN
END SUBROUTINE hdfdinfo

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


SUBROUTINE hdfainfo(sd_id,aindex,name,dtype,nvalues,istat)

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Get the information about an attribute.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2000/10/25
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    sd_id     HDF id of the data set (for HDF file or SDS data set in a file)
!    aindex    Index number for the attribute
!
!  OUTPUT:
!
!    name      Name of the data set
!    dtype     Data type
!    nvalues   Number of values (number of characters if a string)
!    istat     Status returned by sfgainfo
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  INTEGER :: sd_id,aindex
  CHARACTER (LEN=*) :: name
  INTEGER :: dtype,nvalues,istat

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

  INCLUDE 'hdf.f90'    ! HDF4 library include file

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

  INTEGER :: sfgainfo

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

  istat = sfgainfo(sd_id,aindex,name,dtype,nvalues)

  RETURN
END SUBROUTINE hdfainfo

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


SUBROUTINE get_dims_from_hdf(filename,nx,ny,nz,nstyps,ireturn),6

!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Get the grid dimensions from an ARPS hdf file.  (Similar to 
!  get_dims_from_data.)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Gene Bassett
!  2001/04/23
!
!  MODIFICATION HISTORY:
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    filename  File name
!
!  OUTPUT:
!
!    nx,ny,nz  Grid dimensions
!    nstyps    Number of soil types
!    ireturn   Return status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------

  IMPLICIT NONE
  CHARACTER (LEN=*) :: filename
  INTEGER :: nx,ny,nz,nstyps,ireturn

!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------

  INTEGER :: sd_id,istat

!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------

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

  ireturn = 0

  CALL hdfopen(filename,1,sd_id)

  IF (sd_id <= 0) THEN
    ireturn = 1
    RETURN
  ENDIF

  CALL hdfrdi(sd_id,"nx",nx,istat)
  IF (istat > 0) ireturn = 1
  CALL hdfrdi(sd_id,"ny",ny,istat)
  IF (istat > 0) ireturn = 1
  CALL hdfrdi(sd_id,"nz",nz,istat)
  IF (istat > 0) ireturn = 1
  CALL hdfrdi(sd_id,"nstyp",nstyps,istat)
  IF (istat > 0) ireturn = 1

  CALL hdfclose(sd_id,istat)

  RETURN
END SUBROUTINE get_dims_from_hdf