c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                 SUBROUTINE DTADUMP                   ######
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 DTADUMP(nx,ny,nz, houtfmt,nchout,filnam,grdbas,stagrd,
     :           u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,km,
     :           ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,
     :           x,y,z,zp,hterain, j1,j2,j3, tem1,tem2,svitem)
c
c#######################################################################
c
c     PURPOSE:
c
c     Coordinate the history data dump in various data formats by 
c     calling the appropriate history dump routine. 
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     8/23/1992
c
c     MODIFICATION HISTORY:
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       houtfmt  The format flag of history datra dump 
c       nchout   FORTRAN I/O channel number for history data output.
c       filnam   The name of history data dump file
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 x-velocity component (m/s)
c       vbar     Base state y-velocity component (m/s)
c       wbar     Base state vertical 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       tem2     work array.
c       svitem   work array.
c
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 houtfmt           ! Flag of history data dump foramt 
      integer nchout            ! FORTRAN I/O channel number for output
      character filnam*(*)      ! The name of history dump data file
      integer grdbas            ! If this is a grid/base state array 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 x-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state y-velocity (m/s)
      real wbar  (nx,ny,nz)     ! Base state z-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
      real tem2  (nx,ny,nz)     ! Temporary work array
      real svitem(3*nx*ny*nz)   ! Temporary work array

      integer istat, nchout0
c
c#######################################################################
c
c     Routines called:
c
c#######################################################################
c
      external bindump
      external ascdump
      external hdfdump
      external pakdump
      external svidump
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
      return
      IF( houtfmt.eq.0 ) RETURN   ! No data dump 

      IF( houtfmt.eq.1 ) THEN     ! Unformatted binary data dump.
c
c#######################################################################
c
c     History dump of unformatted binary data
c
c#######################################################################
c

        OPEN(unit=nchout,file=filnam,status='new',
     :       form='unformatted',iostat= istat )

        IF( istat.ne.0) GOTO 999

        CALL bindump(nx,ny,nz, nchout , 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)

        CLOSE(unit=nchout)

      ELSEIF( houtfmt.eq.2 ) THEN     ! Formatted ASCII data dump
c
c#######################################################################
c
c     History dump of formatted ASCII data
c
c#######################################################################
c

        OPEN(unit=nchout,file=filnam,status='new',
     :       form='formatted',iostat= istat )

        IF( istat.ne.0) GOTO 999

        CALL ascdump(nx,ny,nz, nchout , 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)

        CLOSE(unit=nchout)

      ELSEIF( houtfmt .eq. 3 ) THEN    ! HDF format dump.
c
c#######################################################################
c
c     History data dump in NCSA HDF file format.
c
c#######################################################################
c

        CALL hdfdump(nx,ny,nz, filnam, 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)

      ELSE IF( houtfmt.eq.4 ) THEN     ! Packed binary data dump.
c
c#######################################################################
c
c     History dump of packed unformatted binary data
c
c#######################################################################
c

        OPEN(unit=nchout,file=filnam,status='new',
     :       form='unformatted',iostat= istat )

        IF( istat.ne.0) GOTO 999

        CALL pakdump(nx,ny,nz, nchout , 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, tem2)

        CLOSE(unit=nchout)

      ELSE IF( houtfmt.eq.5 ) THEN     ! Data dump for Savi3D
c
c#######################################################################
c
c     History dump for Savi3D.
c     Data at all time levels are dumped into a single file.
c
c#######################################################################
c
        nchout0 = 81

        CALL svidump(nx,ny,nz, nchout0,filnam, grdbas, stagrd,
     :       u,v,w,ptprt,pprt,qv,qc,qr,qi,qs,qh,
     :       km,ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar,
     :       x,y,z,zp,hterain, j1,j2,j3, tem1,svitem)

      ENDIF   ! houtfmt switches

      RETURN

