module characteristics
  use common_module
  implicit none
contains
  subroutine characteristic_eigenmatrix(V8,evL,evR,axis)
    ! calculate the characteristics values with eigenmatrix
    ! cf. Powell et al. 1999 JCoPh
    implicit none
    real(DP), intent(IN) :: V8(8)    ! primitive variables
    real(DP), dimension(8,8),intent(OUT) :: evL,evR
    integer,intent(IN) :: axis
    real(DP) :: rho,rhoi,pre,v3(3),b3(3)
    real(DP) :: a,cf,cs,a2,ca2,cb2,cf2,cs2,alpf,alps,bet2,bet3,cf2_cs2
    real(DP) :: case1,case2,case3,sgnB,vn2,Bn2,Bt2
    real(DP) :: rt2,rt2i,rtrho,rtrhoi
    integer :: axis2,axis3,axv(3),axb(3),rw,cl,kdlt(3),kdlt2(3),kdlt3(3)
    axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1
    kdlt=0**abs(axis-(/1,2,3/)) ; kdlt2=0**abs(axis2-(/1,2,3/))
    kdlt3=(/1,1,1/)-kdlt-kdlt2
    axv=(/Vax(axis),Vax(axis2),Vax(axis3)/)
    axb=(/Bax(axis),Bax(axis2),Bax(axis3)/)
    rt2=sqrt(2._DP); rt2i=1._DP/rt2
    rho=V8(Sax(1)); rhoi=1._DP/rho; rtrho=sqrt(rho); rtrhoi=1._DP/rtrho
    v3=(/V8(axv(1)),V8(axv(2)),V8(axv(3))/)
    b3=(/V8(axb(1)),V8(axb(2)),V8(axb(3))/)
    vn2=sum(v3*v3); Bn2=sum(b3*b3); Bt2=b3(2)*b3(2)+b3(3)*b3(3)
    pre=V8(Sax(2)); pre=max(pre,pressure_inf)
    a2=gamma*pre/rho ; ca2=b3(1)*b3(1)/rho ; cb2=Bn2/rho
    case1=1._DP-eq0(Bt2)! Bt2=0 (<--> by=bz=0)
    case2=ge0(a2-cb2)! a>=cb
    case3=lt0(a2-cb2)! a<cb
    cf2=case1*.5_DP*(a2+cb2+sqrt((a2+cb2)*(a2+cb2)-4._DP*a2*ca2))+&
         (1._DP-case1)*(case2*a2+case3*ca2)
    cs2=case1*a2*ca2/cf2+(1._DP-case1)*(case2*ca2+case3*a2)
    a=sqrt(a2); cf=sqrt(cf2); cs=sqrt(cs2)
    cf2_cs2=cf2-cs2
    case2=0._DP**cf2_cs2 ; case3=1._DP-case2
    alpf=case2*rt2i+case3*sqrt(max(a2-cs2,0._DP)/(cf2_cs2+case2))
    alps=case2*rt2i+case3*sqrt(max(cf2-a2,0._DP)/(cf2_cs2+case2))
    case2=1._DP-case1
    bet2=case1*b3(2)/sqrt(Bt2+case2)+case2*rt2i
    bet3=case1*b3(3)/sqrt(Bt2+case2)+case2*rt2i
    sgnB=sign(1._DP,b3(1))
    ! primitive to eigenmatrix
    ! Entropy wave [cf. Eq.(43) of Powell et al. 1999]
    evL(Sax(1),Sax(1))=1._DP ; evL(Sax(1),Sax(2))=-1._DP/a2
    evL(Sax(1),axv)=0._DP ; evL(Sax(1),axb)=0._DP
    ! Leftward Slow-mode wave [cf. Eq.(47)]
    evL(Sax(2),Sax(1))=0._DP
    evL(Sax(2),Sax(2))=alps*.5_DP*rhoi/a2
    evL(Sax(2),axv(1))=-alps*cs*.5_DP/a2
    evL(Sax(2),axv(2))=-alpf*cf*bet2*sgnB*.5_DP/a2
    evL(Sax(2),axv(3))=-alpf*cf*bet3*sgnB*.5_DP/a2
    evL(Sax(2),axb(1))=0._DP
    evL(Sax(2),axb(2))=-alpf*bet2*.5_DP*rtrhoi/a
    evL(Sax(2),axb(3))=-alpf*bet3*.5_DP*rtrhoi/a
    ! Rightward Alfven-mode wave [cf. Eq.(45)]
    evL(axv(1),Sax(1))=0._DP
    evL(axv(1),Sax(2))=0._DP
    evL(axv(1),axv(1))=0._DP
    evL(axv(1),axv(2))=-bet3*rt2i
    evL(axv(1),axv(3))=bet2*rt2i
    evL(axv(1),axb(1))=0._DP
    evL(axv(1),axb(2))=bet3*rt2i*rtrhoi
    evL(axv(1),axb(3))=-bet2*rt2i*rtrhoi
    ! Leftward Alfven-mode wave [cf. Eq.(45)]
    evL(axv(2),Sax(1))=0._DP
    evL(axv(2),Sax(2))=0._DP
    evL(axv(2),axv(1))=0._DP
    evL(axv(2),axv(2))=-bet3*rt2i
    evL(axv(2),axv(3))=bet2*rt2i
    evL(axv(2),axb(1))=0._DP
    evL(axv(2),axb(2))=-bet3*rt2i*rtrhoi
    evL(axv(2),axb(3))=bet2*rt2i*rtrhoi
    ! Rightward Fast-mode wave [cf. Eq.(46)]
    evL(axv(3),Sax(1))=0._DP
    evL(axv(3),Sax(2))=alpf*.5_DP*rhoi/a2
    evL(axv(3),axv(1))=alpf*cf*.5_DP/a2
    evL(axv(3),axv(2))=-alps*cs*bet2*sgnB*.5_DP/a2
    evL(axv(3),axv(3))=-alps*cs*bet3*sgnB*.5_DP/a2
    evL(axv(3),axb(1))=0._DP
    evL(axv(3),axb(2))=alps*bet2*.5_DP*rtrhoi/a
    evL(axv(3),axb(3))=alps*bet3*.5_DP*rtrhoi/a
    ! Magnetic flux [cf. Eq.(44)]
    evL(axb(1),:)=0._DP ; evL(axb(1),axb(1))=1._DP
    ! Leftward Fast-mode wave [cf. Eq.(46)]
    evL(axb(2),Sax(1))=0._DP
    evL(axb(2),Sax(2))=alpf*.5_DP*rhoi/a2
    evL(axb(2),axv(1))=-alpf*cf*.5_DP/a2
    evL(axb(2),axv(2))=alps*cs*bet2*sgnB*.5_DP/a2
    evL(axb(2),axv(3))=alps*cs*bet3*sgnB*.5_DP/a2
    evL(axb(2),axb(1))=0._DP
    evL(axb(2),axb(2))=alps*bet2*.5_DP*rtrhoi/a
    evL(axb(2),axb(3))=alps*bet3*.5_DP*rtrhoi/a
    ! Rightward Slow-mode wave [cf. Eq.(47)]
    evL(axb(3),Sax(1))=0._DP
    evL(axb(3),Sax(2))=alps*.5_DP*rhoi/a2
    evL(axb(3),axv(1))=alps*cs*.5_DP/a2
    evL(axb(3),axv(2))=alpf*cf*bet2*sgnB*.5_DP/a2
    evL(axb(3),axv(3))=alpf*cf*bet3*sgnB*.5_DP/a2
    evL(axb(3),axb(1))=0._DP
    evL(axb(3),axb(2))=-alpf*bet2*.5_DP*rtrhoi/a
    evL(axb(3),axb(3))=-alpf*bet3*.5_DP*rtrhoi/a

    !------------------------------------------------------------!
    !------------------------------------------------------------!
    !------------------------------------------------------------!
    ! primitive to eigenmatrix
    ! Entropy wave [cf. Eq.(43) of Powell et al. 1999]
    evR(:,Sax(1))=0._DP ; evR(Sax(1),Sax(1))=1._DP
    ! Leftward Slow-mode wave [cf. Eq.(47)]
    evR(Sax(1),Sax(2))=rho*alps
    evR(Sax(2),Sax(2))=alps*gamma*pre
    evR(axv(1),Sax(2))=-alps*cs
    evR(axv(2),Sax(2))=-alpf*cf*bet2*sgnB
    evR(axv(3),Sax(2))=-alpf*cf*bet3*sgnB
    evR(axb(1),Sax(2))=0._DP
    evR(axb(2),Sax(2))=-alpf*bet2*rtrho*a
    evR(axb(3),Sax(2))=-alpf*bet3*rtrho*a
    ! Rightward Alfven-mode wave [cf. Eq.(45)]
    evR(Sax(1),axv(1))=0._DP
    evR(Sax(2),axv(1))=0._DP
    evR(axv(1),axv(1))=0._DP
    evR(axv(2),axv(1))=-bet3*rt2i
    evR(axv(3),axv(1))=bet2*rt2i
    evR(axb(1),axv(1))=0._DP
    evR(axb(2),axv(1))=bet3*rtrho*rt2i
    evR(axb(3),axv(1))=-bet2*rtrho*rt2i
    ! Leftward Alfven-mode wave [cf. Eq.(45)]
    evR(Sax(1),axv(2))=0._DP
    evR(Sax(2),axv(2))=0._DP
    evR(axv(1),axv(2))=0._DP
    evR(axv(2),axv(2))=-bet3*rt2i
    evR(axv(3),axv(2))=bet2*rt2i
    evR(axb(1),axv(2))=0._DP
    evR(axb(2),axv(2))=-bet3*rtrho*rt2i
    evR(axb(3),axv(2))=bet2*rtrho*rt2i
    ! Rightward Fast-mode wave [cf. Eq.(46)]
    evR(Sax(1),axv(3))=rho*alpf
    evR(Sax(2),axv(3))=alpf*gamma*pre
    evR(axv(1),axv(3))=alpf*cf
    evR(axv(2),axv(3))=-alps*cs*bet2*sgnB
    evR(axv(3),axv(3))=-alps*cs*bet3*sgnB
    evR(axb(1),axv(3))=0._DP
    evR(axb(2),axv(3))=alps*bet2*rtrho*a
    evR(axb(3),axv(3))=alps*bet3*rtrho*a
    ! Magnetic flux [cf. Eq.(44)]
    evR(:,axb(1))=0._DP ; evR(axb(1),axb(1))=1._DP
    ! Leftward Fast-mode wave [cf. Eq.(46)]
    evR(Sax(1),axb(2))=rho*alpf
    evR(Sax(2),axb(2))=alpf*gamma*pre
    evR(axv(1),axb(2))=-alpf*cf
    evR(axv(2),axb(2))=alps*cs*bet2*sgnB
    evR(axv(3),axb(2))=alps*cs*bet3*sgnB
    evR(axb(1),axb(2))=0._DP
    evR(axb(2),axb(2))=alps*bet2*rtrho*a
    evR(axb(3),axb(2))=alps*bet3*rtrho*a
    ! Rightward Slow-mode wave [cf. Eq.(47)]
    evR(Sax(1),axb(3))=rho*alps
    evR(Sax(2),axb(3))=alps*gamma*pre
    evR(axv(1),axb(3))=alps*cs
    evR(axv(2),axb(3))=alpf*cf*bet2*sgnB
    evR(axv(3),axb(3))=alpf*cf*bet3*sgnB
    evR(axb(1),axb(3))=0._DP
    evR(axb(2),axb(3))=-alpf*bet2*rtrho*a
    evR(axb(3),axb(3))=-alpf*bet3*rtrho*a
  end subroutine characteristic_eigenmatrix

  subroutine jacobian_primitive_conservative(U8,dVdU,dUdV,axis)
    implicit none
    real(DP), intent(IN) :: U8(8)    !conservative variables
    integer,intent(IN) :: axis
    real(DP), dimension(8,8),intent(OUT) :: dUdV,dVdU
    real(DP) :: rho,rhoi,pre,v3(3),b3(3)
    real(DP) :: case1,case2,case3,sgnB,vn2,Bn2,Bt2
    real(DP) :: rt2,rt2i,rtrho,rtrhoi
    integer :: axis2,axis3,axv(3),axb(3),rw,cl
    axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1
    axv=(/Vax(axis),Vax(axis2),Vax(axis3)/)
    axb=(/Bax(axis),Bax(axis2),Bax(axis3)/)
    rho=U8(Sax(1)); rhoi=1._DP/rho
    v3=(/U8(axv(1)),U8(axv(2)),U8(axv(3))/)*rhoi
    b3=(/U8(axb(1)),U8(axb(2)),U8(axb(3))/)
    vn2=sum(v3*v3); Bn2=sum(b3*b3); Bt2=b3(2)*b3(2)+b3(3)*b3(3)

    dUdV=0._DP
    dUdV(Sax(1),Sax(1))=1._DP
    dUdV(Sax(2),Sax(1))=vn2*.5_DP
    dUdV(Sax(2),Sax(2))=gamma1i
    do rw=1,3
       dUdV(Sax(2),axv(rw))=rho*v3(rw)
       dUdV(Sax(2),axb(rw))=b3(rw)
       dUdV(axv(rw),axv(rw))=rho
       dUdV(axb(rw),axb(rw))=1._DP
       dUdV(axv(rw),Sax(1))=v3(rw)
    end do
    dVdU=0._DP
    dVdU(Sax(1),Sax(1))=1._DP
    dVdU(Sax(2),Sax(1))=vn2*.5_DP/gamma1i
    dVdU(Sax(2),Sax(2))=gamma-1._DP
    do rw=1,3
       dVdU(Sax(2),axv(rw))=-v3(rw)/gamma1i
       dVdU(Sax(2),axb(rw))=-b3(rw)/gamma1i
       dVdU(axv(rw),axv(rw))=rhoi
       dVdU(axb(rw),axb(rw))=1._DP
       dVdU(axv(rw),Sax(1))=-v3(rw)*rhoi
    end do
  end subroutine jacobian_primitive_conservative

  subroutine composite_characteristic_eigenmatrix(U8,LdVdU,dUdVR,axis)
    ! calculate the characteristics values with eigenmatrix
    ! cf. Powell et al. 1999 JCoPh
    implicit none
    real(DP),intent(IN) :: U8(8)    !conservative variables
    integer,intent(IN) :: axis
    real(DP),dimension(8,8),intent(OUT) :: LdVdU,dUdVR
    real(DP) :: rho,rhoi,pre,v3(3),b3(3)
    real(DP) :: a,a2i,cf,cs,a2,ca2,cb2,cf2,cs2,alpf,alps,bet2,bet3,cf2_cs2
    real(DP) :: case1,case2,case3,sgnB,vn2,Bn2,Bt2,k1
    real(DP) :: rt2,rt2i,rtrho,rtrhoi
    integer :: axis2,axis3,axv(3),axb(3),s
    axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1
    axv=(/Vax(axis),Vax(axis2),Vax(axis3)/)
    axb=(/Bax(axis),Bax(axis2),Bax(axis3)/)
    rt2=sqrt(2._DP); rt2i=1._DP/rt2
    rho=U8(Sax(1)); rhoi=1._DP/rho; rtrho=sqrt(rho); rtrhoi=1._DP/rtrho
    v3=(/U8(axv(1)),U8(axv(2)),U8(axv(3))/)*rhoi
    b3=(/U8(axb(1)),U8(axb(2)),U8(axb(3))/)
    vn2=sum(v3*v3); Bn2=sum(b3*b3); Bt2=b3(2)*b3(2)+b3(3)*b3(3)
    pre=(U8(Sax(2))-.5_DP*(rho*vn2+Bn2))/gamma1i; pre=max(pre,pressure_inf)
    a2=gamma*pre/rho; a2i=1._DP/a2; ca2=b3(1)*b3(1)/rho; cb2=Bn2/rho
    case1=1._DP-eq0(Bt2)! Bt2=0 (<--> by=bz=0)
    case2=ge0(a2-cb2)! a>=cb
    case3=lt0(a2-cb2)! a<cb
    cf2=case1*.5_DP*(a2+cb2+sqrt((a2+cb2)*(a2+cb2)-4._DP*a2*ca2))+&
         (1._DP-case1)*(case2*a2+case3*ca2)
    cs2=case1*a2*ca2/cf2+(1._DP-case1)*(case2*ca2+case3*a2)
    a=sqrt(a2); cf=sqrt(cf2); cs=sqrt(cs2)
    k1=1._DP-gamma

    cf2_cs2=cf2-cs2
    case2=eq0(cf2_cs2); case3=1._DP-case2
    alpf=case2*rt2i+case3*sqrt(max(a2-cs2,0._DP)/(cf2_cs2+case2))
    alps=case2*rt2i+case3*sqrt(max(cf2-a2,0._DP)/(cf2_cs2+case2))
    case2=1._DP-case1
    bet2=case1*b3(2)/sqrt(Bt2+case2)+case2*rt2i
    bet3=case1*b3(3)/sqrt(Bt2+case2)+case2*rt2i
    sgnB=sign(1._DP,b3(1))

    LdVdU(Sax(1),Sax(1))=1._DP+k1*vn2*.5_DP*a2i
    LdVdU(Sax(1),Sax(2))=k1*a2i
    do s=1,3
       LdVdU(Sax(1),axv(s))=-k1*v3(s)*a2i
       LdVdU(Sax(1),axb(s))=-k1*b3(s)*a2i
    end do

    LdVdU(Sax(2),Sax(1))=(-k1*vn2*alps+2._DP*(alps*cs*v3(1)&
         +sgnB*alpf*cf*(bet2*v3(2)+bet3*v3(3))))*.25_DP*rhoi*a2i
    LdVdU(Sax(2),Sax(2))=-k1*alps*.5_DP*rhoi*a2i
    LdVdU(Sax(2),axv(1))=alps*(k1*v3(1)-cs)*.5_DP*rhoi*a2i
    LdVdU(Sax(2),axv(2))=(k1*alps*v3(2)-sgnB*alpf*cf*bet2)*.5_DP*rhoi*a2i
    LdVdU(Sax(2),axv(3))=(k1*alps*v3(3)-sgnB*alpf*cf*bet3)*.5_DP*rhoi*a2i
    LdVdU(Sax(2),axb)=k1*alps*b3*.5_DP*rhoi*a2i
    LdVdU(Sax(2),axb(2))=LdVdU(Sax(2),axb(2))-alpf*bet2*.5_DP*rtrhoi/a
    LdVdU(Sax(2),axb(3))=LdVdU(Sax(2),axb(3))-alpf*bet3*.5_DP*rtrhoi/a

    LdVdU(axv(1),Sax(1))=(bet3*v3(2)-bet2*v3(3))*rt2i*rhoi
    LdVdU(axv(1),Sax(2))=0._DP
    LdVdU(axv(1),axv(1))=0._DP
    LdVdU(axv(1),axv(2))=-bet3*rt2i*rhoi
    LdVdU(axv(1),axv(3))=bet2*rt2i*rhoi
    LdVdU(axv(1),axb(1))=0._DP
    LdVdU(axv(1),axb(2))=bet3*rt2i*rtrhoi
    LdVdU(axv(1),axb(3))=-bet2*rt2i*rtrhoi

    LdVdU(axv(2),Sax(1))=(bet3*v3(2)-bet2*v3(3))*rt2i*rhoi
    LdVdU(axv(2),Sax(2))=0._DP
    LdVdU(axv(2),axv(1))=0._DP
    LdVdU(axv(2),axv(2))=-bet3*rt2i*rhoi
    LdVdU(axv(2),axv(3))=bet2*rt2i*rhoi
    LdVdU(axv(2),axb(1))=0._DP
    LdVdU(axv(2),axb(2))=-bet3*rt2i*rtrhoi
    LdVdU(axv(2),axb(3))=bet2*rt2i*rtrhoi

    LdVdU(axv(3),Sax(1))=(-k1*vn2*alpf+2._DP*(-alpf*cf*v3(1)&
         +sgnB*alps*cs*(bet2*v3(2)+bet3*v3(3))))*.25_DP*rhoi*a2i
    LdVdU(axv(3),Sax(2))=-k1*alpf*.5_DP*rhoi*a2i
    LdVdU(axv(3),axv(1))=alpf*(k1*v3(1)+cf)*.5_DP*rhoi*a2i
    LdVdU(axv(3),axv(2))=(k1*alpf*v3(2)-sgnB*alps*cs*bet2)*.5_DP*rhoi*a2i
    LdVdU(axv(3),axv(3))=(k1*alpf*v3(3)-sgnB*alps*cs*bet3)*.5_DP*rhoi*a2i
    LdVdU(axv(3),axb)=k1*alpf*b3*.5_DP*rhoi*a2i
    LdVdU(axv(3),axb(2))=LdVdU(axv(3),axb(2))+alps*bet2*.5_DP*rtrhoi/a
    LdVdU(axv(3),axb(3))=LdVdU(axv(3),axb(3))+alps*bet3*.5_DP*rtrhoi/a

    LdVdU(axb(1),Sax)=0._DP; LdVdU(axb(1),axv)=0._DP
    LdVdU(axb(1),axb(1))=1._DP; LdVdU(axb(1),axb(2:3))=0._DP

    LdVdU(axb(2),Sax(1))=-(k1*vn2*alpf+2._DP*(-alpf*cf*v3(1)&
         +sgnB*alps*cs*(bet2*v3(2)+bet3*v3(3))))*.25_DP*rhoi*a2i
    LdVdU(axb(2),Sax(2))=-k1*alpf*.5_DP*rhoi*a2i
    LdVdU(axb(2),axv(1))=alpf*(k1*v3(1)-cf)*.5_DP*rhoi*a2i
    LdVdU(axb(2),axv(2))=(k1*alpf*v3(2)+sgnB*alps*cs*bet2)*.5_DP*rhoi*a2i
    LdVdU(axb(2),axv(3))=(k1*alpf*v3(3)+sgnB*alps*cs*bet3)*.5_DP*rhoi*a2i
    LdVdU(axb(2),axb)=k1*alpf*b3*.5_DP*rhoi*a2i
    LdVdU(axb(2),axb(2))=LdVdU(axb(2),axb(2))+alps*bet2*.5_DP*rtrhoi/a
    LdVdU(axb(2),axb(3))=LdVdU(axb(2),axb(3))+alps*bet3*.5_DP*rtrhoi/a

    LdVdU(axb(3),Sax(1))=-(k1*vn2*alps+2._DP*(alps*cs*v3(1)&
         +sgnB*alpf*cf*(bet2*v3(2)+bet3*v3(3))))*.25_DP*rhoi*a2i
    LdVdU(axb(3),Sax(2))=-k1*alps*.5_DP*rhoi*a2i
    LdVdU(axb(3),axv(1))=alps*(k1*v3(1)+cs)*.5_DP*rhoi*a2i
    LdVdU(axb(3),axv(2))=(k1*alps*v3(2)+sgnB*alpf*cf*bet2)*.5_DP*rhoi*a2i
    LdVdU(axb(3),axv(3))=(k1*alps*v3(3)+sgnB*alpf*cf*bet3)*.5_DP*rhoi*a2i
    LdVdU(axb(3),axb)=k1*alps*b3*.5_DP*rhoi*a2i
    LdVdU(axb(3),axb(2))=LdVdU(axb(3),axb(2))-alpf*bet2*.5_DP*rtrhoi/a
    LdVdU(axb(3),axb(3))=LdVdU(axb(3),axb(3))-alpf*bet3*.5_DP*rtrhoi/a


    dUdVR(Sax(1),Sax(1))=1._DP; dUdVR(Sax(2),Sax(1))=vn2*.5_DP
    dUdVR(axv,Sax(1))=v3; dUdVR(axb,Sax(1))=0._DP

    dUdVR(Sax(1),Sax(2))=rho*alps
    dUdVR(Sax(2),Sax(2))=rho*alps*(vn2*.5_DP-a2/k1-cs*v3(1))&
         -rho*alpf*cf*sgnB*(bet2*v3(2)+bet3*v3(3))&
         -alpf*rtrho*a*(bet2*b3(2)+bet3*b3(3))
    dUdVR(axv(1),Sax(2))=rho*alps*(v3(1)-cs)
    dUdVR(axv(2),Sax(2))=rho*(alps*v3(2)-alpf*cf*bet2*sgnB)
    dUdVR(axv(3),Sax(2))=rho*(alps*v3(3)-alpf*cf*bet3*sgnB)
    dUdVR(axb(1),Sax(2))=0._DP
    dUdVR(axb(2),Sax(2))=-alpf*rtrho*a*bet2
    dUdVR(axb(3),Sax(2))=-alpf*rtrho*a*bet3

    dUdVR(Sax(1),axv(1))=0._DP
    dUdVR(Sax(2),axv(1))=(rho*(bet2*v3(3)-bet3*v3(2))&
         -rtrho*(bet2*b3(3)-bet3*b3(2)))*rt2i
    dUdVR(axv(1),axv(1))=0._DP
    dUdVR(axv(2),axv(1))=-rho*bet3*rt2i
    dUdVR(axv(3),axv(1))=rho*bet2*rt2i
    dUdVR(axb(1),axv(1))=0._DP
    dUdVR(axb(2),axv(1))=rtrho*rt2i*bet3
    dUdVR(axb(3),axv(1))=-rtrho*rt2i*bet2

    dUdVR(Sax(1),axv(2))=0._DP
    dUdVR(Sax(2),axv(2))=(rho*(bet2*v3(3)-bet3*v3(2))&
         +rtrho*(bet2*b3(3)-bet3*b3(2)))*rt2i
    dUdVR(axv(1),axv(2))=0._DP
    dUdVR(axv(2),axv(2))=-rho*bet3*rt2i
    dUdVR(axv(3),axv(2))=rho*bet2*rt2i
    dUdVR(axb(1),axv(2))=0._DP
    dUdVR(axb(2),axv(2))=-rtrho*rt2i*bet3
    dUdVR(axb(3),axv(2))=rtrho*rt2i*bet2

    dUdVR(Sax(1),axv(3))=rho*alpf
    dUdVR(Sax(2),axv(3))=rho*alpf*(vn2*.5_DP-a2/k1+cf*v3(1))&
         -rho*alps*cs*sgnB*(bet2*v3(2)+bet3*v3(3))&
         +alps*rtrho*a*(bet2*b3(2)+bet3*b3(3))
    dUdVR(axv(1),axv(3))=rho*alpf*(v3(1)+cf)
    dUdVR(axv(2),axv(3))=rho*(alpf*v3(2)-alps*cs*bet2*sgnB)
    dUdVR(axv(3),axv(3))=rho*(alpf*v3(3)-alps*cs*bet3*sgnB)
    dUdVR(axb(1),axv(3))=0._DP
    dUdVR(axb(2),axv(3))=alps*rtrho*a*bet2
    dUdVR(axb(3),axv(3))=alps*rtrho*a*bet3

    dUdVR(Sax(1),axb(1))=0._DP; dUdVR(Sax(2),axb(1))=b3(1)
    dUdVR(axv,axb(1))=0._DP; dUdVR(axb(1),axb(1))=1._DP
    dUdVR(axb(2:3),axb(1))=0._DP

    dUdVR(Sax(1),axb(2))=rho*alpf
    dUdVR(Sax(2),axb(2))=rho*alpf*(vn2*.5_DP-a2/k1-cf*v3(1))&
         +rho*alps*cs*sgnB*(bet2*v3(2)+bet3*v3(3))&
         +alps*rtrho*a*(bet2*b3(2)+bet3*b3(3))
    dUdVR(axv(1),axb(2))=rho*alpf*(v3(1)-cf)
    dUdVR(axv(2),axb(2))=rho*(alpf*v3(2)+alps*cs*bet2*sgnB)
    dUdVR(axv(3),axb(2))=rho*(alpf*v3(3)+alps*cs*bet3*sgnB)
    dUdVR(axb(1),axb(2))=0._DP
    dUdVR(axb(2),axb(2))=alps*rtrho*a*bet2
    dUdVR(axb(3),axb(2))=alps*rtrho*a*bet3

    dUdVR(Sax(1),axb(3))=rho*alps
    dUdVR(Sax(2),axb(3))=rho*alps*(vn2*.5_DP-a2/k1+cs*v3(1))&
         +rho*alpf*cf*sgnB*(bet2*v3(2)+bet3*v3(3))&
         -alpf*rtrho*a*(bet2*b3(2)+bet3*b3(3))
    dUdVR(axv(1),axb(3))=rho*alps*(v3(1)+cs)
    dUdVR(axv(2),axb(3))=rho*(alps*v3(2)+alpf*cf*bet2*sgnB)
    dUdVR(axv(3),axb(3))=rho*(alps*v3(3)+alpf*cf*bet3*sgnB)
    dUdVR(axb(1),axb(3))=0._DP
    dUdVR(axb(2),axb(3))=-alpf*rtrho*a*bet2
    dUdVR(axb(3),axb(3))=-alpf*rtrho*a*bet3
  end subroutine composite_characteristic_eigenmatrix
