!$$$$$$ List of SUBROUTINES:
!$$$$$$
!$$$$$$ SUBROUTINE HotStreak(fluid,LCO,CCO,DFN,HSF,m,n,A,t,f,Q,P,DELTAP_F,F_coef,T_in,rho_in,Q_H2O,g,g_c,epsilonD_e, &
!$$$$$$ V_AS,A_fi,e,e_min,e_e,DELTAD,U,w,U25,s,DELTAz,z,phi,D,k_lhf,i_lhf,j_lhf,k_mot,i_mot,k_mfr, &
!$$$$$$ i_mfr,MIN_DELTA,P_lhf,w_lhf,h_lhf,Ts_lhf,T_lhf,q_lhf,T_mot,w_mot,w_mfr,T_mfr):
!$$$$$$ returns the minimum excess superheat/critical heat flux, as well as quantities at limiting locations
!$$$$$$  
!$$$$$$ SUBROUTINE HotSpotHeatFlux(fluid,k,m,n,A,f,Q,P,U,U_bar,U25,z,D_min,T_hs,mu,k_cond,Re_D,Pr,phi,T_Shs,q_flux_hs,h_hs):
!$$$$$$ returns the hot spot head fluxes, surface temperatures, and heat transfer coefficients
!$$$$$$ 
!$$$$$$ SUBROUTINE Burnout(fluid,LCO,m,n,g,g_c,T_hs,T_sat,q_flux_hs,P,D_min,Re_D,Pr,k_cond,DELTA):
!$$$$$$ returns the differences between surface heat fluxes and burnout heat fluxes
!$$$$$$ 
!$$$$$$ SUBROUTINE IncipientBoiling(m,n,U23,T_Shs,T_sat,q_flux_hs,P,DELTA):
!$$$$$$ returns the differences between surface temperatures and incipient boiling temperatures
!$$$$$$ 
!$$$$$$ SUBROUTINE FlowExcursion(fluid,m,n,w_hs,T_hs,T_sat,q_flux_hs,P,D_min,Re_D,Pr,k_cond,DELTA):
!$$$$$$ returns the differences between surface heat fluxes and flow excursion heat fluxes
!$$$$$$ 
!$$$$$$ SUBROUTINE CriticalHeatFlux(fluid,m,n,g_c,w_hs,T_hs,T_sat,q_flux_hs,P,D_min,DELTA):
!$$$$$$ returns the differences between surface heat fluxes and critical heat fluxes
!$$$$$$ 
!$$$$$$ SUBROUTINE NetVaporGen(m,n,w_hs,T_hs,T_sat,q_flux_hs,rho,D_min,DELTA):
!$$$$$$ returns the differences between surface heat fluxes and point of net vapor generation heat fluxes
!$$$$$$ 

