c
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE SVIDUMP                    ######
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 SVIDUMP(nx,ny,nz, nchout,graffn, 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, grafgrid)
c
c#######################################################################
c
c     PURPOSE:
c
c     Output a file in GRAF format for the visualization program Savi3D.
c
c#######################################################################
c
c     AUTHOR: Jason J. Levit
c     07/20/92
c
c     MODIFICATION HISTORY:
c
c     08/02/92 (J. Levit)
c     Added full documentation and performed a clean-up.
c
c     08/04/92 (M. Xue)
c     Subroutine streamlined to conform to data dump format standards.
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     9/18/92 (J. Levit, M. Xue)
c     Added code to produce a configuration file for Savi3D.
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       graffn   Name of the Savi3D GRAF 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                (not used in this routine).
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       wbar     Base state vertial 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       grafgrid Passed from a dummy variable, used to define the grid
c                in Savi3D
c
c     OUTPUT:
c
c       None.
c
c     WORK ARRAY:
c
c       tem1     Work array
c       grafgrid A work arrays of size 3*nx*ny*nz
c
c#######################################################################
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
      character graffn*(*)   ! Name of the Savi3D GRAF file
      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 wbar  (nx,ny,nz)     ! Base state w-velocity (m/s)
      real ptbar (nx,ny,nz)     ! Base state potential temperature (K)
      real pbar  (nx,ny,nz)     ! Base state pressure (Pascal)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real qvbar (nx,ny,nz)     ! Base state water vapor 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  (ny-1,nx-1,nz-1)     ! Temporary work array.
      real grafgrid(3,ny-1,nx-1,nz-1) ! Used by Savi3D to define the
                                      ! data grid.
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c

      integer     ierr          ! Used as an internal error code by Savi3D.
      integer     i,j,k         ! Used by do loops.
      integer     nchout        ! Used to open Savi3D config file.
      integer ishf, jshf, kshf
      character *50 configname  ! Used to create Savi3D config file.
      real conx, cony           ! Used to create Savi3D config file.
      real scale
      real mtokm    ! conversion factor: meters to kilometers
                    ! used for scaling coordinate grid for output
      parameter ( mtokm = 0.001)
      integer gbwrtn  ! See if grid and base state parameter/arrays
                      ! have been written into the data file
      data gbwrtn /0/
      save gbwrtn
      character chtem1*7, chtem2*7 ! Used to create Savi3D config file.
c
c#######################################################################
c
c     Include files:
c
c#######################################################################
c

      include 'globcst.inc'
      include 'grafibm.inc'

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     Data needs to be in meters for Savi3D, so the following
c     statements convert the data from kilometers to meters.
c
c#######################################################################
c
      IF(stagrd.eq.1) THEN
        ishf=1
        jshf=1
        kshf=1
        scale=1.0
      ELSE
        ishf=0
        jshf=0
        kshf=0
        scale=1000.0
      END IF


c
c#######################################################################
c
c      Create the Savi3D configuration file.
c
c#######################################################################
c
        configname=runname(1:6)//'.sviconfig'
        open (unit=nchout,file=configname(1:16),status='unknown',
     :  form='formatted')

        cony=((ny-3)*(1000.0)*(1.0/111111.0))
        conx=((nx-3)*(1000.0)*(1.0/111111.0))

        write(chtem1,'(f7.3)') cony
        write(chtem2,'(f7.3)') conx

        DO 1 i=1,7
          IF ( chtem1(i:i).eq.' ') chtem1(i:i) = '0'
          IF ( chtem2(i:i).eq.' ') chtem2(i:i) = '0'
1       CONTINUE
  
        write (nchout,'(4a)')
     :  'toplat=',chtem1,',btmlat=000.000,lftlon=000.000,rhtlon=',chtem2
        write (nchout,'(a,a)') 'rafs_gridded=',graffn
        write (nchout,'(a)') 'initzstretch=1'

        close (nchout)