end module characteristics
module spatial
  use common_module
  use characteristics
  implicit none
contains
  subroutine interpol(U0,xi,yj,zk,axis,UR,UL)
    implicit none
    real(DP),dimension(xlen,ylen,zlen,8),intent(IN)::U0
    integer,intent(IN) :: xi,yj,zk,axis
    real(DP), intent(OUT) :: UR(8),UL(8)
    real(DP),dimension(8,4)  :: V8s
    real(DP),dimension(8,8) :: evR,evL
    real(DP) :: VR(8),VL(8),dV8(8)
    integer s,is(4),js(4),ks(4),kdlt(3)
    if (udim(axis).eq.0) then
       UR=CT_preprocess(U0,xi,yj,zk)
       UL=UR
       return
    end if
    kdlt=0**abs(axis-(/1,2,3/))
    is=xi+(/-1,0,1,2/)*kdlt(1)
    js=yj+(/-1,0,1,2/)*kdlt(2)
    ks=zk+(/-1,0,1,2/)*kdlt(3)
    do s=1,4
       V8s(:,s)=cnsrvtv2prmtv(CT_preprocess(U0,is(s),js(s),ks(s)),gamma1i)
    end do
    if (powell_et_al_1999 .eq. 1) then
       call characteristic_eigenmatrix(V8s(:,2),evL,evR,axis)
       dV8=v_Mv(evR,kernel(v_Mv(evL,V8s(:,3)-V8s(:,2)),&
                           v_Mv(evL,V8s(:,2)-V8s(:,1))))
       do s=1,8; dV8(s)=dV8(s)*(1._DP-eq0(maxval(abs(V8s(s,:)-V8s(s,2)))))
       end do
       UL=prmtv2cnsrvtv(V8s(:,2)+dV8*.5_DP,gamma1i)

       call characteristic_eigenmatrix(V8s(:,3),evL,evR,axis)
       dV8=v_Mv(evR,kernel(v_Mv(evL,(V8s(:,4)-V8s(:,3))),&
                           v_Mv(evL,(V8s(:,3)-V8s(:,2)))))
       do s=1,8; dV8(s)=dV8(s)*(1._DP-eq0(maxval(abs(V8s(s,:)-V8s(s,3)))))
       end do
       UR=prmtv2cnsrvtv(V8s(:,3)-dV8*.5_DP,gamma1i)
    else
       dV8=kernel(V8s(:,3)-V8s(:,2),V8s(:,2)-V8s(:,1))
       UL=prmtv2cnsrvtv(V8s(:,2)+dV8*.5_DP,gamma1i)
       dV8=kernel(V8s(:,4)-V8s(:,3),V8s(:,3)-V8s(:,2))
       UR=prmtv2cnsrvtv(V8s(:,3)-dV8*.5_DP,gamma1i)       
    end if
  contains
    function kernel(dVR,dVL)
      implicit none
      real(DP), intent(IN),dimension(8) :: dVR,dVL
      real(DP), dimension(8) :: kernel
      integer s
      do s=1,8
         kernel(s)=.5_DP*(sign(1._DP,dVR(s))+sign(1._DP,dVL(s)))*&
              min(2._DP*abs(dVR(s)),.5_DP*abs(dVR(s)+dVL(s)),&
              2._DP*abs(dVL(s)))
      end do
    end function kernel
  end subroutine interpol
  function v_Mv(matrix,vector)
    implicit none
    real(DP),intent(IN) :: matrix(:,:),vector(:)
    integer :: cl
    real(DP) :: v_Mv(size(vector))
    do cl=1,size(vector)
       v_Mv(cl)=sum(matrix(cl,:)*vector)
    end do
  end function V_MV
  function M_MM(matrix1,matrix2)
    implicit none
    real(DP),intent(IN) :: matrix1(:,:),matrix2(:,:)
    real(DP) :: M_MM(int(sqrt(dble(size(matrix1)))),&
         int(sqrt(dble(size(matrix1)))))
    integer :: rw,cl,dim
    dim=int(sqrt(dble(size(matrix1))))
    do rw=1,dim; do cl=1,dim
       M_MM(cl,rw)=sum(matrix1(cl,:)*matrix2(:,rw))
    end do; end do
  end function M_MM
end module spatial