!$$$$$$ 
SUBROUTINE HotStreak(fluid,LCO,CCO,DFN,HSF,m,n,A,t,f,Q,P,DELTAP_F,F_coef,T_in,rho_in,Q_H2O,g,g_c,epsilonD_e, &
V_AS,A_fi,e,e_min,e_e,DELTAD,U,w,U25,s,DELTAz,z,phi,D,k_lhf,i_lhf,j_lhf,k_mot,i_mot,k_mfr, &
i_mfr,MIN_DELTA,P_lhf,w_lhf,h_lhf,Ts_lhf,T_lhf,q_lhf,T_mot,w_mot,w_mfr,T_mfr)
!$$$$$$  
IMPLICIT NONE
CHARACTER*4, INTENT(IN) :: fluid
INTEGER, INTENT(IN) :: LCO,CCO,DFN,HSF,m,n
DOUBLE PRECISION, INTENT(IN) :: A,t,f,Q,P,DELTAP_F,F_coef,T_in,rho_in,Q_H2O,g,g_c,epsilonD_e,V_AS,A_fi
DOUBLE PRECISION, INTENT(IN), DIMENSION(2) :: e,e_min,e_e,DELTAD
DOUBLE PRECISION, INTENT(IN), DIMENSION(25) :: U
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,2) :: w,U25,s
DOUBLE PRECISION, INTENT(IN), DIMENSION(n,2) :: DELTAz,z
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: phi,D
INTEGER, INTENT(OUT) :: k_lhf,i_lhf,j_lhf,k_mot,i_mot,k_mfr,i_mfr
DOUBLE PRECISION, INTENT(OUT) :: MIN_DELTA,P_lhf,w_lhf,h_lhf,Ts_lhf,T_lhf,q_lhf,T_mot,w_mot,w_mfr,T_mfr
INTEGER :: i,j,k
INTEGER, DIMENSION(2) :: i_limit,j_limit,i_mt,i_mw,oldi,oldj,oldi_mw,oldi_mt
DOUBLE PRECISION ::  P_F,SUMMATION213,U_temp,ThermCond,Prandtl,SatTemp
DOUBLE PRECISION, DIMENSION(2) :: MIN_D,oldMIN,oldT_max,oldw_min,P_limit,w_limit,h_limit, &
Ts_limit,T_limit,q_limit,T_max,w_mt,w_min,T_mw
DOUBLE PRECISION, DIMENSION(m) :: w_hso
DOUBLE PRECISION, DIMENSION(m,2) :: w_hs,T_out,U20,U21,U_bar
DOUBLE PRECISION, DIMENSION(m,n) :: T_Shso,q_flux_hso,h_hso,Re_Do,rho_o,mu_o,T_hso
DOUBLE PRECISION, DIMENSION(m,n,2) :: D_min,T_hs,q_flux_hs,rho,mu,Pij,T_sat,h_hs,T_Shs, &
Re_ave,f_fric,Re_D,k_cond,Pr
DOUBLE PRECISION, DIMENSION(3:m-2,2) :: MINi,oldMINi
INTEGER, DIMENSION(3:m-2,2) :: j_i,oldj_i
DOUBLE PRECISION, DIMENSION(m,2:n,2) :: phi_z,D_z,D3_z
DOUBLE PRECISION, DIMENSION(3:m-2,3:n-2,2) :: DELTA

!Product of uncertainty factors, used for heat flux & bulk fluid temperature
U_temp=U(1)*U(2)*U(3)*U(24)
!Pressure at inlet to the fuel assembly
P_F=P-(((7E-8*rho_in*g)/(g_c*144))*V_AS**2)+((rho_in*V_AS**2)/ &
(2*g_c*144*(448.831*A_fi)**2)) !psia

