!========================================================================
! 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
!                      physflux - Compute the physical flux
!
! 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, quadrant problem
! Contains:
! Revision History: Ver. 1.0 Oct. 2000 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,gamma12,gamma3,xQ,yQ
  !  - Initial values in quadrants
  REAL (KIND=qPrec), DIMENSION(4) :: pQ,rhoQ,uQ,vQ  
  
  !=======
  CONTAINS
  !=======
  
    SUBROUTINE setprob      
      INTEGER i
      ! Read gas constant from file
      OPEN(UNIT=7,FILE='setprob.data',STATUS='old',FORM='formatted')      
      READ(7,*) gamma
      gamma1 = gamma - 1.; gamma12=0.5*gamma1; gamma3=3.-gamma;
      READ(7,*) xQ,yQ
      DO i=1,4
        READ(7,*)pQ(i),rhoQ(i),uQ(i),vQ(i)
      END DO
      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)
      ! Interface declarations      
      TYPE (NodeInfo) :: Info  ! Data associated with this grid  
      ! Internal declarations
      INTEGER i,j,iq
      REAL (KIND=xPrec) :: x,y      
      !      
      x=Info%Xlower(1)+Info%dX(1)/2.0
      DO i=1,Info%mX(1)
        y=Info%Xlower(2)+Info%dX(2)/2.0
        DO j=1,Info%mX(2)
          iq=1          
          IF (x <  xQ .AND. y >= yQ) iq = 2
          IF (x <  xQ .AND. y <  yQ) iq = 3
          IF (x >= xQ .AND. y <  yQ) iq = 4      
          Info%q(i,j,1,1,1) = rhoQ(iq)
          Info%q(i,j,1,1,2) = rhoQ(iq)*uQ(iq)
          Info%q(i,j,1,1,3) = rhoQ(iq)*vQ(iq)
          Info%q(i,j,1,1,4) = pQ(iq)/gamma1 + 0.5d0*rhoQ(iq)*(uQ(iq)**2 + vQ(iq)**2)
          y=y+Info%dX(2)
        END DO
        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 physical fluxes for the 2D Euler equations
!
!    q_t + f(q)_x + g(q)_y = 0
!
! -----------------------------------------------------------------------      
    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 NodeInfoDef.f90 for definitions			    
      TYPE (NodeInfo) :: Info  ! Data associated with this grid
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q,f,s
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: A      
      ! Internal declarations
      INTEGER, PARAMETER :: xflux=1, yflux=2
      INTEGER, PARAMETER :: xspeeds=1, yspeeds=2
      INTEGER, PARAMETER :: xJacobian=1, yJacobian=2
      INTEGER, SAVE :: mbc,mx,iError
      INTEGER, POINTER, DIMENSION(:) :: I    
      REAL (KIND=qPrec), SAVE, ALLOCATABLE, DIMENSION (:) :: rho,invrho,p,u,v,u2v2,c,H      
      ! Index range
      I=>Info%I1D      
      SELECT CASE (irequest)
        CASE (Initialize)
          ! Current number of interior cells, ghost cells
          mx=Info%mXnow; mbc=Info%mbc
          ! Allocate some convenient local variables
          ALLOCATE(rho(1-mbc:mx+mbc),invrho(1-mbc:mx+mbc),p(1-mbc:mx+mbc),c(1-mbc:mx+mbc), &
                   u(1-mbc:mx+mbc),v(1-mbc:mx+mbc),u2v2(1-mbc:mx+mbc),H(1-mbc:mx+mbc),STAT=iError)
          IF (iError /=0) THEN
            PRINT *,'Cannot allocate work arrays in physflux'
            STOP
          END IF          
        CASE (Finalize)
          ! Release work arrays      
          DEALLOCATE(rho,invrho,p,c,u,v,u2v2,H,STAT=iError)
          IF (iError /=0) THEN
            PRINT *,'Cannot deallocate work arrays in physflux'
            STOP
          END IF
        CASE (RequestFluxes)
          rho=q(I,1); 
          IF (MINVAL(rho)<=0.) THEN
            PRINT *,'Error in physflux: negative density'
            print *,I
            print *,rho
            print *,indx
            STOP     
          END IF
          invrho=1./rho; u=q(I,2)*invrho; v=q(I,3)*invrho
          u2v2=u**2+v**2; p=gamma1*(q(I,4)-0.5*rho*u2v2); c=gamma*p*invrho
          IF (MINVAL(c)<=0.) THEN
            PRINT *,'Error in physflux: negative or zero sound velocity'	
            STOP
          END IF
          c=SQRT(c)
          SELECT CASE (ixy)
            CASE (xflux) 	    
              f(I,1)=q(I,2)
              f(I,2)=q(I,2)*u+p
              f(I,3)=q(I,2)*v
              f(I,4)=(q(I,4)+p)*u
            CASE (yflux)	    
              f(I,1)=q(I,3)
              f(I,2)=q(I,3)*u
              f(I,3)=q(I,3)*v+p
              f(I,4)=(q(I,4)+p)*v
           END SELECT
        CASE (RequestSpeeds)
           SELECT CASE (ixy)
             CASE (xspeeds)
               s(I,1)=u + c
               s(I,2)=u
               s(I,3)=u
               s(I,4)=u - c
             CASE (yspeeds)
               s(I,1)=v + c
               s(I,2)=v
               s(I,3)=v
               s(I,4)=v - c
             END SELECT
        CASE (RequestJacobians)
          A=0.; H=invrho*(q(I,4)+p)
          SELECT CASE (ixy)	    
            CASE (xJacobian)
              A(I,1,2)=1.
	    
              A(I,2,1)=gamma12*u2v2(I)-u(I)**2; A(I,2,2)=gamma3*u(I); 
              A(I,2,3)=-gamma1*v(I); A(I,2,4)=gamma1

              A(I,3,1)=-u(I)*v(I); A(I,3,2)=v(I); A(I,3,3)=u(I);

              A(I,4,1)=u(I)*(-H(I)+gamma12*u2v2(I)); A(I,4,2)=H(I)-gamma1*u(I)**2
              A(I,4,3)=-gamma1*u(I)*v(I); A(I,4,4)=gamma*u(I)
            CASE (yJacobian)
              A(I,1,3)=1.

              A(I,2,1)=-u(I)*v(I); A(I,2,2)=v(I); A(I,2,3)=u(I);

              A(I,3,1)=gamma12*u2v2(I)-v(I)**2; A(I,3,2)=-gamma1*u(I)
              A(I,3,3)=gamma3*v(I); A(I,3,4)=gamma1
      	    
              A(I,4,1)=v(I)*(-H(I)+gamma12*u2v2(I)); A(I,4,2)=-gamma1*u*v
              A(I,4,3)=H(I)-gamma1*v(I)**2;
              A(I,4,4)=gamma*v(I)
          END SELECT    	        
      END SELECT      
    END SUBROUTINE physflux
              
!*****************
END MODULE problem
!*****************