!========================================================================
! BEARCLAW Boundary Embedded Adaptive Refinement Conservation LAW package
!========================================================================
! (c) Copyright Sorin Mitran, 2000
! Department of Applied Mathematics
! University of Washington
! mitran@amath.washington.edu
! -----------------------------------------------------------------------
! This code may be freely used for educational and research purposes.
! For any other use please contact the author.
! -----------------------------------------------------------------------
! File:             problem.f90
! Purpose:          Definition of problem specific data and procedures
!                   Contains:
!                      setprob  - Initialize problem parameters
!                      qinit    - Initialize field variables
!                      b4step   - Called prior to each time step
!                      setaux   - Initialize auxilliary variable fields
!                      src      - Compute source terms
!                      physflux - Compute the physical flux, related quantities
!                                 Compute wave decomposition
!
! Comments:         This module is problem specific and must be written by the
!                   user. For exampels see the BEARCLAW documentation and also
!                   the examples and applications directories
!                 
!                   The encapsulation of all problem specific data and routines
!                   and data into a module is useful in exposing data that needs
!                   to be shared.
! Application:      Euler equations, shock tube problem
! Contains:
! Revision History: Ver. 1.0 Oct. 2001 Sorin Mitran
! -----------------------------------------------------------------------

!*************
MODULE problem
!*************
  ! User definitions of data structures associated with each node
  USE NodeInfoDef
  IMPLICIT NONE    ! It's safer to require explicit declarations
  SAVE             ! Save module information
  PRIVATE          ! Everything is implicitly private, i.e. accessible only
                   ! to this module. 
  ! These are the public entry points (acessible to other modules)                   
  PUBLIC setprob,qinit,b4step,setaux,src,physflux,problemBC,afterstep,problemIO
  
  ! Problem specific global variables:  
  !  - Problem parameters
  REAL (KIND=qPrec) :: gamma,gamma1
  REAL (KIND=qPrec) :: rhol,rhoul,el,rhor,rhour,er
  
  !
  INTEGER, PARAMETER :: OK=0
  REAL (KIND=qPrec), PARAMETER :: zero=0., half=0.5, one=1., two=2.

  !=======
  CONTAINS
  !=======
  
    SUBROUTINE setprob  
      ! User definitions of data structures associated with each node
      USE NodeInfoDef      
      REAL (KIND=qPrec) :: ul,pl,ur,pr
      ! Set the initial shock tube states
      OPEN(UNIT=7,FILE='setprob.data',STATUS='old',FORM='formatted')  
      READ(7,*) gamma
      gamma1 = gamma - 1.d0
      READ(7,*) rhol,ul,pl
      READ(7,*) rhor,ur,pr
      ! data in left state
      rhoul = rhol*ul
      el = pl/gamma1 + 0.5d0*rhol*ul**2
      ! data in right state
      rhour = rhor*ur
      er = pr/gamma1 + 0.5d0*rhor*ur**2
      CLOSE(7)
    END SUBROUTINE setprob
    
    SUBROUTINE problemBC(Info)
      ! User definitions of data structures associated with each node
      USE NodeInfoDef  
      ! Interface declarations      
      TYPE (NodeInfo) :: Info  ! Data associated with this grid      
      ! Internal declarations
    END SUBROUTINE problemBC

    SUBROUTINE qinit(Info)
      ! User definitions of data structures associated with each node
      USE NodeInfoDef  
      ! Interface declarations      
      TYPE (NodeInfo) :: Info  ! Data associated with this grid  
      ! Internal declarations
      INTEGER i
      REAL (KIND=xPrec), PARAMETER :: xDiscontinuity = 0.5
      REAL (KIND=xPrec) :: x      
      !
      x=Info%Xlower(1)+Info%dX(1)/2.0
      DO i=1,Info%mX(1)
        IF (x<xDiscontinuity) THEN
          !      x y z t m
          Info%q(i,1,1,1,1) = rhol
          Info%q(i,1,1,1,2) = rhoul
          Info%q(i,1,1,1,3) = el
        ELSE
          !      x y z t m
          Info%q(i,1,1,1,1) = rhor
          Info%q(i,1,1,1,2) = rhour
          Info%q(i,1,1,1,3) = er
        END IF
        x=x+Info%dX(1)
      END DO    
    END SUBROUTINE qinit
    
    SUBROUTINE b4step(Info)  
      ! Interface declarations    
      TYPE (NodeInfo) :: Info        
    END SUBROUTINE b4step
    
    SUBROUTINE afterstep(Info)
      ! Interface declarations    
      TYPE (NodeInfo) :: Info
      ! Include operations that need to be carried out after each time step here:
    END SUBROUTINE afterstep
    
    SUBROUTINE problemIO(nframe,tnow,IOrequest,Info,qmax,qmin)
      ! Interface declarations
      INTEGER :: nframe,IOrequest; REAL (KIND=qPrec) :: tnow
      TYPE (NodeInfo), OPTIONAL :: Info      
      REAL (KIND=qPrec), OPTIONAL, DIMENSION(:) :: qmax,qmin
      INTEGER, PARAMETER :: UserBeforeGridIO=-1, UserAfterGridIO=-2
      INTEGER, PARAMETER :: UserIOMinMax=0
      ! This routine provides for output of computation snapshots in user-defined
      ! formats.
      ! IOrequest is a function parameter with values:
      ! -1 UserBeforeGridIO - indicates routine is being called prior to parsing
      !                       of tree of grids. This is where output files
      !                       should be opened
      ! -2 UserAfterGridIO  - indicates routine is being called after parsing
      !                       of tree of grids. This is where output files
      !                       should be closed
      ! -3 UserOutputq      - indicates routine should output the field
      !                       variables in user specified format
      ! 0  UserIOMinMax     - routine is being called after qmin and qmax have
      !                       been determined. These are the min/max field
      !                       variable values. If output of quantities different
      !                       from the field variables is desired, qmin &
      !                       qmax should be defined with respect to those quantities
      ! 1,2 ...             - routine is being called to output field variable
      !                       1,2,... on the current grid in conjunction with
      !                       HDF support routines in BEARIO. Only valid for
      !                       UserHDF output style      
    END SUBROUTINE problemIO
    
    SUBROUTINE setaux(Info)  
      ! Interface declarations    
      TYPE (NodeInfo) :: Info        
    END SUBROUTINE setaux
    
    SUBROUTINE src(Info)
      ! Interface declarations    
      TYPE (NodeInfo) :: Info        
    END SUBROUTINE src
    