999   CONTINUE
      write(6,'(1x,a,a,/1x,i3,a)')
     :  'Error occured when opening file ',filnam,
     :  'using FORTRAN unit ',nchout,' Program stopped in DTADUMP.'
      STOP

      END
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE BINDUMP                    ######
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 BINDUMP(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)
c
c#######################################################################
c
c     PURPOSE:
c
c     Write history data into channel nchanl as formatted binary data.
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.  This is for the
c     convenience of the plotting applications.
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
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 ARRAY:
c
c       tem1     work array.
c
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
c
      integer nchanl            ! FORTRAN I/O channel number for output
      integer grdbas            ! If this is a grid/base state array 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 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 i,j,k,l
      integer n3dary,n3dsclr
      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
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
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
      ENDIF

      write(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim 
c
c#######################################################################
c
c     Write header info
c
c#######################################################################
c
      write(nchanl) fmtver
      write(nchanl) runname
      write(nchanl) nocmnt
      IF( nocmnt.gt.0 ) THEN
        DO 10 l=1,nocmnt
          write(nchanl) cmnt(l)
 10     CONTINUE
      ENDIF
      write(nchanl) curtim,tmunit
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)

        write(nchanl) 1,1,0,mstout,0,0

      ELSE

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

        write(nchanl) grdout,basout,varout,mstout,iceout,trbout

      ENDIF

      write(nchanl) n3dary
      write(nchanl) n3dsclr

      write(nchanl) nx,ny,nz   
c
c#######################################################################
c
c     If grdout=1 or grdbas=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
        write(nchanl) 'x coordinate'
        write(nchanl) 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
        write(nchanl) 'y coordinate'
        write(nchanl) tem1

        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
        write(nchanl) 'z coordinate'
        write(nchanl) tem1

        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
        write(nchanl) 'zp coord    '
        write(nchanl) 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
        write(nchanl) 'ubar        '
        write(nchanl)  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
        write(nchanl) 'vbar        '
        write(nchanl)  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
        write(nchanl) 'wbar        '
        write(nchanl)  tem1 

        write(nchanl) 'ptbar       '
        write(nchanl)  ptbar

        write(nchanl) 'pbar        '
        write(nchanl)  pbar 

        IF(mstout.eq.1) THEN 

          write(nchanl) 'qvbar       '
          write(nchanl)  qvbar

        ENDIF

      ENDIF

      
      IF ( grdbas.eq.1 ) RETURN

c
c#######################################################################
c
c     If varout = 1, Write out uprt, vprt, wprt, ptprt, pprt.
c
c#######################################################################
c
c
c#######################################################################
c
c     Write out u,v and w averaged to grid cell centers (scalar points)
c
c#######################################################################
c
      IF(varout.eq.1) THEN
        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
        write(nchanl) 'uprt        '
        write(nchanl) 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
        write(nchanl) 'vprt        '
        write(nchanl) 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
        write(nchanl) 'wprt        '
        write(nchanl) tem1 
c
c#######################################################################
c
c     Write out scalars 
c
c#######################################################################
c
        write(nchanl) 'ptprt       '
        write(nchanl) ptprt

        write(nchanl) 'pprt        '
        write(nchanl) pprt 

      ENDIF     ! varout
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

        write(nchanl) 'qvprt       '
        write(nchanl)  tem1

        write(nchanl) 'qc          '
        write(nchanl)  qc 
c
        write(nchanl) 'qr          '
        write(nchanl)  qr

        IF(iceout.eq.1) THEN 
c
          write(nchanl) 'qi          '
          write(nchanl)  qi

          write(nchanl) 'qs          '
          write(nchanl)  qs 
c
          write(nchanl) 'qh          '
          write(nchanl)  qh

        ENDIF   !iceout

      ENDIF   !mstout
c
c#######################################################################
c
c     If trbout = 1, write out turbulence parameter, km.
c
c#######################################################################
c
      IF( trbout.eq.1 ) THEN

        write(nchanl) 'km           '
        write(nchanl) km

      ENDIF   ! trbout

      RETURN
      END     
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE ASCDUMP                    ######
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 ASCDUMP(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)
c
c#######################################################################
c
c     PURPOSE:
c
c     Write history data into channel nchanl as formatted ASCII data
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.  This is for the
c     convenience of the plotting applications.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     5/31/92.    
c
c     MODIFICATION HISTORY:
c
c     7/13/92 (K. Brewster)
c     Added runname, comment and version number writing
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 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
c
      integer nchanl            ! FORTRAN I/O channel number for output
      integer grdbas            ! If this is a grid/base state array 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 ASCII Formatted 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 i,j,k,l
      integer n3dary,n3dsclr
      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
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
      ENDIF

      write(6,'(1x,a,f13.3/)') 'Writing history data at time=', curtim 
c
      write(nchanl,901) fmtver
      write(nchanl,901) runname
c
      write(nchanl,'(1x,i4)') nocmnt
      IF(nocmnt.gt.0) THEN
        DO 10 l=1,nocmnt
          write(nchanl,901) cmnt(l)
 10     CONTINUE
      ENDIF