DO k=1,2
  !Minimum possible channel thicknesses
  DO i=1,m
    DO j=1,n
      IF (DFN==2) THEN
        D_min(i,j,k)=D(i,j,k)-e(k)+e_min(k)-DELTAD(k) !mil
        ELSE
          D_min(i,j,k)=D(i,j,k) !mil
      END IF
    END DO
  END DO
  DO i=1,m
    DO j=2,n
      !Product of power densities & axial increments, used in summation below
      phi_z(i,j,k)=((phi(i,j-1,k)+phi(i,j,k))/2)*DELTAz(j,k) !in
      !Product of channel thicknesses and axial increments; gets used in mass flow rate
      D_z(i,j,k)=(DELTAz(j,k)/1000)*((D_min(i,j-1,k)+D_min(i,j,k))/2) !in^2
      D3_z(i,j,k)=(DELTAz(j,k)*1000*12000**2)/((D_min(i,j-1,k)+ &
      D_min(i,j,k))/2)**3 !1/ft^2
    END DO
  END DO
  
  !Hot streak mass flow rates
  CALL MassFlow(fluid,k,m,n,A,t,f,Q,P,z(n,k),DELTAP_F,F_coef,T_in,rho_in, &
  Q_H2O,g,g_c,epsilonD_e,U_temp,e_e,U,w,phi_z,D_z,D3_z,w_hso)
  DO i=1,m
    w_hs(i,k)=w_hso(i) !lb_m/in-s
  END DO
  
  !Hot streak fluid temperature distrubtion
  CALL BulkFluidTemp(fluid,k,m,n,A,f,Q,P,F_coef,T_in,rho_in,Q_H2O,g,g_c, &
  epsilonD_e,U_temp,e_e,U,w_hs,z,D_min,phi_z,D_z,D3_z,T_hso,rho_o,mu_o,Re_Do)
  
  DO i=1,m
    DO j=1,n
      T_hs(i,j,k)=T_hso(i,j) !F
      rho(i,j,k)=rho_o(i,j) !lb_m/ft^3
      mu(i,j,k)=mu_o(i,j) !lb_m/ft-hr
      Re_ave(i,j,k)=Re_Do(i,j)
      k_cond(i,j,k)=ThermCond(fluid,T_hs(i,j,k),P) !Btu/hr-ft-F
      Re_D(i,j,k)=(2*w_hs(i,k)*12*3600)/mu(i,j,k)
      Pr(i,j,k)=Prandtl(fluid,T_hs(i,j,k),P)
    END DO
  END DO
 
  !Local nonbond factors
  DO i=1,m
    IF (k==1) THEN
      U20(i,k)=1.33687-(.35423*s(i,k))+(.14503*s(i,k)**2)-(.01669*s(i,k)**3)
      U21(i,k)=.863686-(.016507*s(i,k))-(.01095*s(i,k)**2)+(.0047976*s(i,k)**3)
      ELSE
        U20(i,k)=1.180171-(.278079*s(i,k))+(.151756*s(i,k)**2)- &
        (.014261*s(i,k)**3)
        U21(i,k)=.881393-(.249204*s(i,k))+(.181639*s(i,k)**2)- &
        (.033932*s(i,k)**3)
    END IF
  END DO
  !Product of local fuel segregation factor and nonbond factor
  DO i=1,m
    IF (HSF==1) THEN
      U_bar(i,k)=1
      ELSE
        IF (CCO==1) THEN
          U_bar(i,k)=U(18)*U20(i,k)
          ELSE
            U_bar(i,k)=U(19)*U21(i,k)
        END IF
    END IF
  END DO
  
  !Hot spot heat fluxes & surface temperatures for hot & cold channels
  CALL HotSpotHeatFlux(fluid,k,m,n,A,f,Q,P,U,U_bar,U25,z,D_min,T_hs,mu, &
  k_cond,Re_D,Pr,phi,T_Shso,q_flux_hso,h_hso)
  DO i=1,m
    DO j=1,n
      T_Shs(i,j,k)=T_Shso(i,j) !F
      q_flux_hs(i,j,k)=q_flux_hso(i,j) !Btu/hr-ft^2
      h_hs(i,j,k)=h_hso(i,j) !Btu/hr-ft^2-F
    END DO
  END DO
  
  DO i=1,m
    !Pressure & saturation temperature at the inlet
    Pij(i,1,k)=P_F+((g/g_c)*((rho(i,1,k)*z(1,k))/12**3))+(((1/(g_c*rho_in))- &
    (.04/(2*g_c*rho_in)))*((w_hs(i,k)*12000)/e_e(k))**2)-((1/(g_c* &
    rho(i,1,k)))*((w_hs(i,k)*12000)/D_min(i,1,k))**2) !psia
    T_sat(i,1,k)=SatTemp(fluid,Pij(i,1,k)) !F
    DO j=2,n
      CALL ColebrookFac(epsilonD_e,Re_ave(i,j,k),F_coef,f_fric(i,j,k))
      !Pressure of fluid at each position
      Pij(i,j,k)=P_F+((g/g_c)*((rho(i,j,k)*z(j,k))/12**3))+(((1/(g_c*rho_in))- &
      (.04/(2*g_c*rho_in)))*((w_hs(i,k)*12000)/e_e(k))**2)-((1/(g_c* &
      rho(i,j,k)))*((w_hs(i,k)*12000)/D_min(i,j,k))**2)-(((f_fric(i,j,k)* &
      w_hs(i,k)**2)/(4*g_c*rho(i,CEILING(j/2.),k)))* &
      SUMMATION213(1,m,2,n,1,2,k,i,2,j,D3_z)) !psia
      !Saturated temperature of fluid at each pressure
      T_sat(i,j,k)=SatTemp(fluid,Pij(i,j,k)) !F
    END DO 
  END DO
END DO

!Limiting criteria
IF (LCO==1.OR.LCO==2) THEN
  CALL Burnout(fluid,LCO,m,n,g,g_c,T_hs,T_sat,q_flux_hs,Pij,D_min,Re_D,Pr, &
  k_cond,DELTA)
END IF

IF (LCO==3) THEN
  CALL IncipientBoiling(m,n,U(23),T_Shs,T_sat,q_flux_hs,Pij,DELTA)
END IF

IF (LCO==4) THEN
  CALL FlowExcursion(fluid,m,n,w_hs,T_hs,T_sat,q_flux_hs,Pij,D_min,Re_D,Pr, &
  k_cond,DELTA)
END IF

