c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE HDFREAD                    ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE hdfread(nx,ny,nz, grdbas, filename, time, x,y,z,zp,
     :           uprt, vprt, wprt, ptprt, pprt,
     :           qvprt, qc, qr, qi, qs, qh, km,
     :           ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar,
     :           ireturn, tem1)

c
c#######################################################################
c
c     PURPOSE:
c
c     Read in data fields from an HDF scientific data set.
c
c     All data read in are located at the grid cell centers.  Variables 
c     u, v and w were averaged to the cell center before being dumped. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     2/15/92.    
c
c     MODIFICATION HISTORY:
c
c     6/08/92  Added full documentation (K. Brewster)
c
c     7/14/92 (K. Brewster)
c     Added runname, comment and version number reading
C
c     8/20/92 (M. Xue)
c     Added reading of array z, and removed checking and copying of first 
c     data array. 
C
c#######################################################################
c
c     INPUT : 
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical
c
c       grdbas   Data read mode flag. 
c                If grdbas=1, only grid and base state arrays will be read 
c                If grdbas=0, all arrays will be read based on data 
c                             parameter setting.  
c       filename  Character variable naming the input HDF file
c
c     OUTPUT:
c
c       time     Time in seconds of data read from "filename"
c
c       x        x coordinate of grid points in physical/comp. space (km)
c       y        y coordinate of grid points in physical/comp. space (km)
c       z        z coordinate of grid points in computational space (km)
c       zp       z coordinate of grid points in physical space (km)
c
c       uprt     x component of perturbation velocity (m/s)
c       vprt     y component of perturbation velocity (m/s)
c       wprt     Vertical component of perturbation velocity in Cartesian
c                coordinates (m/s).
c       ptprt    Perturbation potential temperature (K)
c       pprt     Perturbation pressure (Pascal)
c
c       qvprt    Perturbation water vapor mixing ratio (kg/kg)
c       qc       Cloud water mixing ratio (kg/kg)
c       qr       Rainwater mixing ratio (kg/kg)
c       qi       Cloud ice mixing ratio (kg/kg)
c       qs       Snow mixing ratio (kg/kg)
c       qh       Hail mixing ratio (kg/kg)
c
c       km       Turbulent mixing coefficient (m**2/s)
c
c       ubar     Base state x velocity component (m/s)
c       vbar     Base state y velocity component (m/s)
c       wbar     Base state z velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state air density (kg/m**3)
c       qvbar    Base state water vapor mixing ratio (kg/kg)
c
c       ireturn  Return status indicator
c
c    WORK ARRAY
c
c       tem1
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nx,ny,nz
c
      integer grdbas
      character filename*(*)
c
      real x     (nx,ny,nz)     ! x coord.
      real y     (nx,ny,nz)     ! y coord.
      real z     (nx,ny,nz)     ! z coord.
      real zp    (nx,ny,nz)     ! physical x coord.
c
      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)
c
      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 km    (nx,ny,nz)     ! Turbulent mixing coefficient
c
      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
c
      real tem1  (nx,ny,nz)     ! Temporary work array

      integer ireturn
      real time
c
c#######################################################################
c
c     Parameters describing routine that wrote the gridded data
c
c#######################################################################
c
      character*40 fmtver,fmtverin
      character*10 tmunit
c
c#######################################################################
c
c     Misc. local variables
c
c#######################################################################
c
      integer lchanl
      parameter (lchanl=6)                      ! listing unit number

      integer isizes(3)
      integer i,j,k,iret,ndim,itime
      integer knt,kntscal
      integer kntin,kntscin
      integer nxin,nyin,nzin,nxyz
      integer bgrdin, bbasin, bvarin, bmstin, bicein, btrbin
      character label*256,unit*15,fmt*15,coord*15
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'indtflg.inc'
      include 'globcst.inc'
c
C#######################################################################
C
C     HDF Library routines called
C
C#######################################################################
c
      integer dsgdims,dsgdast,dsgdata
      external dsgdims
      external dsgdast
      external dsgdata
