!========================================================================
! 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, moving mesh, flexible-walled shock tube 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,cstring1,cstring2,mstring1,mstring2, &
                       kappa3,kappa4,kappa5,kappa6,ystring1,ystring2
  !  - Initial values
  REAL (KIND=qPrec) :: rho01,u01,v01,p01,rho02,u02,v02,p02,entropy0
  REAL (KIND=qPrec), DIMENSION(4) :: q01,q02
  !  - Current string shapes, velocities, local force
  INTEGER :: nstring=1000
  REAL (KIND=qPrec) :: dxstring=2.0/1000.0
  REAL (KIND=qPrec), DIMENSION(0:1000), TARGET :: y3,y4,y5,y6,u3,u4,u5,u6, &
                     f3,f4,f5,f6,theta3,theta4,theta5,theta6
  !
  INTEGER, PARAMETER :: OK=0
  REAL (KIND=qPrec), PARAMETER :: zero=0., half=0.5, one=1., two=2.
  
  !=======
  CONTAINS
  !=======
    
    SUBROUTINE get_string(istring,n,x,y,theta,u,f)
      ! For n given abscissas in x return the deflection y, tangent theta,
      ! velocity u and external force of string istring
      INTEGER istring,n
      REAL (KIND=qPrec), DIMENSION(:) :: x,y,theta,u,f
      !
      INTEGER i,ix
      REAL (KIND=qPrec) :: delx
      REAL (KIND=qPrec), DIMENSION(:), POINTER :: ystr,ustr,thetastr,fstr
      ! Select appropriate global variables to obtain values from
      SELECT CASE(istring)
        CASE (3)
          ystr=>y3; ustr=>u3; thetastr=>theta3; fstr=>f3; delx=0.
        CASE (4)
          ystr=>y4; ustr=>u4; thetastr=>theta4; fstr=>f4; delx=2.
        CASE (5)
          ystr=>y5; ustr=>u5; thetastr=>theta5; fstr=>f5; delx=0.
        CASE (6)
          ystr=>y6; ustr=>u6; thetastr=>theta6; fstr=>f6; delx=2.
        CASE DEFAULT
          PRINT *,'Bad code in string, istring=',istring
          STOP
      END SELECT
      DO i=1,n
        ix=(x(i)-delx)/dxstring; ix=MOD(ix,nstring+1)
        y(i)=ystr(ix); theta(i)=thetastr(ix); u(i)=ustr(ix); f(i)=fstr(ix)
      END DO
    END SUBROUTINE get_string

    SUBROUTINE set_string(istring,iset,n,x,y,theta,u,f)
      ! For n given abscissas in x set the deflection y, tangent theta,
      ! velocity u and external force of string istring
      INTEGER istring,iset,n
      REAL (KIND=qPrec), DIMENSION(:) :: x,y,theta,u,f
      !
      INTEGER i,ix,ix1,ix2,jx1,jx2
      REAL (KIND=qPrec) :: delx,dx,dy,du,dtheta,df,xfine
      REAL (KIND=qPrec), DIMENSION(:), POINTER :: ystr,ustr,thetastr,fstr
      ! Select appropriate global variables to set
      SELECT CASE(istring)
        CASE (3)
          ystr=>y3; ustr=>u3; thetastr=>theta3; fstr=>f3; delx=0.
        CASE (4)
          ystr=>y4; ustr=>u4; thetastr=>theta4; fstr=>f4; delx=2.
        CASE (5)
          ystr=>y5; ustr=>u5; thetastr=>theta5; fstr=>f5; delx=0.
        CASE (6)
          ystr=>y6; ustr=>u6; thetastr=>theta6; fstr=>f6; delx=2.
        CASE DEFAULT
          PRINT *,'Bad code in string, istring=',istring
          STOP
      END SELECT
      ! Set string values on fine grid, piecewise linear interpolation
      dx=x(2)-x(1)
      ix1=(x(1)-dx/2.-delx)/dxstring; ix1=MOD(ix1,nstring+1)
      ix2=(x(n)+dx/2.-delx)/dxstring; ix2=MOD(ix2,nstring+1)
      jx1=1; jx2=2;
      SELECT CASE(iset)
        CASE (1:2) ! Set force on string (i.e. gas pressure)
          df=(f(jx2)-f(jx1))/dx
          DO i=ix1,ix2
            xfine=delx+i*dxstring
            IF (xfine>=x(jx2) .AND. jx2<n) THEN
              jx1=jx1+1; jx2=jx2+1
              df=(f(jx2)-f(jx1))/dx
            END IF
            fstr(i)=f(jx1)+df*(xfine-x(jx1))
          END DO
        CASE (3:6) ! Set string geometry, velocity
          dy=(y(jx2)-y(jx1))/dx; dtheta=(theta(jx2)-theta(jx1))/dx
          du=(u(jx2)-u(jx1))/dx
          DO i=ix1,ix2
            xfine=delx+i*dxstring
            IF (xfine>=x(jx2) .AND. jx2<n) THEN
              jx1=jx1+1; jx2=jx2+1
              dy=(y(jx2)-y(jx1))/dx; dtheta=(theta(jx2)-theta(jx1))/dx
              du=(u(jx2)-u(jx1))/dx
            END IF
            ystr(i)=y(jx1)+dy*(xfine-x(jx1))
            thetastr(i)=theta(jx1)+dtheta*(xfine-x(jx1))
            ustr(i)=u(jx1)+du*(xfine-x(jx1))            
          END DO
          IF (iset==3 .OR. iset==4) THEN
            ystr(0)=ystring1; ustr(0)=0.
            ystr(nstring)=ystring1; ustr(nstring)=0.
          ELSE
            ystr(0)=ystring2; ustr(0)=0.
            ystr(nstring)=ystring2; ustr(nstring)=0.
          END IF
      END SELECT
    END SUBROUTINE set_string    
    
    SUBROUTINE setcons(rho,u,v,p,q)
    REAL (KIND=qPrec) :: rho,u,v,p,q(4)
      ! Compute conservative variables from primitive variables
      q(1)=rho; q(2)=rho*u; q(3)=rho*v
      q(4)=p/gamma1+0.5*rho*(u**2+v**2)
    END SUBROUTINE setcons
    
    SUBROUTINE setprob      
      INTEGER i
      REAL (KIND=qPrec) :: x
      ! Read gas constant from file
      OPEN(UNIT=7,FILE='setprob.data',STATUS='old',FORM='formatted')      
      READ(7,*)gamma
      gamma1 = gamma - 1.d0
      READ(7,*)rho01,u01,v01,p01
      CALL setcons(rho01,u01,v01,p01,q01)
      READ(7,*)rho02,u02,v02,p02
      CALL setcons(rho02,u02,v02,p02,q02)      
      entropy0=(LOG(p01)-gamma*LOG(rho01))/gamma1      
      READ(7,*)cstring1,cstring2,mstring1,mstring2
      kappa3=0.5*p01/cstring2**2/mstring2
      kappa4=0.5*p02/cstring2**2/mstring2
      kappa5=0.5*p01/cstring1**2/mstring1
      kappa6=0.5*p02/cstring1**2/mstring1
      READ(7,*)ystring1,ystring2
      CLOSE(7)
      ! Compute initial string shape
      DO i=0,nstring
       x=i*dxstring
       y3(i) = kappa3*x*(x-2.)+ystring1
       theta3(i) = kappa3*2*(x-1.)
       y5(i) = -kappa5*x*(x-2.)+ystring2
       theta5(i) = -kappa5*2*(x-1.)
       x=x+2.
       y4(i) = kappa4*(x-2.)*(x-4.)+ystring1
       theta4(i) = kappa4*2*(x-3.)       
       y6(i) = -kappa6*(x-2.)*(x-4.)+ystring2
       theta6(i) = -kappa6*2*(x-3.)       
      END DO
      ! Set initial string velocity
      u3=0.; u4=0.; u5=0.; u6=0.
    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
      INTEGER ibc,rmbc,r,mbc,mx,my,i,istring,np,iError
      REAL (KIND=qPrec), ALLOCATABLE, DIMENSION(:) :: &
                         xstring,ystring,ustring,thetastring,fstring
      ! Boundary conditions for vibrating string:      
      r=Info%r; mbc=Info%mbc; rmbc=r*mbc; mx=Info%mX(1)
      SELECT CASE(Info%mthbc(Info%bcnow))
        CASE (11) ! left boundary of string
          ! q(1)=theta extended along r*mbc
          ! q(2)=u reflected
          DO ibc=1,rmbc
            Info%q(1-ibc,1,1,1,1) =  Info%q(ibc,1,1,1,1)
            Info%q(1-ibc,1,1,1,2) = -Info%q(ibc,1,1,1,2)
          END DO           
        CASE (12) ! right boundary of string
          ! q(1)=theta extended along r*mbc
          ! q(2)=u reflected
          DO ibc=1,rmbc
            Info%q(mx+ibc,1,1,1,1) =  Info%q(mx+1-ibc,1,1,1,1)
            Info%q(mx+ibc,1,1,1,2) = -Info%q(mx+1-ibc,1,1,1,2)
          END DO
        CASE (13:14) ! fluid-string boundary (bottom of fluid domain)         
          ALLOCATE(xstring(1-rmbc:mx+rmbc),ystring(1-rmbc:mx+rmbc),     &
                   ustring(1-rmbc:mx+rmbc),thetastring(1-rmbc:mx+rmbc), &
                   fstring(1-rmbc:mx+rmbc),STAT=iError)
          IF (iError /= 0) THEN
            PRINT *,'Cannot allocate temporary vectors in problemBC'
            STOP
          END IF                    
          ! Copy all components from domain interior          
          DO ibc=1,rmbc
            Info%q(:,1-ibc,1,1,:) = Info%q(:,ibc,1,1,:)            
          END DO
          ! Compute current abscissas
          DO i=1-rmbc,mx+rmbc
            xstring(i)=Info%Xlower(1)+(i-0.5)*Info%dX(1)
          END DO
          ! Find current string velocity          
          istring=Info%mthbc(Info%bcnow)-10; np=mx+2*rmbc
          CALL get_string(istring,np,xstring,ystring,thetastring,ustring,fstring)
          ! Negate the normal velocity (3rd component):
          DO ibc=1,rmbc
            Info%q(1:mx,1-ibc,1,1,3) = 2*ustring(1:mx)-Info%q(1:mx,ibc,1,1,3)
          END DO
          DEALLOCATE(xstring,ystring,ustring,thetastring,fstring)
        CASE (15:16) ! fluid-string boundary (top of fluid domain)
          my=Info%mX(2)
          ALLOCATE(xstring(1-rmbc:mx+rmbc),ystring(1-rmbc:mx+rmbc),     &
                   ustring(1-rmbc:mx+rmbc),thetastring(1-rmbc:mx+rmbc), &
                   fstring(1-rmbc:mx+rmbc),STAT=iError)
          IF (iError /= 0) THEN
            PRINT *,'Cannot allocate temporary vectors in problemBC'
            STOP
          END IF
          ! Copy all components from domain interior
          DO ibc=1,rmbc
            Info%q(:,my+ibc,1,1,:) = Info%q(:,my+1-ibc,1,1,:)
          END DO
          ! Compute current abscissas
          DO i=1-rmbc,mx+rmbc
            xstring(i)=Info%Xlower(1)+(i-0.5)*Info%dX(1)
          END DO
          ! Find current string velocity          
          istring=Info%mthbc(Info%bcnow)-10; np=mx+2*rmbc
          CALL get_string(istring,np,xstring,ystring,thetastring,ustring,fstring)
          ! Negate the normal velocity (3rd component):
          DO ibc=1,rmbc
            Info%q(:,my+ibc,1,1,3) = 2*ustring(:)-Info%q(:,my+1-ibc,1,1,3)
          END DO
          DEALLOCATE(xstring,ystring,ustring,thetastring,fstring)
      END SELECT
    END SUBROUTINE problemBC

    SUBROUTINE qinit(Info)
      ! Interface declarations      
      TYPE (NodeInfo) :: Info  ! Data associated with this grid  
      ! Internal declarations
      INTEGER i,mx,my
      REAL (KIND=qPrec) :: xcell
      !
      SELECT CASE (Info%nEquationSet)
        CASE (1) ! Gas in left state
          mx=Info%mX(1); my=Info%mX(2)
          Info%q(1:mx,1:my,1,1,1)=q01(1)
          Info%q(1:mx,1:my,1,1,2)=q01(2)
          Info%q(1:mx,1:my,1,1,3)=q01(3)
          Info%q(1:mx,1:my,1,1,4)=q01(4)
        CASE (2) ! Gas in right state
          mx=Info%mX(1); my=Info%mX(2)
          Info%q(1:mx,1:my,1,1,1)=q02(1)
          Info%q(1:mx,1:my,1,1,2)=q02(2)
          Info%q(1:mx,1:my,1,1,3)=q02(3)
          Info%q(1:mx,1:my,1,1,4)=q02(4)
        CASE (3)
          DO i=1,Info%mX(1)
            xcell = Info%Xlower(1) + (i-0.5)*Info%dX(1)        
            Info%q(i,1,1,1,1) = kappa3*2*(xcell-1.)
            Info%q(i,1,1,1,2) = 0.
            Info%q(i,1,1,1,3) = kappa3*xcell*(xcell-2.)+ystring1
          END DO
        CASE (4)
          DO i=1,Info%mX(1)
            xcell = Info%Xlower(1) + (i-0.5)*Info%dX(1)        
            Info%q(i,1,1,1,1) = kappa4*2*(xcell-3.)
            Info%q(i,1,1,1,2) = 0.
            Info%q(i,1,1,1,3) = kappa4*(xcell-2.)*(xcell-4.)+ystring1
          END DO
        CASE (5)
          DO i=1,Info%mX(1)
            xcell = Info%Xlower(1) + (i-0.5)*Info%dX(1)        
            Info%q(i,1,1,1,1) = -kappa5*2*(xcell-1.)
            Info%q(i,1,1,1,2) = 0.
            Info%q(i,1,1,1,3) = -kappa5*xcell*(xcell-2.)+ystring2
          END DO
        CASE (6)
          DO i=1,Info%mX(1)
            xcell = Info%Xlower(1) + (i-0.5)*Info%dX(1)        
            Info%q(i,1,1,1,1) = -kappa6*2*(xcell-3.)
            Info%q(i,1,1,1,2) = 0.
            Info%q(i,1,1,1,3) = -kappa6*(xcell-2.)*(xcell-4.)+ystring2
          END DO
        CASE DEFAULT
          PRINT *,'Bad Info%nEquationSet in qinit'
      END SELECT    
    END SUBROUTINE qinit
    
    SUBROUTINE b4step(Info)  
      ! Interface declarations    
      TYPE (NodeInfo) :: Info
      !
      INTEGER i,j,mbc,r,rmbc,mcsi,meta,iError,i34,i56,np
      REAL (KIND=xPrec) :: t,csi,eta,dcsi,deta,csilower,etalower
      REAL (KIND=qPrec), ALLOCATABLE, DIMENSION(:) :: &
                         xstring,y34,u34,theta34,y56,u56,theta56,fstring
      !
      IF (Info%nEquationSet>2) RETURN ! Nothing to be done on strings
      ! Set auxilliary variables for gas dynamics on moving mesh      
      ! Get information about current grid                           
      mbc=Info%mbc; r=Info%r; rmbc=r*mbc      
      mcsi=Info%mX(1); meta=Info%mX(2) ; dcsi=Info%dX(1); deta=Info%dX(2)     
      csilower=Info%Xlower(1); etalower=Info%Xlower(2)
      ALLOCATE(xstring(1-rmbc:mcsi+rmbc),y34(1-rmbc:mcsi+rmbc),u34(1-rmbc:mcsi+rmbc), &
               theta34(1-rmbc:mcsi+rmbc),y56(1-rmbc:mcsi+rmbc),u56(1-rmbc:mcsi+rmbc), &
               theta56(1-rmbc:mcsi+rmbc),fstring(1-rmbc:mcsi+rmbc),STAT=iError)
      IF (iError /= 0) THEN
        PRINT *,'Cannot allocate temporary vectors in setaux'
        STOP
      END IF
      ! Get current string positions      
      DO i=1-rmbc,mcsi+rmbc
        xstring(i)=csilower+(i-0.5)*dcsi
      END DO
      SELECT CASE(Info%nEquationSet)
        CASE (1)
          i34=3; i56=5
        CASE (2)
          i34=4; i56=6
      END SELECT
      np=mcsi+2*rmbc
      CALL get_string(i34,np,xstring,y34,theta34,u34,fstring)
      CALL get_string(i56,np,xstring,y56,theta56,u56,fstring)
      ! Compute current values of metric coefficients      
      DO i=1-rmbc,mcsi+rmbc
        csi=csilower+(i-0.5)*dcsi
        DO j=1-rmbc,meta+rmbc
          eta=etalower+(j-0.5)*deta                    
          Info%aux(i,j,1,1,3)= &     ! aux(i,j,3) contains Y_csi(i,j)
            0.5*(eta+1)*(theta56(i)-theta34(i))+theta34(i)
          Info%aux(i,j,1,1,4)= &     ! aux(i,j,4) contains Y_eta(i,j)
            0.5*(y56(i)-y34(i))          
          Info%aux(i,j,1,1,6)= &     ! aux(i,j,6) contains Y_t(i,j)
            0.5*(eta+1)*(u56(i)-u34(i))+u34(i)
          ! Jacobian = X_csi Y_eta - X_eta Y_csi
          Info%aux(i,j,1,1,7) = Info%aux(i,j,1,1,1)*Info%aux(i,j,1,1,4) - &
                                Info%aux(i,j,1,1,2)*Info%aux(i,j,1,1,3)
        END DO
      END DO
      DEALLOCATE(xstring,y34,u34,theta34,y56,u56,theta56,fstring)
    END SUBROUTINE b4step

    SUBROUTINE get_pres(mx,j,q,p)
      INTEGER mx,j
      REAL (KIND=qPrec), DIMENSION(:,:,:,:,:) :: q
      REAL (KIND=qPrec), DIMENSION(:) :: p      
      p(1:mx)=gamma1*(q(1:mx,j,1,1,4) -                              &
                      0.5*(q(1:mx,j,1,1,2)**2 + q(1:mx,j,1,1,3)**2)/ &
                           q(1:mx,j,1,1,1))
    END SUBROUTINE get_pres                               
    
    SUBROUTINE afterstep(Info)
      ! Interface declarations    
      TYPE (NodeInfo) :: Info
      ! Include operations that need to be carried out after each time step here:
      INTEGER mcsi,meta,iError,i,istring,iset
      REAL (KIND=qPrec) :: dy,dt
      REAL (KIND=qPrec), PARAMETER :: eta3=-1.0, eta4=-1.0, eta5=1.0, eta6=1.0
      REAL (KIND=qPrec), ALLOCATABLE, DIMENSION(:) :: csi,y,theta,u,f,p
      ! Update global variables that describe fluid-string interaction
      mcsi=Info%mX(1)
      ALLOCATE(csi(mcsi),y(mcsi),theta(mcsi),u(mcsi),f(mcsi),p(mcsi),STAT=iError)
      csi=0.; y=0.; theta=0.; u=0.; f=0.
      IF (iError /= 0) THEN
        PRINT *,'Cannot allocate work arrays in afterstep'
        STOP
      END IF
      DO i=1,mcsi
        csi(i)=Info%Xlower(1)+(i-0.5)*Info%dX(1)
      END DO
      SELECT CASE (Info%nEquationSet)
        CASE (1) ! Fluid domain 1 sets force on strings 3,5
          ! Check whether bottom of fluid domain is adjacent to string
          IF (ABS(Info%Xlower(2)-eta3)<0.25*Info%dX(2)) THEN            
            CALL get_pres(mcsi,1,Info%q,p); f=p
            CALL get_pres(mcsi,1,Info%qold,p)
            ! Average of pressure at current and previous time step.
            ! This implements the trapezoid rule in the integration of the
            ! string source term
            f=0.5*(p+f)/mstring2; f=-f ! On string 3 pressure gives downward force
            istring=3; iset=1
            CALL set_string(istring,iset,mcsi,csi,y,theta,u,f)
          END IF
          ! Check whether top of fluid domain is adjacent to string
          IF (ABS(Info%Xupper(2)-eta5)<0.25*Info%dX(2)) THEN
            meta=Info%mX(2)
            CALL get_pres(mcsi,meta,Info%q,p); f=p
            CALL get_pres(mcsi,meta,Info%qold,p)
            ! Average of pressure at current and previous time step.
            ! This implements the trapezoid rule in the integration of the
            ! string source term
            f=0.5*(p+f)/mstring1
            istring=5; iset=1
            CALL set_string(istring,iset,mcsi,csi,y,theta,u,f)
          END IF
        CASE (2) ! Fluid domain 2 sets force on strings 4,6
          ! Check whether bottom of fluid domain is adjacent to string
          IF (ABS(Info%Xlower(2)-eta4)<0.25*Info%dX(2)) THEN
            istring=4; iset=2
            CALL get_pres(mcsi,1,Info%q,p); f=p
            CALL get_pres(mcsi,1,Info%qold,p)
            ! Average of pressure at current and previous time step.
            ! This implements the trapezoid rule in the integration of the
            ! string source term
            f=0.5*(p+f)/mstring2; f=-f ! On string 4 pressure gives downward force
            CALL set_string(istring,iset,mcsi,csi,y,theta,u,f)
          END IF
          ! Check whether top of fluid domain is adjacent to string
          IF (ABS(Info%Xupper(2)-eta6)<0.25*Info%dX(2)) THEN
            istring=6; iset=2; meta=Info%mX(2)
            CALL get_pres(mcsi,meta,Info%q,p); f=p
            CALL get_pres(mcsi,meta,Info%qold,p)
            ! Average of pressure at current and previous time step.
            ! This implements the trapezoid rule in the integration of the
            ! string source term
            f=0.5*(p+f)/mstring1
            CALL set_string(istring,iset,mcsi,csi,y,theta,u,f)
          END IF
        CASE (3:6)
          istring=Info%nEquationSet; iset=istring; dt=Info%dt
          theta=Info%q(1:mcsi,1,1,1,1); u=Info%q(1:mcsi,1,1,1,2)
          ! Update string position (Trapezoid integration rule)
          Info%q(1:mcsi,1,1,1,3) = Info%q(1:mcsi,1,1,1,3) + 0.5*dt* &
                     (Info%q(1:mcsi,1,1,1,2)+Info%qold(1:mcsi,1,1,1,2))
          y=Info%q(1:mcsi,1,1,1,3)
          CALL set_string(istring,iset,mcsi,csi,y,theta,u,f)
      END SELECT
      DEALLOCATE(csi,y,theta,u,f,p)
    END SUBROUTINE afterstep
    
    SUBROUTINE extend(mx,my,f)
      INTEGER mx,my
      REAL (KIND=qPrec), DIMENSION(0:mx,0:my) :: f,g
      ! From the cell-center data in f(1:mx,1:my) construct cell-node
      ! data in f(0:mx,0:my)      
      f(0,:)=f(1,:)
      f(:,0)=f(:,1)
      RETURN
    END SUBROUTINE extend

    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
      !
      SELECT CASE (Info%nEquationSet)
        CASE (1:2)
          CALL problemIO_gas(nframe,tnow,IOrequest,Info,qmax,qmin)
        CASE (3:6)
          CALL problemIO_string(nframe,tnow,IOrequest,Info,qmax,qmin)
      END SELECT
    END SUBROUTINE problemIO
    
    SUBROUTINE problemIO_gas(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,UserNrOutVars=-4      
      ! 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
      ! -4 UserNrOutVars    - Only for UserHDF output style. routine is being
      !                       called to permit modification of number of
      !                       output variables NrOutVars
      ! 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
      INTEGER, PARAMETER :: NrOutVars=6, NrVars=4, xCoord=5, yCoord=6
      INTEGER :: mx,my,i,j,iError,i34,i56,np
      REAL (KIND=qPrec) :: csi,eta,dcsi,deta,csilower,etalower,x,y
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION(:) :: xstr,y34,y56,ustr,tstr,fstr
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION(:,:) :: rho,u,v,pres, &
                                                              entropy,Mach,rhoE
      mx=Info%mX(1); my=Info%mX(2)                     
      SELECT CASE (IOrequest)        
        CASE (UserNrOutVars) ! Set number of output variables, carry out
                             ! additional initializations
          ! We need to output the current grid coordinates in addition
          ! to the field variables
          Info%NrOutVars=NrOutVars
          ALLOCATE(rho(0:mx,0:my),u(0:mx,0:my),v(0:mx,0:my),pres(0:mx,0:my), &
                   entropy(0:mx,0:my),Mach(0:mx,0:my),rhoE(0:mx,0:my),       &
                   xstr(0:mx),y34(0:mx),y56(0:mx),ustr(0:mx),tstr(0:mx),     &
                   fstr(0:mx),STAT=iError)
          IF (iError/=0) THEN
            PRINT *,'Cannot allocate output buffers in problemIO'
            STOP
          END IF          
          rho(1:mx,1:my)=Info%q(1:mx,1:my,1,1,1);   CALL extend(mx,my,rho)
          u(1:mx,1:my)=Info%q(1:mx,1:my,1,1,2);     CALL extend(mx,my,u)
          v(1:mx,1:my)=Info%q(1:mx,1:my,1,1,3);     CALL extend(mx,my,v)
          rhoE(1:mx,1:my)=Info%q(1:mx,1:my,1,1,4);  CALL extend(mx,my,rhoE)
          u=u/rho; v=v/rho
          entropy=u**2+v**2
          pres=gamma1*(rhoE-0.5*entropy)
          Mach=SQRT( entropy/(gamma*pres) * rho)
          entropy=(LOG(pres)-gamma*LOG(rho))/gamma1 - entropy0
        CASE (UserIOMinMax)
          ! Find minimum and maximum field variable values encountered
          IF (nframe==0) THEN
            qmax(1)=MAXVAL(Mach); qmin(1)=MINVAL(Mach)
            qmax(2)=MAXVAL(u); qmin(2)=MINVAL(u)
            qmax(3)=MAXVAL(v); qmin(3)=MINVAL(v)
            qmax(4)=MAXVAL(entropy); qmin(4)=MINVAL(entropy)
          ELSE
            qmax(1)=MAX(qmax(1),MAXVAL(Mach)); qmin(1)=MIN(qmin(1),MINVAL(Mach))
            qmax(2)=MAX(qmax(2),MAXVAL(u)); qmin(2)=MIN(qmin(2),MINVAL(u))
            qmax(3)=MAX(qmax(3),MAXVAL(v)); qmin(3)=MIN(qmin(3),MINVAL(v))
            qmax(4)=MAX(qmax(4),MAXVAL(entropy))
            qmin(4)=MIN(qmin(4),MINVAL(entropy))         
          END IF
          ! Last request to problemIO. Release buffer space
          DEALLOCATE(rho,u,v,pres,Mach,entropy,rhoE,xstr,y34,y56,ustr,tstr,fstr)
        CASE (1)
          ! Fill the output buffer with the requested field variable          
          Info%qbuf(0:mx,0:my,1,1) = Mach
        CASE (2)
          ! Fill the output buffer with the requested field variable          
          Info%qbuf(0:mx,0:my,1,1) = u
        CASE (3)
          ! Fill the output buffer with the requested field variable          
          Info%qbuf(0:mx,0:my,1,1) = v
        CASE (4)                  
          ! Fill the output buffer with the requested field variable          
          Info%qbuf(0:mx,0:my,1,1) = entropy
        CASE (xCoord)
          ! Fill the output buffer with the x coordinate          
          dcsi=Info%dX(1); csilower=Info%Xlower(1)
          DO i=0,mx
            csi=csilower+i*dcsi
            x=csi
            Info%qbuf(i,0:my,1,1)=x                 
          END DO          
        CASE (yCoord)
          ! Fill the output buffer with the y coordinate
          dcsi=Info%dX(1); deta=Info%dX(2)
          csilower=Info%Xlower(1); etalower=Info%Xlower(2)
          DO i=0,mx
            csi=csilower+i*dcsi
            xstr(i)=csi            
          END DO
          SELECT CASE (Info%nEquationSet)
            CASE (1)
              i34=3; i56=5
            CASE (2)
              i34=4; i56=6
          END SELECT
          np=mx+1
          CALL get_string(i34,np,xstr,y34,tstr,ustr,fstr)
          CALL get_string(i56,np,xstr,y56,tstr,ustr,fstr)
          DO j=0,my
            eta=etalower+j*deta
            tstr=0.5*(eta+1)*(y56-y34)+y34
            Info%qbuf(0:mx,j,1,1) = tstr
          END DO          
      END SELECT 
    END SUBROUTINE problemIO_gas

    SUBROUTINE problemIO_string(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,UserNrOutVars=-4      
      !
      SELECT CASE (IOrequest)        
        CASE (UserNrOutVars) ! Suppress output of string field variables
          Info%NrOutVars=0
      END SELECT
    END SUBROUTINE problemIO_string
    
    SUBROUTINE setaux(Info)  
      ! Interface declarations    
      TYPE (NodeInfo) :: Info
      INTEGER i,j,mbc,r,rmbc,mcsi,meta,iError,i34,i56,np
      REAL (KIND=xPrec) :: t,csi,eta,dcsi,deta,csilower,etalower
      REAL (KIND=qPrec), ALLOCATABLE, DIMENSION(:) :: &
                         xstring,y34,u34,theta34,y56,u56,theta56,fstring
      !
      IF (Info%nEquationSet>2) RETURN
      ! Set auxilliary variables for gas dynamics on moving mesh
      Info%aux=0.                           
      ! Get information about current grid                           
      mbc=Info%mbc; r=Info%r; rmbc=r*mbc      
      mcsi=Info%mX(1); meta=Info%mX(2) ; dcsi=Info%dX(1); deta=Info%dX(2)     
      csilower=Info%Xlower(1); etalower=Info%Xlower(2)
      ALLOCATE(xstring(1-rmbc:mcsi+rmbc),y34(1-rmbc:mcsi+rmbc),u34(1-rmbc:mcsi+rmbc), &
               theta34(1-rmbc:mcsi+rmbc),y56(1-rmbc:mcsi+rmbc),u56(1-rmbc:mcsi+rmbc), &
               theta56(1-rmbc:mcsi+rmbc),fstring(1-rmbc:mcsi+rmbc),STAT=iError)
      IF (iError /= 0) THEN
        PRINT *,'Cannot allocate temporary vectors in setaux'
        STOP
      END IF
      ! Get current string positions      
      DO i=1-rmbc,mcsi+rmbc
        xstring(i)=csilower+(i-0.5)*dcsi
      END DO
      SELECT CASE(Info%nEquationSet)
        CASE (1)
          i34=3; i56=5
        CASE (2)
          i34=4; i56=6
      END SELECT
      np=mcsi+2*rmbc
      CALL get_string(i34,np,xstring,y34,theta34,u34,fstring)
      CALL get_string(i56,np,xstring,y56,theta56,u56,fstring)
      ! Compute initial values of metric coefficients      
      DO i=1-rmbc,mcsi+rmbc
        csi=csilower+(i-0.5)*dcsi
        DO j=1-rmbc,meta+rmbc
          eta=etalower+(j-0.5)*deta          
          Info%aux(i,j,1,1,1)= 1.0   ! aux(i,j,1) contains X_csi(i,j)
          Info%aux(i,j,1,1,2)= 0.0   ! aux(i,j,2) contains X_eta(i,j)          
          Info%aux(i,j,1,1,3)= &     ! aux(i,j,3) contains Y_csi(i,j)
            0.5*(eta+1)*(theta56(i)-theta34(i))+theta34(i)
          Info%aux(i,j,1,1,4)= &     ! aux(i,j,4) contains Y_eta(i,j)
            0.5*(y56(i)-y34(i))
          Info%aux(i,j,1,1,5)= 0.0   ! aux(i,j,5) contains X_t(i,j)
          Info%aux(i,j,1,1,6)= &     ! aux(i,j,6) contains Y_t(i,j)
            0.5*(eta+1)*(u56(i)-u34(i))+u34(i)
          ! Jacobian = X_csi Y_eta - X_eta Y_csi
          Info%aux(i,j,1,1,7) = Info%aux(i,j,1,1,1)*Info%aux(i,j,1,1,4) - &
                                Info%aux(i,j,1,1,2)*Info%aux(i,j,1,1,3)
        END DO
      END DO      
      DEALLOCATE(xstring,y34,u34,theta34,y56,u56,theta56,fstring)
    END SUBROUTINE setaux
    
    SUBROUTINE src(Info)
      ! Interface declarations    
      TYPE (NodeInfo) :: Info
      ! Include code to compute source terms here:
      ! Trapezoidal rule integration of q_t = f 
      INTEGER mx,iError,i,mbc,rmbc,r,np,istring
      REAL (KIND=qPrec) :: dt,dx,dx2
      REAL (KIND=qPrec), ALLOCATABLE, DIMENSION(:) :: x,y,theta,u,f
      IF (Info%nEquationSet<3) RETURN
      ! Source term for string equation      
      mx=Info%mX(1); dt=Info%dt;
      dx=Info%dX(1); dx2=0.5*dx; mbc=Info%mbc; r=Info%r; rmbc=r*mbc
      ALLOCATE(x(1-rmbc:mx+rmbc),y(1-rmbc:mx+rmbc),u(1-rmbc:mx+rmbc), &
               theta(1-rmbc:mx+rmbc),f(1-rmbc:mx+rmbc),STAT=iError)
      IF (iError /=0) THEN
        PRINT *,'Cannot allocate work array in src'
        STOP
      END IF      
      DO i=1-rmbc,mx+rmbc
        x(i) = Info%Xlower(1)+dx2+(i-1)*dx
      END DO
      istring=Info%nEquationSet; np=mx+2*rmbc
      CALL get_string(istring,np,x,y,theta,u,f)
      Info%q(:,1,1,1,2) = Info%q(:,1,1,1,2) + dt*f(:)
      DEALLOCATE(x,y,theta,u,f)
    END SUBROUTINE src

    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
      !
      SELECT CASE(Info%nEquationSet)
        CASE (1:2)
          CALL physflux_gas(ixy,indx,irequest,Info,q,f,s,A)
        CASE (3:6)
          CALL physflux_string(ixy,indx,irequest,Info,q,f,s,A)
      END SELECT
    END SUBROUTINE physflux
    
    SUBROUTINE physflux_gas(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 = .FALSE.
      INTEGER :: i,j,iError
      INTEGER, SAVE :: mbc,mx,nq,mw,NrVars,mwaves
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION(:) :: delta,coef
      ! Cell center quantities
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION (:) :: sqrtrho,pres, &
                         u,v,c,h,rho
      ! Roe averaged quantities
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION (:) :: uR,vR,hR,cR, &
                         rhosq2,g1c2,u2v2,euv,kx,ky,kV,dels,xeta,xcsi, &
                         yeta,ycsi,xt,yt
      ! Shorter local names      
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: Apdq,Amdq,Asdq,BpAsdq,BmAsdq, &
                                                     speed,q1D,aux1D
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: wave      
      REAL (KIND=qPrec) :: rho1,rhou1,en1,p1,c1,rho2,rhou2,en2,p2,c2
      REAL (KIND=qPrec) :: s0,s1,s2,s3,sfract,df(4),rhov1,rhov2
      !   
      ! Current number of interior cells, boundary cells, field variables, waves
      mx=Info%mXnow; mbc=Info%mbc; NrVars=Info%NrVars; mwaves=Info%mwaves
      ! Associate local names with Info components
      q1D=>Info%q1D; aux1D=>Info%aux1D
      Apdq=>Info%Apdq; Amdq=>Info%Amdq; speed=>Info%speed; wave=>Info%wave
      Asdq=>Info%Asdq; BpAsdq=>Info%BpAsdq; BmAsdq=>Info%BmAsdq      
      SELECT CASE (irequest)
        CASE (Initialize)
          ! Allocate local work space
          ALLOCATE(delta(NrVars),coef(mwaves),                                             &
                   sqrtrho(1-mbc:mx+mbc),pres(1-mbc:mx+mbc),rho(1-mbc:mx+mbc),             &
                   u(1-mbc:mx+mbc),v(1-mbc:mx+mbc),h(1-mbc:mx+mbc),c(1-mbc:mx+mbc),        &
                   uR(2-mbc:mx+mbc),vR(2-mbc:mx+mbc),hR(2-mbc:mx+mbc),cR(2-mbc:mx+mbc),    &
                   rhosq2(2-mbc:mx+mbc),g1c2(2-mbc:mx+mbc),u2v2(2-mbc:mx+mbc),             &
                   euv(2-mbc:mx+mbc),kx(2-mbc:mx+mbc),ky(2-mbc:mx+mbc),kV(2-mbc:mx+mbc),   &
                   dels(2-mbc:mx+mbc),xeta(2-mbc:mx+mbc),xcsi(2-mbc:mx+mbc),               &
                   yeta(2-mbc:mx+mbc),ycsi(2-mbc:mx+mbc),xt(2-mbc:mx+mbc),yt(2-mbc:mx+mbc),&
                   STAT=iError)
          IF (iError/=OK) THEN
            PRINT *,'Cannot allocate work space in problem module, physflux'
            STOP
          END IF
        CASE (Finalize)
          ! DeAllocate local work space
          DEALLOCATE( delta,coef,sqrtrho,pres,rho,u,v,h,c,uR,vR,hR,cR,rhosq2,g1c2, &
                      u2v2,euv,kx,ky,kV,dels,xeta,xcsi,yeta,ycsi,xt,yt,STAT=iError)
          IF (iError/=OK) THEN
            PRINT *,'Cannot deallocate 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)
            PRINT *,'indx=',indx
            STOP
          END IF          
          ! Compute cell centered quantities required in Riemann solve & entropy fix
          rho(1-mbc:mx+mbc) = 1./q1D(1-mbc:mx+mbc,1)
          u(1-mbc:mx+mbc) = q1D(1-mbc:mx+mbc,2)*rho(1-mbc:mx+mbc)
          v(1-mbc:mx+mbc) = q1D(1-mbc:mx+mbc,3)*rho(1-mbc:mx+mbc)
          rho(1-mbc:mx+mbc) = q1D(1-mbc:mx+mbc,1)
          sqrtrho(1-mbc:mx+mbc) = SQRT(rho(1-mbc:mx+mbc))
          pres(1-mbc:mx+mbc) = gamma1 * ( q1D(1-mbc:mx+mbc,4) - &
              half*(q1D(1-mbc:mx+mbc,2)*u(1-mbc:mx+mbc) + &
                    q1D(1-mbc:mx+mbc,3)*v(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,4) + 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)*v(1-mbc:mx+mbc)
          f(1-mbc:mx+mbc,4)=f(1-mbc:mx+mbc,1)*h(1-mbc:mx+mbc)                    
        CASE (RequestNormalWaves)
          Apdq=0.; Amdq=0.; speed=0.; wave=0.
          ! Compute Roe averages
          rhosq2(2-mbc:mx+mbc) = one / (sqrtrho(1-mbc:mx+mbc-1)+sqrtrho(2-mbc:mx+mbc))
          uR(2-mbc:mx+mbc) = ( q1D(1-mbc:mx+mbc-1,2)/sqrtrho(1-mbc:mx+mbc-1) + &
                               q1D(2-mbc:mx+mbc,2)/sqrtrho(2-mbc:mx+mbc) )  *  &
                             rhosq2(2-mbc:mx+mbc)
          vR(2-mbc:mx+mbc) = ( q1D(1-mbc:mx+mbc-1,3)/sqrtrho(1-mbc:mx+mbc-1) + &
                               q1D(2-mbc:mx+mbc,3)/sqrtrho(2-mbc:mx+mbc) )  *  &
                             rhosq2(2-mbc:mx+mbc)          
          hR(2-mbc:mx+mbc) =                                                          &
            ( (q1D(1-mbc:mx+mbc-1,4)+pres(1-mbc:mx+mbc-1))/sqrtrho(1-mbc:mx+mbc-1) +   &
              (q1D(2-mbc:mx+mbc,4)+pres(2-mbc:mx+mbc))/sqrtrho(2-mbc:mx+mbc)       ) * &
            rhosq2(2-mbc:mx+mbc)          
          u2v2(2-mbc:mx+mbc)=uR(2-mbc:mx+mbc)**2+vR(2-mbc:mx+mbc)**2
          cR(2-mbc:mx+mbc) = gamma1*(hR(2-mbc:mx+mbc)-half*u2v2(2-mbc:mx+mbc))
          IF (MINVAL(cR(2-mbc:mx+mbc))<=zero) THEN
            PRINT *,'Error: Negative or zero Roe average sound velocity in physflux.'            
            PRINT '(1x,a6,1x,32(f7.4,1x))','c**2= ',cR(2-mbc:mx+mbc)            
            STOP
          END IF
          g1c2(2-mbc:mx+mbc) = gamma1/cR(2-mbc:mx+mbc)
          cR(2-mbc:mx+mbc) = SQRT(cR(2-mbc:mx+mbc))
          euv(2-mbc:mx+mbc) = hR(2-mbc:mx+mbc)-u2v2(2-mbc:mx+mbc)          
          ! Construct the eigenbasis          
          ! Find current projection directions, metric coefficients
          xt(2-mbc:mx+mbc)=aux1D(2-mbc:mx+mbc,5)
          yt(2-mbc:mx+mbc)=aux1D(2-mbc:mx+mbc,6)          
          SELECT CASE(ixy)
            CASE (1)
              xeta(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,2)
              yeta(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,4)
              kx = yeta; ky = -xeta              
              dels = xt*yeta - xeta*yt
            CASE (2)
              xcsi(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,1)
              ycsi(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,3)
              kx = -ycsi; ky = xcsi
              dels = xcsi*yt - xt*ycsi
          END SELECT
          kV = 1./SQRT(kx**2+ky**2); kx=kV*kx; ky=kV*ky
          kV = kx*uR + ky*vR
          ! Backward acoustic
          speed(2-mbc:mx+mbc,1) = kV(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,1) = one
          wave(2-mbc:mx+mbc,2,1) = uR(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc)*kx(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,1) = vR(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc)*ky(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,4,1) = hR(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc)*kV(2-mbc:mx+mbc)
          ! Shear
          speed(2-mbc:mx+mbc,2) = kV(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,2) = zero
          wave(2-mbc:mx+mbc,2,2) = -ky(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,2) = kx(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,4,2) = kx(2-mbc:mx+mbc)*vR(2-mbc:mx+mbc) - &
                                   ky(2-mbc:mx+mbc)*uR(2-mbc:mx+mbc)
          ! Entropy
          speed(2-mbc:mx+mbc,3) = kV(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,3) = one
          wave(2-mbc:mx+mbc,2,3) = uR(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,3) = vR(2-mbc:mx+mbc) 
          wave(2-mbc:mx+mbc,4,3) = half*u2v2(2-mbc:mx+mbc)
          ! Forward acoustic
          speed(2-mbc:mx+mbc,4) = kV(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,4) = one
          wave(2-mbc:mx+mbc,2,4) = uR(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc)*kx(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,4) = vR(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc)*ky(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,4,4) = hR(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc)*kV(2-mbc:mx+mbc)          
          ! Decompose jump in q onto eigenbases
          DO j=2-mbc,mx+mbc
            ! Jump in q
            delta(1:4) = q1D(j,1:4) - q1D(j-1,1:4)
            ! Find coef(1) thru coef(4), the coefficients of the 4 eigenvectors
            coef(1) =   (cR(j)*kV(j) + half*gamma1*u2v2(j)) * delta(1)   &
                      - (cR(j)*kx(j)+gamma1*uR(j)) * delta(2)            &
                      - (cR(j)*ky(j)+gamma1*vR(j)) * delta(3)            &
                      + gamma1 * delta(4)
            coef(1) = half*coef(1)/cR(j)**2
            coef(2) = (ky(j)*uR(j)-kx(j)*vR(j))*delta(1) - &
                      ky(j)*delta(2) + kx(j)*delta(3)
            coef(3) =   (cR(j)**2 - half*gamma1*u2v2(j)) * delta(1)       &
                      + gamma1*uR(j)*delta(2) + gamma1*vR(j)*delta(3)     &
                      - gamma1*delta(4)
            coef(3) = coef(3)/cR(j)**2
            coef(4) =   (-cR(j)*kV(j) + half*gamma1*u2v2(j)) * delta(1)   &
                      + (cR(j)*kx(j)-gamma1*uR(j)) * delta(2)             &
                      + (cR(j)*ky(j)-gamma1*vR(j)) * delta(3)             &
                      + gamma1 * delta(4)
            coef(4) = half*coef(4)/cR(j)**2
            DO mw=1,mwaves
              wave(j,1:NrVars,mw) = coef(mw)*wave(j,1:NrVars,mw)              
            END DO
          END DO
          ! Compute Godunov flux f0
          IF (efix) THEN
            ! With entropy fix
            ! Compute flux differences amdq and apdq.
            ! First compute amdq as sum of s*wave for left going waves.
            ! Incorporate entropy fix by adding a modified fraction of wave
            ! if s should change sign.
            DO j=2-mbc,mx+mbc
              ! 1-wave. Check for fully supersonic case
              s0=u(j-1)-c(j-1)
              IF ( s0>=zero .AND. speed(j,1) > zero ) THEN
                ! Everything is right-going
                Amdq(j,1:NrVars) = zero  
                CYCLE
              END IF
              ! u-c to right of 1-wave
              rho1 = q1D(j-1,1) + wave(j,1,1)
              rhou1 = q1D(j-1,2) + wave(j,2,1)
              rhov1 = q1D(j-1,3) + wave(j,3,1)
              en1 = q1D(j-1,4) + wave(j,4,1)
              p1 = gamma1*(en1  - 0.5*(rhou1**2+rhov1**2)/rho1)
              c1 = SQRT(gamma*p1/rho1)
              s1 = rhou1/rho1 - c1
              IF ( s0<zero .AND. s1>0 ) THEN
                ! Transonic rarefaction in the 1-wave                
                sfract = s0 * (s1-speed(j,1)) / (s1-s0)
              ELSE IF (speed(j,1) < zero) THEN
                ! Left-going 1-wave
                sfract=speed(j,1)
              ELSE
                ! Right-going 1-wave
                sfract=zero  ! Never should reach this instruction
              END IF
              Amdq(j,1:NrVars) = sfract*wave(j,1:NrVars,1)
              
              ! Check 2,3-wave (contact discontinuity)
              IF (speed(j,2) >= zero) CYCLE ! 2- and 3- wave is right-going
              Amdq(j,1:NrVars) = Amdq(j,1:NrVars) + speed(j,2)*wave(j,1:NrVars,2) + &
                                 speed(j,3)*wave(j,1:NrVars,3)
              
              ! Check 4-wave. Compute u+c to left of 4-wave
              rho2 = q1D(j,1) - wave(j,1,4)
              rhou2 = q1D(j,2) - wave(j,2,4)
              rhov2 = q1D(j,3) - wave(j,3,4)
              en2 = q1D(j,4) - wave(j,3,4)
              p2 = gamma1*(en2 - 0.5d0*(rhou2**2+rhov2**2)/rho2)
              c2 = SQRT(gamma*p2/rho2)
              s2 = rhou2/rho2 + c2; s3=u(j)+c(j)
              IF ( s2 < zero .AND. s3 > zero ) THEN
                ! Transonic rarefaction in the 4-wave
                sfract = s2 * (s3-speed(j,4)) / (s3-s2)
              ELSE IF (speed(j,4) < zero) THEN
                ! Left-going 4-wave
                sfract = speed(j,4)
              ELSE
                ! Right-going 4-wave
                CYCLE
              END IF
              Amdq(j,1:NrVars) = Amdq(j,1:NrVars) + sfract*wave(j,1:NrVars,4)
              
            END DO

            ! Compute the right-going flux-differences
            DO j=2-mbc,mx+mbc
              df=zero
              DO mw=1,mwaves
                df(1:NrVars) = df(1:NrVars) + speed(j,mw) * wave(j,1:NrVars,mw)
              END DO
              Apdq(j,1:NrVars) = df(1:NrVars) - Amdq(j,1:NrVars)
            END DO
          !===
          ELSE
          !===
            ! No entropy fix
            ! amdq = SUM s*wave   over left-going waves
            ! apdq = SUM s*wave   over right-going waves            
            DO mw=1,mwaves              
              DO j=2-mbc,mx+mbc
                IF (speed(j,mw) < zero) THEN
                  Amdq(j,1:NrVars) = Amdq(j,1:NrVars) +  &
                                     speed(j,mw) * wave(j,1:NrVars,mw)
                ELSE
                  Apdq(j,1:NrVars) = Apdq(j,1:NrVars) +  &
                                     speed(j,mw) * wave(j,1:NrVars,mw)
                END IF
              END DO              
            END DO            
          END IF
        CASE (RequestTransverseWaves)
          BpAsdq=0.; BmAsdq=0.; speed=0.; wave=0.
          ! Construct the eigenbasis          
          ! Find current projection directions, metric coefficients
          xt(2-mbc:mx+mbc)=aux1D(2-mbc:mx+mbc,5)
          yt(2-mbc:mx+mbc)=aux1D(2-mbc:mx+mbc,6)          
          SELECT CASE(ixy)
            CASE (1)
              xcsi(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,1)
              ycsi(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,3)
              kx = -ycsi; ky = xcsi
              dels = xcsi*yt - xt*ycsi              
            CASE (2)
              xeta(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,2)
              yeta(2-mbc:mx+mbc) = aux1D(2-mbc:mx+mbc,4)
              kx = yeta; ky = -xeta              
              dels = xt*yeta - xeta*yt
          END SELECT
          kV = 1./SQRT(kx**2+ky**2); kx=kV*kx; ky=kV*ky
          kV = kx*uR + ky*vR
          ! Backward acoustic
          speed(2-mbc:mx+mbc,1) = kV(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,1) = one
          wave(2-mbc:mx+mbc,2,1) = uR(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc)*kx(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,1) = vR(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc)*ky(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,4,1) = hR(2-mbc:mx+mbc) - cR(2-mbc:mx+mbc)*kV(2-mbc:mx+mbc)
          ! Shear
          speed(2-mbc:mx+mbc,2) = kV(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,2) = zero
          wave(2-mbc:mx+mbc,2,2) = -ky(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,2) = kx(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,4,2) = kx(2-mbc:mx+mbc)*vR(2-mbc:mx+mbc) - &
                                   ky(2-mbc:mx+mbc)*uR(2-mbc:mx+mbc)
          ! Entropy
          speed(2-mbc:mx+mbc,3) = kV(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,3) = one
          wave(2-mbc:mx+mbc,2,3) = uR(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,3) = vR(2-mbc:mx+mbc) 
          wave(2-mbc:mx+mbc,4,3) = half*u2v2(2-mbc:mx+mbc)
          ! Forward acoustic
          speed(2-mbc:mx+mbc,4) = kV(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc) - dels(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,1,4) = one
          wave(2-mbc:mx+mbc,2,4) = uR(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc)*kx(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,3,4) = vR(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc)*ky(2-mbc:mx+mbc)
          wave(2-mbc:mx+mbc,4,4) = hR(2-mbc:mx+mbc) + cR(2-mbc:mx+mbc)*kV(2-mbc:mx+mbc)          
          ! Decompose jump in q onto eigenbases
          DO j=2-mbc,mx+mbc
            ! Jump in q
            delta(1:4) = Asdq(j,1:4)
            ! Find coef(1) thru coef(4), the coefficients of the 4 eigenvectors
            coef(1) =   (cR(j)*kV(j) + half*gamma1*u2v2(j)) * delta(1)   &
                      - (cR(j)*kx(j)+gamma1*uR(j)) * delta(2)            &
                      - (cR(j)*ky(j)+gamma1*vR(j)) * delta(3)            &
                      + gamma1 * delta(4)
            coef(1) = half*coef(1)/cR(j)**2
            coef(2) = (ky(j)*uR(j)-kx(j)*vR(j))*delta(1) - &
                      ky(j)*delta(2) + kx(j)*delta(3)
            coef(3) =   (cR(j)**2 - half*gamma1*u2v2(j)) * delta(1)       &
                      + gamma1*uR(j)*delta(2) + gamma1*vR(j)*delta(3)     &
                      - gamma1*delta(4)
            coef(3) = coef(3)/cR(j)**2
            coef(4) =   (-cR(j)*kV(j) + half*gamma1*u2v2(j)) * delta(1)   &
                      + (cR(j)*kx(j)-gamma1*uR(j)) * delta(2)             &
                      + (cR(j)*ky(j)-gamma1*vR(j)) * delta(3)             &
                      + gamma1 * delta(4)
            coef(4) = half*coef(4)/cR(j)**2
            DO mw=1,mwaves
              wave(j,1:NrVars,mw) = coef(mw)*wave(j,1:NrVars,mw)              
            END DO
          END DO
          ! Compute the transverse fluctuations            
          DO mw=1,mwaves              
              DO j=2-mbc,mx+mbc
                IF (speed(j,mw) < zero) THEN
                  BmAsdq(j,1:NrVars) = BmAsdq(j,1:NrVars) +  &
                                     speed(j,mw) * wave(j,1:NrVars,mw)
                ELSE
                  BpAsdq(j,1:NrVars) = BpAsdq(j,1:NrVars) +  &
                                     speed(j,mw) * wave(j,1:NrVars,mw)
                END IF
              END DO              
            END DO            
        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_gas

    SUBROUTINE physflux_string(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
      REAL (KIND=qPrec) :: invc1,c1
      REAL (KIND=qPrec), ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta,coef
      ! Shorter local names for the wave propagation quantities
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:) :: Apdq,Amdq,speed,aux1D
      REAL (KIND=qPrec), POINTER, DIMENSION (:,:,:) :: wave
      INTEGER mbc,mx,nq,n,iError
      INTEGER, POINTER, DIMENSION(:) :: I
      ! Current number of interior cells, boundary cells
      mx=Info%mXnow   ! This is the number of cells in the current 1D slice
                      ! This may be different from any Info%mX(n) from bear.data
                      ! Always use this procedure to set the number of cells
      mbc=Info%mbc    ! This is the number ghost cells added at each boundary
      nq=Info%NrVars  ! This is the number of components of q
      ! Index range - contains indices where a request should be answered
      I=>Info%I1D     ! I is set to point to the same locations as Info%I1D
      ! Associate local names with Info components
      Apdq=>Info%Apdq; Amdq=>Info%Amdq; speed=>Info%speed; wave=>Info%wave
      Apdq=0.; Amdq=0.; speed=0.; wave=0.
      aux1D=>Info%aux1D
      SELECT CASE (irequest)
        CASE (Initialize)
          ! You can carry out problem specific initialization here.
          ! Examples: allocation of dynamic arrays, computation of
          !           quantities that are needed in successive calls
          !           to physflux routine. 
          ! Remember to add the SAVE attribute to quantities that you 
          ! want to be available in successive physflux calls
          ALLOCATE (delta(2-mbc:mx+mbc,2),coef(2-mbc:mx+mbc,2),STAT=iError)
          IF (iError /=0) THEN
            PRINT *,'Cannot allocate work array in physflux'
            STOP
          END IF
        CASE (Finalize)
          ! Problem specific clean-up, e.g. deallocation of arrays
          DEALLOCATE (delta,coef,STAT=iError)
          IF (iError/=0) THEN
            PRINT *,'Cannot allocate work space in problem module, physflux'
            STOP
          END IF
        CASE (RequestFluxes)
          ! Used in solver modules: ClassicBEAR, CompactBEAR
          ! Write code that computes f(1-mbc:mx+mbc,1:nq), the physical fluxes
          ! at the cell centers. This may also be written as f(I,1:nq) with I=>Info%I1D
          ! f(I,1:nq) = 
        CASE (RequestSpeeds)
          ! Used in solver modules: ClassicBEAR, CompactBEAR
          ! Write code that computes s(1-mbc:mx+mbc,1:nq), the physical speeds
          ! at the cell centers. This may also be written as s(I,1:nq) with I=>Info%I1D
          ! s(I,1:nq) =           
        CASE (RequestJacobians)
          ! Used in solver modules: ClassicBEAR, CompactBEAR
          ! Write code that computes A(1-mbc:mx+mbc,1:nq,1:nq), the flux Jacobians
          ! at the cell centers. This may also be written as A(I,1:nq,1:nq) with I=>Info%I1D
          ! A(I,1:nq,1:nq) =                     
        CASE (RequestNormalWaves)
          ! Specific to WAVEBEAR solver, the wave propagation algorithms.
          ! Compute the solution to the normal Riemann problem, see User Manual for details
          SELECT CASE(Info%nEquationSet)
            CASE (3:4)
              c1=cstring2
            CASE (5:6)
              c1=cstring1
          END SELECT
          invc1=1./c1
          delta(I,:)=q(I,:)-q(I-1,:); delta(I,2)=invc1*delta(I,2)
          coef(I,1)=0.5*(delta(I,1) + delta(I,2))
          coef(I,2)=0.5*(delta(I,1) - delta(I,2))
          ! Waves:
          speed(I,1) = -c1
          wave(I,1,1) = coef(I,1)
          wave(I,2,1) = coef(I,1)*c1

          speed(I,2) =  c1
          wave(I,1,2) = coef(I,2)
          wave(I,2,2) = -coef(I,2)*c1
          ! Fluctuations:          
          DO n=1,2  ! Use wave propagation update for theta, u only
            Amdq(I,n) = speed(I,1)*wave(I,n,1)
            Apdq(I,n) = speed(I,2)*wave(I,n,2)
          END DO
          Amdq(I,3)=0.; Apdq(I,3)=0.  
        CASE (RequestTransverseWaves)
          ! Specific to WAVEBEAR solver, the wave propagation algorithms.
          ! Compute the transverse correction waves, see User Manual for details
        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_string
              
!*****************
END MODULE problem
!*****************