!========================================================================
! 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:             classicbear.f90
! Purpose:          Coding of some classic schemes as the time stepping
!                   routine for BEARCLAW
! Contains:
! Revision History: Ver. 1.0 Aug. 2001 Sorin Mitran
! -----------------------------------------------------------------------
! Comments:
!   1) This is an example implementation of the BEARez TimeStep function.
!      The basic task a TimeStep function must perform is to provide an
!      updated Info%qnew reflecting field values at t+dt, starting
!      from the field values at time t given in Info%qold.
! -----------------------------------------------------------------------

!************
MODULE scheme
!************
  ! User definitions of data structures associated with each node
  USE NodeInfoDef
  ! Problem specific routines (setprob,qinit,b4step,src,physflux)
  USE problem  
  
  IMPLICIT NONE
  SAVE
  PRIVATE		! Everything is implicitly PRIVATE  
  PUBLIC TimeStep,AllocSchemeFields,DeAllocSchemeFields,ReadSchemeRootData
  
  ! Global variables
  INTEGER, PARAMETER :: Centered=0,LaxFriedrichs=1,Upwind=2,LaxWendroff=3,Richtmyer=4
  INTEGER, PARAMETER :: BeamWarming=5,Fromm=6,minmod=7,superbee=8,MC=9,vanLeer=10
  INTEGER, PARAMETER, DIMENSION(0:10) :: AlgStages = (/ 1,1,1,1,2,1,1,1,1,1,1 /)
  
  REAL (KIND=xPrec), PARAMETER :: zero=0.0
  
  ! Error codes 
  PUBLIC  err_BadCFL
  INTEGER, PARAMETER :: err_OK = 0
  INTEGER, PARAMETER :: err_BadCFL = 1001  ! Error codes > 1000 are considered warnings  
    
  !-------         
  CONTAINS  
  !-------
    
  INTEGER FUNCTION TimeStep(Info,TimeStepParam)
    ! Interface declarations
    TYPE (NodeInfo) :: Info
    TYPE (FuncParam) :: TimeStepParam    
    !
    INTEGER iError
    TimeStep=err_OK
    ! Check for inactive nodes that have had their field released
    ! to conserve memory but are not yet eliminated from the tree
    ! (Elimination from tree occurs after SynchLevels)
    IF (.NOT. Info%FieldsAllocated) RETURN    
    IF (Info%tnow<TimeStepParam%tend) THEN
      Info%dt=TimeStepParam%tend-TimeStepParam%tstart    
      ! Save current field variables
      Info%qold=Info%q
      ! Advance the grid in time
      iError=advance(Info)      
      TimeStep=iError
      Info%NewGrid=.FALSE.; Info%ValidData=.FALSE.            
    ELSE
      ! Shouldn't ever see this message
      PRINT *,'Grid at level ',Info%level,' already at tnow=',Info%tnow,' >= tend=',TimeStepParam%tend
    END IF
  END FUNCTION TimeStep
  
  INTEGER FUNCTION advance(Info)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info    
    ! Internal declarations
    INTEGER nstage
    REAL (KIND=xPrec) :: dt    
    ! Save grid time step 
    dt=Info%dt    
    !=================
    ! Take a time step    
    !=================
    ! Notes: 1) Boundary conditions have been set by the AMR procedure using the
    !           ghost cell approach. The ghost cell region is mbc*r cells wide
    !           with r the refinement ratio and mbc the standard CLAWPACK
    !           ghost cell region width.
    !        2) r time steps will be taken between boundary condition updates                
    !  --------------------------------------------------------
    !  1. Call user-supplied routine which might set aux arrays
    !     for this time step, for example.
    !  --------------------------------------------------------
    CALL b4step(Info)        
    !  ----------------------------------------
    !  2. Take one step on the conservation law
    !  ----------------------------------------    
    Info%NrStages=AlgStages(Info%method(8))
    DO nstage=1,Info%NrStages     ! Some algorithms require intermediate stages
      Info%nstage=nstage
      CALL beforestage(Info)      
      CALL step(Info)
      CALL afterstage(Info)       ! Multistage algorithms typically require some
                                  ! storage gymnastics
    END DO
    ! ---------------------------------------------------
    !  3. Suggest time step to be first tried on next call
    ! ---------------------------------------------------    
    !IF (Info%cfl>0.d0) Info%dt=dt*Info%cflv(2)/Info%cfl    
    !  ---------------------------------------------------
    !  4. Check to see if the Courant number was too large
    !  ---------------------------------------------------
    IF (Info%cfl .gt. Info%cflv(1)) THEN
      advance = err_BadCFL   ! Signal the bad CFL number error
      print *,'Bad cfl =',Info%cfl      
      !RETURN                  ! Let AMR figure it out ...
    END IF    
    !  --------------------------------------------
    !  5. Time step is accepted. Apply source terms
    !  --------------------------------------------    
    CALL src(Info)          
    advance=err_OK
  END FUNCTION advance
  
  SUBROUTINE beforestage(Info)
    ! Interface declarations
    TYPE (NodeInfo) :: Info
    
  END SUBROUTINE beforestage
  
  ! Carry out a time step
  SUBROUTINE step(Info)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info    
    ! Internal declarations    
    INTEGER j,iError
    INTEGER n,n1,indx(MaxDims),i1,i2,i3,i4
    INTEGER ilow(MaxDims),ihigh(MaxDims),inlow,inhigh
    INTEGER nlow(MaxDims),nhigh(MaxDims)
    INTEGER mbc,meqn,nDim,nstep,nsteps    
    !
    Info%cfl=zero
    nDim=Info%nDim; nstep=Info%AMRStep; nsteps=Info%AMRSteps;    
    ! Figure out what the "interior" domain is in this AMR subiteration
    mbc=Info%mbc
    ilow(1:nDim)=1-(nsteps-nstep)*mbc
    ihigh(1:nDim)=Info%mX(1:nDim)+(nsteps-nstep)*mbc    
    ! Perform sweeps along coordinate directions
    DO n=1,nDim    ! This is the dimension in the edge normal direction
      ! Loop over the transverse dimensions 
      ! BEARCLAW works up to 4D, so we have up to 3 transverse dimensions.
      ! Loops over the maximum 4 dimensions are coded; loops corresponding
      ! to the normal dimension or to dimensions not included in this problem
      ! are only executed once.
      !
      ! Figure out over which dimensions we're looping:
      ! 1. First assume that we're in 1D and all outer loops have just one pass
      nlow=1; nhigh=1
      ! 2. Now modify loop limits to reflect other dimensions
      DO n1=1,nDim              ! loop limits for dimensions>nDim remain 1,1
        IF (n1 /= n) THEN       ! loop limits for normal dimension remain 1,1
          nlow(n1)=ilow(n1)-1   ! loop limits for active transverse dimensions
          nhigh(n1)=ihigh(n1)+1 ! correspond to interior plus one row of ghost cells
        END IF
      END DO
      Info%mXnow=ihigh(n)-ilow(n)+1; inlow=ilow(n)-mbc; inhigh=ihigh(n)+mbc
      ! Set up an index array with current interior domain + ghost cells along 1D slice
      ! We need to allocate this with the current interior/ghost domain sizes so that
      ! the user routines (in problem module) have the current index range to work with
      ! We're redeclaring one of the components of Info that is common to all BEARCLAW
      ! applications, so we must first release previous memory.
      DEALLOCATE(Info%I1D); NULLIFY(Info%I1D)
      ALLOCATE(Info%I1D(1-mbc:Info%mXnow+mbc),STAT=iError)
      IF (iError /= 0) THEN
        PRINT *,'Cannot allocate index work array in step'
      STOP
      END IF
      Info%I1D=0
      DO j=1-mbc,Info%mXnow+mbc
        Info%I1D(j)=j
      END DO
      DO i1=nlow(1),nhigh(1)
        indx(1)=i1
        DO i2=nlow(2),nhigh(2)
          indx(2)=i2
          DO i3=nlow(3),nhigh(3)
            indx(3)=i3
            DO i4=nlow(4),nhigh(4)              
              indx(4)=i4              
              ! Compute the fluxes              
              CALL numflux(n,indx,inlow,inhigh,Info)                  
              ! Compute qnew
              CALL updateq(n,indx,inlow,inhigh,Info)
            END DO
          END DO
        END DO
      END DO      
    END DO    
  END SUBROUTINE step
  
  SUBROUTINE afterstage(Info)
    ! Interface declarations
    TYPE (NodeInfo) :: Info
    ! Internal declarations
    INTEGER nq
    SELECT CASE(Info%method(8))
      CASE (Richtmyer)
        IF (Info%nstage==1) THEN
          ! A half step has been taken and stored in Info%q. Place these
          ! values in Info%qold and put the initial values in Info%q
          ! Use the output buffer as temporary storage          
          DO nq=1,Info%NrVars
            Info%qbuf(:,:,:,:)=Info%qold(:,:,:,:,nq)
            Info%qold(:,:,:,:,nq)=Info%q(:,:,:,:,nq)
            Info%q(:,:,:,:,nq)=Info%qbuf(:,:,:,:)
          END DO
        END IF
      CASE DEFAULT
    END SELECT
  END SUBROUTINE afterstage
  
  ! Numerical flux routine:
  !  Compute the numerical flux along dimension n
  !  Slice is at indices indx
  SUBROUTINE numflux(n,indx,inlow,inhigh,Info)
    ! Interface declarations
    INTEGER n,indx(MaxDims),inlow,inhigh
    TYPE (NodeInfo) :: Info
    ! Internal declarations
    INTEGER mbc,mx,m,j,iError,nq,i1,i2,i,indx2(MaxDims),nDim
    INTEGER, ALLOCATABLE, DIMENSION(:) :: Ip1,Im1
    REAL (KIND=qPrec) :: fact,dtdx,dtdx2
    REAL (KIND=qPrec), POINTER, DIMENSION (:) :: dtdx1D
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q1D,fql,fqr,sl,sr,numfql,numfqr,qcenter,qedge
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: fq,s,work
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: Ar,Al,A,Acenter,Aedge
    ! Associate local names with fields of grid data structure
    q1D=>Info%q1D; dtdx1D=>Info%dtdx1D
    fql=>Info%fql; numfql=>Info%numfql; sl=>Info%sl; Al=>Info%Aql
    fqr=>Info%fqr; numfqr=>Info%numfqr; sr=>Info%sr; Ar=>Info%Aqr    
    ! More suggestive names for algorithms which don't differentiate
    ! between right and left edges of a cell
    fq=>Info%fqr; s=>Info%sr; A=>Info%Aqr;
    ! Clean out memory contents
    fql=zero; numfql=zero; sl=zero; Al=zero
    fqr=zero; numfqr=zero; sr=zero; Ar=zero
    ! Current number of interior cellsm ghost cells, field variables
    mx=Info%mXnow; mbc=Info%mbc; nq=Info%NrVars    
    ! Algorithms are programmed for uniform spacing
    dtdx=Info%dt/Info%dX(n); dtdx2=0.5*dtdx
    ! Useful index range vectors
    ALLOCATE(Ip1(nq),Im1(nq),STAT=iError)
    IF (iError /= 0) THEN
      PRINT *,'Cannot allocate index range arrays in numflux'
      STOP
    END IF
    Im1=SPREAD(-1,1,nq); Ip1=SPREAD(1,1,nq)        
    SELECT CASE (Info%method(8))
      CASE (centered)
        CALL extract1D(n,indx,inlow,inhigh,Info,q1D)
        CALL physflux(n,indx,RequestFluxes,Info,q1D,fq,s,A)
        numfqr = 0.5*(fq + EOSHIFT(fq,Ip1))
        numfql = 0.5*(EOSHIFT(fq,Im1) + fq)
      CASE (LaxFriedrichs)            
        fact=1.0/2.0**Info%nDim/dtdx
        CALL extract1D(n,indx,inlow,inhigh,Info,q1D)
        CALL physflux(n,indx,RequestFluxes,Info,q1D,fq,s,A)		
        numfqr = 0.5*(fq + EOSHIFT(fq,Ip1)) - fact*(EOSHIFT(q1D,Ip1) - q1D)
        numfql = 0.5*(EOSHIFT(fq,Im1) + fq) - fact*(q1D - EOSHIFT(q1D,Im1))	
      CASE (Upwind)
        CALL extract1D(n,indx,inlow,inhigh,Info,q1D)
        CALL physflux(n,indx,RequestFluxes,Info,q1D,fq,s,A)      	
        CALL physflux(n,indx,RequestSpeeds,Info,q1D,fq,s,A)
	      ! Numerical flux at left edge of cell
	      WHERE (s<zero)
	        numfql = fq
	      ELSEWHERE
          numfql = EOSHIFT(fq,Im1)
        END WHERE
        ! Numerical flux at right edge of cell			       	
        WHERE (s>zero)
          numfqr = fq
        ELSEWHERE
          numfqr = EOSHIFT(fq,Ip1)
        END WHERE
      CASE (LaxWendroff)
        ! More descriptive names
        Acenter=>Info%Aql; Aedge=>Info%Aqr; qcenter=>Info%q1D; 
        qedge=>Info%sr ! We won't need the speeds, so use the memory space for
	                     ! evaluating cell edge values
        work=>Info%sl			
        qedge=0.; Aedge=0.
        qcenter=0.; Acenter=0.
        work=0.
        ! Get the physical fluxes at the cell centers
        CALL extract1D(n,indx,inlow,inhigh,Info,qcenter)
        CALL physflux(n,indx,RequestFluxes,Info,qcenter,fq,s,Acenter)	
        ! Get the flux Jacobians at the cell centers
        CALL physflux(n,indx,RequestJacobians,Info,qcenter,fq,s,Acenter)	
        numfqr=0.5*(fq+EOSHIFT(fq,Ip1))
        numfql=0.5*(EOSHIFT(fq,Im1)+fq)
        ! Evaluate cell edge field values by arithmetic average
        qedge=0.5*(qcenter+EOSHIFT(qcenter,Ip1))
        CALL physflux(n,indx,RequestJacobians,Info,qedge,fq,s,Aedge)
        fact=0.5*dtdx
        work=EOSHIFT(fq,Ip1)-fq
        DO i=2-mbc,mx+mbc-1
          work(i,1:nq)=fact*MATMUL(Aedge(i,1:nq,1:nq),work(i,1:nq))	  
        END DO
        numfqr=numfqr-work
        numfql=numfql-EOSHIFT(work,Im1)
        ! Include cross-derivative terms	
        DO nDim=1,Info%nDim
          indx2=indx
          IF (nDim==n) CYCLE
          fact=0.125*Info%dt/Info%dX(nDim)          
          indx2(nDim)=indx(nDim)+1
          CALL extract1D(n,indx2,inlow,inhigh,Info,qcenter)
          CALL physflux(n,indx,RequestFluxes,Info,qcenter,fq,s,Acenter)
          work=fq
          indx2(nDim)=indx(nDim)-1
          CALL extract1D(n,indx2,inlow,inhigh,Info,qcenter)
          CALL physflux(n,indx,RequestFluxes,Info,qcenter,fq,s,Acenter)
          work=work-fq
          DO i=2-mbc,mx+mbc-1
            work(i,1:nq)=fact*MATMUL(Acenter(i,1:nq,1:nq),work(i,1:nq))	  
          END DO	  
          numfqr=numfqr-work-EOSHIFT(work,Ip1)
          numfql=numfql-EOSHIFT(work,Im1)-work
        END DO
      CASE (Richtmyer)
        SELECT CASE(Info%nstage)
          CASE (1) ! First stage of algorithm - take a half Lax-Friedrichs step
            fact=1.0/2.0**Info%nDim/dtdx
            CALL extract1D(n,indx,inlow,inhigh,Info,q1D)
            CALL physflux(n,indx,RequestFluxes,Info,q1D,fq,s,A)		
            numfqr = 0.5*(fq + EOSHIFT(fq,Ip1)) - fact*(EOSHIFT(q1D,Ip1) - q1D)
            numfql = 0.5*(EOSHIFT(fq,Im1) + fq) - fact*(q1D - EOSHIFT(q1D,Im1))	
          CASE (2) ! Second stage of algorithm - evaluate fluxes at half space,time step
            
        END SELECT  
      CASE (BeamWarming)
      
      CASE (Fromm)
      CASE (minmod)
      CASE (superbee)
      CASE (MC)
      CASE (vanLeer)
      CASE DEFAULT
        PRINT *,'Bad method(8)=',Info%method(8),'. No such scheme implemented.'
        STOP
    END SELECT
    DEALLOCATE(Im1,Ip1)
  END SUBROUTINE numflux

  ! Update field variable routine  
  SUBROUTINE updateq(n,indx,inlow,inhigh,Info)
    ! Interface declarations    
    INTEGER n,indx(4),inlow,inhigh
    TYPE (NodeInfo) :: Info            
    ! Internal declarations  
    INTEGER mbc,mcapa,i1,i2,i3,i4,i1Dlow,i1Dhigh,mX,nq
    REAL (KIND=qPrec), POINTER, DIMENSION (:) :: dtdx1D
    mcapa=Info%method(6); mbc=Info%mbc; nq=Info%NrVars
    dtdx1D=>Info%dtdx1D
    ! Link between standard 1-mbc:mx+mbc index space and AMR index space
    i1Dlow=inlow+mbc; i1Dhigh=inhigh-mbc    
    i1=indx(1); i2=indx(2); i3=indx(3); i4=indx(4); mX=Info%mXnow
    SELECT CASE (n)
      CASE (1)        
        Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) = Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) - &
           SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq))
      CASE (2)        
        Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) = Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) - &
	         SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq))
      CASE (3)
        Info%q(i1,i2,i1Dlow:i1Dhigh,i4,1:nq) = Info%q(i1,i2,i1Dlow:i1Dhigh,i4,1:nq) - &
	         SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq))
      CASE (4)
        Info%q(i1,i2,i3,i1Dlow:i1Dhigh,1:nq) = Info%q(i1,i2,i3,i1Dlow:i1Dhigh,1:nq) - &
	         SPREAD(dtdx1D(1:mX),1,nq)*(Info%numfqr(1:mX,1:nq)-Info%numfql(1:mX,1:nq))
    END SELECT      
  END SUBROUTINE updateq
  
  ! Extract a 1D slice of data
  SUBROUTINE extract1D(n,indx,inlow,inhigh,Info,q1D)
    ! Interface declarations    
    INTEGER n,indx(4),inlow,inhigh
    TYPE (NodeInfo) :: Info        
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q1D
    ! Internal declarations
    INTEGER i1Dlow,i1Dhigh,i1,i2,i3,i4    
    ! Transform 1D slice of index space to standard 1-mbc:mx+mbc
    i1Dlow=1-Info%mbc; i1Dhigh=inhigh-inlow+i1Dlow    
    ! Copy a 1D slice of data along dimension n from qold to q1D    
    i1=indx(1); i2=indx(2); i3=indx(3); i4=indx(4)
    SELECT CASE (n)
      CASE (1)
        q1D(i1Dlow:i1Dhigh,:)=Info%qold(inlow:inhigh,i2,i3,i4,:)        
      CASE (2)
        q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,inlow:inhigh,i3,i4,:)
      CASE (3)
        q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,i2,inlow:inhigh,i4,:)
      CASE (4)
        q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,i2,i3,inlow:inhigh,:)
    END SELECT
    ! Compute step ratio
    Info%dtdx1D(i1Dlow:i1Dhigh)=Info%dt/Info%dX(n)
  END SUBROUTINE extract1D
  
  ! Allocate fields associated with a particular numerical scheme
  SUBROUTINE AllocSchemeFields(Info)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info  
  END SUBROUTINE AllocSchemeFields
  
  ! Release fields associated with a particular numerical scheme
  SUBROUTINE DeAllocSchemeFields(Info)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info    
  END SUBROUTINE DeAllocSchemeFields
  
  ! Read scheme-specific root level data
  SUBROUTINE ReadSchemeRootData(RootInfo)
    ! Interface declarations
    TYPE (NodeInfo) :: RootInfo
    ! Internal declarations
    INTEGER i,iError
    ! Carry out user-specified initialization of root level grid from input file
    ! None required.      
    ! Input parameters for clawpack routines      
    DO i=1,8
      READ(55,*) RootInfo%method(i)
    END DO
    RootInfo%maux=RootInfo%method(7)      
  END SUBROUTINE ReadSchemeRootData            
  ! End of module CONTAINS clause

!****************
END MODULE scheme
!****************