c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE HDFDUMP                    ######
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 HDFDUMP(nx,ny,nz, filename , grdbas, stagrd,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :           ubar,vbar,ptbar,pbar,rhobar,qvbar,
     :           x,y,z,zp,hterain, j1,j2,j3, tem1)
c
c#######################################################################
c
c     PURPOSE:
c
c     Produces a history data file "filename" in NCSA HDF file format by 
c     calling HDF library subroutines. All output data are located at the 
c     grid cell centers. For example, u, v and w are all averaged to the 
c     cell center before being dumped. Note: HDF uses a special location 
c     for the coordinate fields which are stored as 1-dimensional vectors 
c     in the HDF "scale vectors". The coordinates are converted from meters 
c     to km. This is for compatibility with HDF-package applications like 
c     Xdataslice. If grdout is set to "2", the coordinate fields are also 
c     dumped as 3 dimensional data fields which have been converted from 
c     meters to kilometers.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     3/10/92.    
c
c     MODIFICATION HISTORY:
c
c     6/06/92 (M. Xue)
c     Added full documentation.
c
c     7/13/92 (K. Brewster)
c     Added runname, comment and version number writing
c     Restructured to dump coordinate info as an HDF scale field
c     instead of as an HDF data field.  General overhaul.
c
c     8/23/92 (M. Xue)
c     Modify to perform the dumping of both base and t-dependent arrays
c     and added control on grid staggering.
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       filename File name of history dump data.
c       grdbas   Flag indicating if this is a call for the data dump
c                of grid and base state arrays only. If so, grdbas=1.
c       stagrd   Indicating if the input variables are on a staggered
c                grid. If not (stgred=0), all variables are located at
c                the grid cell center, and the units of the coordinates
c                are in km insted of meter otherwise.
c
c       u        x component of velocity at a given time level (m/s)
c       v        y component of velocity at a given time level (m/s)  
c       w        Vertical component of Cartesian velocity at a given 
c                time level (m/s)   
c       ptprt    Perturbation potential temperature at a given time 
c                level (K)
c       pprt     Perturbation pressure at  a given time level (Pascal)
c       qv       Water vapor specific humidity at a given time level (kg/kg)
c       qc       Cloud water mixing ratio at a given time level (kg/kg)
c       qr       Rainwater mixing ratio at a given time level (kg/kg)
c       qi       Cloud ice mixing ratio at a given time level (kg/kg)
c       qs       Snow mixing ratio at a given time level (kg/kg)
c       qh       Hail mixing ratio at a given time level (kg/kg)
c
c       km       Turbulent mixing coefficient (m**2/s)
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       ptbar    Base state potential temperature (K)
c       pbar     Base state pressure (Pascal)
c       rhobar   Base state density (kg/m**3)
c       qvbar    Base state water vapor specific humidity (kg/kg)
c
c       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        z coordinate of grid points in computational space (m)
c       zp       Vertical coordinate of grid points in physical space (m)
c       hterain  The height of terrain (m)
c
c       j1       Coordinate transform Jacobian -d(zp)/dx
c       j2       Coordinate transform Jacobian -d(zp)/dy
c       j3       Coordinate transform Jacobian  d(zp)/dz
c
c     OUTPUT:
c
c       None.
c
c     WORK ARRAY:
c
c       tem1     work array.
c
c#######################################################################
c
c     The following parameters are passed into this subroutine through 
c     a common block in globcst.inc, and they control which 
c     variables are output.
c
c     grdout =0 or 1. If grdout=0, grid variables are not dumped. 
c     basout =0 or 1. If basout=0, base state variables are not dumped.
c     varout =0 or 1. If varout=0, model perturbation variables are not dumped. 
c     mstout =0 or 1. If mstout=0, water variables are not dumped. 
c     iceout =0 or 1. If iceout=0, dumping of qi, qs and qh is skipped.
c     trbout =0 or 1. If trbout=0, dumping of turbulence param km is skipped.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nx,ny,nz          ! Number of grid points in 3 directions

      character filename*(*)
      integer grdbas            ! If this is a grid/base state dump
      integer stagrd            ! If input arrays are on a staggered grid
c
      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)
c
      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)
c
      real km    (nx,ny,nz)     ! The turbulent mixing coefficient for
                                ! momentum. ( m**2/s )
