anchor.diff:     use IAPWS97, only : cowat
   
     implicit none
  -  
  +
     real*8 p,h, energyscale
     real*8 T,dw, hx
     integer ierr
  -  
  +
     real*8 tx0, tx, dt, txd
     real*8 hw, hxd, dwmol
     logical succ
     integer itr
  -  
  +
     tx0 = T ! get initial guess
  -  
  +
     if (tx0 <1D0) then ! obtain initial guess
  -    tx0 = 100D0 ! or use constant cp 
  -  endif 
  -
  -  itr= 0  
  +    tx0 = 100D0 ! or use constant cp
  +  endif
  +
  +  itr= 0
     tx=tx0
     dt=-1D-7
  -  
  -  do 
  +
  +  do
      succ= cowat(tx, p, dw, hw)
      hx= hw * energyscale
  -!   if(dabs((hx-h)/h)<= 1D-14) exit 
  -    if(dabs((hx-h))<= 1D-14) exit 
  +!   if(dabs((hx-h)/h)<= 1D-14) exit
  +    if(dabs((hx-h))<= 1D-14) exit
      txd = tx + dt
      succ = cowat(txd, p, dw, hw)
      hxd= hw * energyscale
  
positive.code:       stop
      endif
     enddo
     t=tx  
    End Subroutine steameos_PH
    
  !*************************************************************
  
  
   
  
  !*************************************************************
   
  !*************************************************************
  ! Subroutine water_steam_prop_PH
  !-------------------------------------------------------------  
    subroutine water_steam_prop_PH(p, h, T, Sw, &
      Den,Denw, Dens, &
      hw, hs, dDendh, dDendp,  dhwdh,dhsdh,&
      dTdh, dswdh, ierror, dhwdp, dhsdp,dTdp)
     use IAPWS97, only : cowat, supst, tsat
    implicit none      
    
    real*8, parameter::  fmwh2o = 18.01534d0
    real*8, parameter::  ps350c = 16.529D6, hw350c=1670.9D0
    real*8, parameter::  energyscale=1D-6 ! MJ/Kg
    real*8, parameter::  epi_p = 0.1D-0, epi_h = 0.1D0 * energyscale
   
  
    real*8, intent(inout)  :: p, h ! Pscal, KJ/Kg, C (initial guess)
    real*8, intent(out) :: T, Sw, hw,hs,Den,Denw,Dens,  dTdh, dDendp,dDendh, dhwdh,dhsdh, dswdh, dhwdp,dhsdp
    integer,intent(out) :: ierror  ! ierr = 10  critical point   
   
    real*8 Ts 
    real*8 dw, dwmol
    real*8 dg, dgmol, dgp,dgt,hg,hgp,hgt
    integer ierr
  
negative.code:   real*8 dg0, dg1, hg0,hg1, h1
    real*8 t0,t1, d0,d1
    real*8 dTdp
    logical succ
  
  
  ! determine phase condition
    succ= TSAT( p, Ts)
  ! if(h<1D-1) h=1D-1
   if(p> ps350c )then
     if(h<=hw350c)iphase=1
     if(Ts>370D0) iphase=2
    else
      succ= cowat(Ts,p, dw, hw)
      hw=hw * energyscale
      succ= supst(Ts, p, dg,hg)
      hg=hg* energyscale
  
      if(h >= hg)then
        iphase = 2
      else if(h <= hw)then
        iphase=1
      else
        iphase=3
      endif
    endif
  
  
    select case(iphase)
    case(1) ! water only
      T= Ts-10D0
  