IF (LCO==5) THEN
  CALL CriticalHeatFlux(fluid,m,n,g_c,w_hs,T_hs,T_sat,q_flux_hs,Pij, &
  D_min,DELTA)
END IF

IF (LCO==6) THEN
  CALL NetVaporGen(m,n,w_hs,T_hs,T_sat,q_flux_hs,rho,D_min,DELTA)
END IF 

IF (LCO==7) THEN
  DO k=1,2
    DO i=3,m-2
      DO j=3,m-2
        DELTA(i,j,k)=T_sat(i,j,k)-T_hs(i,j,k)
      END DO
    END DO
  END DO
END IF

DO k=1,2  
  !Limiting heat flux
  DO i=3,m-2
    MINi(i,k)=DELTA(i,3,k)
    j_i(i,k)=3
  END DO
  DO i=3,m-2
    DO j=4,n-2
      oldMINi(i,k)=MINi(i,k)
      oldj_i(i,k)=j_i(i,k)
      IF (DELTA(i,j,k)<MINi(i,k)) THEN
        MINi(i,k)=DELTA(i,j,k); j_i(i,k)=j
        ELSE
          MINi(i,k)=oldMINi(i,k)
          j_i(i,k)=oldj_i(i,k)
      END IF
    END DO
  END DO
  MIN_D(k)=MINi(3,k)
  j_limit(k)=j_i(3,k)
  i_limit(k)=3
  DO i=4,m-2
    oldMIN(k)=MIN_D(k)
    oldj(k)=j_limit(k)
    oldi(k)=i_limit(k)
    IF (MINi(i,k)<oldMIN(k)) THEN
      MIN_D(k)=MINi(i,k); i_limit(k)=i; j_limit(k)=j_i(i,k)
      ElSE
        MIN_D(k)=oldMIN(k)
        j_limit(k)=oldj(k)
        i_limit(k)=oldi(k)
    END IF
  END DO
  q_limit(k)=q_flux_hs(i_limit(k),j_limit(k),k) !Btu/hr-ft^2
  T_limit(k)=T_hs(i_limit(k),j_limit(k),k) !F
  Ts_limit(k)=T_Shs(i_limit(k),j_limit(k),k) !F
  h_limit(k)=h_hs(i_limit(k),j_limit(k),k) !Btu/hr-ft^2-F
  w_limit(k)=w_hs(i_limit(k),k) !lb_m/in-s
  P_limit(k)=Pij(i_limit(k),j_limit(k),k) !psia

  !Outlet bulk fluid temperatures
  DO i=3,m-2
    T_out(i,k)=T_hs(i,n,k)
  END DO
  !Maximum hot streak outlet bulk fluid temperature
  T_max(k)=T_out(3,k)
  i_mt(k)=3
  DO i=4,m-2
    oldT_Max(k)=T_max(k)
    oldi_mt(k)=i_mt(k)
    IF (T_out(i,k)>oldT_max(k)) THEN
      T_max(k)=T_out(i,k) !F
      i_mt(k)=i
      ELSE
        T_max(k)=oldT_max(k) !F
        i_mt(k)=oldi_mt(k)
    END IF
  END DO
  w_mt(k)=w_hs(i_mt(k),k) !lb_m/in-s
  
  !Minimum hot streak mass flow rate
  w_min(k)=w_hs(3,k)
  i_mw(k)=3
  DO i=4,m-2
    oldw_min(k)=w_min(k)
    oldi_mw(k)=i_mw(k)
    IF (w_hs(i,k)<oldw_min(k)) THEN
      w_min(k)=w_hs(i,k) !lb_m/in-s
      i_mw(k)=i
      ELSE
        w_min(k)=oldw_min(k) !lb_m/in-s
        i_mw(k)=oldi_mw(k)
    END IF
  END DO
  T_mw(k)=T_hs(i_mw(k),n,k) !F
    
END DO

!Overall limiting heat flux
IF (MIN_D(1)<MIN_D(2)) THEN
  MIN_DELTA=MIN_D(1)
  k_lhf=1; i_lhf=i_limit(1); j_lhf=j_limit(1);
  q_lhf=q_limit(1); T_lhf=T_limit(1); Ts_lhf=Ts_limit(1)
  h_lhf=h_limit(1); w_lhf=w_limit(1); P_lhf=P_limit(1)
  ELSE 
    MIN_DELTA=MIN_D(2)
    k_lhf=2; i_lhf=i_limit(2); j_lhf=j_limit(2);
    q_lhf=q_limit(2); T_lhf=T_limit(2); Ts_lhf=Ts_limit(2)
    h_lhf=h_limit(2); w_lhf=w_limit(2); P_lhf=P_limit(2)
