anchor.diff:    real*8 dg0, dg1, hg0,hg1, h1
     real*8 t0,t1, d0,d1
     real*8 dTdp
  -  logical succ 
  - 
  -
  -! determine phase condition    
  +  logical succ
  +
  +
  +! determine phase condition
     ierror=0
     succ= TSAT( p, Ts)
  - 
  - if(p> ps350c )then 
  +
  + if(p> ps350c )then
      if(h<=hw350c)iphase=1
  -   if(Ts>370) iphase=2    
  +   if(Ts>370) 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 
  - 
  -     
  +    endif
  +  endif
  +
  +
     select case(iphase)
     case(1) ! water only
       T= Ts-10D0
  
positive.code:    dswdh=(sw1-sw0)/delh
  
    end select
    T=T0
                
   end subroutine water_steam_prop_PH
   
  !*************************************************************
   
  !*************************************************************
  ! Subroutine water_steam_prop_PH
  !-------------------------------------------------------------   
   subroutine water_steam_prop_PH_noderiv(p, h, T, Sw, &
      Den,Denw, Dens, hw, hs,visw,viss,ierror)
     use IAPWS97, only : cowat, supst, tsat, visc
    implicit none      
    
    real*8, parameter::  fmwh2o = 18.01534d0
    real*8, parameter::  epi_p = 1D0, epi_h=1D-3
    real*8, parameter::  ps350c = 16.529D6, hw350c=1670.9D0
    real*8, parameter::  energyscale=1D-6 ! MJ/Kg
   
  
    real*8, intent(in)  :: p, h ! Pscal, KJ/Kg, C (initial guess)
    real*8, intent(out) :: T, Sw, hw,hs,Den,Denw,Dens, viss,visw
    integer,intent(out) :: ierror  ! ierr = 10  critical point   
   
    real*8 Ts 
    real*8 dw  
    real*8 dg, hg
    integer ierr
    integer :: iphase = 0
  
negative.code:    succ=supst(t1, p+delp, dg,hg)
     hw=hw*energyscale; hg=hg*energyscale
     sw1 = 1D0/(1D0 -dw/dg*(hw-h)/(hg-h))
     d1=sw1*dw+(1D0-sw1)*dg
     dDendp=(d1-d0)/delp
  
  !+++++++++++++++++++++++++++
  
     dhwdp=(hw-hw0)/delp
     dhsdp=(hg-hg0)/delp
  
  