c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE PAKDUMP                    ######
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 PAKDUMP(nx,ny,nz, nchanl, 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, intem)
c
c#######################################################################
c
c     IMPORTANT NOTICE:
c
c     This version is designed for the IBM RISC/6000 and VAX computers
c     (that write IEEE standard integer*4 binary files).  
c     The files are generally not portable between different machine types.
c
c     PURPOSE:
c
c     Write history data into direct access channel nchanl,
c     as packed integers.
c
c     All output data are located at the grid cell centers, therefore
c     u, v and w are averaged to the cell center before being dumped. 
c
c     Note: coordinate fields are dumped as 3 dimensional fields which
c     have been converted from meters to kilometers.  The coordinate
c     values are also averaged from the u position on the staggered grid
c     to the scalar position.  This is done for the convenience of the 
c     plotting applications.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     6/25/92.    
c
c     MODIFICATION HISTORY:
c
c     7/14/92 (K. Brewster)
c     Added runname, comment and version number writing
c
c     7/25/92 (K. Brewster)
c     Replaced use of integer*2 with packing into integer*4's
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       nchanl   Fortran I/O channel number for history data output.
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 ARRAYS:
c
c       tem1       work array.
c       intem      integer work array.
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c
      integer nchanl            ! Fortran I/O channel number for output
      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
      integer intem (nx*ny*nz)  ! Temporary integer work array