END IF

!Overall maximum hot streak outlet bulk water temperature
IF (T_max(1)>T_max(2)) THEN
  k_mot=1; i_mot=i_mt(1); T_mot=T_max(1); w_mot=w_mt(1)
  ELSE
    k_mot=2; i_mot=i_mt(2); T_mot=T_max(2); w_mot=w_mt(2)
END IF

!Overall minimum flow rate
IF (w_min(1)<w_min(2)) THEN
  k_mfr=1; i_mfr=i_mw(1); T_mfr=T_mw(1); w_mfr=w_min(1)
  ELSE  
    k_mfr=2; i_mfr=i_mw(2); T_mfr=T_mw(2); w_mfr=w_min(2)
END IF

END SUBROUTINE

!$$$$$$
SUBROUTINE HotSpotHeatFlux(fluid,k,m,n,A,f,Q,P,U,U_bar,U25,z,D_min,T_hs,mu,k_cond,Re_D,Pr,phi,T_Shs,q_flux_hs,h_hs)
!$$$$$$ 
IMPLICIT NONE
CHARACTER*4, INTENT(IN) :: fluid
INTEGER, INTENT(IN) :: k,m,n
DOUBLE PRECISION, INTENT(IN) :: A,f,Q,P
DOUBLE PRECISION, INTENT(IN), DIMENSION(25) :: U
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,2) :: U_bar,U25
DOUBLE PRECISION, INTENT(IN), DIMENSION(n,2) :: z
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: D_min,T_hs,mu,k_cond,Re_D,Pr,phi
DOUBLE PRECISION, INTENT(OUT), DIMENSION(m,n) :: T_Shs,q_flux_hs,h_hs
INTEGER :: i,j
DOUBLE PRECISION :: Visc
DOUBLE PRECISION, DIMENSION(m,n) ::mu_S,Nu,DIFF_T,oldT_Shs

DO i=1,m
  DO j=1,n
    !Initial guess for surface temperature
    T_Shs(i,j)=T_hs(i,j,k) !F
  END DO
END DO

!Assume differences between older & newer values of surface temperature are initially 1
DIFF_T=1 !F
    
DO i=1,m
  DO j=1,n
    !Loop will keep running until surface temperature converges to within 0.1F
    DO WHILE (ABS(DIFF_T(i,j))>=.1)
      !Previous values of surface temperature are kept
      oldT_Shs(i,j)=T_Shs(i,j) !F
      !Viscosity of fluid at T_s(i,j)
      mu_S(i,j)=Visc(fluid,oldT_Shs(i,j),P) !lb_m/ft-hr
      !Nusselt number using a modified Hausen correlation
      IF (z(j,k)>0) THEN
        Nu(i,j)=.0235*(Re_D(i,j,k)**.8-230)*((1.8*Pr(i,j,k)**.3)-.8)*(1+((1./3.)* &
        ((2*D_min(i,j,k))/(1000*z(j,k)))**(2./3.)))*(mu(i,j,k)/mu_S(i,j))**.14
        ELSE
          Nu(i,j)=.0235*(Re_D(i,j,k)**.8-230)*((1.8*Pr(i,j,k)**.3)-.8)*(mu(i,j,k)/ &
          mu_S(i,j))**.14
      END IF
      !Heat transfer coefficient
      h_hs(i,j)=(U(8)*k_cond(i,j,k)*Nu(i,j)*12000)/(2*D_min(i,j,k)) !Btu/hr-ft^2-F
      !Fuel plate surface temperature; U25 is only used for j=29
      IF (j==n-2) THEN
        T_Shs(i,j)=T_hs(i,j,k)+((3.412E6*U(1)*U(2)*U(3)*U25(i,k)*((Q*f)/A)* &
        phi(i,j,k)*(1+((U_bar(i,k)-1)*(h_hs(i,j)/15000))))/h_hs(i,j)) !F
        ELSE
          T_Shs(i,j)=T_hs(i,j,k)+((3.412E6*U(1)*U(2)*U(3)*((Q*f)/A)*phi(i,j,k)* &
          (1+((U_bar(i,k)-1)*(h_hs(i,j)/15000))))/h_hs(i,j)) !F
      END IF
      !Difference between previous & current value of surface temperature
      DIFF_T(i,j)=T_Shs(i,j)-oldT_Shs(i,j)
    END DO
  END DO