c
      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)       ! The height of the terrain.

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

      real tem1  (nx,ny,nz)     ! Temporary work array
c
c#######################################################################
c
c     Parameters describing this routine
c
c#######################################################################
c
      character*40 fmtver
      parameter (fmtver='003.00 HDF Coded Data')
      character*10 tmunit
      parameter (tmunit='seconds   ')
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      real mtokm    ! conversion factor: meters to kilometers
                    ! used for scaling coordinate grid for output
      parameter ( mtokm=0.001)
      integer iret,itime,i,j,k
      integer n3dary,n3dsclr
      integer isizes(3)
      integer dssdims,dssdast,dsadata,dspdata
      character label*256,header*256
      integer ishf, jshf, kshf
      real scale
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external kntary
      external hdfgdmp
      external edgfill
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

c
c
c#######################################################################
c
c     Set the parameters depending on the input data attributes.
c
c     For stagrd=0 the data arrays are located at the C-grid cell center
c         stagrd=1 the data arrays are staggered on a C-grid.
c
c#######################################################################
c
      IF(stagrd.eq.1) THEN
        ishf=1
        jshf=1
        kshf=1
        scale=mtokm
      ELSE
        ishf=0
        jshf=0
        kshf=0
        scale=1.0
      END IF

      isizes(1) = nx
      isizes(2) = ny
      isizes(3) = nz

      iret = dssdims(3,isizes)

      write(6,'(1x,a,f13.3,a,a/)') 
     :     'Writing HDF data at time=', curtim,' into file ',filename  

c#######################################################################
c
c     Pack header info into the allowable character strings
c
c#######################################################################
c

c
c#######################################################################
c
c     Check the switch settings to see how many 3-D arrays and
c     3-D scalars will be dumped.
c
c#######################################################################
c
      IF( grdbas.eq.1 ) THEN

        CALL kntary(fmtver,1,1,0,mstout,0,0, n3dary,n3dsclr)

        itime = nint(curtim)
        write(header,'(5x,i6,a10,a21,6(i1),2(i2),a40,a80,a80)') 
     :                itime,tmunit,fmtver,
     :                1,1,0,mstout,0,0,
     :                n3dary,n3dsclr,runname,cmnt(1),cmnt(2)

      ELSE

        CALL kntary(fmtver,grdout,basout,varout,mstout,iceout,trbout,
     :              n3dary,n3dsclr)

        itime = nint(curtim)
        write(header,'(5x,i6,a10,a21,6(i1),2(i2),a40,a80,a80)') 
     :                itime,tmunit,fmtver,
     :                grdout,basout,varout,mstout,iceout,trbout,
     :                n3dary,n3dsclr,runname,cmnt(1),cmnt(2)

      ENDIF

c     write(label,'(a5,a251)') 'u    ',header(6:256)
c     iret = dssdast(label,'(m/s)','E15.9',' ')

c
c#######################################################################
c
c     Write grid coordinate scales. This is to be used by software such
c     as x-image that reads HDF files. The same information will be 
c     written in toerms of three-D arrays.
c
c#######################################################################
c
      CALL hdfgdmp(nx,ny,nz, x,y,z, scale, ishf, jshf, kshf, tem1)
