con2primHydro Subroutine

  pure subroutine con2primHydro(lfac,xi,d,sqrs,tau)
    
    double precision, intent(out) :: xi,lfac
    double precision, intent(in)  :: d,sqrs,tau

    ! .. local ..
    integer          :: ni,niiter,nit,n2it,ni3
    double precision :: pcurrent,pnew
    double precision :: er,er1,ff,df,dp,v2
    double precision :: pmin,lfac2inv,pLabs,pRabs,pprev
    double precision :: xicurrent,h,dhdp
    double precision :: oldff1,oldff2,Nff
    double precision :: pleft,pright
    double precision :: rho
    !---------------------------------------------------------------------

    ! left and right brackets for p-range
    pmin=dsqrt(sqrs)/(one-dmaxvel)-tau-d
    pLabs=max(srhd_small_pressure,pmin)
    pRabs=1.0d99
    ! start value from input
    pcurrent=pLabs

    er1=one
    pprev=pcurrent

    ! Fudge Parameters
    oldff1=1.0d7  ! High number
    oldff2=1.0d9  ! High number bigger then oldff1
    n2it = 0
    nit  = 0

    LoopNR:  do ni=1,maxitnr
       nit = nit + 1
       !=== Relax NR iteration accuracy=======!
       if(nit>maxitnr/4)then
          ! mix pressure value for convergence
          pcurrent=half*(pcurrent+pprev)
          ! relax accuracy requirement
          er1=10.0d0*er1
          nit = nit - maxitnr/10
       endif
       !=======================================!

       niiter=ni
       xicurrent=tau+d+pcurrent

       v2=sqrs/xicurrent**2
       lfac2inv=one - v2
          lfac=one/dsqrt(lfac2inv)

       !== calculation done using the EOS ==!
       rho=d*dsqrt(lfac2inv)
       h = rho + pcurrent*phys_gamma/(phys_gamma-1.0d0)
       dhdp = phys_gamma/(phys_gamma-1.0d0) + d/dsqrt(lfac2inv)*sqrs/xicurrent**3
       !=======================================!
       ff=-xicurrent*lfac2inv + h
       df=- two*sqrs/xicurrent**2  + dhdp - lfac2inv

       if (ff*df==zero) then
          if (ff==zero) then
             exit ! zero found
          endif
       else
          pnew=pcurrent-ff/df
          if (ff*df>zero) then
             ! pressure iterate has decreased
             ! restrict to left
             pnew=max(pnew,pLabs)
          else  ! ff*df<0
             ! pressure iterate has increased
             ! restrict to right
             pnew=min(pnew,pRabs)
          endif
       endif

       !===============================================!
       dp=pcurrent-pnew
       er=two*dabs(dp)/(pnew+pcurrent)
       if(((er<tolernr*er1).or.(dabs(dp)<absaccnr))) exit LoopNR
       !===============================================!

       ! For very small values of pressure, NR algorithm is not efficient to
       ! find root, use Euler algorithm to find precise value of pressure
       if((dabs(oldff2-ff) < 1.0d-8 .or. niiter >= maxitnr-maxitnr/20).and.&
            ff * oldff1 < zero    .and.  dabs(ff)>absaccnr)then

          n2it=n2it+1
          if(n2it<=3) pcurrent=half*(pnew+pcurrent)
          if(n2it>3)then
             pright =pcurrent
             pleft=pprev
             pcurrent=half*(pleft+pright)
             Dicho:  do ni3=1,maxitnr
                !===================!
                xicurrent=tau+d+pcurrent
                v2=sqrs/xicurrent**2
                lfac2inv=one - v2

                   lfac=one/dsqrt(lfac2inv)

                !== calculation done using the EOS ==!
                rho=d*dsqrt(lfac2inv)
                h = rho + pnew*phys_gamma/(phys_gamma-1.0d0)
                Nff=-xicurrent*lfac2inv + h
                !=======================================!
                !==== Iterate ====!
                if(ff * Nff < zero)then
                   pleft=pcurrent
                else
                   pright=pcurrent
                endif

                pcurrent=half*(pleft+pright)
                !==================!

                !=== The iteration converged ===!
                if(2.0d0*dabs(pleft-pright)/(pleft+pright)< absaccnr &
                     .or. dabs(ff)<absaccnr)then
                   pnew=pcurrent
                   exit LoopNR
                endif
                !==============================!

                !==============================!

                !=== conserve the last value of Nff ===!
                ff=Nff
                !======================================!
             enddo    Dicho
          endif

       else
          !====== There is no problems, continue the NR iteration ======!
          pprev=pcurrent
          pcurrent=pnew
          !=============================================================!
       endif


       !=== keep the values of the 2 last ff ===!
       oldff2=oldff1
       oldff1=ff
       !========================================!
    enddo LoopNR

    !------------------------------!
    xi=tau+d+pcurrent
    v2=sqrs/xicurrent**2
    lfac2inv=one - v2
       lfac=one/dsqrt(lfac2inv)

  end subroutine con2primHydro