! -----------------------------------------------------------------------
!
! Compute the Riemann problem solution for the 1D Euler equations
!
!    q_t + f(q)_x = 0
!
! Routine should return wave decomposition resulting from solving Riemann
! problem at interfaces j=2-mbc:mxnow+mbc, between left states in
! q1D(1-mbc:mxnow+mbc-1) and right states in q1D(2-mbc:mxnow+mbc)
! -----------------------------------------------------------------------
    SUBROUTINE physflux(ixy,indx,irequest,Info,q,f,s,A)
      ! User definitions of data structures associated with each node
      USE NodeInfoDef  
      ! Interface declarations    
      INTEGER ixy           ! The dimension along which to solve the normal Riemann problem    
      INTEGER indx(MaxDims) ! Indices of this slice, indx(ixy) has no significance, other
                            ! indices give position of this 1D slice in index space
      INTEGER irequest      ! Request code, see NodeInfoGlobal.f90 for definitions
                            ! Additional request codes may be defined in solver modules
      TYPE (NodeInfo) :: Info  ! Data associated with this grid
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q,f,s
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: A
      ! Internal declarations      
      LOGICAL, PARAMETER :: efix = .TRUE.
      INTEGER mbc,mx,nq,n,NrVars,iError,j
      ! Shorter local names
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q1D
      ! Some additional local fields
      ! Cell center quantities
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION (:) :: sqrtrho, &
           pres,rhoRsq2,delta,coef,rho,u,c,h      
      !
      ! Current number of interior cells, boundary cells, field variables, waves
      mx=Info%mXnow; mbc=Info%mbc; NrVars=Info%NrVars
      ! Associate local names with Info components
      q1D=>Info%q1D      
      SELECT CASE (irequest)
        CASE (Initialize)
          ! Allocate local work space
          ALLOCATE(sqrtrho(1-mbc:mx+mbc),pres(1-mbc:mx+mbc),              &
                   h(1-mbc:mx+mbc),c(1-mbc:mx+mbc),rho(1-mbc:mx+mbc),     &
                   u(1-mbc:mx+mbc),     &
                   STAT=iError)
          IF (iError/=OK) THEN
            PRINT *,'Cannot allocate work space in problem module, physflux'
            STOP
          END IF
        CASE (Finalize)
          ! Release local work space
          DEALLOCATE(sqrtrho,pres,h,c,rho,u,STAT=iError)
          IF (iError/=OK) THEN
            PRINT *,'Cannot allocate work space in problem module, physflux'
            STOP
          END IF
        CASE (RequestFluxes)
          IF (MINVAL(q1D(1-mbc:mx+mbc,1))<=zero) THEN
            PRINT *,'Error: Negative or zero density in physflux.'
            PRINT '(1x,a6,1x,32(f7.4,1x))','rho=  ',q1D(1-mbc:mx+mbc,1)
            STOP
          END IF
          ! Compute cell centered quantities
          rho(1-mbc:mx+mbc) = q1D(1-mbc:mx+mbc,1)
          u(1-mbc:mx+mbc) = q1D(1-mbc:mx+mbc,2)/rho(1-mbc:mx+mbc)
          sqrtrho(1-mbc:mx+mbc) = SQRT(rho(1-mbc:mx+mbc))
          pres(1-mbc:mx+mbc) = gamma1 *                    &
            ( q1D(1-mbc:mx+mbc,3) - half*q1D(1-mbc:mx+mbc,2)*u(1-mbc:mx+mbc) )
          c(1-mbc:mx+mbc) = gamma*pres(1-mbc:mx+mbc)/rho(1-mbc:mx+mbc)
          IF (MINVAL(c(1-mbc:mx+mbc))<=zero) THEN
            PRINT *,'Error: Non-physical centered sound velocity in physflux.'
            PRINT '(1x,a6,1x,32(f7.4,1x))','c**2 =  ',c(1-mbc:mx+mbc)
            STOP
          END IF
          c(1-mbc:mx+mbc) = SQRT(c(1-mbc:mx+mbc))
          h(1-mbc:mx+mbc) = (q1D(1-mbc:mx+mbc,3) + pres(1-mbc:mx+mbc)) / &
                            rho(1-mbc:mx+mbc)
          ! Compute the physical fluxes at cell centers
          f(1-mbc:mx+mbc,1)=q1D(1-mbc:mx+mbc,2)
          f(1-mbc:mx+mbc,2)=f(1-mbc:mx+mbc,1)*u(1-mbc:mx+mbc)+pres(1-mbc:mx+mbc)
          f(1-mbc:mx+mbc,3)=f(1-mbc:mx+mbc,1)*h(1-mbc:mx+mbc)        
        CASE (RequestSpeeds)
          ! Not coded. Fill in later
        CASE (RequestJacobians)
          ! Not coded. Fill in later
        CASE DEFAULT          
          PRINT *,'Error: the solver module Scheme.f90 emitted a request that'
          PRINT *,'is not implemented in the problem module problem.f90'
          STOP
      END SELECT      
    END SUBROUTINE physflux
              
!*****************
END MODULE problem
!*****************