
c No bloody filtering (dhp)



c
      subroutine do_ppm_de0_3dles_gam
     1       (xl, ro_pen, pp_pen, u1_pen, u2_pen, u3_pen, sk_pen,

     1             rhonu,    pnu,   u1nu,   u2nu,   u3nu, sknu,
     1        eosgam, dt, smlrho, smallp, smallu,
     1        smalle, courmx, n_real_zones_in, nbdy_set, nrows, ky, kz)
c
c  **************************************************************************
c  *                          THIS CODE IS do_ppm_de0_3dles_gam.             *
c  *          THIS CODE WAS WRITTEN BY: WOODWARD RESEARCH GROUP. (c) 2001   *
c  *                          ALL RIGHTS RESERVED.                          *
c  **************************************************************************
c
c
c
c
C  WHAT              This routine performs 1 pass of a 3D direct Eulerian 
c                    hydrodynamic update,  with artificial diffusion.
c
c  ASSUMPTIONS       A uniform, Cartesian grid is assumed.
c                    A gamma-law equation of state is used:
C                         p(i) = ( eosgam - 1.0) * rho(I) * e_internal(i)
C                             where eosgam [gamma] is the ratio of specific
C                             heats.
C
C  BOUNDARIES       The necessary number of boundary zones to set (NBDY_SET)
c                   is 5, three for interpolation, one for integration, and
c                   Our Fiend Diffuse takes one.  It will be further assumed
c                   that at the outset the 2*NBDY_SET boundary zones [often
c                   referred to as fake zones] contain valid values
c                   appropriate to the user's chosen boundary conditions.
C                   This is the only means of signalling boundary conditions
C                   to this routine The user signals his understanding and
c                   acknowledges this by passing to Do_PPM_DE0_1DGAMMA
c                   the number of boundary zones set on each side.  This
c                   routine verifies that sufficient fake zones have been
c                   provided before continuing.  If not,  an advisory is
c                   displayed and the program terminated.
C
C  DIMENSIONS
c      input        It will be assumed that the dimensions of all zone-averaged
c                   type quantities { e.g. RHO, P, U1, U2, U3 RHONU, PNU, 
c                   etc, ]  are
C                             [1-NBDY_SET:N_REAL_ZONES_IN+NBDY_SET].
C                   The dimensions of all [and I mean ALL] interface type
C                   quantities [e.g. xl, dvoll, dmassl, pavl, uavl, anything
C                   referenced to and interface, usually the lefthand zone
c                   interface] are
C                              [1-NBDY_SET:N_REAL_ZONES_IN+NBDY_SET+1].
C                   Note the "+1" please.  There may be a quiz later.
C

C      output       At the completion of this entire routine, new zone-averaged
C                   quantities [e.g. rhonu, pnu, u1nu, u2nu, u3nu] will have 
c                   updated values on the range
C                                [1:N_REAL_ZONES_IN],
C                   and any returned interface type thingies will have valid
C                  values on the range [There are none at this time]:
C                                [1:N_REAL_ZONES_IN+1].
c
C
c
C  DATA
C     input        Input arrays are:
C
C                   xl:  the locations of the lefthand zone interfaces,
c                        defining the computational grid.  Please realize
c                        that a numerical grid of n zones would require n+1 xl
c                        values (interface locations.)
c               ro_pen:  The zone-averaged densities (3-D pencil array)
c               pp_pen:  The zone-averaged pressures (3-D pencil array)
c               u1_pen:  The zone-averaged long-vels (3-D pencil array)
c               u2_pen:  The zone-averaged tran-vels (3-D pencil array)
c               u3_pen:  The zone-averaged tran-vels (3-D pencil array)
c               sk_pen:  The SGS kinetirc energy (3-D pencil array)
c
C                  Input scalers are:
C
C                   eosgam: the ratio of specific heats (gamma) for a
c                           gamma-law gas.
C                       dt:   the timestep.
C                   smlrho:  a trivial value for density
C                   smallp:  a trivial value for pressure
C                   smallu:  a trivial value for velocity
C                   smalle:  a trivial value for nergy
C                            the  above 4 values are generally about
C                                      1.0e-06 for single precision and
C                                      1.0d-08 for dourble precision.
C          n_real_zones_in:  The number zones (NOT counting fake zones) to be
c                               processed.  Upon completion of the timestep,
c                               this many zones will have valid entries in the
c                               output arrays.
C                 nbdy_set:  Number of boundary zones on either end set.  This
c                            is a check to ensure sufficient values have been
c                            set by the calling routine.  The number of
c                            boundary zones required on either side must be
c                            at least equat to the required values: 5
c
c                    nrows:  internal transverse dimension of pencil
c                       ky:  transverse offest of update strip
c                       kz:  transverse offest of update strip
c
C
C     output       Output arrays are:

C                      rhonu:  The new zone-averaged densities at the new
c                              timestep.
C                        pnu:  The new zone-averaged pressures at the new
c                              timestep.
C                       u1nu:  The new zone-averaged longitudinal-velocities 
c                              at the new timestep.
C                       u2nu:  The new zone-averaged transverse-velocities 
c                              at the new timestep.
C                       u3nu:  The new zone-averaged transverse-velocities 
c                              at the new timestep.
C                       sknu:  The new zone-averaged SGS K.E.
c                              at the new timestep.
C
C                  Output scaler is:
C
C                     courmx:  The maximum Courant number found by this call
C                              to this subroutine.
C
C
C  SPECIAL NOTES
c
C    directional splitting 
c
c                   PPM performs multi-dimensional hydrodynamics by directional 
c                   splitting, in which a one-dimensional operator is applied
c                   in alteration, e. g. a sequence of passes x-y-z-z-y-x.  
c                   During the hydrodynamic step, the transverse velocities are 
c                   advected passively.  Therefore the following proceedure 
c                   will be followed:
c
c                      1)  Find the the fluxes of the active hydrodynamical 
c                          quantities using do_ppm98_deul0_flux_gamma_1d.  The 
c                          active hydro. quantities are: rho, pressure, 
c                          longitudinal velocity (u1).
c
c                      2)  using these fluxes to update the zone averages of 
c                          density (rhonu), longitudinal velocity (u1nu), and 
c                          a partial update of the total energy (e_totnu) 
c                          using ppm98_deul_update_1d.
c
c                          The portion of the total energy updated is the 
c                          internal energy + the longitudinal kinetic energy.  
c                          That is  we will write the total energy as
c
c                             ET(I) = EI(I) + EKL(I) + EKT(I), 
c                          where
c
c                             ET(I) = total specific energy in zone I
c
c                             EI(I) = specific internal energy in zone I
c
c                             EKL(I) = total specific longitudinal kinetic 
c                                      energy in zone I  
c                                      ekl(i) = 0.5 * u1(i)*u1(i)
c            
c                             EKT(I) = total specific transverse kinetic 
c                                     energy in zone I  
c                                     ekt(i) = 0.5 * (u2(i)*u2(i) + u3(i)*u3(I))
c
c                      3)  From the fluxes at the zone interfaces combined 
c                          with  the old and new zone masses, it is possible 
c                          to find the new zone averages of the transverse 
c                          velocity (UTNU).  The fluxes of the transverse 
c                          velocity are then used to determine the transverse 
c                          kinetic energy portion of the total energy, which 
c                          is then added to the earlier partial ET update of 
c                          step2.
c    Velocities
c
c                   It is assumed that the velocities u1, u2, and u3 are from
c                   a right hand coordinate system, hence:
c                            dir1 = dir2 cross dir3,
c                            dir2 = dir3 cross dir1,
c                            dir3 = dir1 cross dir2,
c                   dir# is the unit vector in the # direction.
c
c
c
C
C
C
C
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C
C these permit switching btw double and single precision

C

























c



C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C*************************************************************************
C
c
C
c
c
C  
      