c
      IF( grdout.eq.1  .or. grdbas.eq.1 ) THEN

        DO 100 k=1,nz-1
        DO 100 j=1,ny-1
        DO 100 i=1,nx-1
          tem1(i,j,k) = scale*(x(i)+x(i+ishf))*0.5
 100    CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        write(label,'(a5,a251)') 'x    ',header(6:256)
        iret = dssdast(label,'(km)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

        DO 110 k=1,nz-1
        DO 110 j=1,ny-1
        DO 110 i=1,nx-1
          tem1(i,j,k) = scale*(y(j)+y(j+jshf))*0.5
 110    CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        write(label,'(a5,a251)') 'y    ',header(6:256)
        iret = dssdast(label,'(km)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

        DO 115 k=1,nz-1
        DO 115 j=1,ny-1
        DO 115 i=1,nx-1
          tem1(i,j,k) = scale*(z(k)+z(k+kshf))*0.5
 115    CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        write(label,'(a5,a251)') 'z    ',header(6:256)
        iret = dssdast(label,'(km)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )
  
        DO 120 k=1,nz-1
        DO 120 j=1,ny-1
        DO 120 i=1,nx-1
          tem1(i,j,k) = mtokm*(zp(i,j,k)+zp(i,j,k+kshf))*0.5
 120    CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        write(label,'(a5,a251)') 'zp   ',header(6:256)
        iret = dssdast(label,'(km)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

      ENDIF    ! grdout
c
c#######################################################################
c
c     If basout=1, write out base state variables.
c
c#######################################################################
c
      IF( basout.eq.1  .or. grdbas.eq.1 ) THEN 

        DO 550 k=1,nz-1
        DO 550 j=1,ny-1
        DO 550 i=1,nx-1
          tem1(i,j,k) = (ubar(i,j,k)+ubar(i+ishf,j,k))*0.5
550     CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'ubar ',header(6:256)
        iret = dssdast(label,'(m/s)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

        DO 600 k=1,nz-1
        DO 600 j=1,ny-1
        DO 600 i=1,nx-1
          tem1(i,j,k) = (vbar(i,j,k)+vbar(i,j+jshf,k))*0.5
600     CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'vbar ',header(6:256)
        iret = dssdast(label,'(m/s)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

        DO 700 k=1,nz-1
        DO 700 j=1,ny-1
        DO 700 i=1,nx-1
          tem1(i,j,k) = 0.0 
700     CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'wbar ',header(6:256)
        iret = dssdast(label,'(m/s)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

        CALL edgfill(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'ptbar',header(6:256)
        iret = dssdast(label,'(degrees K)','E15.9',' ')
        iret = dsadata(filename,3,isizes, ptbar )

        CALL edgfill(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'pbar ',header(6:256)
        iret = dssdast(label,'(Pascals)','E15.9',' ')
        iret = dsadata(filename,3,isizes, pbar )

        IF( mstout.eq.1 ) THEN 

          CALL edgfill(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          write(label,'(a5,a251)') 'Qvbar',header(6:256)
          iret = dssdast(label,'(kg/kg)','E15.9',' ')
          iret = dsadata(filename,3,isizes, qvbar )

        ENDIF   ! mstout

      ENDIF   ! basout

      
      IF ( grdbas.eq.1 ) RETURN

c
c
c#######################################################################
c
c     If varout = 1, Write out uprt, vprt, wprt, ptprt, pprt.
c
c#######################################################################
c
      IF( varout.eq.1 ) THEN
c
c#######################################################################
c
c     Write out u,v and w averaged to grid cell centers (scalar points)
c
c#######################################################################
c
        DO 200 k=1,nz-1
        DO 200 j=1,ny-1
        DO 200 i=1,nx-1
          tem1(i,j,k)=(u(i,j,k)+u(i+ishf,j,k)
     :                -ubar(i,j,k)-ubar(i+ishf,j,k))*.5
200     CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        write(label,'(a5,a251)') 'uprt ',header(6:256)
        iret = dssdast(label,'(m/s)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

        DO 300 k=1,nz-1
        DO 300 j=1,ny-1
        DO 300 i=1,nx-1
          tem1(i,j,k)=(v(i,j,k)+v(i,j+jshf,k)
     :                -vbar(i,j,k)-vbar(i,j+jshf,k))*.5
300     CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        write(label,'(a5,a251)') 'vprt ',header(6:256)
        iret = dssdast(label,'(m/s)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )

        DO 400 k=1,nz-1
        DO 400 j=1,ny-1
        DO 400 i=1,nx-1
          tem1(i,j,k)=(w(i,j,k)+w(i,j,k+kshf))*0.5
400     CONTINUE
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)

        write(label,'(a5,a251)') 'wprt ',header(6:256)
        iret = dssdast(label,'(m/s)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1 )
c
c#######################################################################
c
c     Write out scalars 
c
c#######################################################################
c

        CALL edgfill(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'ptprt',header(6:256)
        iret = dssdast(label,'(degrees K)','E15.9',' ')
        iret = dsadata(filename,3,isizes, ptprt)

        CALL edgfill(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'pprt ',header(6:256)
        iret = dssdast(label,'(Pascals)','E15.9',' ')
        iret = dsadata(filename,3,isizes, pprt )

      ENDIF
c
      IF(mstout.eq.1) THEN 

        DO 500 k=1,nz-1
        DO 500 j=1,ny-1
        DO 500 i=1,nx-1
          tem1(i,j,k)=qv(i,j,k)-qvbar(i,j,k)
500     CONTINUE

        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'Qvprt',header(6:256)
        iret = dssdast(label,'(kg/kg)','E15.9',' ')
        iret = dsadata(filename,3,isizes, tem1)
c
        CALL edgfill(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'Qc   ',header(6:256)
        iret = dssdast(label,'(kg/kg)','E15.9',' ')
        iret = dsadata(filename,3,isizes, qc )

        CALL edgfill(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'Qr   ',header(6:256)
        iret = dssdast(label,'(kg/kg)','E15.9',' ')
        iret = dsadata(filename,3,isizes, qr )

        IF(iceout.eq.1) THEN 
c
          CALL edgfill(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          write(label,'(a5,a251)') 'Qi   ',header(6:256)
          iret = dssdast(label,'(kg/kg)','E15.9',' ')
          iret = dsadata(filename,3,isizes, qi )

          CALL edgfill(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          write(label,'(a5,a251)') 'Qs   ',header(6:256)
          iret = dssdast(label,'(kg/kg)','E15.9',' ')
          iret = dsadata(filename,3,isizes, qs )

          CALL edgfill(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          write(label,'(a5,a251)') 'Qh   ',header(6:256)
          iret = dssdast(label,'(kg/kg)','E15.9',' ')
          iret = dsadata(filename,3,isizes, qh )

        ENDIF

      ENDIF

      IF( trbout.eq.1 ) THEN

        CALL edgfill(km,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        write(label,'(a5,a251)') 'Km   ',header(6:256)
        iret = dssdast(label,'(m**2/s)','E15.9',' ')
        iret = dsadata(filename,3,isizes, km )

      ENDIF   ! trbout

      RETURN
      END     
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE HDFGDMP                    ######
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 HDFGDMP(nx,ny,nz, x,y,z, scale,ishf,jshf,kshf, tem1)
c
c#######################################################################
c
c     PURPOSE:
c
c     Set up one-dimensional arrays for storage of the coordinate values as 
c     HDF scale vectors.  Call appropriate HDF routines to write the scale 
c     vectors. Note: HDF can only handle coordinates on a regular grid.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/21/92.    
c
c     MODIFICATION HISTORY:
c       8/23/92 (M. Xue)
c       Arguments scale,ishf,jshf,kshf added.
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       x        x coordinate of grid points in physical/comp. space (m)
c       y        y coordinate of grid points in physical/comp. space (m)
c       z        Vertical coordinate of grid points in computational 
c                space (m)
c       scale    Multiplication factor for x, y and z coordinates
c       ishf     I-index shift parameter for coordinate averaging
c       jshf     J-index shift parameter for coordinate averaging
c       kshf     K-index shift parameter for coordinate averaging
c
c     OUTPUT:
c
c       None.
c
c     WORK ARRAY:
c
c       tem1     work array, declared herein as one-dimensional
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
      integer nx,ny,nz
      real x(nx)
      real y(ny)
      real z(nz)
      real scale 
      integer ishf, jshf, kshf
      real tem1(nx*ny*nz)
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i,j,k
      integer iret
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      integer dssdisc
      integer dssdist
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      iret = dssdist(1, 'x','(km)','F10.3')
      iret = dssdist(2, 'y','(km)','F10.3')
      iret = dssdist(3, 'z','(km)','F10.3')
c
      DO 100 i=1,nx-1
        tem1(i) = (x(i)+x(i+ishf))*scale
100   CONTINUE
      iret = dssdisc(1, nx ,tem1)

      DO 110 j=1,ny-1
        tem1(j) = (y(j)+y(j+jshf))*scale
110   CONTINUE
      iret = dssdisc(2, ny ,tem1)
c
c#######################################################################
c
c     We use the computational coordinate z as the grid coordinate scale
c     in the vertical direction.
c
c#######################################################################
c
      DO 120 k=1,nz-1
        tem1(k) = (z(k)+z(k+kshf))*scale
120   CONTINUE
      iret = dssdisc(3, nz ,tem1)
c
      RETURN
      END