END DO

!Hot spot heat flux
DO i=1,m
  DO j=1,n
    IF (j==n-2) THEN
      q_flux_hs(i,j)=3.412E6*U(1)*U(2)*U(3)*U25(i,k)*((Q*f)/A)*phi(i,j,k)* &
      (1+((U_bar(i,k)-1)*(h_hs(i,j)/15000))) !Btu/hr-ft^2
      ELSE
        q_flux_hs(i,j)=3.412E6*U(1)*U(2)*U(3)*((Q*f)/A)*phi(i,j,k)*(1+ &
        ((U_bar(i,k)-1)*(h_hs(i,j)/15000))) !Btu/hr-ft^2
    END IF
  END DO
END DO

END SUBROUTINE

!$$$$$$ 
SUBROUTINE Burnout(fluid,LCO,m,n,g,g_c,T_hs,T_sat,q_flux_hs,P,D_min,Re_D,Pr,k_cond,DELTA)
!$$$$$$
IMPLICIT NONE  
CHARACTER*4, INTENT(IN) :: fluid
INTEGER, INTENT(IN) :: LCO,m,n
DOUBLE PRECISION, INTENT(IN) :: g,g_c
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: T_hs,T_sat,q_flux_hs,P,D_min,Re_D,Pr,k_cond
DOUBLE PRECISION, INTENT(OUT), DIMENSION(3:m-2,3:n-2,2) :: DELTA
INTEGER :: i,j,k
DOUBLE PRECISION :: K_prime,C,K_boboil,x_f,x_g,SatDens,HeatVap,SurfTens,SatCp
DOUBLE PRECISION, DIMENSION(3:m-2,3:n-2,2) :: rho_f,rho_g,h_fg,sigma,cp_f,X,DTW,q_flux_bonb, &
q_flux_boboil,q_flux_bo

!Constants used in burnout correlation
IF (LCO==1) THEN
  K_prime=.019; C=.8; K_boboil=.1
  ELSE
    K_prime=.023; C=1; K_boboil=.18
END IF
!Static equilibrium quality of liquid & vapor
x_f=0.; x_g=1.

DO k=1,2
  DO i=3,m-2
    DO j=3,n-2
      !Saturation properties of fluid
      rho_f(i,j,k)=SatDens(fluid,P(i,j,k),x_f) !lb_m/ft^3
      rho_g(i,j,k)=SatDens(fluid,P(i,j,k),x_g) !lb_m/ft^3
      h_fg(i,j,k)=HeatVap(fluid,P(i,j,k)) !Btu/lb_m
      sigma(i,j,k)=SurfTens(fluid,T_sat(i,j,k)) !lb_f/ft
      cp_f(i,j,k)=SatCp(fluid,P(i,j,k),x_f) !Btu/lb_m-F
      !Superheat
      X(i,j,k)=T_sat(i,j,k)/1000
      IF (T_sat(i,j,k)<254) THEN
        DTW(i,j,k)=-127.04+(929.94*X(i,j,k))+(89.266*X(i,j,k)**2)- &
        (3714.4*X(i,j,k)**3) !F
        ELSE
          DTW(i,j,k)=165.51-(700.94*X(i,j,k))+(1259.4*X(i,j,k)**2)- &
          (837.21*X(i,j,k)**3) !F
      END IF
      !Nonboiling part of burnout heat flux
      q_flux_bonb(i,j,k)=K_prime*((k_cond(i,j,k)*12000)/(2*D_min(i,j,k)))* &
      Re_D(i,j,k)**.8*Pr(i,j,k)**(1./3.)*(DTW(i,j,k)+T_sat(i,j,k)- &
      T_hs(i,j,k)) !Btu/hr-ft^2
      !Boiling term
      q_flux_boboil(i,j,k)=K_boboil*3600*h_fg(i,j,k)*sqrt(rho_g(i,j,k))* &
      (sigma(i,j,k)*g_c*g*(rho_f(i,j,k)-rho_g(i,j,k)))**.25*(1+(C* &
      (rho_f(i,j,k)/rho_g(i,j,k))**.75*((cp_f(i,j,k)*(T_sat(i,j,k)- &
      T_hs(i,j,k)))/(9.8*h_fg(i,j,k))))) !Btu/hr-ft^2
      !Burnout heat flux
      q_flux_bo(i,j,k)=q_flux_bonb(i,j,k)+q_flux_boboil(i,j,k) !Btu/hr-ft^2
      !Difference between wall heat flux & corresponding burnout heat flux
      DELTA(i,j,k)=q_flux_bo(i,j,k)-q_flux_hs(i,j,k)
    END DO
  END DO
