module temporal
  use boundary_module
  use spatial
  use source_term_module
  use ct_module
  use heat_conduction_module
  implicit none
contains
  subroutine ticktock(U0)
    implicit none
    real(DP),dimension(xlen,ylen,zlen,8),intent(INOUT) :: U0
    real(DP),dimension(xlen,ylen,zlen,8) :: dU
    real(DP) :: t0,dt
    t0=tax(mod(n-1,taxlen)+1); dt=CFL_ideal(U0)
    tax(mod(n,taxlen)+1)=t0+dt*cflf
    if (n .eq. 1) tax(taxlen+1)=tax(2)-tax(1)

    if (heat_conduction.eq.1) then
       call heat_conduction_RKL2(U0,dU,t0,dt*cflf*.5_DP)
       U0=bc(U0,dU,tax(mod(n,taxlen)+1))
!       call rk_tvd(U0,dU,dt,torder)
!       U0=bc(U0,dU,tax(mod(n,taxlen)+1))
       call heat_conduction_RKL2(U0,dU,t0+dt*cflf*.5_DP,dt*cflf*.5_DP)
       U0=bc(U0,dU,tax(mod(n,taxlen)+1))

    else
       call rk_tvd(U0,dU,dt,torder)
       U0=bc(U0,dU,tax(mod(n,taxlen)+1))
    end if
  end subroutine ticktock

  function CFL_ideal(U0)
    implicit none
    real(DP),dimension(xlen,ylen,zlen,8),intent(IN)::U0
    real(DP) :: CFL_ideal
    real(DP) :: V8(8),rhoi,cs2,ca2,cf
    real(DP) :: pc(3),hc(3),dcell_c(2),SLR(2)
    real(DP) :: pfR(3),pfL(3),hfR(3),hfL(3),sfR,sfL
    real(DP) :: eps=p_yocto*p_yocto
    real(DP) :: dt_a(udim(1)+udim(2)+udim(3))
    integer :: ir(2),jr(2),kr(2),szi,szj,szk,s(3),kdlt(3),cnt,cnt1
    integer :: axis
    CFL_ideal=p_giga
    ir=xr+(/-1,1/)*udim(1); jr=yr+(/-1,1/)*udim(2)
    kr=zr+(/-1,1/)*udim(3)
    szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1
    !$omp parallel do private(s,cnt1,pc,hc,axis,kdlt,pfR,hfR,pfL,hfL,&
    !$omp dcell_c,V8,rhoi,cs2,ca2,cf,SLR,dt_a),reduction(min:CFL_ideal)
    do cnt=1,szi*szj*szk
       call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),s)
       call cell_center(s(3),s(2),s(1),pc,hc)
       V8=cnsrvtv2prmtv(CT_preprocess(U0,s(3),s(2),s(1)),gamma1i)
       rhoi=1._DP/V8(Sax(1)); cs2=gamma*V8(Sax(2))*rhoi
       ca2=sum(V8(Bax)*V8(Bax))*rhoi
       do cnt1=1,sum(udim)
          axis=udir(cnt1); kdlt=0**abs(axis-(/3,2,1/))
          call cell_face(s(3),s(2),s(1),axis, 1,pfR,hfR,sfR)
          call cell_face(s(3),s(2),s(1),axis,-1,pfL,hfL,sfL)
          dcell_c(1)=hfL(axis)*abs(pc(axis)-pfL(axis))*2._DP
          dcell_c(2)=hfR(axis)*abs(pfR(axis)-pc(axis))*2._DP
          cf=abs((cs2+ca2)**2-4._DP*cs2*V8(Bax(axis))*V8(Bax(axis))*rhoi)
          cf=sqrt(.5_DP*(cs2+ca2+sqrt(cf)))
          SLR=(/-min(-eps,V8(Vax(axis))-cf),max(eps,V8(Vax(axis))+cf)/)

          dt_a(cnt1)=minval(dcell_c/SLR)
       end do
       CFL_ideal=min(CFL_ideal,1._DP/sum(1._DP/dt_a))
    end do
    !$omp end parallel do
  contains
    subroutine loopcounter(count,sizes,org,subs)
      implicit none
      integer,intent(IN) :: count,sizes(:),org(:)
      integer,intent(OUT) :: subs(:)
      integer :: dim,s,quo(size(subs))
      dim=size(sizes)
      quo=count
      do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1
      end do
      do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s)
      end do
    end subroutine loopcounter
  end function CFL_ideal
  function Euler(U0)
    !$ use omp_lib
    implicit none
    real(DP),dimension(xlen,ylen,zlen,8),intent(IN)::U0
    real(DP),dimension(xlen,ylen,zlen,8)::Euler
    real(DP),dimension(xlen,ylen,zlen,8,3)::Fxyz
    real(DP),dimension(xlen,ylen,zlen,3,2)::SLRs
    real(DP),dimension(8) :: F8,UR,UL,FR,FL,dF8
    real(DP) :: SLR(2)
    real(DP) :: dcell_c,dcell_f,dcell_e2,dcell_e3,sfR,sfL,dvc
    real(DP),dimension(3) :: peR,peL,pfR,pfL,heR2,heL2,heR3,heL3,hfR,hfL,hc,tmp
    integer axis,axis2,axis3,kdlt(3)
    integer :: ir(2),jr(2),kr(2)
    integer :: cnt,cnt1,szi,szj,szk,s(3),sb(3)
    Fxyz=0._DP; SLRs=0._DP; Euler=0._DP
    !$omp parallel
    do axis=1,3
       kdlt=0**abs(axis-(/1,2,3/))
       ir=xr+(/-1-kdlt(1),1/)*udim(1)
       jr=yr+(/-1-kdlt(2),1/)*udim(2)
       kr=zr+(/-1-kdlt(3),1/)*udim(3)
       szi=ir(2)-ir(1)+1; szj=jr(2)-jr(1)+1; szk=kr(2)-kr(1)+1
       !$omp barrier
       !$omp do private(s,UR,UL,F8,SLR)
       do cnt=1,szi*szj*szk
          call loopcounter(cnt,(/szk,szj,szi/),(/kr(1),jr(1),ir(1)/),s)
          call interpol(U0,s(3),s(2),s(1),axis,UR,UL)
          call hlld_flux(UR,UL,F8,axis,SLR)
          Fxyz(s(3),s(2),s(1),:,axis)=F8
          SLRs(s(3),s(2),s(1),axis,:)=SLR
       end do
       !$omp end do
    end do
    !$omp end parallel
    call CT(U0,Fxyz,SLRs)
    szi=xr(2)-xr(1)+1; szj=yr(2)-yr(1)+1; szk=zr(2)-zr(1)+1
    !$omp parallel do shared(Euler,Fxyz),default(firstprivate)
    do cnt=1,szi*szj*szk
       call loopcounter(cnt,(/szk,szj,szi/),(/zr(1),yr(1),xr(1)/),s)
       hc=scale_factors(s(3),s(2),s(1))
       dF8=0._DP
       do cnt1=1,sum(udim)
          axis=udir(cnt1); kdlt=0**abs(axis-(/1,2,3/))
          axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1
          sb=s-(/kdlt(3),kdlt(2),kdlt(1)/)
          call cell_face(s(3),s(2),s(1),axis, 1,pfR,hfR,sfR)
          call cell_face(s(3),s(2),s(1),axis,-1,pfL,hfL,sfL)
          
          dcell_f=pfR(axis)-pfL(axis)
          call cell_edge(s(3),s(2),s(1),axis,axis2,peR,heR2)
          call cell_edge(sb(3),sb(2),sb(1),axis,axis2,peL,heL2)
          call cell_edge(s(3),s(2),s(1),axis,axis3,peR,heR3)
          call cell_edge(sb(3),sb(2),sb(1),axis,axis3,peL,heL3)
          FR=Fxyz(s(3),s(2),s(1),:,axis)
          FL=Fxyz(sb(3),sb(2),sb(1),:,axis)
          FR(Sax)=FR(Sax)*sfR
          FR(Vax)=FR(Vax)*sfR*hfR
          FR(Bax(axis2))=FR(Bax(axis2))*heR2(axis3)
          FR(Bax(axis3))=FR(Bax(axis3))*heR3(axis2)

          FL(Sax)=FL(Sax)*sfL
          FL(Vax)=FL(Vax)*sfL*hfL
          FL(Bax(axis2))=FL(Bax(axis2))*heL2(axis3)
          FL(Bax(axis3))=FL(Bax(axis3))*heL3(axis2)


          dF8(Sax)=dF8(Sax)+(FR(Sax)-FL(Sax))/dcell_f
          dF8(Vax)=dF8(Vax)+(FR(Vax)-FL(Vax))/dcell_f
          dF8(Bax(axis2))=dF8(Bax(axis2))&
               +(FR(Bax(axis2))-FL(Bax(axis2)))/dcell_f
          dF8(Bax(axis3))=dF8(Bax(axis3))&
               +(FR(Bax(axis3))-FL(Bax(axis3)))/dcell_f
       end do
       dvc=hc(1)*hc(2)*hc(3)
       dF8(Sax)=dF8(Sax)/dvc
       dF8(Vax)=dF8(Vax)/hc/dvc
       do axis=1,3
          call cell_face(s(3),s(2),s(1),axis, 1,pfR,hfR,sfR)
          dF8(Bax(axis))=dF8(Bax(axis))/sfR
       end do
       Euler(s(3),s(2),s(1),:)=dF8
    end do
    !$omp end parallel do
    call source_term(U0,Fxyz,Euler)
  end function Euler

  subroutine rk_tvd_3rd(U0,dU,dt)
    !$use omp_lib
    implicit none
    real(DP),intent(IN) :: dt
    real(DP),dimension(xlen,ylen,zlen,8),intent(INOUT) :: U0
    real(DP),dimension(xlen,ylen,zlen,8) :: rk_U,rk_F,dU
    integer :: szi,szj,szk,sza,cnt,s(4)
    rkstep=1; rk_F=Euler(U0)
    szi=xlen; szj=ylen; szk=zlen; sza=8
    !$omp parallel do private(s),default(shared)
    do cnt=1,szi*szj*szk*sza
       call loopcounter(cnt,(/sza,szk,szj,szi/),(/1,1,1,1/),s)
       dU(s(4),s(3),s(2),s(1))=-rk_F(s(4),s(3),s(2),s(1))*dt*cflf
    end do
    !$omp end parallel do
    rk_U=bc(U0,dU,tax(mod(n-1,taxlen)+1)+dt*cflf)
    rkstep=2; rk_F=Euler(rk_U)
    !$omp parallel do private(s),default(shared)
    do cnt=1,szi*szj*szk*sza
       call loopcounter(cnt,(/sza,szk,szj,szi/),(/1,1,1,1/),s)
       dU(s(4),s(3),s(2),s(1))=.25_DP*(rk_U(s(4),s(3),s(2),s(1))&
            -U0(s(4),s(3),s(2),s(1))-rk_F(s(4),s(3),s(2),s(1))*dt*cflf)
    end do
    !$omp end parallel do
    rk_U=bc(U0,dU,tax(mod(n-1,taxlen)+1)+dt*cflf*.5_DP)

    rkstep=3; rk_F=Euler(rk_U)
    !$omp parallel do private(s),default(shared)
    do cnt=1,szi*szj*szk*sza
       call loopcounter(cnt,(/sza,szk,szj,szi/),(/1,1,1,1/),s)
       dU(s(4),s(3),s(2),s(1))&
            =(rk_U(s(4),s(3),s(2),s(1))-U0(s(4),s(3),s(2),s(1))&
            -rk_F(s(4),s(3),s(2),s(1))*dt*cflf)*2._DP/3._DP
    end do
    !$omp end parallel do
  contains
    subroutine loopcounter(count,sizes,org,subs)
      implicit none
      integer,intent(IN) :: count,sizes(:),org(:)
      integer,intent(OUT) :: subs(:)
      integer :: dim,s,quo(size(subs))
      dim=size(sizes)
      quo=count
      do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1
      end do
      do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s)
      end do
    end subroutine loopcounter
  end subroutine rk_tvd_3rd
  subroutine rk_tvd(U0,dU,dt,norder)
    implicit none
    real(DP),intent(IN) :: dt
    real(DP),dimension(xlen,ylen,zlen,8),intent(INOUT) :: U0
    real(DP),dimension(xlen,ylen,zlen,8) :: F0,dU
    real(DP),allocatable,dimension(:,:,:,:,:) :: rk_U,rk_F
    integer,intent(IN)::norder
    integer rkstep1,rk_c(norder)
    integer :: szi,szj,szk,sza,cnt,s(4)
    real(DP),allocatable,dimension(:)::rk_a,rk_b
    if (norder.eq.3) then
       call rk_tvd_3rd(U0,dU,dt)
       return
    end if
    allocate(rk_a(norder*(norder+1)/2))
    allocate(rk_b(norder*(norder+1)/2))
    allocate(rk_U(xlen,ylen,zlen,8,norder))
    allocate(rk_F(xlen,ylen,zlen,8,norder))
    select case (norder)
    case (1)
       rk_a=1._DP-1._DP
       rk_b=1._DP; rk_c=1._DP
    case (2)
       rk_a=(/1._DP,.5_DP,.5_DP/)-(/1._DP,1._DP,0._DP/)
       rk_b=(/1._DP,.0_DP,.5_DP/)
       rk_c=(/1._DP,1._DP/)
    case (3)
       rk_a=(/1._DP,.75_DP,.25_DP,1._DP/3._DP,0._DP,2._DP/3._DP/)-&
            (/1._DP, 1._DP, 0._DP,      1._DP,0._DP,0._DP/)
       rk_b=(/1._DP,0._DP,.25_DP,0._DP,0._DP,2._DP/3._DP/)
       rk_c=(/1._DP,.5_DP,1._DP/)
    case (4)
       rk_a=(/1._DP,.5_DP, .5_DP,1._DP/9._DP,2._DP/9._DP,2._DP/3._DP,&
            0._DP, 1._DP/3._DP,1._DP/3._DP,1._DP/3._DP/)-&
            (/1._DP,1._DP,0._DP,       1._DP,0._DP,0._DP,&
            1._DP,0._DP,0._DP,0._DP/)
       rk_b=(/.5_DP,-.25_DP, .5_DP,-1._DP/9._Dp, -1._DP/3._DP, 1._DP,&
            0._DP, 1._DP/6._DP, 0._DP, 1._DP/6._DP/)
       rk_c=(/1._DP,.5_DP,.5_DP,1._DP/)
    end select

    szi=xlen; szj=ylen; szk=zlen; sza=8
    rk_U(:,:,:,:,1)=U0
    if(norder.ge.2)then
       do rkstep=1,norder-1
          dU=0._DP
          rk_F(:,:,:,:,rkstep)=Euler(rk_U(:,:,:,:,rkstep))
          do rkstep1=1,rkstep
             !$omp parallel do private(s),default(shared)
             do cnt=1,szi*szj*szk*sza
                call loopcounter(cnt,(/sza,szk,szj,szi/),(/1,1,1,1/),s)
                dU(s(4),s(3),s(2),s(1))=dU(s(4),s(3),s(2),s(1))&
                     +rk_a(rkstep*(rkstep-1)/2+rkstep1)&
                     *rk_U(s(4),s(3),s(2),s(1),rkstep1)&
                     -rk_b(rkstep*(rkstep-1)/2+rkstep1)&
                     *rk_F(s(4),s(3),s(2),s(1),rkstep1)*dt*cflf
             end do
             !$omp end parallel do
          end do
          rk_U(:,:,:,:,rkstep+1)=&
               bc(U0,dU,tax(mod(n-1,taxlen)+1)+dt*cflf*rk_c(rkstep))
       end do
    end if
    rkstep=norder
    rk_F(:,:,:,:,rkstep)=Euler(rk_U(:,:,:,:,rkstep))
    dU=0._DP
    do rkstep1=1,rkstep
       !$omp parallel do private(s),default(shared)
       do cnt=1,szi*szj*szk*sza
          call loopcounter(cnt,(/sza,szk,szj,szi/),(/1,1,1,1/),s)
          dU(s(4),s(3),s(2),s(1))=dU(s(4),s(3),s(2),s(1))&
               +rk_a(rkstep*(rkstep-1)/2+rkstep1)&
               *rk_U(s(4),s(3),s(2),s(1),rkstep1)&
               -rk_b(rkstep*(rkstep-1)/2+rkstep1)&
               *rk_F(s(4),s(3),s(2),s(1),rkstep1)*dt*cflf
       end do
       !$omp end parallel do
    end do
  contains
    subroutine loopcounter(count,sizes,org,subs)
      implicit none
      integer,intent(IN) :: count,sizes(:),org(:)
      integer,intent(OUT) :: subs(:)
      integer :: dim,s,quo(size(subs))
      dim=size(sizes)
      quo=count
      do s=2,dim; quo(s)=(quo(s-1)-1)/sizes(s-1)+1
      end do
      do s=1,dim; subs(s)=mod(quo(s)-1,sizes(s))+org(s)
      end do
    end subroutine loopcounter
  end subroutine rk_tvd

  subroutine hlld_flux(UR,UL,F,axis,SLR)
    implicit none
    real(DP), intent(IN)  :: UR(8),UL(8)
    integer , intent(IN)  :: axis
    real(DP), intent(OUT) :: F(8)
    real(DP), intent(OUT) :: SLR(2)
    real(DP), parameter :: eps=1.D-40
    real(DP)              :: BR(3),BL(3),vR(3),vL(3),pTR,pTL,Bp
    real(DP)              :: BR1(3),BL1(3),vR1(3),vL1(3),pT1
    real(DP)              :: BR2(3),BL2(3),vR2(3),vL2(3),pT2
    real(DP)              :: SR,SL,SM,CsR,CsL,CaR,CaL,CfR,CfL
    real(DP)              :: SR1,SL1,SR2,SL2
    real(DP)              :: factor,denomi,rrhoR,rrhoL,sgnBp
    real(DP)              :: UR1(8),UL1(8),UR2(8),UL2(8)
    real(DP)              :: BnR2,BnL2,vnR2,vnL2,pR,pL
    real(DP)              :: SRmvR,SLmvL,SMmvR,SMmvL,SRmSM,SLmSM
    real(DP)              :: case1,case2
    integer :: axis2,axis3,axv(3),axb(3),nthrd=1
    F=0._DP
    axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1
    axv=(/Vax(axis),Vax(axis2),Vax(axis3)/)
    axb=(/Bax(axis),Bax(axis2),Bax(axis3)/)

    Bp=(UL(axb(1))+UR(axb(1)))*.5_DP
    vR=(/UR(axv(1)),UR(axv(2)),UR(axv(3))/)/UR(Sax(1))
    BR=(/Bp,UR(axb(2)),UR(axb(3))/)
    vL=(/UL(axv(1)),UL(axv(2)),UL(axv(3))/)/UL(Sax(1))
    BL=(/Bp,UL(axb(2)),UL(axb(3))/)
    vnR2=sum(vR*vR); vnL2=sum(vL*vL)
    BnR2=sum(BR*BR); BnL2=sum(BL*BL)

    pR=(UR(Sax(2))-.5_DP*(UR(Sax(1))*vnR2+BnR2))/gamma1i
    pL=(UL(Sax(2))-.5_DP*(UL(Sax(1))*vnL2+BnL2))/gamma1i
    pR=abs(pR); pL=abs(pL)
    CsR=sqrt(gamma*pR/UR(Sax(1))); CsL=sqrt(gamma*pL/UL(Sax(1)))
    CaR=sqrt(BnR2/UR(Sax(1))); CaL=sqrt(BnL2/UL(Sax(1)))
    CfR=CsR*CsR+CaR*CaR; CfL=CsL*CsL+CaL*CaL
    CfR=sqrt(.5_DP*(CfR+sqrt(CfR*CfR-4._DP*(CsR*Bp)*(CsR*Bp)/UR(Sax(1)))))
    CfL=sqrt(.5_DP*(CfL+sqrt(CfL*CfL-4._DP*(CsL*Bp)*(CsL*Bp)/UL(Sax(1)))))
    SR=max(0._DP,max(vL(1),vR(1))+max(CfL,CfR))
    SL=min(0._DP,min(vL(1),vR(1))-max(CfL,CfR))