c
      write(nchanl,'(1x,e16.8,1x,a10)') curtim,tmunit
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)

        write(nchanl,'(6(1x,i4))') 1,1,0,mstout,0,0 

      ELSE

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

        write(nchanl,'(6(1x,i4))') 
     :        grdout,basout,varout,mstout,iceout,trbout

      ENDIF

      write(nchanl,'(1x,i12)')   n3dary
      write(nchanl,'(1x,i12)')   n3dsclr

      write(nchanl,'(1x,3i12 )') nx,ny,nz   
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
        write(nchanl,'(1x,a)') 'x coordinate'
        write(nchanl,900) 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
        write(nchanl,'(1x,a)') 'y coordinate'
        write(nchanl,900) tem1

        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
        write(nchanl,'(1x,a)') 'z coordinate'
        write(nchanl,900) tem1

        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
        write(nchanl,'(1x,a)') 'zp coordinat'
        write(nchanl,900) tem1

      ENDIF
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
        write(nchanl,'(1x,a)') 'ubar        '
        write(nchanl,900) 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
        write(nchanl,'(1x,a)') 'vbar        '
        write(nchanl,900) 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
        write(nchanl,'(1x,a)') 'wbar        '
        write(nchanl,900) tem1 

        write(nchanl,'(1x,a)') 'ptbar       '
        write(nchanl,900) ptbar

        write(nchanl,'(1x,a)') 'pbar        '
        write(nchanl,900) pbar 

        IF(mstout.eq.1) THEN 

          write(nchanl,'(1x,a)') 'qvbar       '
          write(nchanl,900) qvbar

        ENDIF

      ENDIF

      
      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
        write(nchanl,'(1x,a)') 'uprt        '
        write(nchanl,900) 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
        write(nchanl,'(1x,a)') 'vprt        '
        write(nchanl,900) 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
        write(nchanl,'(1x,a)') 'wprt        '
        write(nchanl,900) tem1 
c
c#######################################################################
c
c       Write out scalars 
c
c#######################################################################
c

        write(nchanl,'(1x,a)') 'ptprt       '
        write(nchanl,900) ptprt

        write(nchanl,'(1x,a)') 'pprt        '
        write(nchanl,900) pprt 
      ENDIF
c
c#######################################################################
c
c     If mstout = 1, Write out moisture variables
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
        write(nchanl,'(1x,a)') 'qvprt       '
        write(nchanl,900) tem1

        write(nchanl,'(1x,a)') 'qc          '
        write(nchanl,900) qc 
c
        write(nchanl,'(1x,a)') 'qr          '
        write(nchanl,900) qr

        IF( iceout.eq.1 ) THEN 
c
          write(nchanl,'(1x,a)') 'qi          '
          write(nchanl,900) qi

          write(nchanl,'(1x,a)') 'qs          '
          write(nchanl,900) qs 
c
          write(nchanl,'(1x,a)') 'qh          '
          write(nchanl,900) qh

        ENDIF    ! iceout

      ENDIF  ! mstout

      IF( trbout.eq.1 ) THEN

        write(nchanl,'(1x,a)') 'km          '
        write(nchanl,900) km

      ENDIF   ! trbout

900   format(1x,8e16.10)
901   format(1x,a80)

      RETURN
      END     
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE KNTARY                     ######
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 KNTARY(fmtver,
     :            grdout,basout,varout,mstout,iceout,trbout,
     :            n3dary,n3dsclr)