END DO

END SUBROUTINE

!$$$$$$ 
SUBROUTINE IncipientBoiling(m,n,U23,T_Shs,T_sat,q_flux_hs,P,DELTA)
!$$$$$$
IMPLICIT NONE  
INTEGER, INTENT(IN) :: m,n
DOUBLE PRECISION, INTENT(IN) :: U23
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: T_Shs,T_sat,q_flux_hs,P
DOUBLE PRECISION, INTENT(OUT), DIMENSION(3:m-2,3:n-2,2) :: DELTA
INTEGER :: i,j,k
DOUBLE PRECISION, DIMENSION(3:m-2,3:n-2,2) :: T_Sib

DO k=1,2
  DO i=3,m-2
    DO j=3,n-2
      !Incipient boiling surface temperature
      T_Sib(i,j,k)=T_sat(i,j,k)+(U23*(q_flux_hs(i,j,k)/(15.6* &
      P(i,j,k)**1.156))**(P(i,j,k)**.0234/2.3)) !F
      DELTA(i,j,k)=T_Sib(i,j,k)-T_Shs(i,j,k)
    END DO
  END DO
END DO

END SUBROUTINE

!$$$$$$ 
SUBROUTINE FlowExcursion(fluid,m,n,w_hs,T_hs,T_sat,q_flux_hs,P,D_min,Re_D,Pr,k_cond,DELTA)
!$$$$$$
IMPLICIT NONE  
CHARACTER*4, INTENT(IN) :: fluid
INTEGER, INTENT(IN) :: m,n
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,2) :: w_hs
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: T_hs,T_sat,q_flux_hs,P,D_min,Re_D,Pr,k_cond
DOUBLE PRECISION, INTENT(OUT), DIMENSION(3:m-2,3:n-2,2) :: DELTA
INTEGER :: i,j,k
DOUBLE PRECISION :: x_f,SatCp
DOUBLE PRECISION, DIMENSION(3:m-2,3:n-2,2) :: cp_f,eta_sub,G_flux,Pe,q_flux_fe

!Static equilibrium quality of liquid
x_f=0.

DO k=1,2
  DO i=3,m-2
    DO j=3,n-2
      !Saturated liquid heat capacity
      cp_f(i,j,k)=SatCp(fluid,P(i,j,k),x_f) !Btu/lb_m-F
      !Subcooling correction factor
      eta_sub(i,j,k)=.55+(20.178/(T_sat(i,j,k)-T_hs(i,j,k)))
      !Mass flux
      G_flux=(w_hs(i,k)*12*12000*3600)/D_min(i,j,k) !lb_m/ft^2-hr
      !Peclet number
      Pe(i,j,k)=Re_D(i,j,k)*Pr(i,j,k)
      !Flow excursion heat flux
      IF (Pe(i,j,k)<=7E4) THEN
        q_flux_fe(i,j,k)=455*eta_sub(i,j,k)*((k_cond(i,j,k)*12000)/ &
        D_min(i,j,k))*(T_sat(i,j,k)-T_hs(i,j,k)) !Btu/hr-ft^2
        ELSE
          q_flux_fe(i,j,k)=.0065*eta_sub(i,j,k)*G_flux(i,j,k)*cp_f(i,j,k)* &
          (T_sat(i,j,k)-T_hs(i,j,k)) !Btu/hr-ft^2
      END IF
      DELTA(i,j,k)=q_flux_fe(i,j,k)-q_flux_hs(i,j,k)
    END DO
  END DO
END DO

END SUBROUTINE