!    SR=max(vR(1)+CfR,vL(1)+CfL); SL=min(vR(1)-CfR,vL(1)-CfL)
    SLR=(/SL,SR/)
    SRmvR=SR-vR(1); SLmvL=SL-vL(1)
    pTR=pR+BnR2*.5_DP; pTL=pL+BnL2*.5_DP
    case1=ge0(SL); case2=le0(SR)
    if (case1+case2 .eq. 1._DP) then
       F=U2flux(UL,Bp,pTL,axis)*case1+U2flux(UR,Bp,pTR,axis)*case2
       return
    end if

    SM=(SRmvR*UR(Sax(1))*vR(1)-SLmvL*UL(Sax(1))*vL(1)-pTR+pTL)&
         /(SRmvR*UR(Sax(1))-SLmvL*UL(Sax(1)))
    SMmvR=SM-vR(1); SMmvL=SM-vL(1)
    SRmSM=SR-SM; SLmSM=SL-SM

    UR1(Sax(1))=UR(Sax(1))*SRmvR/SRmSM; UL1(Sax(1))=UL(Sax(1))*SLmvL/SLmSM
    SR1=SM+abs(BR(1))/sqrt(UR1(Sax(1))); SL1=SM-abs(BL(1))/sqrt(UL1(Sax(1)))

    vR1(1)=SM; BR1(1)=BR(1)
    denomi=UR(Sax(1))*SRmvR*SRmSM-BR(1)*BR(1)
    case1=gt0(abs(denomi)-eps)
    denomi=case1/(denomi+(1._DP-case1))
    ! NOTE +
    ! if (abs(denomi) .gt. eps) then
    !    case1=1._DP i.e. denomi=1/denomi else case1=0._DP i.e. denomi=0._DP
    !-
    vR1(2:3)=vR(2:3)-BR(2:3)*BR(1)*SMmvR*denomi
    BR1(2:3)=BR(2:3)*((UR(Sax(1))*SRmvR*SRmvR-BR(1)*BR(1))*denomi+1._DP-case1)
    UR1(axv(1))=vR1(1); UR1(axv(2))=vR1(2); UR1(axv(3))=vR1(3)
    UR1(Vax)=UR1(Vax)*UR1(Sax(1))
    UR1(axb(1))=BR1(1); UR1(axb(2))=BR1(2); UR1(axb(3))=BR1(3)

    vL1(1)=SM; BL1(1)=BL(1)
    denomi=UL(Sax(1))*SLmvL*SLmSM-BL(1)*BL(1)
    case1=gt0(abs(denomi)-eps)
    denomi=case1/(denomi+(1._DP-case1))
    ! NOTE +
    ! if (abs(denomi) .gt. eps) then
    !    case1=1._DP i.e. denomi=1/denomi else case1=0._DP i.e. denomi=0._DP
    !-
    vL1(2:3)=vL(2:3)-BL(2:3)*BL(1)*SMmvL*denomi
    BL1(2:3)=BL(2:3)*((UL(Sax(1))*SLmvL*SLmvL-BL(1)*BL(1))*denomi+1._DP-case1)
    UL1(axv(1))=vL1(1); UL1(axv(2))=vL1(2); UL1(axv(3))=vL1(3)
    UL1(Vax)=UL1(Vax)*UL1(Sax(1))
    UL1(axb(1))=BL1(1); UL1(axb(2))=BL1(2); UL1(axb(3))=BL1(3)

    pT1=(UR(Sax(1))*pTL*SRmvR-UL(Sax(1))*pTR*SLmvL&
         +UR(Sax(1))*UL(Sax(1))*SRmvR*SLmvL*(vR(1)-vL(1)))&
         /(UR(Sax(1))*SRmvR-UL(Sax(1))*SLmvL)
    UR1(Sax(2))=(SRmvR*UR(Sax(2))-pTR*vR(1)+pT1*SM&
         +BR(1)*(sum(vR*BR)-sum(vR1*BR1)))/SRmSM
    UL1(Sax(2))=(SLmvL*UL(Sax(2))-pTL*vL(1)+pT1*SM&
         +BL(1)*(sum(vL*BL)-sum(vL1*BL1)))/SLmSM

    case1=lt0(SL)*gt0(SL1); case2=gt0(SR)*lt0(SR1)
    if (case1+case2 .eq. 1._DP) then
       F=U2flux(UL1,Bp,pT1,axis)*case1+U2flux(UR1,Bp,pT1,axis)*case2; return
    else
       if (abs(case1)+abs(case2) .ne. 0._DP) then
          print*,'% Case Error @ hlld_module.f90'
          print*,n,rkstep,axis
          print'("UR",8e15.7)',UR
          print'("UL",8e15.7)',UL
          stop
       end if
    end if
    sgnBp=sign(1._DP,Bp)

    UR2(Sax(1))=UR1(Sax(1)) ; UL2(Sax(1))=UL1(Sax(1))
    rrhoR=sqrt(UR2(Sax(1))) ; rrhoL=sqrt(UL2(Sax(1)))
    factor=1._DP/(rrhoR+rrhoL)
    vR2(1)=SM
    vR2(2:3)=(rrhoR*vR1(2:3)+rrhoL*vL1(2:3)+(BR1(2:3)-BL1(2:3))*sgnBp)*factor
    vL2=vR2
    UR2(axv(1))=vR2(1); UR2(axv(2))=vR2(2); UR2(axv(3))=vR2(3)
    UR2(Vax)=UR2(Vax)*UR2(Sax(1))
    UL2(axv(1))=vL2(1); UL2(axv(2))=vL2(2); UL2(axv(3))=vL2(3)
    UL2(Vax)=UL2(Vax)*UL2(Sax(1))

    BR2(1)=BR1(1)
    BR2(2:3)=(rrhoL*BR1(2:3)+rrhoR*BL1(2:3)&
         +rrhoR*rrhoL*(vR1(2:3)-vL1(2:3))*sgnBp)*factor
    BL2=BR2

    UR2(axb(1))=BR2(1); UR2(axb(2))=BR2(2); UR2(axb(3))=BR2(3)
    UL2(axb(1))=BL2(1); UL2(axb(2))=BL2(2); UL2(axb(3))=BL2(3)

    pT2=pT1
    UR2(Sax(2))=UR1(Sax(2))+rrhoR*(sum(vR1*BR1)-sum(vR2*BR2))*sgnBp
    UL2(Sax(2))=UL1(Sax(2))-rrhoL*(sum(vL1*BL1)-sum(vL2*BL2))*sgnBp
    case1=le0(SL1)*gt0(SM); case2=ge0(SR1)*le0(SM)
    if (case1+case2 .eq. 1) then
       F=U2flux(UL2,Bp,pT2,axis)*case1+U2flux(UR2,Bp,pT2,axis)*case2
       return
    else
       if (HLLDerror .eq. 0) then
          print'("  Axis...",i5,"  Time...",i5)',axis,n
          print'("UR   ",8e12.3)',UR
          print'("UL   ",8e12.3)',UL
          print'("pR,CfR,SR,SR1,SM   ",5e12.3)',pR,CfR,SR,SR1,SM
       end if
       if (HLLDerror .le. 100) then
          write (2,'   ("UR   ",8e12.3)') UR
          write (2,'   ("UL   ",8e12.3)') UL
          write (2,'   ("pR,CfR,SR,SR1,SM   ",5e12.3)') pR,CfR,SR,SR1,SM
          HLLDerror=HLLDerror+1
       endif
    end if
  contains
    function U2flux(U8,Bp,pT,axis)
      implicit none
      real(DP),intent(IN) :: U8(8),Bp,pT
      real(DP)  :: U2flux(8)
      integer :: axis,axis2,axis3,axv(3),axb(3)
      axis2=mod(axis,3)+1 ; axis3=mod(axis2,3)+1
      axv=(/Vax(axis),Vax(axis2),Vax(axis3)/)
      axb=(/Bax(axis),Bax(axis2),Bax(axis3)/)

      U2flux(Sax(1))=U8(axv(1))

      U2flux(axv(1))=U8(axv(1))*U8(axv(1))/U8(Sax(1))+pT-Bp*Bp
      U2flux(axv(2))=U8(axv(1))*U8(axv(2))/U8(Sax(1))-Bp*U8(axb(2))
      U2flux(axv(3))=U8(axv(1))*U8(axv(3))/U8(Sax(1))-Bp*U8(axb(3))

      U2flux(axb(1))=0._DP
      U2flux(axb(2))=(U8(axv(1))*U8(axb(2))-U8(axv(2))*Bp)/U8(Sax(1))
      U2flux(axb(3))=(U8(axv(1))*U8(axb(3))-U8(axv(3))*Bp)/U8(Sax(1))

      U2flux(Sax(2))=U8(axv(1))*Bp+U8(axv(2))*U8(axb(2))+U8(axv(3))*U8(axb(3))
      U2flux(Sax(2))=U2flux(Sax(2))*Bp
      U2flux(Sax(2))=((U8(Sax(2))+pT)*U8(axv(1))-U2flux(Sax(2)))/U8(Sax(1))
    end function U2flux
  end subroutine hlld_flux

end module temporal