c
c#######################################################################
c
c     Parameters describing this routine
c
c#######################################################################
c
      character*40 fmtver
      parameter (fmtver='003.00 Packed Binary 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 ihdlen,lablen
      parameter (ihdlen=5000, lablen=40)
      integer ihead(ihdlen)
      integer label(lablen)
      integer i,j,k,l
      integer nxyz
      integer ilen,illen
      integer n3dary,n3dsclr
      real amax,amin,ref
      integer iscale,nbits,ihend,iend
      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 mkhead
      external a3dmax
      external edgfill
      external packdat
      external mklabel
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

      nxyz = nx*ny*nz
      iend=(nxyz+1)/2

      write(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim 
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)

      ELSE

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

      ENDIF

c
c#######################################################################
c
c     Write header info
c
c#######################################################################
c
      CALL mkhead(nx,ny,nz,ihdlen,ilen,
     :            fmtver,tmunit,n3dary,n3dsclr,ihead)
      write(nchanl) ihead
c
c#######################################################################
c
c     If grdout=1, write out grid variables
c
c#######################################################################
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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'x coordinate',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'y coordinate',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

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

        CALL a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'z coordinate',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        DO 130 k=1,nz-1
        DO 130 j=1,ny-1
        DO 130 i=1,nx-1
          tem1(i,j,k) = scale*(zp(i,j,k)+zp(i,j,k+kshf))*0.5
 130    CONTINUE

        CALL a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'zp coordinat',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

       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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'ubar        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'vbar        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'wbar        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        CALL a3dmax0(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :               amax,amin)
        CALL edgfill(ptbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, ptbar, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'ptbar       ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        CALL a3dmax0(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :               amax,amin)
        CALL edgfill(pbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, pbar, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'pbar        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        IF(mstout.eq.1) THEN 

          CALL a3dmax0(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :                 amax,amin)
          CALL edgfill(qvbar,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          CALL packdat(nxyz, qvbar, intem,amin,amax,ref,iscale,nbits)
          CALL mklabel(lablen,illen,
     :               'qvbar       ',amin,amax,ref,iscale,nbits,label)
          write(nchanl) label
          write(nchanl) (intem(i),i=1,iend)

        ENDIF
      ENDIF

      
      IF ( grdbas.eq.1 ) RETURN

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))*0.5  
200     CONTINUE

        CALL a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'uprt        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'vprt        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'wprt        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)
c
c#######################################################################
c
c     Write out scalars 
c
c#######################################################################
c
        CALL a3dmax0(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,
     :               amax,amin)
        CALL edgfill(ptprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, ptprt, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'ptprt       ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend) 

        CALL a3dmax0(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(pprt,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, pprt, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'pprt        ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

      ENDIF
c
c#######################################################################
c
c     If mstout = 1, write out moisture scalars 
c
c#######################################################################
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 a3dmax0(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(tem1,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, tem1, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'qvprt       ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        CALL a3dmax0(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(qc,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, qc, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'qc          ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)
c
        CALL a3dmax0(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(qr,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, qr, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'qr          ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

        IF(iceout.eq.1) THEN 
c
          CALL a3dmax0(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
          CALL edgfill(qi,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          CALL packdat(nxyz, qi, intem,amin,amax,ref,iscale,nbits)
          CALL mklabel(lablen,illen,
     :               'qi          ',amin,amax,ref,iscale,nbits,label)
          write(nchanl) label
          write(nchanl) (intem(i),i=1,iend)

          CALL a3dmax0(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
          CALL edgfill(qs,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          CALL packdat(nxyz, qs, intem,amin,amax,ref,iscale,nbits)
          CALL mklabel(lablen,illen,
     :               'qs          ',amin,amax,ref,iscale,nbits,label)
          write(nchanl) label
          write(nchanl) (intem(i),i=1,iend)
c
          CALL a3dmax0(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
          CALL edgfill(qh,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
          CALL packdat(nxyz, qh, intem,amin,amax,ref,iscale,nbits)
          CALL mklabel(lablen,illen,
     :               'qh          ',amin,amax,ref,iscale,nbits,label)
          write(nchanl) label
          write(nchanl) (intem(i),i=1,iend)

        ENDIF  ! iceout
      ENDIF  ! mstout

      IF( trbout.eq.1 ) THEN

        CALL a3dmax0(km,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1,amax,amin)
        CALL edgfill(km,1,nx,1,nx-1,1,ny,1,ny-1,1,nz,1,nz-1)
        CALL packdat(nxyz, km, intem,amin,amax,ref,iscale,nbits)
        CALL mklabel(lablen,illen,
     :               'km          ',amin,amax,ref,iscale,nbits,label)

        write(nchanl) label
        write(nchanl) (intem(i),i=1,iend)

      ENDIF   ! trbout
c
      RETURN
      END     
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE PACKDAT                    ######
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 PACKDAT(nxyz, a, ivec,amin,amax,ref,iscale,nbits)
c
c#######################################################################
c
c     IMPORTANT NOTICE:
c
c     This version is designed for the IBM RISC/6000 and VAX computers
c     (that write IEEE standard integer*4 binary files).  
c     The files are generally not portable between different machine types.
c
c     PURPOSE:
c
c     Two to one packing of real data vector "a" into integer 
c     vector "ivec".
c
c     Structure of integer element ivec(k=i/2) for integer*4 ivec:
c      
c     |-| | | | | | | | | | | | | | | |1| | | | | | | | | | | | | | | |
c      1          15 bits              1              15 bits
c      s|<-----ifix(scaled a(i-1))--->|s|<-----ifix(scaled a(i))----->|
c      i                               i
c      g                               g
c      n of (a(i-1)-ref)               n of (a(i)-ref), 1 if neg, 0 if positive
c
c#######################################################################
c
c     AUTHOR: Keith Brewster, based on Rich Carpenter's
c                             subroutine GRIBER
c     6/25/92.    
c
c     MODIFICATION HISTORY:
c
c     8/28/92.   (Keith Brewster)
c     Beefed up documentation.
c
c#######################################################################
c
c     INPUT:
c
c       nxyz     Length of data vector
c
c       a        Real data vector
c
c       amax     Maximum value within vector a
c       amin     Minimum value within vector a
c
c     OUTPUT:
c
c       ivec     Packed data vector
c
c       ref      Reference value which is subtracted from all the
c                elements in vector "a" in the packing process.
c
c       iscale   Scaling constant i2=(a-ref)*2**(-iscale)
c
c       nbits    Number of bits in the output integer
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nxyz
c
      real a(nxyz)
      integer ivec(nxyz)
c
      real amax
      real amin
c
      real ref
      integer iscale
      integer nbits
c
c#######################################################################
c
c     Safety check variables
c     These variables are only used for the verification of the 
c     packing result.  If this "safety check" is removed to improve
c     efficiency, these variables are not needed.
c
c#######################################################################
c
      real delmx,delta,delta1,delta2,rf,recx1,recx2
      integer irec,int1,int2,isgn2
      logical warned
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer kbits
      parameter (kbits=16)
      integer ishift,isgnmk,ihfsgn
      parameter (ishift=2**(kbits),
     :           isgnmk=2**(kbits-1),
     :           ihfsgn=2**(kbits-2))
      integer i,k,iend
      integer intem1,intem2,intem
      real range,rscale,rslog
      real d1,d2,factor
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Initialize safety check variables
c
c#######################################################################
c
      delmx=0.
      warned=.false.
c
c#######################################################################
c
c     Set constants and scale parameters
c
c#######################################################################
c
      nbits=kbits
      iend=(nxyz+1)/2
      range=amax-amin
      ref=0.5*(amax+amin)
c
      IF( range.NE.0. ) THEN
        rscale=range/((2.0**kbits)-4.0)
        rslog=log(rscale)/log(2.)
        IF( rslog.gt.0. ) rslog=rslog+1.0
        iscale=rslog
        factor=2.0**(     -iscale)
      ELSE
        iscale=0
        factor=1.
      ENDIF
c
      rf = 2.0**iscale
c
c#######################################################################
c
c     Subtract the reference,
c     scale the result, and transfer to integer array
c
c#######################################################################
c
c
      DO 20 i=2,nxyz,2
        k=i/2
c
        d1 = a(i-1)-ref
        intem1 = NINT(factor*d1)
        intem1 = ishift * intem1
        d2 = a(  i)-ref
        intem2 = NINT(factor*d2)
        intem  = abs(intem1) + abs(intem2)
        intem  = intem+(1-isign(1,intem2))*ihfsgn
        intem  = isign(intem,intem1)
        ivec(k) = intem
c
c#######################################################################
c
c     Temporary safety check -- see if the variables can
c     be recovered successfully.
c
c#######################################################################
c
        irec=intem/ishift
        int1=irec*ishift
        recx1 = ref + (rf*irec)
c
        int2=abs(intem) - abs(int1)
        isgn2=int2/isgnmk
        int2=int2-(isgn2*isgnmk)
        irec=(1-2*isgn2)*int2
        recx2 = ref + (rf*irec)
c
        delta1=abs(a(i-1)-recx1)
        delmx=amax1(delmx,delta1)
        delta2=abs(  a(i)-recx2)
c
        delmx=amax1(delmx,delta2)

        IF( .not.warned .and. (delta1.gt.10. .or. delta2.gt.10.))THEN
          print *, ' delta1= ',delta1,'  delta2= ',delta2
          print *, ' ref = ',ref
          print *, '  rscale = ',rscale
          print *, '   rslog = ',rslog
          print *, '  iscale = ',iscale
          print *, '  factor = ',factor
          print *, ' i= ',i-1,i,'  k= ',k
          print *, ' a(i-1)= ',a(i-1),' intem1= ',intem1
          print *, ' d1 = ',d1
          print *, '   a(i)= ',  a(i),' intem2= ',intem2
          print *, ' d2 = ',d2
          print *, ' intem = ',intem
          warned=.true.
        ENDIF
c
c#######################################################################
c
c     End of temporary safety check 
c
c#######################################################################
c
 20   CONTINUE
c
c#######################################################################
c
c     If nxyz is odd, there is one last variable to convert.
c     Do it.
c
c#######################################################################
c
      IF(k.lt.iend) THEN
        intem = NINT(factor*(a(nxyz)-ref))
        intem = ishift * intem
        ivec(iend) = intem
c
c#######################################################################
c
c     Temporary safety check -- see if the last variable can
c     be recovered successfully.
c
c#######################################################################
c

        irec=intem/ishift
        recx1 = ref + (rf*irec)
c
        delta=abs(a(nxyz)-recx1)
        delmx=amax1(delmx,delta)
c
c#######################################################################
c
c     End of temporary safety check 
c
c#######################################################################
c
      END IF
      write(6,'(a,E14.6,a,E14.6)')
     :        ' Field min= ',amin,'   Field max= ',amax
      write(6,'(a,E14.6)')
     :        ' Maximum difference due to compression= ',delmx
      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE MKHEAD                     ######
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 MKHEAD(nx,ny,nz,ihdlen,ilen,
     :                  fmtver,tmunit,n3dary,n3dsclr,ihead)
c
c#######################################################################
c
c     PURPOSE:
c
c     Take information needed for packed data file header and convert
c     into integer*4 vector which will be written as the first
c     record in the packed data file.
c
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/25/92.    
c
c     MODIFICATION HISTORY:
c
c     8/28/92 (K. Brewster)
c     Added full documentation.
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       ihdlen   Length of header vector.
c
c       fmtver   Character string containing description of data format
c                and version number.
c
c       tmunit   Character string containing description of time units.
c
c       n3dary   Number of 3d arrays of data to be dumped.
c
c       n3dsclr  Number of 3d scalar arrays to be dumped.
c
c     OUTPUT:
c
c       ilen     Number of elements of header vector filled with info.
c
c       ihead    Integer vector containing the header information
c                translated into integers.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz
      integer ihdlen,ilen
      character*40 fmtver
      character*10 tmunit
      integer n3dary,n3dsclr
      integer ihead(ihdlen)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer iclen,jlen,i,l
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c
      include 'globcst.inc'
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external trnchar
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
c     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      ilen=len(fmtver)

      CALL trnchar(fmtver,1,ilen,ihead)

      jlen=len(runname)

      CALL trnchar(runname,(ilen+1),jlen,ihead)

      ilen=ilen+jlen
      ihead(ilen+1) = nocmnt
      ilen=ilen+1
      jlen=len(cmnt(1))

      IF( nocmnt.gt.0 ) THEN

        iclen=ilen
        DO 10 l=1,nocmnt

          CALL trnchar(cmnt(l),(iclen+1),jlen,ihead)
          iclen=iclen+jlen

 10     CONTINUE

      ENDIF

      ilen=ilen+(50*jlen)
      ihead(ilen+1) = curtim
      ilen=ilen+1
      jlen=len(tmunit)

      CALL trnchar(runname,(ilen+1),jlen,ihead)

      ilen=ilen+jlen
      ihead(ilen+1) = grdout
      ihead(ilen+2) = basout
      ihead(ilen+3) = varout
      ihead(ilen+4) = mstout
      ihead(ilen+5) = iceout
      ihead(ilen+6) = trbout
      ihead(ilen+7) = n3dary
      ihead(ilen+8) = n3dsclr
      ihead(ilen+9) = nx
      ihead(ilen+10) = ny
      ihead(ilen+11) = nz
      ilen=ilen+11

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE MKLABEL                    ######
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 MKLABEL(lablen,illen,
     :               chlabel,amin,amax,ref,iscale,nbits,label)
c
c#######################################################################
c
c     PURPOSE:
c
c     Take information needed for data label and convert
c     into integer*4 vector which will be written in the data
c     file just before the data array.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/25/92.    
c
c     MODIFICATION HISTORY:
c
c     8/28/92 (K. Brewster)
c     Added full documentation.
c
c#######################################################################
c
c     INPUT:
c
c       lablen    Length of the label vector.
c       chlabel   Character string decribing the variable in the array
c                 about to be dumped.
c
c                 **The values of the following input variables*******
c                 **must be the same as those used to pack the data***
c
c       amin      Minimum value of the array about to be dumped.
c       amax      Maximum value of array about to be dumped.
c       ref       Reference value used in the data packing.
c       iscale    Base-2 exponent used in the data packing.
c       nbits     Number of bits used to compress the data.
c
c     OUTPUT:
c
c       illen     Number of elements of the label vector filled with info.
c       label     Vector containing translated label info.
c
c#######################################################################
c
c     Variable declarations
c
c#######################################################################
c
      implicit none
c
      integer lablen,illen
      character*12 chlabel
      real amin,amax,ref
      integer iscale,nbits
      integer label(lablen)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer ilen
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external trnchar
      external trnreal
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      ilen=len(chlabel)

      CALL trnchar(chlabel,1,ilen,label)
      CALL trnreal(amin,label(ilen+1),label(ilen+2))
      CALL trnreal(amax,label(ilen+3),label(ilen+4))
      CALL trnreal( ref,label(ilen+5),label(ilen+6))

      label(ilen+7)=iscale
      label(ilen+8)=nbits
      illen=ilen+8
c
      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TRNCHAR                    ######
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 TRNCHAR(string,istrt,ilen,ivect)
c
c
c#######################################################################
c
c     PURPOSE:
c
c     Translates a character string into an integer vector.
c     Each character in the string is represented by its ASCII value.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/25/92.    
c
c     MODIFICATION HISTORY:
c
c     8/28/92 (K. Brewster)
c     Added full documentation.
c
c#######################################################################
c
c     INPUT:
c
c       string     Character string to be translated.
c       istrt      Position in the integer vector to place the first
c                  translated character of the character string.
c       ilen       Length of the character string.
c
c     OUTPUT:
c
c       ivect      Integer vector to be filled with translated string.
c
c#######################################################################
c
c     Variable declarations
c
c#######################################################################
c
      implicit none
      integer istrt,ilen
      character string*(*)
      integer ivect(*)
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i,ism1,iloc
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c

      ism1=istrt-1

      DO 100 i=1,ilen
        iloc=ism1 + i
        ivect(iloc)=ichar(string(i:i))
 100  CONTINUE

      RETURN
      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE TRNREAL                    ######
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 TRNREAL(rvar,int1,int2)
c
c#######################################################################
c
c     PURPOSE:
c
c     Translates a real variable, rvar, into two integers, int1 and int2.
c     The real variable can be expressed as    rvar=int2 * 2**int1.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/25/92.    
c
c     MODIFICATION HISTORY:
c
c     8/28/92 (K. Brewster)
c     Added full documentation.
c
c#######################################################################
c
c     INPUT:
c
c       rvar     Real number to be translated.
c
c     OUTPUT:
c
c       int1     Integer describing the base-2 exponent of rvar.
c       int2     Integer describing the mantissa of rvar.
c
c#######################################################################
c
      implicit none
c
      real rvar
      integer int1,int2
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer kbits
      parameter (kbits=32)
      real rscale,rslog,factor,rtemp
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF( rvar.ne.0. ) THEN
        rscale=rvar/((2.0**(kbits-2))-2.0)
        rslog=log(abs(rscale))/log(2.)
        IF( rslog.gt.0. ) rslog=rslog+1.0
        int1=rslog
        factor=2.0**(-int1)
        rtemp=factor*rvar
        int2=nint(rtemp)
      ELSE
        int1=0
        int2=0        
      ENDIF

      RETURN
      END