!$$$$$$ 
SUBROUTINE CriticalHeatFlux(fluid,m,n,g_c,w_hs,T_hs,T_sat,q_flux_hs,P,D_min,DELTA)
!$$$$$$
IMPLICIT NONE  
CHARACTER*4, INTENT(IN) :: fluid
INTEGER, INTENT(IN) :: m,n
DOUBLE PRECISION, INTENT(IN) :: g_c
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,2) :: w_hs
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: T_hs,T_sat,q_flux_hs,P,D_min
DOUBLE PRECISION, INTENT(OUT), DIMENSION(3:m-2,3:n-2,2) :: DELTA
INTEGER :: i,j,k
DOUBLE PRECISION :: C_1,C_2,C_3,C_4,C_5,x_f,x_g,Enthalpy,SatHeat,HeatVap,SatDens,SurfTens
DOUBLE PRECISION, DIMENSION(3:m-2,2) :: h_exit,h_f_exit,h_fg_exit,x_exit
DOUBLE PRECISION, DIMENSION(3:m-2,3:n-2,2) :: rho_f,rho_g,h_fg,sigma,G_flux,q_flux_crit

!Constants used in critical heat flux correlation
C_1=.0722; C_2=-.312; C_3=-.644; C_4=.9; C_5=.724
!Static equilibrium quality of liquid & vapor
x_f=0.; x_g=1.

DO k=1,2
  DO i=3,m-2
    !Outlet properties
    h_exit(i,k)=Enthalpy(fluid,T_hs(i,n,k),P(i,n,k)) !Btu/lb_m
    h_f_exit(i,k)=SatHeat(fluid,P(i,n,k),x_f) !Btu/lb_m
    h_fg_exit(i,k)=HeatVap(fluid,P(i,n,k)) !Btu/lb_m
    x_exit(i,k)=(h_exit(i,k)-h_f_exit(i,k))/h_fg_exit(i,k)
    DO j=3,n-2
      !Saturation properties
      rho_f(i,j,k)=SatDens(fluid,P(i,j,k),x_f) !lb_m/ft^3
      rho_g(i,j,k)=SatDens(fluid,P(i,j,k),x_g) !lb_m/ft^3
      h_fg(i,j,k)=HeatVap(fluid,P(i,j,k)) !Btu/lb_m
      sigma(i,j,k)=SurfTens(fluid,T_sat(i,j,k)) !lb_f/ft
      !Mass flux
      G_flux=(w_hs(i,k)*12*12000*3600)/D_min(i,j,k) !lb_m/ft^2-hr
      !Critical heat flux
      q_flux_crit(i,j,k)=C_1*G_flux(i,j,k)*h_fg(i,j,k)*((G_flux(i,j,k)**2* &
      2*D_min(i,j,k))/(12000.*3600.**2*g_c*rho_f(i,j,k)*sigma(i,j,k)))**C_2* &
      (rho_f(i,j,k)/rho_g(i,j,k))**C_3*(1-(C_4*(rho_f(i,j,k)/ &
      rho_g(i,j,k))**C_5*x_exit(i,k))) !Btu/hr-ft^2
      DELTA(i,j,k)=q_flux_crit(i,j,k)-q_flux_hs(i,j,k)
    END DO
  END DO
END DO

END SUBROUTINE

!$$$$$$ 
SUBROUTINE NetVaporGen(m,n,w_hs,T_hs,T_sat,q_flux_hs,rho,D_min,DELTA)
!$$$$$$
IMPLICIT NONE  
INTEGER, INTENT(IN) :: m,n
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,2) :: w_hs
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: T_hs,T_sat,q_flux_hs,rho,D_min
DOUBLE PRECISION, INTENT(OUT), DIMENSION(3:m-2,3:n-2,2) :: DELTA
INTEGER :: i,j,k
DOUBLE PRECISION :: K_nvg
DOUBLE PRECISION, DIMENSION(3:m-2,3:n-2,2) :: v_hs,q_flux_nvg

!Constant used in net vapor generation correlation
K_nvg=1.3167E-4

DO k=1,2
  DO i=3,m-2
    DO j=3,n-2
      !Velocity of hot streak fluid
      v_hs(i,j,k)=(w_hs(i,k)*12*12000)/(rho(i,j,k)*D_min(i,j,k)) !ft/s
      !Net vapor generation heat flux
      q_flux_nvg(i,j,k)=((T_sat(i,j,k)-T_hs(i,j,k))* &
      sqrt(v_hs(i,j,k)))/K_nvg !Btu/hr-ft^2
      DELTA(i,j,k)=q_flux_nvg(i,j,k)-q_flux_hs(i,j,k)
    END DO
  END DO
END DO

END SUBROUTINE