c


      dimension       xl(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^               rho(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^                 p(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^                u1(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^                u2(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^                u3(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^                sk(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^             rhonu(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^               pnu(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^              u1nu(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^              u2nu(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^              u3nu(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^              sknu(1-nbdy_set:n_real_zones_in + nbdy_set)
c      character*100 msg
      integer msg
      parameter(msg=0)

c   note that ALL interface thingies are dimensioned 1.0e+00 more.
c
C F77 \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ \/ 
      dimension         courno(1-nbdy_set:n_real_zones_in + nbdy_set),   
     ^              difuse_vel(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^                  unsmad(1-nbdy_set:n_real_zones_in + nbdy_set),       
     ^                  unsmup(1-nbdy_set:n_real_zones_in + nbdy_set),   
     ^                unsmooth(1-nbdy_set:n_real_zones_in + nbdy_set),   
     ^                  shockd(1-nbdy_set:n_real_zones_in + nbdy_set),       
     ^                 e_totnu(1-nbdy_set:n_real_zones_in + nbdy_set),       
     ^                     dxl(1-nbdy_set:n_real_zones_in + nbdy_set),     
     ^                    dvol(1-nbdy_set:n_real_zones_in + nbdy_set),     
     ^                   dvoll(1-nbdy_set:n_real_zones_in + nbdy_set),     
     ^                      dx(1-nbdy_set:n_real_zones_in + nbdy_set),     
     ^                   dxinv(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^                    ceul(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^                    cesq(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^                    clagl(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^                 e_total(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^                  stpadb(1-nbdy_set:n_real_zones_in + nbdy_set) 
      dimension       rhol(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     ^                  pl(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     ^                  ul(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     1                rhor(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     ^                  pr(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     ^                  ur(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     2                drho(1-nbdy_set:n_real_zones_in+nbdy_set),
     ^                  dp(1-nbdy_set:n_real_zones_in+nbdy_set),
     ^                  du(1-nbdy_set:n_real_zones_in+nbdy_set),
     3                rho6(1-nbdy_set:n_real_zones_in+nbdy_set),
     ^                  p6(1-nbdy_set:n_real_zones_in+nbdy_set),
     ^                  u6(1-nbdy_set:n_real_zones_in+nbdy_set) 
      dimension      ullrp(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     ^               pllrp(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     ^               urlrm(1-nbdy_set:n_real_zones_in+nbdy_set+1),
     1               prlrm(1-nbdy_set:n_real_zones_in+nbdy_set+1) 

      dimension         rhoavl(1-nbdy_set:n_real_zones_in + nbdy_set+1),    
     ^                    pavl(1-nbdy_set:n_real_zones_in + nbdy_set+1),     
     ^                    uavl(1-nbdy_set:n_real_zones_in + nbdy_set+1),  
     ^                 etotavl(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^                  difusl(1-nbdy_set:n_real_zones_in + nbdy_set+1),   
     ^                  dmassl(1-nbdy_set:n_real_zones_in + nbdy_set+1),    
     ^                  dmasll(1-nbdy_set:n_real_zones_in + nbdy_set+1),    
     ^                  dmasrl(1-nbdy_set:n_real_zones_in + nbdy_set+1),    
     ^                  dmom1l(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^                    dmnu(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^                   dmnew(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^                      dm(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^                  detotl(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^                 dvolinv(1-nbdy_set:n_real_zones_in + nbdy_set+1),
     ^                   dminv(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^                 dmnuinv(1-nbdy_set:n_real_zones_in + nbdy_set) 
      dimension  smalldu(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^           smallsk(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^            dvolfl(1-nbdy_set:n_real_zones_in + nbdy_set+1), 
     ^          unsmutu2(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^          unsmutu3(1-nbdy_set:n_real_zones_in + nbdy_set),
     ^          unsmutuk(1-nbdy_set:n_real_zones_in + nbdy_set)
      dimension  steepu2(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^           steepu3(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^           steepuk(1-nbdy_set:n_real_zones_in + nbdy_set), 
     ^            dmom2l(1-nbdy_set:n_real_zones_in + nbdy_set+1), 
     ^            dmom3l(1-nbdy_set:n_real_zones_in + nbdy_set+1), 
     ^            drhokl(1-nbdy_set:n_real_zones_in + nbdy_set+1), 
     ^           dektl(1-nbdy_set:n_real_zones_in + nbdy_set) 
      dimension       niterl(1-nbdy_set:n_real_zones_in+nbdy_set+1) 
C F77 /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ /\ 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


      dimension ro_pen(1-nbdy_set:n_real_zones_in + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set  )
      dimension pp_pen(1-nbdy_set:n_real_zones_in + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set  )
      dimension u1_pen(1-nbdy_set:n_real_zones_in + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set  )
      dimension u2_pen(1-nbdy_set:n_real_zones_in + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set  )
      dimension u3_pen(1-nbdy_set:n_real_zones_in + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set  )
      dimension sk_pen(1-nbdy_set:n_real_zones_in + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set,
     1                 1-nbdy_set:nrows           + nbdy_set  )

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     integer ifirst
c     common /dhp/ ifirst
c
c     data ifirst/1/
c     if(ifirst .eq. 1) then
c       ifirst = 0
c       write (6,*) 'nbdy_set = ', nbdy_set
c       write (6,*) 'u1(i), u1_pen(i,ky,kz): ky, kz = ', ky, kz
c       do i = 1-nbdy_set, 10
c         write (6,*) u1(i), u1_pen(i,ky,kz)
c       enddo
c     endif
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      small = 1.0e-08


c
c   We will create and use the following "K" indices as pointers to
c     the beginning and end of the working array.  A count KDO will
c     also be maintained of the number of zones in the array to 
c     process.  Note that all of these shall refer to ZONE counts.
c     Quantities which depend upon lefthand zone interfaces, will
c     require 1.0e+00 additional reference on the right.
c     the K*1's are appropriate whenever starting with interpolating

      smlfrc = 5.0e-03

      courmx = 0.0

      n_required =  5
      if ( nbdy_set .lt. n_required ) then
        write(msg,6666) n_required
6666    format(1x,  i3,' fake zones at each end of the array')
        write(msg,6667)  nbdy_set,n_required,n_real_zones_in
6667    format('The calling program passed',i6,' boundary zones,',
     1       'but ',i3, ' zones are required.')
        stop
      endif
c
c
      N = n_real_zones_in
      KSTART1 = 1-NBDY_SET
      KEND1 = N_REAL_ZONES_IN + NBDY_SET
      KDO1 = KEND1 - KSTART1 +1
      K = KSTART1
      nzones = KDO1
c   
866     format(1h ,i4,1p,9(1x,e10.3))

c      compute some essential quantities:
c
       do i=KSTART1, KEND1
          rho(i) = ro_pen(i,ky,kz)
          p(i)   = pp_pen(i,ky,kz)
          u1(i)  = u1_pen(i,ky,kz)
          u2(i)  = u2_pen(i,ky,kz)
          u3(i)  = u3_pen(i,ky,kz)
          sk(i)  = sk_pen(i,ky,kz)

	  dx(i) = xl(i+1) - xl(i)
          dxinv(i) = 1.0 / dx(i)
          dvol(i) = dx(i)
          dvolinv(i) = dxinv(i)
c         these last 2.0e+00 are true only in cartesian coords.
          dm(i) = rho(i) * dvol(i)
          dminv(i) = 1.0 / dm(i)
          cesq(i) = (eosgam * p(i) / rho(i) )
          ceul(i) = sqrt (cesq(i))
          smalldu(i) = small * ceul(i)
          smallsk(i) = smalle
       enddo

c  Get the fluxes of the active hydrodynamical quantities:
c    These are mass, longitudinal momentum, and energy.

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c   Generate momentum fluxs due to SGS turbulence in opperator split
c   fasion.   Use imperical LES model based on velocity gradients of
c   the resolved field.
c   1) Evaluate all 9 velocity gradients
c   2) Evaluate model for SGS stresses from velocity gradients
c   3) Add these SGS stresses to momentum fluxs below
c
c   Open questions:
c   1) What is correct way to modify energy equation?
c      a) follow an SGS field of k?
c      b) dump energy into heat?
c      c) not follow SGS energy at all?
c   2) What is atom for velocity gradients 
c      a) is there a proper filter for resolved fields?
c      b) make for isotropic evaluation of gradients?
c   3) correct filter width?
c      a) Should it be ~8dx in wavelength?
c      b) what shape? (Gaussian or implicity with PPM or what?)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc



c            
c    This requires the following steps:
c
c
c                   1) Interpolate, with monotonicity constraints the
c                          pressure and longitudinal velocity by constructing
c                          interpolations of the +- Riemann invarients.
c                          (Returns  UNSMUP)
c                   2) Interpolate, with monotonicity constraints and
c                          discontinuity detection, the density by
c                          constructing interpolations of the  0 Riemann
c                          invarient (streamline, or adiabatic invarient.)
c                          (Returns  UNSMAD, STPADB)
c                   3) Averaging the piecewise parabolic interpolations for
c                          density, velocity, and pressure, find the
c                          averages in the appropriate domains of dependence
c                          for a Riemann problem at each zone interface.
c                   4) Solve Riemann shock tube problems at each lefthand
c                          zone interface to obtain appropriate time-averaged
c                          values for density, pressure, and longitudinal
c                          momentum at these interfaces.
c                          (Returns RHOAVL, PAVL, UAVL.)
c                   5) From these time averages, construct interface fluxes
c                           (in the Eulerian sense) of mass, momentum, and
c                           total energy.
c                           (Returns DMASSL, DMOML, and DETOTL.)
c                    6) Determine the Courant condition (stability) for
c                           each _INTERFACE_.  Simply put, in order for
c                           the method to be stable, the various Domains
c                           of dependence for the Riemann problems
c                           cannont exceed the width of either the zone
c                           to the left or the zone to the right.  This
c                           make the courant number array COURNO appear
c                           to be an interface-like quantity.  Its valid
c                           range shall be [K:K+NZONES]

c############################################################################
c
c
c    Now find the spatially-averaged values in the left and right
c    domain of dependence for the lefthand zone interface for each zone.
c
c    PPM98_EULER0_DOMAIN_AVG calculates the left and right states for
c       a Riemann problem to each lefthand zone interface in the
c       following steps:
c                   1) Interpolate, with monotonicity constraints the
c                          pressure and longitudinal velocity by constructing
c                          interpolations of the +- Riemann invarients.
c                          (Returns  UL, UR, DU, U6, PL, PR, DP, P6,  UNSMUP)
c                   2) Interpolate, with monotonicity constraints and
c                          discontinuity detection, the density by
c                          constructing interpolations of the  0 Riemann
c                          invarient (streamline, or adiabatic invarient.)
c                          (Returns  UNSMAD, STPADB, RHOL, RHOR, DRHO, RHO6)
c                   3) Using the piecewise parabolic interpolations for
c                          velocity, and pressure, find the averages in the
c                          appropriate domains of dependence
cc
c
c   Note:
c     The average of density over the r0 DOD cannot
c     be found until the Riemann problem has been
c     solved.  Therefore, the density interplotions
c     are returned only.


      N = n_real_zones_in
      KSTART = 1-nbdy_set
      KEND = N+nbdy_set
      KDO = N + 2*nbdy_set
c
      nnbdy_set = 4
      K = KSTART
      N = KDO - 2*nnbdy_set

      call   ppm98_Euler0_domain_avg
     ^       (  dx(K),   rho(K),      p(K),     u1(K),   ceul(K),
     ^       ullrp(K), pllrp(K),  urlrm(K),  prlrm(K),  clagl(K),
     ^        rhol(K),  rhor(K),   drho(K),   rho6(K),     ul(K),
     ^          ur(K),    du(K),     u6(K),     pl(K),     pr(K),
     ^          dp(K),    p6(K), unsmad(K), stpadb(K), unsmup(K),
     ^      courno(K),       dt,        N,  nnbdy_set)

c

c
c   Update the zone counts
      kstart2 = kstart+nnbdy_set
      kend2 = kend-nnbdy_set
      kdo2 = kdo - 2*nnbdy_set

c   If 1.0e+00 wished to examine the domain averages 1.0e+00 could :
c
c       print *,' i, ullrp(i), pllrp(i), urlrm(i), prlrm(i)'
c      do i=kstart2,kend2+1
c       write(6,66) i, ullrp(i), pllrp(i), urlrm(i), prlrm(i)
c      enddo
c       print *,' i, rhol(i), rho(i), rhor(i), drho(i), rho6(i)'
c      do i=kstart2-1,kend2+1
c       write(6,66) i, rhol(i), rho(i), rhor(i), drho(i), rho6(i)
c      enddo
c
c
c ***********************************************************************
c *                                                                     *
c *   at this point, modifications to pllrp, prlrm to include body      *
c *   forces, or to ullrp, urlrm to include geometrical effects         *
c *    may be done. (See Colella and Woodward 1984).  By modifying      *
c *    the domain averages in this fashon, we may use a single          *
c *    Riemann solver.                                                  *
c *                                                                     *
c ***********************************************************************
c
c
c  Solve the Riemann problems and get the interface fluxes.
c
c
c
c    PPM98_E_RIEMANN_GAMMA uses the domain averages of pressure and
c       velocity to solve Riemann's shocktube problem at each zone
c       interface, to provide time-averaged values of pressure
c       and velocity (returns PAVL and UAVL) at the zone interfaces.
c       The density pressure interpolations are then used to find
c       the time-averaged density (returns RHOAVL).  From these
c       values, advected fluxes of volume, mass, momentum, and total
c       energy are computed. (Returns DVOLL, DMASSL, DMOML, and
c       DETOT1L)
c
      nnbdy = 1

c      KSTART2, KEND2, KDO2 refer to the dod interface values
c
c      we require nnbdy additional zone on the left and right
c      of the density and pressure interpolations.  We will use
c      these as our points of reference.  Hence we adjust these \c
c      K values

       KSTART_RIEMANN = KSTART2-nnbdy
       KEND_RIEMANN = KEND2+nnbdy
       KDO_RIEMANN = KDO2+2*nnbdy
c
c      NITERL is an integer array specifying the number of non-linear iterationsc        to take at each lefthand zone interface location.  In smooth flow, or if
c        1.0e+00 intends to rely upon the later addition of additional dissapation,
c        iterations are necessary. (NITERL(I)=0 for all i).  However, situationsc        may arise inwich the user may wish to take additional  iterations.  Generally
c        only in the most extreme situations would more than 1 additional iteration
c        be required.  The results from shock detection may be used to identify
c        such zones.  Here, we rely upon additional dissapation.
c


      K=KSTART_RIEMANN
      N = KDO_RIEMANN - 2*nnbdy

      do i=k,k+n+2*nnbdy+1
         niterl(i) = 0
      enddo

      call   PPM98_E_RIEMANN_GAMMA
     ^       (  dx(K),  ullrp(K),   pllrp(K), urlrm(K), prlrm(K),
     ^       clagl(K),   rhol(K),     rho(K),  rhor(K),  drho(K),
     ^        rho6(K),     pl(K),       p(K),    pr(K),    dp(K),
     ^          p6(K), rhoavl(K),    pavl(K),  uavl(K),  dvoll(K),
     ^      dmassl(K), dmom1l(K),  detotl(K),   eosgam,    smlrho,
     ^         smallp,        dt,  niterl(K), n, nnbdy)

c
c   Update the zone counts
c
        kstart3 = kstart2
        kend3 = kend2
        kdo3 = kdo2
c
c   If 1.0e+00 wished to examine the fluxes 1.0e+00 could:
c   (Note, these are interface values, hence the +1)
c      print *,' i, rhoavl(i), pavl(i), uavl(i)'
c     do i=kstart3,kend3+1
c      write(6,66) i, rhoavl(i), pavl(i), uavl(i)
c     enddo
c
c      print *,' i, dvoll(i), dmassl(i),dmom1l(i), detotl(i)'
c     do i=kstart3,kend3+1
c      write(6,66) i, dvoll(i), dmassl(i), dmom1l(i), detotl(i)
c     enddo
c
c


c
c
c ***************************************************************************
* *                                                                         *
c *   The conservation laws may now be applied to obtain new zone-averaged  *
c *   values of density, pressure (or energy), and velocity.                *
c *                                                                         *
c ***************************************************************************
c


666     format(1h ,i4,1p,9(1x,e10.3))
66      format(1h ,i4,1p,9(1x,e10.3))

c   the K*2 will be appropriate for working with fluxes (before diffuse)


c   find the courant number from the above step.  When a routine is
c   called with an array such as COURNO, the old values of COURNO are
c   overwritten.  Thus the Courant number should be determined after each
c   such routine, or different arrays used.
c
c   To be stable when a flux at xl(i), the left hand edge of zone I, is
c   computed, the courant numbers in zone I-1 and I need to be <= 1.0
c    Thus we examine Courant numbers over the KSTART2-1,KEND2+1


c   
      courmx = courno(KSTART2 -1)
      do i=KSTART2 , KEND2 +1
        courmx = max (courmx, courno(i))
      enddo

c     in order to update the zones, which is done with *conserved quantities*
c     (mass, momenta, energy), we need to compute the old zone averaged total
c     energy, which was not provided as initial data. How rude.  It is
c     computed only over the range required.


      gamm1i = 1.0 / ( eosgam - 1.0)
      do i=KSTART2 ,KEND2 
	 ekl = 0.5 * u1(I)*u1(I)
	 eint = gamm1i * p(i) / rho(i)
	 e_total(i) = eint + ekl
      enddo

c

c
c   PPM98_DEUL_UPDATE_CART  upates the zone averages by conserviatvely
c      differencing the fluxes found above, assuming cartesian
c      coordinates and no body forces.  Non-cartesian geometries
c      may be treated by including the area terms into the interface
c      fluxes BEFORE this routine is called.  Body forces would
c      entail and additional term added to the new zone-averages of
c      velocity and total energy AFTER this routine is called.
c
c
      K = KSTART2
      nzones = KDO2

      call ppm98_deul_update_cart
     ^     (     dx(K),   dvol(K),      dm(K),    rho(K),     u1(K),
     ^      e_total(K), dmassl(K),  dmom1l(K),   pavl(K), detotl(K),
     ^        rhonu(K),   u1nu(K), e_totnu(K),   dmnu(K),    smlrho,
     ^          smallp,    smallu,     smalle,        dt,    nzones)


c     no fake zones necessary in the above routine, so the "K" count is
c     unchanged.  The updated zone-averages are valid on the range
c     kstart2,kend2, and could be examined by:
c
c      print *,' Finished with the hydro step.'
c      print *,'    i,  rhonu(i),   u1nu(i), e_totnu(i), pnu(i)'
c      do i=kstart2,kend2
c           ekk = 0.5 * (u1nu(i)**2)
c           eii = e_totnu(i) - ekk
c           peenu = (eosgam-1.0) * rhonu(i)*eii
c         write(6,666) i, rhonu(i), u1nu(i), e_totnu(i), peenu
c      enddo


c   Prepare for advecting the transverse velocities U2, U3:
c

c
c     First, compute the advected volume fractions.  Historicaly, PPM
c     has used a spatial fraction: dxl/dx,  when logically it would
c     seem that a mass fraction should be used.  I got no satisfactory
c     explaination from Woodward, as usual, other than "that's how I did
c     it."  For now, I follow this tradition.  I expect to change it.  I 
c     also agree with Woodward that the difference will probably be 
c     negligable.  I like consistancy. And in this formulation, advected
c     mass fractions are no more expensive than spatial or volume fractions.
c     For know, and in Cartesian coordinates (spatial = volume) we use
c     dvoll/dvol

         do i=KSTART2,KEND2+1

             if( dvoll(i) .le. 0.0 ) then
                dvolfl(i) = -dvoll(i) * dvolinv(i)
             else
                dvolfl(i) = dvoll(i) * dvolinv(i-1)
             endif
                dvolfl(i) = min (1.0e+00, dvolfl(i))
          enddo
         do i=KSTART2,KEND2
             if ( dmnu(i) .ne. 0.0 ) then
                dmnuinv(i) = 1.0 / dmnu(i)
             else
                dmnuinv(i) = 0.0
             endif
          enddo
c
c   Now passively advect the transverse volocity, using mass weightings:

c     nbdy_za is the number of fake zones of zone-averaged data
      nbdy_za = 4

c     nbdy_zl is the number of fake zones of zone-interface data
      nbdy_zl = 1
      NZONES =  KDO1 -2*nbdy_za


c
c   PPM98_ADVECT_PASSIVE0 performs the entire advection step.  Its
c
c            steps are:
c                   1) Interpolate, with discontinuity detection and
c                          monotonicity constraints. (Returns UNSMUTUT
c                          and STEEPUT)
c                   2) Integrate the piecewise parabolic interpolations
c                          over the appropriate fraction of each zone
c                          indicated by DVOLFL. (As noted above, logically
c                          it would appear that this should be DMASFL, the 
c                          mass fraction.)
c                   3) Construct weighted fluxes from step 2) using the
c                          provided interface weightings, here dmassl,
c                          the mass advected at the lefthand zone interface.
c                          (Returns DMOMTL)
c                   4) Update the zone averages by conservatively differencing
c                           the interface fluxes from step 3) and using the
c                           old zone mass (DM) and the inverse of the new
c                           zone volume (DMNUINV) as weightings.  
c                           Eulerian scheme the old and new zone volumes are
c                           the same and DVOLINV(i) = 1.0/DVOL(i).) (Returns
c                           PAQNU).
c
c                     These individual steps are also available as separate
c                           modules.

c
c      Thus the following routine requires various kinds of input:
c             1) zone averages of UT for interpolation.  These will
c                require as many boundary cells as were required for
c                the call to DO_PPM98_DEUL0_FLUX_GAMMA.
c             2) Flux factors DVOLFL, DMASSL,  required on a smaller
c                range than the zone averages in 1.  This is because 
c                they will be applied only at interfaces where  valid 
c                interpolations exist to be integrated over the indicated
c                advected zone fractions, and the interpolation routines 
c                require 3 cells usually.  The number of interface boundary 
c                zone is passed here as a check.  At this point, it should
c                usually be 1.
c             3) the zone update weightings (DM, DMNUINV=1/DMNU), which
c                require no zone boundary zones.  They need be valid only 
c                on the range 1:nzones, the same as the updated 
c                transverse velocity UTNU.
c
c      All these arrays begin with the same offset K, for sanity.  This
c      ensures that when refferencing array elements of "I", all such
c      references are to the same physical zone.  Chaos wold result otherwise.

c     nbdy_za is the number of fake zones of zone-averaged data
      nbdy_za = 4

c     nbdy_zl is the number of fake zones of zone-interface data
      nbdy_zl = 1
      NZONES =  KDO1 -2*nbdy_za


      K = KSTART1

            call  ppm98_advect_passive0
     ^             (dvolfl(K),    dmassl(K),    dm(K),
     ^             dmnuinv(K), u2(K), smalldu(K),
     ^             unsmutu2(K),   steepu2(K),     dmom2l(K),
     ^             u2nu(K),      nzones,      nbdy_za )

            call  ppm98_advect_passive0
     ^             (dvolfl(K),    dmassl(K),    dm(K),
     ^             dmnuinv(K), u3(K), smalldu(K),
     ^             unsmutu3(K),   steepu3(K),     dmom3l(K),
     ^             u3nu(K),      nzones,      nbdy_za )

            call  ppm98_advect_passive0
     ^             (dvolfl(K),    dmassl(K),    dm(K),
     ^             dmnuinv(K), sk(K), smallsk(K),
     ^             unsmutuk(K),   steepuk(K),     drhokl(K),
     ^             sknu(K),      nzones,      nbdy_za )

ccc
ccc   Update the transverse kinetic energy and add it to the new total E

      do i=KSTART2 ,KEND2+1
         if (dmassl(i) .eq. 0.0) then
	    dektl(i) = 0.0
         else
	    dektl(i) = 0.5 * (dmom2l(I)*dmom2l(I) + 
     ^                                dmom3l(i)*dmom3l(I)) / dmassl(i)
         endif
      enddo

      do i=KSTART2 ,KEND2 
	 ekt = 0.5 * (u2(I)*u2(I) + u3(i)*u3(I))
         ektnu = ( ekt * dm(i) + dektl(i) - dektl(i+1))*dmnuinv(i)
	 e_totnu(i) =  e_totnu(i) + ektnu
      enddo

C  **********************************************************
C  *                                                        *
C  *  Done with the hydro step, begin diffusion step        *
C  *                                                        *
C  *********************************************************

c
c   The computation of the artifical diffusion is somewhat expensive, and
c   it is only applied in shocks.  Clearly, a shock (particularly shocks
c   requiring additional dissapation) occur in "unsmooth" regions. of
c   the data.  Thus the unsmoothness results from the above interpolations
c   may be used to reduct the additional work caused by diffuse.  Thus
c   we will construct a quantity UNSMOOTH which is the combination of
c   UNSMUTUT and UNSMUP.  Only where UNSMOOTH is non-0.0e+00 will the 
c   diffusion calculation be performed.  PLEASE NOTE THE INCLUSION OF
c   THE TRANSVERS VELOCITY UNSMOOTHNESS MEASURE.  The difusion program
c   attemps to determine the orientation of the shock (if any) relative
c   to the grid.  Without the transverse velocity unsmoothness measure
c   we would miss shock which were all or substantially parallel to
c   the direction of this particular pass, and no addition diffusion
c   would be added here.  This would be undesirable.

      do i=KSTART2 ,KEND2 
         unsmooth(i) =  unsmutu2(i) + unsmutu3(i) + unsmup(i)
      enddo


c
c      zone centerred difusion velocities are required on the range
c      [0,N_REAL_ZONES_IN+1].  From these, diffusive interfaces fluxes
c      for zone on the range [1,N_REAL_ZONES_IN+1] can be computed.  After
c      differencing these fluxes, we will have final zone-averaged
c      values on the range [1:N_REAL_ZONES_IN].  This discussion of
c      zone counts assumes that the number of fake zones sent to this
c      routine is exactly equal to the minimum requirement for this
c      program.  In practice, situations may exist, where the
c      number of fake zones sent may be larger than required.  Therefore
c      all counts and offset calculations are always done relative to
c      the initial values.
c
c      BUT REALIZE THAT PPM98_EUL_DIFFUSE_1D OPERATES ON THE INITIAL
C      DATA.
c
c      Therefore, in the call to difuse below the starting location K is
c      K = KSTART2 - nbdy_df (rather than just KSTART2).
c


      KSTARTDF = KSTART2  
      KENDDF = KEND2  
      KDODF = KENDDF - KSTARTDF + 1 
c
c



c     STEP 1:  Get zone-centered difusion velocities
      
      K=KSTARTDF
      nzones = KDODF 

c   PPM98_DIFUSE0_GAMMA_3D  obtains zone centered Eulerian diffusion velocities.c                  (It's like sausage, you don't really want to know it's done)
c
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c u2_pen(o,ky-2, kz)    <-->  u2zbb
c u2_pen(o,ky-1, kz)    <-->  u2zb
c u2_pen(o,ky  , kz)    <-->  u2
c u2_pen(o,ky+1, kz)    <-->  u2zt
c u2_pen(o,ky+2, kz)    <-->  u2ztt
c u2_pen(o,ky  , kz-2)  <-->  u2zmm
c u2_pen(o,ky  , kz+2)  <-->  u2zpp
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      call  ppm98_eul_difuse0_gamma_3d (

     1  u1_pen(k  ,ky-2,kz  ), u1_pen(k-2,ky  ,kz  ), u1_pen(k-1,ky  ,kz  ),
     1  u1_pen(k,  ky  ,kz  ), u1_pen(k+1,ky  ,kz  ), u1_pen(k+2,ky  ,kz  ),
     1  u1_pen(k,  ky+2,kz  ), u1_pen(k,  ky  ,kz-2), u1_pen(k,  ky  ,kz+2),

     1  u2_pen(k  ,ky-2,kz  ), u2_pen(k  ,ky-1,kz  ), u2_pen(k-2,ky  ,kz  ),
     1  u2_pen(k,  ky  ,kz  ), u2_pen(k+2,ky  ,kz  ), u2_pen(k,  ky+1,kz  ),
     1  u2_pen(k,  ky+2,kz  ), u2_pen(k,  ky  ,kz-2), u2_pen(k,  ky  ,kz+2),

     1  u3_pen(k  ,ky-2,kz  ), u3_pen(k-2,ky  ,kz  ), u3_pen(k,  ky  ,kz  ),
     1  u3_pen(k+2,ky  ,kz  ), u3_pen(k,  ky+2,kz  ), u3_pen(k,  ky  ,kz-2),
     1  u3_pen(k,  ky  ,kz-1), u3_pen(k,  ky  ,kz+1), u3_pen(k,  ky  ,kz+2),

     1  pp_pen(k  ,ky-2,kz  ), pp_pen(k-2,ky  ,kz  ), pp_pen(k  ,ky  ,kz  ),
     1  pp_pen(k+2,ky  ,kz  ), pp_pen(k  ,ky+2,kz  ),
     1  pp_pen(k  ,ky  ,kz-2), pp_pen(k  ,ky  ,kz+2),

     1  ro_pen(k  ,ky-2,kz  ), ro_pen(k-2,ky  ,kz  ), ro_pen(k  ,ky  ,kz  ),
     1  ro_pen(k+2,ky  ,kz  ), ro_pen(k  ,ky+2,kz  ),
     1  ro_pen(k  ,ky  ,kz-2), ro_pen(k  ,ky  ,kz+2),

     1       unsmooth(k), difuse_vel(k),  shockd(k), eosgam, smallu, nzones)

cc     call  ppm98_eul_difuse0_gamma_3d (
cc    ^                u1zbb(K),   u1(k-2),  u1(k-1),   u1(K), u1(k+1), u1(k+2), 
cc    ^     u1ztt(k),  u1zmm(k),  u1zpp(k), u2zbb(K), u2zb(k),   u2(K-2), u2(k), 
cc    ^     u2(k+2),  u2zt(k), u2ztt(K),  u2zmm(k), u2zpp(k), u3zbb(K), u3(k-2),
cc    ^        u3(K), u3(k+2), u3ztt(k),  u3zmm(K),  u3zm(k), u3zp(k), u3zpp(k),
cc    ^     pzbb(K),    p(k-2), p(K), p(k+2), pztt(k), pzmm(K), pzpp(k), 
cc    ^    rhozbb(K), rho(k-2), rho(K), rho(k+2), rhoztt(k),rhozmm(K),rhozpp(k),
cc    ^       unsmooth(K), difuse_vel(K),  shockd(K), eosgam, smallu,  
cc    ^               nzones)

c   update the zone counts


      KSTARTDF1 = KSTARTDF
      KENDDF1 = KENDDF
      KDODF1 = KDODF 


c   At this point, we should have 1.0e+00 extra diffusion velocity
c   at each end of the working array.  These are ZONE-CENTERED
c   velocities, from which INTERFACE centered fluxes will
c   be constructed.  This construction requires this extra fake zone.
c
c      STEP 2: Apply diffuseion
c
c              Get the diffusive zone interface fluxes
c              from the zone centerred difusion velocities
c              Like the hydro step above, this is divided into
c              longitudinal and transverse parts.

      nbdy_difuse_fluxes = 1
      K = KSTARTDF1
      nzones = KDODF1 - 2*nbdy_difuse_fluxes



      nbdy = 1
      K = KSTARTDF1
      nzones = KDODF1 - 2*nbdy


c get the longitudinal hydrodynamic fluxes

      call  ppm98_difuseflux_hydro_cart
     ^    (      dx(K), difuse_vel(K),    
     ^        rhonu(K),       u1nu(K), e_totnu(K), 
     ^        dvoll(K),     dmasll(K),  dmasrl(K), 
     ^       dmassl(K),     dmom1l(K),  detotl(K), 
     ^       courno(K),            dt,     nzones,      
     ^        nbdy )
c
c    get the transverse momentum diffusion fluxes
c
c   PPM98_DIFUSEFLUX_PASSIVE uses the the leftward  and rightward
c     diffusion volumes (MASLL, MASRL) computed above to compute
c     lefthand zone interface mass weighted interface difusion momenta dmomtl.


               call  ppm98_difuseflux_passive
     ^         (u2nu(K), dmasll(K), dmasrl(K), dmom2l(K), nzones,  nbdy)

               call  ppm98_difuseflux_passive
     ^         (u3nu(K), dmasll(K), dmasrl(K), dmom3l(K), nzones,  nbdy)

               call  ppm98_difuseflux_passive
     ^         (sknu(K), dmasll(K), dmasrl(K), drhokl(K), nzones,  nbdy)



      KSTARTDF = KSTARTDF1 + nbdy
      KENDDF = KENDDF1 - nbdy
 
      courdf2 = courno(KSTARTDF)
      do i=KSTARTDF+1,KENDDF
         courdf2 = max (courdf2,courno(i))
      enddo


c      Next apply the diffusive fluxes
c
c       (Final update:, whew!)
c

c     We occasionally get underflows by allowing absurdly small
c     velocities to be computed.  Therefore, we will set such
c     velocities to 0.0e+00 below.
c
      ssmalu = small * smallu
      ssmlu2 = ssmalu * ssmalu
c
      gamma1 = eosgam - 1.0


ccccc      write (6,*) 'i, u1nu(i), u2nu(i), u3nu(i)'

      do 9100   i = KSTARTDF,KENDDF
c
c         save old dmnu to cherish later:
          dmnew(i) = dmnu(i)
ccccccc          if ((dvoll(i) + dvoll(i+1)) .ne. 0.0e+00)   then
c            dmnew(i) = rhonu(i) * dvol(i)
             dmnu(i) = dmnew(i) + dmassl(i) - dmassl(i+1)
             rhonu(i) = dmnu(i) * dvolinv(i)
             dmnui = 1.0e+00 / dmnu(i)
             u1nu(i) = (u1nu(i) * dmnew(i)  +  dmom1l(i) - dmom1l(i+1)) * dmnui
             u2nu(i) = (u2nu(i) * dmnew(i)  +  dmom2l(i) - dmom2l(i+1)) * dmnui
             u3nu(i) = (u3nu(i) * dmnew(i)  +  dmom3l(i) - dmom3l(i+1)) * dmnui
             sknu(i) = (sknu(i) * dmnew(i)  +  drhokl(i) - drhokl(i+1)) * dmnui
             e_totnu(i) = (e_totnu(i) * dmnew(i) + 
     1                                         detotl(i) - detotl(i+1)) * dmnui
ccccccccc         endif
c_dhp    why?  Because - bke
	 if ((u1nu(i) * u1nu(i)) -ssmlu2 .lt. 0)   u1nu(i) = 0.0
         if ((u2nu(i) * u2nu(i)) -ssmlu2 .lt. 0)   u2nu(i) = 0.0
         if ((u3nu(i) * u3nu(i)) -ssmlu2 .lt. 0)   u3nu(i) = 0.0
c_dhp    why?
         eknu = 0.5 * (u1nu(i) * u1nu(i)  +  u2nu(i) * u2nu(i) + 
     ^                                       u3nu(i) * u3nu(i)  )
         e_intnu = max (smalle, (e_totnu(i) - eknu))
         pnu(i) = gamma1 * rhonu(i) * e_intnu
c          print *,' rpu=',i,rhonu(i),pnu(i),u1nu(i)

ccccc
cccccc for non-self consistent energy-momentum fluxs -- but at least no gross effects.
ccccc         u1old = u1nu(i)
ccccc         u2old = u2nu(i)
ccccc         u3old = u3nu(i)

c         u1nu(i) = (u1nu(i) * dmnew(i) + TAU_L11(i) - TAU_L11(i+1)) * dmnui
c         u2nu(i) = (u2nu(i) * dmnew(i) + TAU_L12(i) - TAU_L12(i+1)) * dmnui
c         u3nu(i) = (u3nu(i) * dmnew(i) + TAU_L13(i) - TAU_L13(i+1)) * dmnui

ccccc         write (6,999) i, u1old, u1nu(i), u2old, u2nu(i), u3old, u3nu(i)
ccccc999      format ( i5, 3('   ', 2f10.5))

9100  continue
c
ccccc      if(dmnu(1) .gt. 0.) stop
c
c

      courmx = max (courmx, courdf2 )


CF90      deallocate (       niterl, cesq, clagl)
CF90       deallocate ( rhol,     pl,     ul,
CF90     1              rhor,     pr,     ur,
CF90     2              drho,     dp,     du,
CF90     3              rho6,     p6,     u6)
CF90       deallocate (ullrp,  pllrp,  urlrm, prlrm)
CF90      deallocate       (courno,   difuse_vel,
CF90     ^      unsmad,  unsmup,   stpadb,   ceul,
CF90     ^      rhoavl,    pavl,     uavl,  etotavl,
CF90     ^     e_totnu,       dx,    dxinv,  
CF90     ^      difusl,  dmassl,   dmom1l,  detotl,
CF90     ^      dmom2l,   dmom3l, 
CF90     ^      dmasll, dmasrl,     shockd)
CF90      deallocate (dxl,   dvol, dvoll,   unsmooth)
CF90      deallocate (dm,   dmnu, dmnew )
CF90      deallocate (dvolinv, dminv, dmnuinv)
CF90      deallocate ( smalldu, dvolfl, unsmutu2, unsmutu3 )
CF90      deallocate ( steepu2, steepu3, dektl )

      return
      end
c
c
 