c
c#######################################################################
c
c     Setup for the GRAF interface.  This is only done at time=0,
c     and calls are made to Savi3D subroutines to define the type of
c     grid which the data are being written to.
c
c     Also, the scalar variables are defined for Savi3D here.
c
c#######################################################################
c

      IF ( gbwrtn.eq.0 ) THEN

        CALL GrafOpen(graffn,GrafNEWONLY,'3-D Model','ARPS 3.0',grafh)

        CALL GrafDefnGrid(grafh,0,'ARPS 3.0','3D',
     :       ny-1, nx-1, nz-1,
     :       0.0D0, 0.0D0, 0.0D0,
     :       GrafGridMeters,0,gridh)
 
        DO 1000 k=1,nz-1
        DO 1000 i=1,nx-1
        DO 1000 j=1,ny-1
          grafgrid(1,j,i,k)=scale*(x(i)+x(i+ishf))*0.5
          grafgrid(2,j,i,k)=scale*(y(j)+y(j+jshf))*0.5
          grafgrid(3,j,i,k)=scale*(zp(i,j,k)+zp(i,j,k+kshf))*0.5
1000    CONTINUE 

        CALL GrafWriteGridArray(gridh,grafgrid,ierr)

        IF (varout .eq. 1) THEN

          CALL GrafDefnScalar(gridh,'u','m/s',
     :          'X-velocity total wind',
     :          GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'v','m/s',
     :          'Y-velocity total wind',
     :          GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'w','m/s',
     :          'Z-velocity total wind',
     :          GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'uprt','m/s',
     :          'X-velocity perturbation',
     :          GrafSAVEMAXMIN,GrafTIMEVAR,ierr)
   
          CALL GrafDefnScalar(gridh,'vprt','m/s',
     :          'Y-velocity perturbation',
     :          GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'wprt','m/s',
     :          'Z-velocity perturbation',
     :          GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'ptprt','deg K',
     :         'Perturbation Potential Temperature', 
     :         GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'pprt','mb',
     :          'Perturbation Pressure',
     :          GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

        ENDIF

        IF (mstout .eq. 1) THEN

          CALL GrafDefnScalar(gridh,'qvprt','g/kg',
     :         'Water Vapor Mixing Ratio Perturbation',
     :         GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'qc','g/kg',
     :         'Cloud Water Mixing Ratio',   
     :         GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'qr','g/kg',
     :         'Rain Water Mixing Ratio',
     :         GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          CALL GrafDefnScalar(gridh,'qw','g/kg',
     :         'Total Water Mixing Ratio',
     :         GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

          IF (iceout .eq. 1) THEN

            CALL GrafDefnScalar(gridh,'qi','g/kg',
     :           'Cloud Ice Mixing Ratio',
     :           GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

            CALL GrafDefnScalar(gridh,'qs','g/kg',
     :           'Snow Mixing Ratio',
     :           GrafSAVEMAXMIN,GrafTIMEVAR,ierr)

            CALL GrafDefnScalar(gridh,'qh','g/kg',
     :           'Hail Mixing Ratio',
     :           GrafSAVEMAXMIN,GrafTIMEVAR,ierr)


          ENDIF

        ENDIF

c
c#######################################################################
c
c     The definitions for the time invariant scalars for Savi3D.
c     They are only written once to the Savi3D GRAF file (time=0).
c
c#######################################################################
c

        IF (basout .eq. 1) THEN

        CALL GrafDefnScalar(gridh,'ubar','m/s',
     :       'Base State X Wind Velocity',
     :       GrafSAVEMAXMIN,GrafTIMEINVAR,ierr)

        CALL GrafDefnScalar(gridh,'vbar','m/s',
     :       'Base State Y Wind Velocity',
     :       GrafSAVEMAXMIN,GrafTIMEINVAR,ierr)

        CALL GrafDefnScalar(gridh,'wbar','m/s',
     :       'Base State Z Wind Velocity',
     :       GrafSAVEMAXMIN,GrafTIMEINVAR,ierr)

        CALL GrafDefnScalar(gridh,'ptbar','deg K',
     :       'Base State Potential Temperature',
     :        GrafSAVEMAXMIN,GrafTIMEINVAR,ierr)

        CALL GrafDefnScalar(gridh,'pbar','mb',
     :       'Base State Pressure',
     :       GrafSAVEMAXMIN,GrafTIMEINVAR,ierr)

        CALL GrafDefnScalar(gridh,'rhobar','kg/m**3',
     :       'Base State Air Density',
     :       GrafSAVEMAXMIN,GrafTIMEINVAR,ierr)

          IF (mstout .eq. 1) THEN

            CALL GrafDefnScalar(gridh,'qvbar','g/kg',
     :           'Base State Water Vapor Mixing Ratio', 
     :           GrafSAVEMAXMIN,GrafTIMEINVAR,ierr)

          ENDIF

        ENDIF

c
c
c#######################################################################
c
c     Vector definitions for Savi3D.  Also, define the starting date,
c     time, and time step for Savi3D.
c
c#######################################################################
c
        IF (varout .eq. 1) THEN

          CALL GrafDefnVect(gridh, 'Pertubation_Wind',
     :        'uprt','vprt', 'wprt', ierr)

          CALL GrafDefnVect(gridh, 'Total_Wind',
     :        'u', 'v', 'w', ierr)

        ENDIF

        IF (basout .eq. 1) THEN

          CALL GrafDefnVect(gridh, 'basewind',
     :        'ubar', 'vbar', 'wbar', ierr)

        ENDIF

        CALL GrafTimeStart(grafh,1970,1,1,00,00,00)
        CALL GrafTimeStep(grafh,0,0,0,nint(thisdmp))

      ENDIF

c#######################################################################
c
c     Write the base state time invariant scalars to the Savi3D file.
c     Data is written the first time this routine is called. 
c
c#######################################################################
c

      CALL GrafNewFrame(grafh)

      IF ( gbwrtn .eq. 0 ) THEN

        IF (basout .eq. 1 ) THEN

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

          CALL GrafWriteScalarArray(gridh,'ubar',tem1,ierr)

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

          CALL GrafWriteScalarArray(gridh,'vbar',tem1,ierr)

          DO 1300 k=1,nz-1
          DO 1300 i=1,nx-1
          DO 1300 j=1,ny-1
c           tem1(j,i,k)=((wbar(i,j,k)+wbar(i,j,k+kshf))*0.5)
            tem1(j,i,k)=0.0
1300      CONTINUE

          CALL GrafWriteScalarArray(gridh,'wbar',tem1,ierr)

          DO 1400 k=1,nz-1
          DO 1400 i=1,nx-1
          DO 1400 j=1,ny-1
            tem1(j,i,k)=ptbar(i,j,k)
1400      CONTINUE

          CALL GrafWriteScalarArray(gridh,'ptbar',tem1,ierr)

          DO 1500 k=1,nz-1
          DO 1500 i=1,nx-1
          DO 1500 j=1,ny-1
            tem1(j,i,k)=(pbar(i,j,k)/100)
1500      CONTINUE

          CALL GrafWriteScalarArray(gridh,'pbar',tem1,ierr)

          DO 1600 k=1,nz-1
          DO 1600 i=1,nx-1
          DO 1600 j=1,ny-1
            tem1(j,i,k)=rhobar(i,j,k)
1600      CONTINUE

          CALL GrafWriteScalarArray(gridh,'rhobar',tem1,ierr)

          IF (mstout .eq. 1) THEN

            DO 1700 k=1,nz-1
            DO 1700 i=1,nx-1
            DO 1700 j=1,ny-1
              tem1(j,i,k)=qvbar(i,j,k)*1000.0
1700        CONTINUE

            CALL GrafWriteScalarArray(gridh,'qvbar',tem1,ierr)

          ENDIF

        ENDIF

      ENDIF

c
c#######################################################################
c
c     Write the scalar arrays to the Savi3D file, and mark the end
c     of the frame to the file.
c
c#######################################################################
c

      IF (varout .eq. 1) THEN

        DO 1800 k=1,nz-1
        DO 1800 i=1,nx-1
        DO 1800 j=1,ny-1
          tem1(j,i,k)=(u(i,j,k)+u(i+ishf,j,k))*0.5
1800    CONTINUE

        CALL GrafWriteScalarArray(gridh,'u',tem1,ierr)

        DO 1900 k=1,nz-1
        DO 1900 i=1,nx-1
        DO 1900 j=1,ny-1
          tem1(j,i,k)=(v(i,j,k)+v(i,j+jshf,k))*0.5
1900    CONTINUE

        CALL GrafWriteScalarArray(gridh,'v',tem1,ierr)

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

        CALL GrafWriteScalarArray(gridh,'w',tem1,ierr)

        DO 2100 k=1,nz-1
        DO 2100 i=1,nx-1
        DO 2100 j=1,ny-1
          tem1(j,i,k)=(u(i,j,k)+u(i+ishf,j,k)
     :                -ubar(i,j,k)-ubar(i+ishf,j,k))*.5
2100    CONTINUE

        CALL GrafWriteScalarArray(gridh,'uprt',tem1,ierr)

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

        CALL GrafWriteScalarArray(gridh,'vprt',tem1,ierr)

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

        CALL GrafWriteScalarArray(gridh,'wprt',tem1,ierr)

        DO 2400 k=1,nz-1
        DO 2400 i=1,nx-1
        DO 2400 j=1,ny-1
          tem1(j,i,k)=ptprt(i,j,k)
2400    CONTINUE

        CALL GrafWriteScalarArray(gridh,'ptprt',tem1,ierr)

        DO 2500 k=1,nz-1
        DO 2500 i=1,nx-1
        DO 2500 j=1,ny-1
          tem1(j,i,k)=pprt(i,j,k)*0.01
2500    CONTINUE

        CALL GrafWriteScalarArray(gridh,'pprt',tem1,ierr)

      ENDIF

      IF (mstout .eq. 1) THEN

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

        CALL GrafWriteScalarArray(gridh,'qvprt',tem1,ierr)

        DO 2700 k=1,nz-1
        DO 2700 i=1,nx-1
        DO 2700 j=1,ny-1
          tem1(j,i,k)=qc(i,j,k)*1000.0
2700    CONTINUE

        CALL GrafWriteScalarArray(gridh,'qc',tem1,ierr)

        DO 2800 k=1,nz-1
        DO 2800 i=1,nx-1
        DO 2800 j=1,ny-1
          tem1(j,i,k)=qr(i,j,k)*1000.0
2800    CONTINUE

        CALL GrafWriteScalarArray(gridh,'qr',tem1,ierr)

        DO 2900 k=1,nz-1
        DO 2900 i=1,nx-1
        DO 2900 j=1,ny-1
          tem1(j,i,k)=(qr(i,j,k)+qc(i,j,k))*1000.0
2900    CONTINUE

        CALL GrafWriteScalarArray(gridh,'qw',tem1,ierr)

        IF (iceout .eq. 1) THEN

          DO 3000 k=1,nz-1
          DO 3000 i=1,nx-1
          DO 3000 j=1,ny-1
            tem1(j,i,k)=qi(i,j,k)*1000.0
3000      CONTINUE

          CALL GrafWriteScalarArray(gridh,'qi',tem1,ierr)

          DO 3100 k=1,nz-1
          DO 3100 i=1,nx-1
          DO 3100 j=1,ny-1
            tem1(j,i,k)=qs(i,j,k)*1000.0
3100      CONTINUE

          CALL GrafWriteScalarArray(gridh,'qs',tem1,ierr)

          DO 3200 k=1,nz-1
          DO 3200 i=1,nx-1
          DO 3200 j=1,ny-1
            tem1(j,i,k)=qh(i,j,k)*1000.0
3200      CONTINUE

          CALL GrafWriteScalarArray(gridh,'qh',tem1,ierr)

        ENDIF

      ENDIF

4100  CALL GrafEndFrame(grafh,GrafUPDATE)

      gbwrtn = 1 ! This routine has been called once. 

      RETURN
      END