c
c#######################################################################
c
c     PURPOSE:
c
c     Count the number of 3-D arrays and 3-D scalars that should be written
c     into a history file. This number depends on the switches that have 
c     been set and the particular version of the code. For the current 
c     version of the code, the format version is not consulted.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/14/92
c
c     MODIFICATION HISTORY:
c
c#######################################################################
c
c     INPUT:
c
c     fmtver  Format/Version number a character string describing the
c             format and version number of the grid writing routine
c
c     grdout =0 or 1. If grdout=0, grid variables are not written. 
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 written. 
c     mstout =0 or 1. If mstout=0, water variables are not written. 
c     iceout =0 or 1. If iceout=0, writing of qi, qs and qh is skipped.
c     trbout =0 or 1. If trbout=0, writing of turbulence param km is skipped.
c
c     OUTPUT:
c
c     n3dary     Total number of 3-D arrays written into history file
c     n3dsclr    Total number of 3-D scalars written into history file
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      character*40 fmtver
      integer grdout,basout,varout,mstout,iceout,trbout
      integer n3dary,n3dsclr
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      n3dary=0
      n3dsclr=0
      IF( grdout.eq.1 ) n3dary=n3dary+4
      IF( basout.eq.1 ) THEN
        n3dary=n3dary+7
        n3dsclr=n3dsclr+4
      ENDIF
      IF( varout.eq.1 ) THEN
        n3dary=n3dary+5
        n3dsclr=n3dsclr+2
      ENDIF
      IF( mstout.eq.1 ) THEN
        n3dary=n3dary+3
        n3dsclr=n3dsclr+3
        IF( iceout.eq.1 ) THEN
          n3dary=n3dary+3
          n3dsclr=n3dsclr+4
        ENDIF
      ENDIF
      IF( trbout.eq.1 ) THEN
        n3dary=n3dary+1
        n3dsclr=n3dsclr+1
      ENDIF

      RETURN
      END
c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE EDGFILL                    ######
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 EDGFILL(a,nx1,nx2,ibgn,iend,ny1,ny2,jbgn,jend,
     :                     nz1,nz2,kbgn,kend)
c
c#######################################################################
c
c     PURPOSE:
c
c     Fill in the edges of a data array from the valid interior grid
c     points so that data occupies each storage. This is done primarily
c     for the benefit of the compression routine which must have a valid
c     value in each array location.
c
c#######################################################################
c
c     AUTHOR: Keith Brewster
c     7/14/92.
c
c#######################################################################
c
c     INPUT:
c
c       a        Array to have edges filled in
c
c       nx1:nx2  Dimensioned range of the  first index of array "a"
c       ny1:ny2  Dimensioned range of the second index of array "a"
c       nz1:nz2  Dimensioned range of the  third index of array "a"
c
c       ibgn,iend  Valued range of the  first index of array "a"
c       jgbn,jend  Valued range of the second index of array "a"
c       kbgn,kend  Valued range of the  third index of array "a"
c
      implicit none
c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      integer nx1,nx2,ny1,ny2,nz1,nz2
      real  a(nx1:nx2,ny1:ny2,nz1:nz2)
      integer ibgn,iend,jbgn,jend,kbgn,kend
c
c#######################################################################
c
c     Misc. local variables: 
c
c#######################################################################
c
      integer i,j,k
c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
c     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
c
c#######################################################################
c
c     Fill in top and bottom grid levels, as needed.
c
c#######################################################################
c
      IF( kbgn.gt.nz1 ) THEN
        DO 10 k=nz1,(kbgn-1)
        DO 10 j=jbgn,jend
        DO 10 i=ibgn,iend
          a(i,j,k)=a(i,j,kbgn)
 10     CONTINUE
      ENDIF
c
      IF( kend.lt.nz2 ) THEN
        DO 20 k=(kend+1),nz2
        DO 20 j=jbgn,jend
        DO 20 i=ibgn,iend
          a(i,j,k)=a(i,j,kend)
 20     CONTINUE
      ENDIF
c
c#######################################################################
c
c     Fill in north and south edges, as needed.
c
c#######################################################################
c
      IF( jbgn.gt.ny1 ) THEN
        DO 30 k=nz1,nz2
        DO 30 j=ny1,(jbgn-1)
        DO 30 i=ibgn,iend
          a(i,j,k)=a(i,jbgn,k)
 30     CONTINUE
      ENDIF
c
      IF( jend.lt.ny2 ) THEN
        DO 40 k=nz1,nz2
        DO 40 j=(jend+1),ny2
        DO 40 i=ibgn,iend
          a(i,j,k)=a(i,jend,k)
 40     CONTINUE
      ENDIF
c
c#######################################################################
c
c     Fill in east and west edges, as needed.
c
c#######################################################################
c
      IF( ibgn.gt.nx1 ) THEN
        DO 50 k=nz1,nz2
        DO 50 j=ny1,ny2
        DO 50 i=nx1,(ibgn-1)
          a(i,j,k)=a(ibgn,j,k)
 50     CONTINUE
      ENDIF
c
      IF( iend.lt.nx2 ) THEN
        DO 60 k=nz1,nz2
        DO 60 j=ny1,ny2
        DO 60 i=(iend+1),nx2
          a(i,j,k)=a(iend,j,k)
 60     CONTINUE
      ENDIF
c
      RETURN
      END