c
c#######################################################################
c
c     Other routines called:
c
c#######################################################################
c
      external kntary
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      nxyz=nx*ny*nz
c
C#######################################################################
c
c     Get dimensions of data in HDF file and check against
c        the dimensions passed to HDFREAD
c
C#######################################################################
c
      iret = dsgdims(filename,ndim,isizes,3)
      IF( iret.ne.0) THEN 
        write(6,'(a)') ' End of data reached. '
        ireturn = 1
        RETURN
      ENDIF

      nxin = isizes(1)
      nyin = isizes(2)
      nzin = isizes(3)
      print *, '  nxin= ',nxin,'  nyin= ',nyin,'  nzin= ',nzin

      IF( nxin.ne.nx .or. nyin.ne.ny .or. nzin.ne.nz ) THEN 
        write(6,'(1x,a)') 
     :       ' Dimensions in HDFREAD inconsistent with data.'
        write(6,'(1x,a,3I15)') ' In the data were: ', nxin, nyin, nzin
        write(6,'(1x,a,3I15)') ' In program  were: ', nx, ny, nz
        write(6,'(1x,a)') 
     :       ' Program aborted in HDFREAD.'
        STOP
      ENDIF
c
c#######################################################################
c
c     Read header info from the label for the first data set. 
c
c#######################################################################
c
      iret = dsgdast(label,unit,fmt,coord)

      IF( grdbas.eq.1 ) THEN ! Read grid and base state arrays

        read(label,'(5x,i6,a10,a21,6(i1),2(i2),a40,a80,a80)') 
     :      itime,tmunit,fmtverin,
     :      bgrdin,bbasin,bvarin,bmstin,bicein,btrbin,
     :      kntin,kntscin,runname,cmnt(1),cmnt(2)

        time=float(itime)
        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 kntary(fmtver,bgrdin,bbasin,bvarin,bmstin,bicein,btrbin,
     :              knt,kntscal)

      ELSE ! Normal data reading

        read(label,'(5x,i6,a10,a21,6(i1),2(i2),a40,a80,a80)') 
     :      itime,tmunit,fmtverin,
     :      grdin,basin,varin,mstin,icein,trbin,
     :      kntin,kntscin,runname,cmnt(1),cmnt(2)

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

        CALL kntary(fmtver,grdin,basin,varin,mstin,icein,trbin,
     :              knt,kntscal)
c
c     Check if fmtverin=fmtver, kntin=knt, kntscal=kntscin
c     If not, there is inconsistency in the data.
c
      ENDIF
c
      DO 10 i=1,2
        write(6,'(a,i2,a,a)') ' cmnt(',i,')= ',cmnt(i)
 10   CONTINUE
c
c#######################################################################
c
c     Establish the expected array count based on the switch settings
c
c#######################################################################
c
      CALL kntary(fmtver,grdin ,basin ,varin ,mstin ,icein ,trbin ,
     :            knt,kntscal)
c
C#######################################################################
c
c     Read in the coordinate arrays 
c
C#######################################################################
c
      IF( grdin.eq.1 .or. grdbas.eq.1 ) THEN

        iret = dsgdata(filename,3,isizes,x)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' x. '

        iret = dsgdata(filename,3,isizes,y)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' y. '

        iret = dsgdata(filename,3,isizes,z)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' z. '

        iret = dsgdata(filename,3,isizes,zp)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' zp.'

      END IF
c
C#######################################################################
c
c     Read in base state fields
c
C######################################################################
c
      IF( basin.eq.1 .or. grdbas.eq.1 ) THEN 

