module steamtables
       use datafile_type
       implicit none
       
       type(datafile), private :: sub, sat, super;

       interface init
          module procedure init_steam
       end interface
       
       interface T_sat
          module procedure m_T_sat
       end interface
       interface P_sat
          module procedure m_P_sat
       end interface
       interface v_sat
          module procedure m_v_sat
       end interface
       interface u_sat
          module procedure m_u_sat
       end interface
       interface h_sat
          module procedure m_h_sat
       end interface
       interface s_sat
          module procedure m_s_sat
       end interface
       interface Cv_sat
          module procedure m_Cv_sat
       end interface
       interface Cp_sat
          module procedure m_Cp_sat
       end interface
       interface v_sup
          module procedure m_v_sup
       end interface
       interface u_sup
          module procedure m_u_sup
       end interface
       interface h_sup
          module procedure m_h_sup
       end interface
       interface s_sup
          module procedure m_s_sup
       end interface
       interface v_sub
          module procedure m_v_sub
       end interface
       interface u_sub
          module procedure m_u_sub
       end interface
       interface h_sub
          module procedure m_h_sub
       end interface
       interface s_sub
          module procedure m_s_sub
       end interface
       
contains
       
       real function interp_sup(T, P, col)
               real, intent(in) :: T, P
               integer, intent(in) :: col

               integer :: i, P_low_index, P_high_index, &
                    T_low_index, T_high_index
               real    ::  P_low, P_high, P_factor, T_low_low, T_low_high, T_high_low, &
                    T_high_high, T_high_factor, T_low_factor, lowval, highval

               if (.not. isgood(super)) stop "*** superheat tables not initialized ***"

               i = 1 
               do while (P>get(super,i,2) .and. i<numrows(super))
                  i= i + 1
               end do

               if (T<get(super,i,1) .or. i==1) stop "*** invalid pressure or temperature ***"

               P_high_index = i

               if (col > numcols(super)) stop "*** invalid column ***"
               
               do while(get(super,i-2,1)<get(super,i-1,1))
                  i = i-1
               end do

               P_low_index = i-1
               P_low = get(super, P_low_index, 2)

               P_factor = (P-P_low)/(P_high -P_low)

               i = P_low_index

               if(T>get(super,P_high_index-1,1)) stop "*** invalid temperature ***"

               do while(T>get(super,i,1)) 
                  i = i+1
               end do

               if(i<P_high_index) then
                  T_low_index = i-1
                  T_low_low = get(super,T_low_index,1)
                  T_low_high = get(super,T_low_index+1,1)
                  T_low_factor = (T-T_low_low)/(T_low_high-T_low_low)

                  lowval = get(super,T_low_index,col) + &
                       T_low_factor* (get(super,T_low_index+1,col) - &
                       get(super,T_low_index,col))
               else
                  lowval = interp_sat(2, P_low, 1.0, (col-2)/2)
               end if

               i = P_high_index

               do while(get(super,i+1,1)>get(super,i,1)) 
                  i = i+1
               end do

               if (get(super,i,1)<T) stop "*** invalid temperature ***"

               do while(get(super,i,1)>T .and. i>=P_high_index)
                  i = i -1
               end do

               if (i>=P_high_index) then
                  T_high_index = i-1
                  T_high_low = get(super,T_high_index,1)
                  T_high_high = get(super,T_high_index+1,1)
                  T_high_factor = (T-T_high_low)/(T_high_high-T_high_low)

                  highval = get(super,T_high_index,col) + &
                       T_high_factor*(get(super,T_high_index+1,col)- &
                       get(super,T_high_index,col))
               else
                  highval = interp_sat(2, P_high, 1.0, (col-2)/2)
               end if

               interp_sup = lowval + P_factor*(highval-lowval)
	end function interp_sup
	
	function interp_sat(col1, val1, x, col3)
		integer, intent(in) :: col1, col3
		real, intent(in) :: val1, x
		real :: interp_sat, f, g, factor
		integer :: i
		
		if(get(sat,numrows(sat),col1)<val1) stop "*** supercritical state in saturated interpolator ***"
		
		if(x<0.0 .or. x>1.0) stop "*** invalid quality ***"
		
		i = 1
		
		do while(get(sat,i,col1)<val1)
			i = i + 1
		end do
		factor = (val1-get(sat, i-1, col1))/(get(sat, i, col1)-get(sat, i-1, col1))
		
         f = get(sat,i-1,col3)+factor*(get(sat,i,col3)-get(sat,i-1,col3))
		
		if(col3>2) then
         	g = get(sat,i-1,col3+1)+factor*(get(sat,i,col3+1)-get(sat,i-1,col3+1))
            interp_sat = f + x * (g - f)
		else
			interp_sat = f
		end if
	end function
	
	function interp_sub(T, P, col)
		integer, intent(in) :: col
		real, intent(in) :: T, P
		real :: interp_sub
		
		integer :: i, P_low_index, P_high_index, T_low_index, T_high_index
		real :: P_high, P_low, P_factor, T_low_low, T_low_high, T_high_low, &
					T_high_high, T_high_factor, T_low_factor, lowval, highval
		
		if (.not. isgood(sub)) stop "*** subcooled tables not initialized ***"
         
		if (T_sat(P)<=T) stop "*** saturated or superheated request to subcooled interpolator ***"
		i = 1 
		do while (P>get(sub,i,2) .and. i<numrows(sub))
			i= i + 1
		end do
		if (P>get(sub,i,2) .or. i==1) stop "*** invalid pressure ***"
		
		P_high_index = i
		P_high = get(sub, P_high_index,2)
		
		do while(get(sub,i-2,1)<get(sub,i-1,1))
			i = i-1
		end do
		
		P_low_index = i-1
		P_low = get(sub, P_low_index, 2)
		
		P_factor = (P-P_low)/(P_high -P_low)
		
		i = P_low_index
		
		if(T<get(sub,P_low_index,1)) stop "*** invalid temperature ***"
		
		do while(T>get(sub,i,1)) 
			i = i+1
		end do 

		if(i<P_high_index) then
			T_low_index = i-1
			T_low_low = get(sub,T_low_index,1)
			T_low_high = get(sub,T_low_index+1,1)
			T_low_factor = (T-T_low_low)/(T_low_high-T_low_low)
		
			lowval = get(sub,T_low_index,col) + &
					T_low_factor* (get(sub,T_low_index+1,col) - &
					get(sub,T_low_index,col))
		else
			lowval = interp_sat(2, P_low, 1.0, (col-2)/2)
		end if
		
		i = P_high_index
		
		if (get(sub,P_high_index,1)>T) stop "*** invalid temperature ***"

         do while(get(sub,i+1,1)>get(sub,i,1) .and. T>get(sub,i,1)) 
			i = i+1
		end do
		
		if (T<=get(sub,i,1)) then
			T_high_index = i-1
			T_high_low = get(sub,T_high_index,1)
			T_high_high = get(sub,T_high_index+1,1)
			T_high_factor = (T-T_high_low)/(T_high_high-T_high_low)

			highval = get(sub,T_high_index,col) + &
					T_high_factor*(get(sub,T_high_index+1,col)- &
					get(sub,T_high_index,col))
		else
			highval = interp_sat(2, P_high, 1.0, (col-2)/2)
		end if
		
		interp_sub = lowval + P_factor*(highval-lowval)
	end function interp_sub
	
	subroutine init_steam(satfile, superfile, subfile)
		character, intent(in), optional :: satfile*(len(satfile)), superfile*(len(superfile)), subfile*(len(subfile))
		character*80 :: fname_sat, fname_super, fname_sub
		
		if(.not. present(satfile)) then
			fname_sat = "/ncsu/apwingo/lib/steam/sat_steam.txt"
		else
			fname_sat = satfile
		end if
		
		if(.not. present(superfile)) then
			fname_super = "/ncsu/apwingo/lib/steam/superheat.txt"
		else
			fname_super = superfile
		end if
		
		if(.not. present(subfile)) then
			fname_sub = "/ncsu/apwingo/lib/steam/subcooled.txt"
		else
			fname_sub = subfile
		end if
		
		call open(sat, fname_sat)
		call open(super, fname_super)
		call open(sub, fname_sub)
	end subroutine init_steam
		
	function m_T_sat(P)
		real, intent(in) :: P
		real :: m_T_sat
		
		m_T_sat = interp_sat(2, P, 0.0, 1)
	end function m_T_sat
	
	function m_P_sat(T)
		real, intent(in) :: T
		real :: m_P_sat
		
		m_P_sat = interp_sat(1, T, 0.0, 2)
	end function m_P_sat
	
	function m_v_sat(T, P, x)
		real, intent(in), optional :: T,P
		real, intent(in) :: x
		
		real :: m_v_sat, val
		integer :: col
		
		if(present(T).and.present(P)) stop "*** invalid arguments ***"
		
		if(.not.(present(T).or.present(P))) stop "*** invalid arguments ***"
		
		if(present(T)) then
			col = 1
			val = T
		else
			col = 2
			val = P
		end if
		
		m_v_sat = interp_sat(col, val, x, 3)
	end function m_v_sat
	
	function m_u_sat(T, P, x)
		real, intent(in), optional :: T,P
		real, intent(in) :: x
		
		real :: m_u_sat, val
		integer :: col
		
		if(present(T).and.present(P)) stop "*** invalid arguments ***"
		
		if(.not.(present(T).or.present(P))) stop "*** invalid arguments ***"
		
		if(present(T)) then
			col = 1
			val = T
		else
			col = 2
			val = P
		end if
		
		m_u_sat = interp_sat(col, val, x, 5)
	end function m_u_sat
	
	function m_h_sat(T, P, x)
		real, intent(in), optional :: T,P
		real, intent(in) :: x
		
		real :: m_h_sat, val
		integer :: col
		
		if(present(T).and.present(P)) stop "*** invalid arguments ***"
		
		if(.not.(present(T).or.present(P))) stop "*** invalid arguments ***"
		
		if(present(T)) then
			col = 1
			val = T
		else
			col = 2
			val = P
		end if
		
		m_h_sat = interp_sat(col, val, x,7)
	end function m_h_sat
 	
	function m_s_sat(T, P, x)
		real, intent(in), optional :: T,P
		real, intent(in), optional :: x
		
		real :: m_s_sat, val
		integer :: col
		
		if(.not. present(x)) stop "*** invalid arguments ***"
		
		if(present(T).and.present(P)) stop "*** invalid arguments ***"
		
		if(.not.(present(T).or.present(P))) stop "*** invalid arguments ***"
		
		if(present(T)) then
			col = 1
			val = T
		else
			col = 2
			val = P
		end if
		
		m_s_sat = interp_sat(col, val, x, 9)
	end function m_s_sat
	
	function m_Cv_sat(T, P, x)
		real, intent(in), optional :: T,P
		real, intent(in), optional :: x
		
		real :: m_Cv_sat, val
		integer :: col
		
		if(.not. present(x)) stop "*** invalid arguments ***"
		
		if(present(T).and.present(P)) stop "*** invalid arguments ***"
		
		if(.not.(present(T).or.present(P))) stop "*** invalid arguments ***"
		
		if(present(T)) then
			col = 1
			val = T
		else
			col = 2
			val = P
		end if
		
		m_Cv_sat = interp_sat(col, val, x, 11)
	end function m_Cv_sat
	
	function m_Cp_sat(T, P, x)
		real, intent(in), optional :: T,P
		real, intent(in), optional :: x
		
		real :: m_Cp_sat, val
		integer :: col
		
		if(.not. present(x)) stop "*** invalid arguments ***"
		
		if(present(T).and.present(P)) stop "*** invalid arguments ***"
		
		if(.not.(present(T).or.present(P))) stop "*** invalid arguments ***"
		
		if(present(T)) then
			col = 1
			val = T
		else
			col = 2
			val = P
		end if
		
		m_Cp_sat = interp_sat(col, val, x, 13)
	end function m_Cp_sat
	
	function m_v_sup(T, P)
		real, intent(in) :: T,P
		real m_v_sup
		
		m_v_sup = interp_sup(T, P, 3)
	end function m_v_sup

	function m_u_sup(T, P)
		real, intent(in) :: T,P
		real m_u_sup
		
		m_u_sup = interp_sup(T, P, 4)
	end function m_u_sup
	
	function m_h_sup(T, P)
		real, intent(in) :: T,P
		real m_h_sup
		
		m_h_sup = interp_sup(T, P, 5)
	end function m_h_sup
	
	function m_s_sup(T, P, h)
		real, intent(in),optional :: T,P,h
		real m_s_sup, T_, P_, h_
		
		if(.not. present(T)) then
			P_ = P
			h_ = h
			T_ = T_sat(P_)
			do while (h_sup(T_, P_)<h_)
				T_ = T_ + 0.05
			end do
		else if (.not. present(P)) then
			h_ = h
			T_ = T
			P_ = P_sat(T_)
			do while (h_sup(T_, P_)<h_)
				P_ = P_ + 0.05
			end do
		else
			T_ = T
			P_ = P
		end if
		
		m_s_sup = interp_sup(T_, P_, 6)
	end function m_s_sup
	
	function m_v_sub(T, P)
		real, intent(in) :: T,P
		real m_v_sub
		
		m_v_sub = interp_sub(T, P, 3)
	end function m_v_sub

	function m_u_sub(T, P)
		real, intent(in) :: T,P
		real m_u_sub
		
		m_u_sub = interp_sub(T, P, 4)
	end function m_u_sub
	
	function m_h_sub(T, P)
		real, intent(in) :: T,P
		real m_h_sub
		
		m_h_sub = interp_sub(T, P, 5)
	end function m_h_sub
	
	function m_s_sub(T, P, h)
		real, intent(in),optional :: T,P,h
		real m_s_sub, T_, P_, h_
		
		if(.not. present(T)) then
			P_ = P
			h_ = h
			T_ = T_sat(P_)
			do while (h_sub(T_, P_)<h_)
				T_ = T_ + 0.05
			end do
		else if (.not. present(P)) then
			h_ = h
			T_ = T
			P_ = P_sat(T_)
			do while (h_sub(T_, P_)<h_)
				P_ = P_ + 0.05
			end do
		else
			T_ = T
			P_ = P
		end if
		
		m_s_sub = interp_sub(T_, P_, 6)
	end function m_s_sub
end module steamtables