!========================================================================
! 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:             clawbear.f90
! Purpose:          Coding of wave propagation algorithms as time stepping
!                   routine for BEARCLAW
! Contains:
! Revision History: Ver. 1.0 Oct. 2000 Sorin Mitran
!                   Ver. 2.0 Nov. 2001 SM
! -----------------------------------------------------------------------

!************
MODULE scheme
!************
  ! User definitions of data structures associated with each node
  USE NodeInfoDef
  ! Problem specific routines (setprob,qinit,b4step,src,rpn,rpt)
  USE problem
  ! Conservative fix-up routines. Provided as a service to all solver modules
  USE fixup
  IMPLICIT NONE
  SAVE
  PRIVATE		! Everything is implicitly PRIVATE
  PUBLIC   advance,AllocSchemeFields,DeAllocSchemeFields,ReadSchemeRootData
  
  ! Global variables  
  INTEGER, PARAMETER :: StrangSplitting=2   ! A notation to improve readability
  INTEGER, PARAMETER :: ZeroOrder=0, FirstOrder=1, SecondOrder=2, ThirdOrder=3
  INTEGER, PARAMETER :: DimSplitFirstOrder=-1, DimSplitSecondOrder=-2
  INTEGER, PARAMETER :: MinModLimiter=1, SuperBeeLimiter=2, vanLeerLimiter=3, &
                        MonotonizedCenterLimiter=4, BeamWarmingLimiter=5
                        
  REAL (KIND=xPrec), PARAMETER :: zero=0.0, half=0.5, one=1.0, two=2.0
  
  ! Error codes 
  PUBLIC  err_BadCFL  
  INTEGER, PARAMETER :: ROOT_LEVEL = 0, err_OK=0
  INTEGER, PARAMETER :: err_BadCFL = 1001  ! Error codes > 1000 are considered warnings  
  
  !-------         
  CONTAINS  
  !-------
  
  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. Evaluate source terms in first half of Strang split
    !  ------------------------------------------------------
    IF (Info%method(5) == StrangSplitting) THEN  ! Strang splitting for source term?
      ! Yes. Compute source terms over a half time step
      Info%dt=dt/2.0
      CALL src(Info)
      Info%dt=dt
    END IF       
    !  ----------------------------------------
    !  3. Take one step on the conservation law
    !  ----------------------------------------    
    CALL step(Info)    
    ! ---------------------------------------------------
    !  4. Suggest time step to be first tried on next call
    ! ---------------------------------------------------    
    !IF (Info%cfl>0.d0) Info%dt=dt*Info%cflv(2)/Info%cfl    
    !  ---------------------------------------------------
    !  5. 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    
    !  --------------------------------------------
    !  6. Time step is accepted. Apply source terms
    !  --------------------------------------------    
    SELECT CASE (Info%method(5))
      CASE (1)
        ! Source terms over a full time step        
        CALL src(Info)
      CASE (StrangSplitting)
        ! Source terms over a second half time step for Strang splitting:        
        Info%dt=half*dt
        CALL src(Info)
        Info%dt=dt
      CASE DEFAULT
    END SELECT
    !  -----------------------------------------------
    !  7. Routine called after finishing the time step
    !  ----------------------------------------------- 
    CALL afterstep(Info)
    advance=err_OK
  END FUNCTION advance
  
  ! Carry out a time step
  SUBROUTINE step(Info)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info    
    ! Internal declarations    
    INTEGER n,n1,indx(MaxDims),i1,i2,i3,i4,iError,j
    INTEGER intlow(MaxDims),inthigh(MaxDims),extlow,exthigh
    INTEGER nlow(MaxDims),nhigh(MaxDims)
    INTEGER mbc,meqn,nDim,nstep,nsteps    
    !    
    Info%cfl=zero; nDim=Info%nDim
    nstep=Info%AMRStep; nsteps=Info%AMRSteps; mbc=Info%mbc
    ! Figure out what the domain is in this AMR subiteration
    ! ilow:ihigh is the index span of the current interior cells        
    intlow(1:nDim)=1-(nsteps-nstep)*mbc
    inthigh(1:nDim)=Info%mX(1:nDim)+(nsteps-nstep)*mbc
    ! Perform sweeps along coordinate directions
    DO n=1,nDim    ! This is the dimension along which we'll be solving
                   ! normal Riemann problems      
      ! 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
          IF (Info%method(3) >= 0) THEN
            ! Unsplit version            
            nlow(n1)=intlow(n1)-1   ! loop limits for active transverse dimensions
            nhigh(n1)=inthigh(n1)+1 ! correspond to interior plus one row of ghost cells
          ELSE
            ! Dimensional splitting (fractional steps)      
            nlow(n1)=intlow(n1)-mbc   ! loop limits for active transverse dimensions
            nhigh(n1)=inthigh(n1)+mbc ! correspond to interior plus mbc rows of ghost cells
          END IF
        END IF
      END DO
      ! Current number of cells      
      Info%mXnow=inthigh(n)-intlow(n)+1
      ! extlow:exthigh is the index span of ghost & interior cells
      ! the 'extended' domain
      extlow=intlow(n)-mbc; exthigh=inthigh(n)+mbc
      ! Set up an index array with the interfaces where Riemann problems
      ! will be solved. We're redeclaring one of the components of Info that
      ! is common to all BEARCLAW applications, so we must first release
      ! previous memory.
      IF (ASSOCIATED(Info%I1D)) THEN
        DEALLOCATE(Info%I1D); NULLIFY(Info%I1D)
      END IF
      ALLOCATE(Info%I1D(2-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=2-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
              Info%extlow=extlow; Info%exthigh=exthigh
              ! Compute the fluxes              
              CALL numflux(n,indx,extlow,exthigh,Info)                  
              ! Compute qnew              
              CALL updateq(n,indx,extlow,exthigh,Info)
            END DO
          END DO
        END DO
      END DO      
    END DO    
  END SUBROUTINE step
  
  ! Numerical flux routine:
  !  Compute the numerical flux along dimension n
  !  1D Slice is at indices indx
  SUBROUTINE numflux(n,indx,extlow,exthigh,Info)
    ! Interface declarations
    INTEGER n,indx(MaxDims),extlow,exthigh
    TYPE (NodeInfo) :: Info
    ! Internal declarations
    INTEGER j,mw,mbc,mx,nq,mwaves,NrVars,iError,nDim,indx2(MaxDims)
    INTEGER ydir,zdir
    INTEGER, ALLOCATABLE, DIMENSION (:) :: I    
    REAL (KIND=qPrec) :: fact,cfl
    REAL (KIND=qPrec), POINTER, DIMENSION (:) :: dtdx1D,dtdx1DEdge
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q1D,numfql,numfqr,fq,s
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: speed,Ftilde
    REAL (KIND=qPrec), ALLOCATABLE, DIMENSION (:,:) :: work_s
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: wave,A    
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:,:,:) :: tflux
    !---------------------------------------------------------------------
    ! Associate local names with fields of grid data structure
    ! 1. Current 1D slice & time step
    q1D=>Info%q1D; dtdx1D=>Info%dtdx1D; dtdx1DEdge=>Info%dtdx1DEdge
    ! 2. Cell center Jacobian, flux, speeds
    !    Not used here, but needed for compatibility of CALLs to physflux
    A=>Info%Aql; fq=>Info%fql; s=>Info%sl
    ! 3. Wave propagation local names
    speed=>Info%speed; wave=>Info%wave; Ftilde=>Info%Ftilde
    numfql=>Info%numfql; numfqr=>Info%numfqr; tflux=>Info%tflux
    !---------------------------------------------------------------------
    ! Clean out memory contents
    q1D=zero; dtdx1D=zero; dtdx1DEdge = zero;
    A=zero; fq=zero; s=zero; speed=zero; wave=zero
    Ftilde=zero; numfql=zero; numfqr=zero; tflux=zero
    ! Current number of interior cellsm ghost cells, field variables
    mX=Info%mXnow; mbc=Info%mbc; NrVars=Info%NrVars; mwaves=Info%mwaves
    ! Allocate and initialize a vector with the interior indices 1:mx
    ! + 1 cell to the right
    ALLOCATE(I(1:mx+1),work_s(1:mx+1,1:mwaves),STAT=iError)
    IF (iError /= 0) THEN
      PRINT *,'Error: cannot allocate work vectors in clawbear'
      STOP
    END IF
    DO j=1,mx+1
     I(j)=j
    END DO
    work_s=zero    
    !---------------------------------------------------------------------
    ! Wave propagation scheme
    ! 0. Initialize user physflux routine
    CALL physflux(n,indx,Initialize,Info,q1D,fq,s,A)
    numfql=0.; numfqr=0.
    ! 1. Extract data along the current 1D slice
    CALL extract1D(n,indx,extlow,exthigh,Info,q1D)
    ! 2. Call user routine to compute cell-centered physical flux
    CALL physflux(n,indx,RequestFluxes,Info,q1D,fq,s,A)    
    numfql(1:mX,:) = fq(1:mX,:)
    numfqr(1:mX,:) = fq(1:mX,:)
    ! 3. Call user routine to solve normal Riemann problem
    CALL physflux(n,indx,RequestNormalWaves,Info,q1D,fq,s,A)
    ! 4. Compute the first order (Godunov) numerical flux    
    numfql(1:mX,:) = numfql(1:mX,:)-Info%Apdq(1:mX,:)
    numfqr(1:mX,:) = numfqr(1:mX,:)+Info%Amdq(2:mX+1,:)
    ! 5. Compute CFL
    DO mw=1,mwaves
      cfl=MAXVAL( dtdx1D(1:mX) * ABS(speed(1:mX,mw)) )      
      Info%cfl=MAX(Info%cfl,cfl)
    END DO
    ! 6. Add higher order correction terms
    SELECT CASE (Info%method(2))
      CASE (ZeroOrder,FirstOrder)
      CASE (SecondOrder)
        IF (Info%limiter) CALL ApplyLimiter(Info)        
        dtdX1DEdge(I) = half*(dtdx1D(I-1)+dtdx1D(I))
        Ftilde=zero        
        DO mw=1,mwaves
          work_s(I,mw) = dtdx1DEdge(I)*ABS(speed(I,mw))
          work_s(I,mw) = ABS(speed(I,mw))*(one-work_s(I,mw))
          DO nq=1,NrVars
            Ftilde(I,nq) = Ftilde(I,nq) + work_s(I,mw)*wave(I,nq,mw)
          END DO
        END DO        
        numfql(I,1:NrVars) = numfql(I,1:NrVars)+half*Ftilde(I,1:NrVars)
        numfqr(I,1:NrVars) = numfqr(I,1:NrVars)+half*Ftilde(I+1,1:NrVars)        
        IF (Info%method(3)==2) THEN
          Info%Amdq(I,1:NrVars)=Info%Amdq(I,1:NrVars)+Ftilde(I,1:NrVars)
          Info%Apdq(I,1:NrVars)=Info%Apdq(I,1:NrVars)-Ftilde(I,1:NrVars)
        END IF        
      !CASE (ThirdOrder)
      !  IF (Info%limiter) CALL ApplyLimiter(Info)   
      CASE DEFAULT
        PRINT *,'Normal order of precision, method(2)=', &
                 Info%method(2),'>3 not implemented'
        STOP
    END SELECT
    ! 7. Add transverse corrections
    IF ((Info%method(3)==1 .OR. Info%method(3)==2) .AND. Info%nDim>1) THEN
      SELECT CASE (Info%nDim)
        CASE (2)
          SELECT CASE (n)
            CASE (1)
              ydir=2
            CASE (2)
              ydir=1            
          END SELECT
          Info%tdir=ydir; Info%Asdq=Info%Amdq; Info%ndir=-1
          CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
          DO nq=1,NrVars
            tflux(I-1,-1,1,1,nq) = tflux(I-1,-1,1,1,nq) - &
                                   0.5*dtdx1D(I)*Info%BmAsdq(I,nq)
            tflux(I-1, 1,1,1,nq) = tflux(I-1, 1,1,1,nq) - &
                                   0.5*dtdx1D(I)*Info%BpAsdq(I,nq)
          END DO
          Info%Asdq=Info%Apdq; Info%ndir=1
          CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
          DO nq=1,NrVars
            tflux(I,-1,1,1,nq) = tflux(I,-1,1,1,nq) - &
                                0.5*dtdx1D(I)*Info%BmAsdq(I,nq)
            tflux(I, 1,1,1,nq) = tflux(I, 1,1,1,nq) - &
                                0.5*dtdx1D(I)*Info%BpAsdq(I,nq)
          END DO                               
        CASE (3)
          ! Figure out current orientation
          SELECT CASE (n)
            CASE (1)
              ydir=2; zdir=3
            CASE (2)
              ydir=3; zdir=1
            CASE (3)
              ydir=1; zdir=2
          END SELECT
          ! As=Am branch
          Info%Asdq=Info%Amdq
          ! Split along y
          Info%tdir=ydir
          CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! BmAmdq in BmAsdq, BpAmdq in BpAsdq
            Info%BsAsdq=Info%BpAsdq; Info%Asdq=Info%BmAsdq
            Info%tdir=zdir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! CmBmAmdq in BmAsdq, CpBmAmdq in BpAsdq
            
            Info%Asdq=Info%BsAsdq  ! BpAmdq
            Info%tdir=zdir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! CmBpAmdq in BmAsdq, CpBpAmdq in BpAsdq
          
          Info%tdir=zdir
          CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! CmAmdq in BmAsdq, CpAmdq in BpAsdq
            Info%BsAsdq=Info%BpAsdq; Info%Asdq=Info%BmAsdq
            Info%tdir=ydir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! BmCmAmdq in BmAsdq, BpCmAmdq in BpAsdq
            
            Info%Asdq=Info%BsAsdq  ! BpAmdq
            Info%tdir=ydir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! BmCpAmdq in BmAsdq, BpCpAmdq in BpAsdq
            
          ! As=Ap branch
          Info%Asdq=Info%Apdq
          ! Split along y
          Info%tdir=ydir
          CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! BmApdq in BmAsdq, BpApdq in BpAsdq
            Info%BsAsdq=Info%BpAsdq; Info%Asdq=Info%BmAsdq
            Info%tdir=zdir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! CmBmApdq in BmAsdq, CpBmApdq in BpAsdq
            
            Info%Asdq=Info%BsAsdq  ! BpApdq
            Info%tdir=zdir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! CmBpApdq in BmAsdq, CpBpApdq in BpAsdq
          
          Info%tdir=zdir
          CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! CmApdq in BmAsdq, CpApdq in BpAsdq
            Info%BsAsdq=Info%BpAsdq; Info%Asdq=Info%BmAsdq
            Info%tdir=ydir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! BmCmApdq in BmAsdq, BpCmApdq in BpAsdq
            
            Info%Asdq=Info%BsAsdq  ! BpApdq
            Info%tdir=ydir
            CALL physflux(n,indx,RequestTransverseWaves,Info,q1D,fq,s,A)
            ! BmCpApdq in BmAsdq, BpCpApdq in BpAsdq
          
          
      END SELECT
    END IF
    ! 8. Finalize user physflux routine
    CALL physflux(n,indx,Finalize,Info,q1D,fq,s,A)
    ! 9. Clean up    
    DEALLOCATE(I,work_s,STAT=iError)
    IF (iError /= 0) THEN
      PRINT *,'Error in deallocate of work vectors in clawbear/numflux'
      STOP
    END IF
  END SUBROUTINE numflux
  
  ! Update field variable routine  
  SUBROUTINE updateq(n,indx,extlow,exthigh,Info)
    ! Interface declarations    
    INTEGER n,indx(4),extlow,exthigh
    TYPE (NodeInfo) :: Info            
    ! Internal declarations  
    INTEGER mbc,i1Dlow,i1Dhigh,mX,nq,i1,i2,i3,i4,nmbc
    REAL (KIND=qPrec), POINTER, DIMENSION (:) :: dtdx1D
    REAL (KIND=qPrec) :: dtdx,dtdy,dtdz
    ! Update flux accumulators at boundaries of this grid
    CALL UpdateFixups(n,indx,extlow,exthigh,Info)
    ! Set fluxes at interfaces of this grid with its children to zero
    CALL CoarseFineFlux(n,indx,extlow,exthigh,Info)
    ! Local names
    mbc=Info%mbc; nq=Info%NrVars
    dtdx1D=>Info%dtdx1D
    ! Link between standard 1-mbc:mx+mbc index space and AMR index space
    i1Dlow=extlow+mbc; i1Dhigh=exthigh-mbc    
    i1=indx(1); i2=indx(2); i3=indx(3); i4=indx(4); mX=Info%mXnow
    ! Update values on current 1D slice, current grid
    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
    IF (Info%method(3)>0 .AND. Info%nDim==2) THEN
      ! Add transverse corrections
      nmbc=(Info%nsteps-Info%nnow+1)*mbc
      SELECT CASE (n)
        CASE (1)
          dtdy=Info%dt/Info%dX(2)          
          Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) = Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) &
             - dtdy*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))
          IF (i2+1 <= Info%mX(2)+nmbc) THEN
            Info%q(i1Dlow:i1Dhigh,i2+1,i3,i4,1:nq) = Info%q(i1Dlow:i1Dhigh,i2+1,i3,i4,1:nq) &
               + dtdy*Info%tflux(1:mX,1,1,1,1:nq)
          END IF
          IF (i2-1 >= 1-nmbc) THEN
            Info%q(i1Dlow:i1Dhigh,i2-1,i3,i4,1:nq) = Info%q(i1Dlow:i1Dhigh,i2-1,i3,i4,1:nq) &
               - dtdy*Info%tflux(1:mX,-1,1,1,1:nq)
          END IF
        CASE (2)        
          dtdx=Info%dt/Info%dX(1)
          Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) = Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) &
		         - dtdx*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))
		      IF (i1+1 <= Info%mX(1)+nmbc) THEN
  		      Info%q(i1+1,i1Dlow:i1Dhigh,i3,i4,1:nq) = Info%q(i1+1,i1Dlow:i1Dhigh,i3,i4,1:nq) &
	  	         + dtdx*Info%tflux(1:mX,1,1,1,1:nq)
	  	    END IF
	  	    IF (i1-1 >= 1-nmbc) THEN
  		      Info%q(i1-1,i1Dlow:i1Dhigh,i3,i4,1:nq) = Info%q(i1-1,i1Dlow:i1Dhigh,i3,i4,1:nq) &
	  	         - dtdx*Info%tflux(1:mX,-1,1,1,1:nq)
	  	    END IF
      END SELECT  
    END IF
    IF (Info%method(3)>0 .AND. Info%nDim==3) THEN
      ! Add transverse corrections
      SELECT CASE (n)
        CASE (1)
          Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) = Info%q(i1Dlow:i1Dhigh,i2,i3,i4,1:nq) &
             - dtdy*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))      &
             - dtdz*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))
        CASE (2)        
          Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) = Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) &
             - dtdy*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))      &
             - dtdz*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))
        CASE (3)
          Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) = Info%q(i1,i1Dlow:i1Dhigh,i3,i4,1:nq) &
		         - dtdy*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))      &
             - dtdz*(Info%tflux(1:mX,1,1,1,1:nq) - Info%tflux(1:mX,-1,1,1,1:nq))        
      END SELECT  
    END IF
  END SUBROUTINE updateq
  
  ! Extract a 1D slice of data
  SUBROUTINE extract1D(n,indx,extlow,exthigh,Info,q1D)
    ! Interface declarations    
    INTEGER n,indx(MaxDims),extlow,exthigh
    TYPE (NodeInfo) :: Info        
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: q1D
    ! Internal declarations
    INTEGER i1Dlow,i1Dhigh,i1,i2,i3,i4
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: aux1D
    q1D=0.; Info%ChildSign1D=0; aux1D=>Info%aux1D
    ! Transform 1D slice of index space to standard 1-mbc:mx+mbc
    i1Dlow=1-Info%mbc; i1Dhigh=Info%mXnow+Info%mbc    
    ! 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(extlow:exthigh,i2,i3,i4,:)        
      CASE (2)
        q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,extlow:exthigh,i3,i4,:)        
      CASE (3)
        q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,i2,extlow:exthigh,i4,:)        
      CASE (4)
        q1D(i1Dlow:i1Dhigh,:)=Info%qold(i1,i2,i3,extlow:exthigh,:)
    END SELECT
    IF (Info%maux>0) THEN
      SELECT CASE (n)
        CASE (1)
          aux1D(i1Dlow:i1Dhigh,:)=Info%aux(extlow:exthigh,i2,i3,i4,:)
        CASE (2)
          aux1D(i1Dlow:i1Dhigh,:)=Info%aux(i1,extlow:exthigh,i3,i4,:)
        CASE (3)
          aux1D(i1Dlow:i1Dhigh,:)=Info%aux(i1,i2,extlow:exthigh,i4,:)
        CASE (4)
          aux1D(i1Dlow:i1Dhigh,:)=Info%aux(i1,i2,i3,extlow:exthigh,:)
      END SELECT
    END IF
    ! Compute step ratio
    Info%dtdx1D(i1Dlow:i1Dhigh)=Info%dt/Info%dX(n)
    ! Call fix-up routine to identify coarse/fine boundaries
    CALL CoarseFine1D(n,indx,extlow,exthigh,Info)  
  END SUBROUTINE extract1D             

  SUBROUTINE ApplyLimiter(Info)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info    
    ! Internal declarations
    INTEGER mw,mwaves,i,mx
    REAL (KIND=qPrec) :: wnorm2,dotl,dotr,wlimit
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: speed
    REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: wave    
    !---------------------------------------------------------------------
    ! Associate local names with fields of grid data structure    
    speed=>Info%speed; wave=>Info%wave
    mwaves=Info%mwaves; mX=Info%mXnow
    !---------------------------------------------------------------------       
    DO mw=1,mwaves
      dotr=zero
      DO i=0,mX+1    
        dotl=dotr
        wnorm2=DOT_PRODUCT(wave(i,:,mw),wave(i,:,mw))
        dotr=DOT_PRODUCT(wave(i+1,:,mw),wave(i,:,mw))        
        IF (i/=0 .AND. wnorm2/=zero) THEN        
          IF (speed(i,mw) > zero) THEN
            wlimit=philim(wnorm2,dotl,Info%mthlim(mw))
          ELSE
            wlimit=philim(wnorm2,dotr,Info%mthlim(mw))
          END IF
          wave(i,:,mw)=wlimit * wave(i,:,mw)
        END IF
      END DO      
    END DO
  END SUBROUTINE ApplyLimiter
  
  REAL (KIND=qPrec) FUNCTION philim(a,b,limiter)
    REAL (KIND=qPrec) a,b
    INTEGER limiter
    REAL (KIND=qPrec) r,c
    r=b/a
    SELECT CASE(limiter)
      CASE (MinModLimiter)
        philim = MAX(zero, MIN(one, r))
      CASE (SuperBeeLimiter)
        philim = MAX(zero,MIN(one, two*r), MIN(two, r))
      CASE (vanLeerLimiter)
        philim = (r + ABS(r)) / (one + ABS(r))
      CASE (MonotonizedCenterLimiter)
        c = half*(one + r)
        philim = MAX(zero,MIN(c,two,two*r))
      CASE (BeamWarmingLimiter)
        philim = r
      CASE DEFAULT
        PRINT *,'Limiter nr. ',limiter,' not implemented'
        STOP
    END SELECT  
  END FUNCTION philim

  ! Allocate fields associated with a particular numerical scheme
  SUBROUTINE AllocSchemeFields(Info,Parent)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info,Parent
    ! Internal declarations    
    INTEGER n,maxm,iError,bcbuf(MaxDims),auxbuf(MaxDims),mbc,rmbc
    INTEGER tDims(3,2)
    maxm = MAXVAL(Info%mX); mbc=Info%mbc; rmbc=mbc*Info%r      
    DO n=1,MaxDims
      IF (Info%mX(n) < 1) Info%mX(n)=1
      IF (Info%mX(n)==1) THEN
        bcbuf(n)=0; auxbuf(n)=0
      ELSE
        bcbuf(n)=rmbc; auxbuf(n)=mbc       
      END IF
    END DO    
    ! Copy over components from parent
    IF (Info%level>ROOT_LEVEL) THEN
      Info%limiter=Parent%limiter; Info%meqn=Parent%meqn; Info%mwaves=Parent%mwaves
      Info%method=Parent%method; Info%mthlim=Parent%mthlim
    END IF
    ! Allocate space for intermediate time field variable arrays, if needed
    IF (Info%method(5)<2) THEN
      ! Only need qold
      Info%qhalf => Info%qold 
    ELSE
      ALLOCATE(Info%qhalf(1-bcbuf(1):Info%mX(1)+bcbuf(1),  &
                          1-bcbuf(2):Info%mX(2)+bcbuf(2),  &
                          1-bcbuf(3):Info%mX(3)+bcbuf(3),  &
                          1-bcbuf(4):Info%mX(4)+bcbuf(4),  &      
                          Info%meqn),                      &
               STAT=iError)
      Info%qhalf=zero
      IF (iError /= 0) THEN
        PRINT *,'Error in allocation of qhalf array in AllocSchemeFields'
        STOP
      END IF
    END IF
    ! Compute the number of transverse dimensions required
    tDims=1
    IF (Info%method(3)>0 .AND. Info%nDim>1) THEN
      tDims(1:Info%nDim-1,1)=-1; tDims(1:Info%nDim-1,2)=1;
    END IF
    ! Work space allocation            
    ALLOCATE(Info%dtdx1DEdge(1-mbc:maxm+2*rmbc),                     &
             Info%Apdq(1-mbc:maxm+2*rmbc,Info%NrVars),               &
             Info%Amdq(1-mbc:maxm+2*rmbc,Info%NrVars),               &
             Info%Ftilde(1-mbc:maxm+2*rmbc,Info%NrVars),             &
             Info%speed(1-mbc:maxm+2*rmbc,Info%mwaves),              &             
             Info%Asdq(1-mbc:maxm+2*rmbc,Info%NrVars),               &
             Info%BpAsdq(1-mbc:maxm+2*rmbc,Info%NrVars),             &
             Info%BmAsdq(1-mbc:maxm+2*rmbc,Info%NrVars),             &
             Info%BsAsdq(1-mbc:maxm+2*rmbc,Info%NrVars),             &
             Info%wave(1-mbc:maxm+2*rmbc,Info%NrVars,Info%mwaves),   &
             Info%tflux(1-mbc:maxm+2*rmbc,tDims(1,1):tDims(1,2),     &
                        tDims(2,1):tDims(2,2),tDims(3,1):tDims(3,2), &
                        Info%NrVars),                                &
             STAT=iError)             
    IF (iError /= 0) THEN
      PRINT *,'Error in allocation of work arrays in AllocSchemeFields'
      STOP
    ELSE      
      Info%dtdX1DEdge=zero; Info%Apdq=zero; Info%Amdq=zero
      Info%Ftilde=zero; Info%speed=zero; 
      Info%Asdq=zero; Info%BpAsdq=zero; Info%BmAsdq=zero
      Info%wave=zero;
    END IF  
  END SUBROUTINE AllocSchemeFields
  
  ! Release fields associated with a particular numerical scheme
  SUBROUTINE DeAllocSchemeFields(Info)
    ! Interface declarations    
    TYPE (NodeInfo) :: Info
    ! Internal declarations
    INTEGER iError
    IF (Info%method(5)>=2) THEN
      IF (ASSOCIATED(Info%qhalf)) DEALLOCATE(Info%qhalf)
      NULLIFY(Info%qhalf)
    END IF
    IF (ASSOCIATED(Info%dtdX1DEdge)) DEALLOCATE(Info%dtdX1DEdge)
    IF (ASSOCIATED(Info%Apdq)) DEALLOCATE(Info%Apdq)
    IF (ASSOCIATED(Info%Amdq)) DEALLOCATE(Info%Amdq)
    IF (ASSOCIATED(Info%Ftilde)) DEALLOCATE(Info%Ftilde)
    IF (ASSOCIATED(Info%speed)) DEALLOCATE(Info%speed)
    IF (ASSOCIATED(Info%Asdq)) DEALLOCATE(Info%Asdq)
    IF (ASSOCIATED(Info%BpAsdq)) DEALLOCATE(Info%BpAsdq)
    IF (ASSOCIATED(Info%BmAsdq)) DEALLOCATE(Info%BmAsdq)
    IF (ASSOCIATED(Info%BmAsdq)) DEALLOCATE(Info%BsAsdq)
    IF (ASSOCIATED(Info%wave)) DEALLOCATE(Info%wave)
    IF (ASSOCIATED(Info%tflux)) DEALLOCATE(Info%tflux)
    NULLIFY(Info%dtdX1DEdge,Info%Apdq,Info%Amdq,Info%Ftilde,Info%speed, &
            Info%Asdq,Info%BpAsdq,Info%BmAsdq,Info%BsAsdq,Info%wave,Info%tflux)
  END SUBROUTINE DeAllocSchemeFields
  
  SUBROUTINE ReadSchemeRootData(RootInfo)
    ! Interface declarations
    TYPE (NodeInfo) :: RootInfo
    ! Internal declarations
    INTEGER i,iError
    ! Carry out CLAWBEAR specific initialization of root level grid from input file
    RootInfo%meqn=RootInfo%NrVars    
    ! Input parameters for clawpack routines
    DO i=1,7
      READ(55,*) RootInfo%method(i)
    END DO    
    READ(55,*) RootInfo%mwaves
    IF (RootInfo%mwaves>RootInfo%MaxWaves) THEN
      PRINT *,'Too many waves. Increase parameter in NodeInfoTypeWaveBEAR.f90'
      STOP
    END IF    
    RootInfo%mthlim=0    
    READ(55,*)(RootInfo%mthlim(i),i=1,RootInfo%mwaves)    
    IF (RootInfo%method(6)<0) THEN      
      PRINT *,'BEARCLAW error: method(6)=mcapa must be positive or zero'
      STOP
    ENDIF
    IF (RootInfo%method(6)>RootInfo%method(7)) THEN      
      PRINT *,'BEARCLAW error: method(6) > method(7)'
      STOP
    ENDIF      
    ! Flag limiters    
    RootInfo%limiter = .FALSE.
    DO i=1,RootInfo%mwaves
     IF (RootInfo%mthlim(i)>0) RootInfo%limiter = .TRUE.
    END DO    
  END SUBROUTINE ReadSchemeRootData
  
  ! End of module CONTAINS clause

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