c
C#######################################################################
c
c     Read in u,v and w averaged to grid cell centers (scalar points).
c
C######################################################################
c
        iret = dsgdata(filename,3,isizes,ubar)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' ubar.'

        iret = dsgdata(filename,3,isizes,vbar)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' vbar.'

        iret = dsgdata(filename,3,isizes,wbar)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' wbar.'

        iret = dsgdata(filename,3,isizes,ptbar)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' ptbar.'

        iret = dsgdata(filename,3,isizes,pbar)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' pbar.'

        IF( mstin.eq.1) THEN 
          iret = dsgdata(filename,3,isizes,qvbar)
          iret = dsgdast(label,unit,fmt,coord)
          write(lchanl,910) label(1:6),' qvbar.'
        ENDIF

      ENDIF

      IF( grdbas.eq.1 ) GOTO 930 

      IF( varin.eq.1 ) THEN

        iret = dsgdata(filename,3,isizes,uprt)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' uprt.'

        iret = dsgdata(filename,3,isizes,vprt)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' vprt.'

        iret = dsgdata(filename,3,isizes,wprt)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' wprt.'

        iret = dsgdata(filename,3,isizes,ptprt)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' ptprt.'

        iret = dsgdata(filename,3,isizes,pprt)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' pprt.'

      END IF

      IF( mstin.eq.1 ) THEN

        iret = dsgdata(filename,3,isizes,qvprt)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' qvprt.'

        iret = dsgdata(filename,3,isizes,qc)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' qc.'

        iret = dsgdata(filename,3,isizes,qr)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' qr.'

        IF( icein.eq.1 ) THEN 

          iret = dsgdata(filename,3,isizes,qi)
          iret = dsgdast(label,unit,fmt,coord)
          write(lchanl,910) label(1:6),' qi.'

          iret = dsgdata(filename,3,isizes,qs)
          iret = dsgdast(label,unit,fmt,coord)
          write(lchanl,910) label(1:6),' qs.'

          iret = dsgdata(filename,3,isizes,qh)
          iret = dsgdast(label,unit,fmt,coord)
          write(lchanl,910) label(1:6),' qh.' 

        ENDIF  ! icein
      ENDIF  ! mstin

      IF( trbin.eq.1 ) THEN

        iret = dsgdata(filename,3,isizes,km)
        iret = dsgdast(label,unit,fmt,coord)
        write(lchanl,910) label(1:6),' km.'

      ENDIF

910   FORMAT(1x,'Field ',a,' was read in array ',a)
c
C#######################################################################
c
c     Friendly exit message
c
C######################################################################
c
930   CONTINUE

      write(lchanl,'(/a,F10.3,a/)')
     : ' Data at time=', time/60,' (min) were successfully read.'

      ireturn = 0

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE HDFSCRD                    ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c

      SUBROUTINE hdfscrd(nx,ny,nz,nxyz, x,y,z, tem1)
c
c#######################################################################
c
c     PURPOSE:
c
c     Convert HDF "scale" information into ARPS coordinate arrays.
c     This is sufficient when ARPS is on a rectangular grid, but for
c     an adaptive grid the grid should be written (and subsequently
c     read) as a data field.  That data field would then overwrite
c     the coordinates stored by means of this routine.
c
c     AUTHOR: Keith Brewster
c     7/14/92.    
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
      integer nx,ny,nz,nxyz
c
      real x(nx,ny,nz)
      real y(nx,ny,nz)
      real z(nx,ny,nz)
      real tem1(nxyz)
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i,j,k
      integer iret
      character*256 xlabel,ylabel,zlabel
      character*256 xunit,yunit,zunit
      character*256 xfmt,yfmt,zfmt
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      integer dsgdisc
      integer dsgdist
      external dsgdisc
      external dsgdist
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      iret = dsgdist(1,xlabel,xunit,xfmt)
      iret = dsgdist(2,ylabel,yunit,yfmt)
      iret = dsgdist(3,zlabel,zunit,zfmt)
c
      iret = dsgdisc(1,nx,tem1)
      DO 100 k=1,nz
      DO 100 j=1,ny
      DO 100 i=1,nx
        x(i,j,k) = tem1(i)
100   CONTINUE

      iret = dsgdisc(2,ny,tem1)
      DO 110 k=1,nz
      DO 110 j=1,ny
      DO 110 i=1,nx
        y(i,j,k) = tem1(j)
110   CONTINUE

      iret = dsgdisc(3,nz,tem1)
      DO 120 k=1,nz
      DO 120 j=1,ny
      DO 120 i=1,nx
        z(i,j,k) = tem1(k)
120   CONTINUE
c
      RETURN
      END
