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:   real*8 dg0, dg1, hg0,hg1, h1
    real*8 t0,t1, d0,d1
    real*8 dTdp
    logical succ 
   
  
  ! determine phase condition    
    ierror=0
    succ= TSAT( p, Ts)
   
   if(p> ps350c )then 
     if(h<=hw350c)iphase=1
     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 
   
       
    select case(iphase)
    case(1) ! water only
      T= Ts-10D0
  
negative.code:     Den= sw*dw+(1D0-sw)*dg
      Denw=dw
      Dens=dg
      hs=hg
      delp = epi_p
      delh = epi_h
      call viss_noderiv1(Dens,T, viss)
      call viss_noderiv1(Denw,T, visw)
    !  print *, 'eos: ',p,h,hs,hw, sw
     end select
     if(h<0 .or.h> 3D3) print *,'h ph: ', p,t,h,hw,hg,iphase
  
  
  
   end subroutine water_steam_prop_PH_noderiv
  
  
  