!$$$$$$ List of SUBROUTINES:
!$$$$$$
!$$$$$$ SUBROUTINE PartIIDeflections(fluid,l,m,n,noti,A,t,f,Q,P,DELTAP_F,F_coef,T_in,rho_in,Q_H2O,g,g_c,epsilonD_e,k_oxide,delta_VR, &
!$$$$$$ alpha_SP,lambda_g,DELTAP_nwfile,e_n,e_w,e_ne,e_we,DELTAD,xi,theta,U,w_A,s,DELTAs,DELTAz,z,T_MA,T_SPii,T_SPio, &
!$$$$$$ T_SPoi,T_SPoo,U4,U5,phi,psi_nH2,psi_nC2,psi_wH2,psi_wC2,T_SnH1,T_SnC1,T_SwH1,T_SwC1,w_n,w_w,D_n,D_w, &
!$$$$$$ T_SnH,T_SnC,T_SwH,T_SwC,psi_nH1,psi_nC1,psi_wH1,psi_wC1):
!$$$$$$ returns the narrow and wide coolant channel thicknesses, mass flow rates, and hot/cold surface temperatures
!$$$$$$ 

!$$$$$$ 
SUBROUTINE PartIIDeflections(fluid,l,m,n,noti,xi,A,t,f,Q,P,DELTAP_F,F_coef,T_in,rho_in,Q_H2O,g,g_c,epsilonD_e,k_oxide,delta_VR, &
alpha_SP,lambda_g,DELTAP_nwfile,e_n,e_w,e_ne,e_we,DELTAD,theta,U,w_A,s,DELTAs,DELTAz,z,T_MA,T_SPii,T_SPio, &
T_SPoi,T_SPoo,U4,U5,phi,psi_nH2,psi_nC2,psi_wH2,psi_wC2,T_SnH1,T_SnC1,T_SwH1,T_SwC1,w_n,w_w,D_n,D_w, &
T_SnH,T_SnC,T_SwH,T_SwC,psi_nH1,psi_nC1,psi_wH1,psi_wC1)
!$$$$$$ 
IMPLICIT NONE
CHARACTER*4, INTENT(IN) :: fluid
INTEGER, INTENT(IN) :: l,m,n,noti
INTEGER, INTENT(IN), DIMENSION(4) :: xi
DOUBLE PRECISION, INTENT(IN) :: A,t,f,Q,P,DELTAP_F,F_coef,T_in,rho_in,Q_H2O,g,g_c,epsilonD_e,k_oxide,delta_VR,alpha_SP
DOUBLE PRECISION, INTENT(IN), DIMENSION(2) :: lambda_g,DELTAP_nwfile,e_n,e_w,e_ne,e_we,DELTAD
DOUBLE PRECISION, INTENT(IN), DIMENSION(noti) :: theta
DOUBLE PRECISION, INTENT(IN), DIMENSION(25) :: U
DOUBLE PRECISION, INTENT(IN), DIMENSION(n) :: T_SPii,T_SPio,T_SPoi,T_SPoo
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,2) :: w_A,s
DOUBLE PRECISION, INTENT(IN), DIMENSION(m+1,2) :: DELTAs
DOUBLE PRECISION, INTENT(IN), DIMENSION(n,2) :: DELTAz,z,T_MA
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n) :: U4,U5
DOUBLE PRECISION, INTENT(IN), DIMENSION(m,n,2) :: phi,psi_nH2,psi_nC2,psi_wH2,psi_wC2,T_SnH1,T_SnC1,T_SwH1,T_SwC1
DOUBLE PRECISION, INTENT(OUT), DIMENSION(m,2) :: w_n,w_w
DOUBLE PRECISION, INTENT(OUT), DIMENSION(m,n,2) :: D_n,D_w,T_SnH,T_SnC,T_SwH,T_SwC,psi_nH1,psi_nC1,psi_wH1,psi_wC1
INTEGER :: i,j,k
DOUBLE PRECISION :: pi,U_temp,SUMMATION12,SUMMATION123,MAX_T
DOUBLE PRECISION, DIMENSION(2) :: T_MnwA_H,T_MnwA_C,ER_MnwA_H,ER_MnwA_C,DELTAP_nw
DOUBLE PRECISION, DIMENSION(5) :: DIFF_T
DOUBLE PRECISION, DIMENSION(m) :: w_no,w_wo
DOUBLE PRECISION, DIMENSION(n) :: DIFF_TMnnH,DIFF_TMnnC,DIFF_TMnwH,DIFF_TMnwC,DIFF_TMww
DOUBLE PRECISION, DIMENSION(m,2) :: DELTAP_nwi,rho_ex_w,rho_ex_n,delta_DELTA_PH,delta_DELTA_PC,oldw_n,oldw_w, &
DELTAP_inlet,DELTAP_exit
DOUBLE PRECISION, DIMENSION(n,2) :: ER_MnnH,ER_MnnC,ER_MnwH,ER_MnwC,ER_Mww,deltaT_MnnH,deltaT_MnnC, &
deltaT_MnwH,deltaT_MnwC,deltaT_Mww,oldT_MnnH,oldT_MnnC,oldT_MnwH,oldT_MnwC,oldT_Mww,T_MnnH, &
T_MnnC,T_MnwH,T_MnwC,T_Mww,sinz
DOUBLE PRECISION, DIMENSION(m,n) :: T_SnHo,T_SnCo,T_SwHo,T_SwCo,mu_no,mu_wo,T_wo,T_no,rho_no,rho_wo, &
Re_Dno,Re_Dwo
DOUBLE PRECISION, DIMENSION(m,n,2) :: q_flux_H,q_flux_C,T_MnH,T_MnC,T_MwH,T_MwC,T_n,T_w,delta_onH, &
delta_onC,delta_owH,delta_owC,X_nH,X_nC,X_wH,X_wC,T_MnnHi,T_MnnCi,T_MnwHi,T_MnwCi,T_Mwwi,delta_nnH, &
delta_nnC,delta_nwH,delta_nwC,delta_ww,delta_MnnH,delta_MnnC,delta_MnwH,delta_MnwC,delta_Mww, &
delta_Tnw_TAC,delta_Tnw_TAH,delta_Tnn_TAC,delta_Tnn_TAH,delta_Tww_TA,D_1n,D_1w,mu_n,mu_w
DOUBLE PRECISION, DIMENSION(m,2:n,2) :: phi_z,D_nz,D_wz,D3_nz,D3_wz
DOUBLE PRECISION, DIMENSION(4:m-2,n,2) :: T_MnnSH,T_MnnSC,T_MnwSH,T_MnwSC,T_MwwS
DOUBLE PRECISION, DIMENSION(4:n-2,2) :: T_MnwZH,T_MnwZC
DOUBLE PRECISION, DIMENSION(4:m-2,2) :: DPS

pi=3.14159265359
U_temp=U(1)*U(2)*U(3)

DO k=1,2
  !Product of power densities & axial increments, used in summation for bulk fluid temperatures
  DO i=1,m
    DO j=2,n
      phi_z(i,j,k)=((U4(i,j)+U5(i,j))/2)*((phi(i,j-1,k)+phi(i,j,k))/2)* &
      DELTAz(j,k) !in
    END DO
  END DO
  !Heat flux at each position for hot & cold plates
  DO i=1,m
    DO j=1,n
      q_flux_H(i,j,k)=3.412E6*U_temp*U4(i,j)* &
      ((Q*f)/A)*phi(i,j,k) !Btu/hr-ft^2
      q_flux_C(i,j,k)=3.412E6*U_temp*U5(i,j)* &
      ((Q*f)/A)*phi(i,j,k) !Btu/hr-ft^2
    END DO
  END DO
    
  !For the first iteration of the DO loop, assume the following quantities for fuel plates between channels:
  !Assume initial values for the mass flow rates are those calculated from Part I
  DO i=1,m
    w_n(i,k)=w_A(i,k) !lb_m/in-s
    w_w(i,k)=w_A(i,k) !lb_m/in-s
  END DO
  !Assume initial values for average metal temperatures down the length of a fuel plate
  DO j=1,n
    T_MnnH(j,k)=T_MA(j,k); T_MnnC(j,k)=T_MA(j,k) !F
    T_MnwH(j,k)=T_MA(j,k); T_MnwC(j,k)=T_MA(j,k) !F
    T_Mww(j,k)=T_MA(j,k) !F
  END DO
  DO i=1,m
    DO j=1,n
      !Assume initial values for metal temperatures of the fuel plate
      T_MnnHi(i,j,k)=T_MA(j,k);  T_MnnCi(i,j,k)=T_MA(j,k) !F
      T_MnwHi(i,j,k)=T_MA(j,k);  T_MnwCi(i,j,k)=T_MA(j,k) !F
      T_Mwwi(i,j,k)=T_MA(j,k) !F
      !Assume initial values for the oxide film thicknesses
      X_nH(i,j,k)=0; X_nC(i,j,k)=0 !mil
      X_wH(i,j,k)=0; X_wC(i,j,k)=0 !mil
    END DO
  END DO
  !Initial guess for the average pressure difference across the
  !fuel plate between a narrow and wide coolant channel
  DELTAP_nw(k)=DELTAP_nwfile(k)
  !Initial guess for the maximum difference between older & newer values of average metal temperatures
  MAX_T=1
  
  !This loop will keep running until the maximum average metal temperature difference between older & 
  !newer values, which includes narrow-narrow hot, narrow-narrow cold, narrow-wide hot, narrow-wide 
  !cold, & wide-wide, is within 0.1 F. This ensures that all of the average metal temperatures converge
  !to within 0.1 F.
  DO WHILE(MAX_T>=.1)
    !Keeping the previous values of the mass flow rate
    DO i=1,m
      oldw_n(i,k)=w_n(i,k) !lb_m/in-s
      oldw_w(i,k)=w_w(i,k) !lb_m/in-s
    END DO
    !Keeping the previous values of the average metal temperatures
    DO j=1,n
      oldT_MnnH(j,k)=T_MnnH(j,k); oldT_MnnC(j,k)=T_MnnC(j,k) !F
      oldT_MnwH(j,k)=T_MnwH(j,k); oldT_MnwC(j,k)=T_MnwC(j,k) !F
      oldT_Mww(j,k)=T_Mww(j,k) !F
    END DO
    !Average value of T_Mnw, hot & cold
    DO j=4,n-2
      T_MnwZH(j,k)=((oldT_MnwH(j-1,k)+oldT_MnwH(j,k))/2)*DELTAz(j,k) !F-in
      T_MnwZC(j,k)=((oldT_MnwC(j-1,k)+oldT_MnwC(j,k))/2)*DELTAz(j,k) !F-in
    END DO
    T_MnwA_H(k)=SUMMATION12(4,n-2,1,2,k,4,n-2,T_MnwZH)/ &
    SUMMATION12(1,n,1,2,k,4,n-2,DELTAz) !F
    T_MnwA_C(k)=SUMMATION12(4,n-2,1,2,k,4,n-2,T_MnwZC)/ &
    SUMMATION12(1,n,1,2,k,4,n-2,DELTAz) !F
    !Elastic Modulus ratio at T=T_MnwA
    ER_MnwA_H(k)=(-1.624E-6*T_MnwA_H(k)**2)+(4.719E-4*T_MnwA_H(k))+.9737
    ER_MnwA_C(k)=(-1.624E-6*T_MnwA_C(k)**2)+(4.719E-4*T_MnwA_C(k))+.9737
    !Deflection for hot and cold plates due to differential pressure
    DO i=1,m
      IF (k==1) THEN
        delta_DELTA_PH(i,k)=((U(10)*DELTAP_nw(k))/ER_MnwA_H(k))*((2.8372E-2* &
        s(i,k)**5)-(.20491*s(i,k)**4)+(.27529*s(i,k)**3)+(.57806*s(i,k)**2)- &
        (.89329*s(i,k))) !mil
        delta_DELTA_PC(i,k)=((U(10)*DELTAP_nw(k))/ER_MnwA_C(k))*((2.8372E-2* &
        s(i,k)**5)-(.20491*s(i,k)**4)+(.27529*s(i,k)**3)+(.57806*s(i,k)**2)- &
        (.89329*s(i,k))) !mil
        ELSE
          delta_DELTA_PH(i,k)=((U(10)*DELTAP_nw(k))/ER_MnwA_H(k))*((3.0799E-2* &
          s(i,k)**5)-(.197*s(i,k)**4)+(.23697*s(i,k)**3)+(.42645*s(i,k)**2)- &
          (.59068*s(i,k))) !mil
          delta_DELTA_PC(i,k)=((U(10)*DELTAP_nw(k))/ER_MnwA_C(k))*((3.0799E-2* &
          s(i,k)**5)-(.197*s(i,k)**4)+(.23697*s(i,k)**3)+(.42645*s(i,k)**2)- &
          (.59068*s(i,k))) !mil
      END IF
    END DO  
    !Deflection due to temperature difference between an individual plate and that of an average fuel plate
    !For the fuel plate between two narrow coolant channels
    DO j=1,n
      ER_MnnH(j,k)=(-1.624E-6*oldT_MnnH(j,k)**2)+(4.719E-4*oldT_MnnH(j,k))+.9737
      ER_MnnC(j,k)=(-1.624E-6*oldT_MnnC(j,k)**2)+(4.719E-4*oldT_MnnC(j,k))+.9737
    END DO
    DO i=1,m
      DO j=1,n
        IF (k==1) THEN
          delta_Tnn_TAH(i,j,k)=((U(11)*(oldT_MnnH(j,k)-T_MA(j,k)))/ &
          ER_MnnH(j,k))*((-3.2084E-4*s(i,k)**5)+(4.3177E-3*s(i,k)**4)- &
          (.018077*s(i,k)**3)+(.011087*s(i,k)**2)+(.043908*s(i,k))) !mil
          delta_Tnn_TAC(i,j,k)=((U(11)*(oldT_MnnC(j,k)-T_MA(j,k)))/ &
          ER_MnnC(j,k))*((-3.2084E-4*s(i,k)**5)+(4.3177E-3*s(i,k)**4)- &
          (.018077*s(i,k)**3)+(.011087*s(i,k)**2)+(.043908*s(i,k))) !mil
          ELSE
            delta_Tnn_TAH(i,j,k)=((U(11)*(oldT_MnnH(j,k)-T_MA(j,k)))/ &
            ER_MnnH(j,k))*((-1.8518E-4*s(i,k)**5)+(4.5322E-3*s(i,k)**4)- &
            (.02105*s(i,k)**3)-(6.0732E-4*s(i,k)**2)+(.083033*s(i,k))) !mil
            delta_Tnn_TAC(i,j,k)=((U(11)*(oldT_MnnC(j,k)-T_MA(j,k)))/ &
            ER_MnnC(j,k))*((-1.8518E-4*s(i,k)**5)+(4.5322E-3*s(i,k)**4)- &
            (.02105*s(i,k)**3)-(6.0732E-4*s(i,k)**2)+(.083033*s(i,k))) !mil
        END IF
      END DO
    END DO
    !For the fuel plate between a narrow coolant channel and a wide coolant channel
    DO j=1,n
      ER_MnwH(j,k)=(-1.624E-6*oldT_MnwH(j,k)**2)+(4.719E-4*oldT_MnwH(j,k))+.9737
      ER_MnwC(j,k)=(-1.624E-6*oldT_MnwC(j,k)**2)+(4.719E-4*oldT_MnwC(j,k))+.9737
    END DO
    DO i=1,m
      DO j=1,n
        IF (k==1) THEN
          delta_Tnw_TAH(i,j,k)=((U(11)*(oldT_MnwH(j,k)-T_MA(j,k)))/ &
          ER_MnwH(j,k))*((-3.2084E-4*s(i,k)**5)+(4.3177E-3*s(i,k)**4)- &
          (.018077*s(i,k)**3)+(.011087*s(i,k)**2)+(.043908*s(i,k))) !mil
          delta_Tnw_TAC(i,j,k)=((U(11)*(oldT_MnwC(j,k)-T_MA(j,k)))/ &
          ER_MnwC(j,k))*((-3.2084E-4*s(i,k)**5)+(4.3177E-3*s(i,k)**4)- &
          (.018077*s(i,k)**3)+(.011087*s(i,k)**2)+(.043908*s(i,k))) !mil
          ELSE
            delta_Tnw_TAH(i,j,k)=((U(11)*(oldT_MnwH(j,k)-T_MA(j,k)))/ &
            ER_MnwH(j,k))*((-1.8518E-4*s(i,k)**5)+(4.5322E-3*s(i,k)**4)- &
            (.02105*s(i,k)**3)-(6.0732E-4*s(i,k)**2)+(.083033*s(i,k))) !mil
            delta_Tnw_TAC(i,j,k)=((U(11)*(oldT_MnwC(j,k)-T_MA(j,k)))/ &
            ER_MnwC(j,k))*((-1.8518E-4*s(i,k)**5)+(4.5322E-3*s(i,k)**4)- &
            (.02105*s(i,k)**3)-(6.0732E-4*s(i,k)**2)+(.083033*s(i,k))) !mil
        END IF
      END DO
    END DO
    !For the fuel plate between two wide coolant channels
    DO j=1,n
      ER_Mww(j,k)=(-1.624E-6*oldT_Mww(j,k)**2)+(4.719E-4*oldT_Mww(j,k))+.9737
    END DO
    DO i=1,m
      DO j=1,n
        IF(k==1) THEN
          delta_Tww_TA(i,j,k)=((U(11)*(oldT_Mww(j,k)-T_MA(j,k)))/ &
          ER_Mww(j,k))*((-3.2084E-4*s(i,k)**5)+(4.3177E-3*s(i,k)**4)- &
          (.018077*s(i,k)**3)+(.011087*s(i,k)**2)+(.043908*s(i,k))) !mil
          ELSE
            delta_Tww_TA(i,j,k)=((U(11)*(oldT_Mww(j,k)-T_MA(j,k)))/ &
            ER_Mww(j,k))*((-1.8518E-4*s(i,k)**5)+(4.5322E-3*s(i,k)**4)- &
            (.02105*s(i,k)**3)-(6.0732E-4*s(i,k)**2)+(.083033*s(i,k))) !mil
        END IF
      END DO
    END DO
    !Decrease in coolant channel thickness because of the increase in fuel plate thickness due to heating
    DO i=1,m
      DO j=1,n
        !For fuel plate between two narrow coolant channels
        delta_nnH(i,j,k)=U(12)*.5*alpha_sp*t*(T_MnnHi(i,j,k)-70) !mil
        delta_nnC(i,j,k)=U(12)*.5*alpha_sp*t*(T_MnnCi(i,j,k)-70) !mil
        !For fuel plates between a wide coolant channel and a narrow coolant channel
        delta_nwH(i,j,k)=U(12)*.5*alpha_sp*t*(T_MnwHi(i,j,k)-70) !mil
        delta_nwC(i,j,k)=U(12)*.5*alpha_sp*t*(T_MnwCi(i,j,k)-70) !mil
        !For fuel plates between two wide coolant channels
        delta_ww(i,j,k)=U(12)*.5*alpha_sp*t*(T_Mwwi(i,j,k)-70) !mil
      END DO
    END DO

    !Increase in fuel plate thickness due to oxide formation
    DO i=1,m
      DO j=1,n
        IF (X_nH(i,j,k)<=3) THEN
          delta_onH(i,j,k)=.2026*X_nH(i,j,k) !mil
          ELSE
            delta_onH(i,j,k)=.6078 !mil
        END IF
        IF (X_nC(i,j,k)<=3) THEN
          delta_onC(i,j,k)=.2026*X_nC(i,j,k) !mil
          ELSE
            delta_onC(i,j,k)=.6078 !mil
        END IF
        IF (X_wH(i,j,k)<=3) THEN
          delta_owH(i,j,k)=.2026*X_wH(i,j,k) !mil
          ELSE
            delta_owH(i,j,k)=.6078 !mil
        END IF
        IF (X_wC(i,j,k)<=3) THEN
          delta_owC(i,j,k)=.2026*X_wC(i,j,k) !mil
          ELSE
            delta_owC(i,j,k)=.6078 !mil
        END IF
      END DO
    END DO
    !Average difference between fuel plate temperature and side plate temperature
    DO j=1,n
      !For plates between two narrow coolant channels
      IF (k==1) THEN
        DELTAT_MnnH(j,k)=oldT_MnnH(j,k)-((T_SPii(j)+T_SPio(j))/2) !F
        DELTAT_MnnC(j,k)=oldT_MnnC(j,k)-((T_SPii(j)+T_SPio(j))/2) !F
        ELSE
          DELTAT_MnnH(j,k)=oldT_MnnH(j,k)-((T_SPoi(j)+T_SPoo(j))/2) !F
          DELTAT_MnnC(j,k)=oldT_MnnC(j,k)-((T_SPoi(j)+T_SPoo(j))/2) !F
      END IF
      !For plates between a narrow and wide coolant channel
      IF (k==1) THEN
        DELTAT_MnwH(j,k)=oldT_MnwH(j,k)-((T_SPii(j)+T_SPio(j))/2) !F
        DELTAT_MnwC(j,k)=oldT_MnwC(j,k)-((T_SPii(j)+T_SPio(j))/2) !F
        ELSE
          DELTAT_MnwH(j,k)=oldT_MnwH(j,k)-((T_SPoi(j)+T_SPoo(j))/2) !F
          DELTAT_MnwC(j,k)=oldT_MnwC(j,k)-((T_SPoi(j)+T_SPoo(j))/2) !F
      END IF
      !For plates between two wide coolant channels
      IF (k==1) THEN
        DELTAT_Mww(j,k)=oldT_Mww(j,k)-((T_SPii(j)+T_SPio(j))/2) !F
        ELSE
          DELTAT_Mww(j,k)=oldT_Mww(j,k)-((T_SPoi(j)+T_SPoo(j))/2) !F
      END IF
    END DO
    !Deflection of fuel plate due to the temperature difference between fuel plate and side plate
    DO j=1,n
      IF (z(j,k)<6) THEN
        sinz(j,k)=sin((pi*z(j,k))/12.)
        ELSE IF (z(j,k)<18) THEN
          sinz(j,k)=1.
          ELSE 
            sinz(j,k)=sin((pi*(24.-z(j,k)))/12.)
      END IF
    END DO
    DO i=1,m
      DO j=1,n
        IF(DELTAT_MnnH(j,k)<=0) THEN
          delta_MnnH(i,j,k)=0
          ELSE
            delta_MnnH(i,j,k)=U(14)*.0864*DELTAT_MnnH(j,k)*(sin((pi*s(i,k))/ &
            SUMMATION12(1,m+1,1,2,k,2,m,DELTAs)))*sinz(j,k) !mil
        END IF
        IF(DELTAT_MnnC(j,k)<=0) THEN
          delta_MnnC(i,j,k)=0
          ELSE
            delta_MnnC(i,j,k)=U(14)*.0864*DELTAT_MnnC(j,k)*(sin((pi*s(i,k))/ &
            SUMMATION12(1,m+1,1,2,k,2,m,DELTAs)))*sinz(j,k) !mil
        END IF
        IF(DELTAT_MnwH(j,k)<=0) THEN
          delta_MnwH(i,j,k)=0
          ELSE
            delta_MnwH(i,j,k)=U(14)*.0864*DELTAT_MnwH(j,k)*(sin((pi*s(i,k))/ &
            SUMMATION12(1,m+1,1,2,k,2,m,DELTAs)))*sinz(j,k) !mil
        END IF
        IF(DELTAT_MnwC(j,k)<=0) THEN
          delta_MnwC(i,j,k)=0
          ELSE
            delta_MnwC(i,j,k)=U(14)*.0864*DELTAT_MnwC(j,k)*(sin((pi*s(i,k))/ &
            SUMMATION12(1,m+1,1,2,k,2,m,DELTAs)))*sinz(j,k) !mil
        END IF
        IF(DELTAT_Mww(j,k)<=0) THEN
          delta_Mww(i,j,k)=0
          ELSE
            delta_Mww(i,j,k)=U(14)*.0864*DELTAT_Mww(j,k)*(sin((pi*s(i,k))/ &
            SUMMATION12(1,m+1,1,2,k,2,m,DELTAs)))*sinz(j,k) !mil
        END IF
      END DO
    END DO
    !Thicknesses of the coolant channels
    DO i=1,m
      DO j=1,n
        !For the narrow channel
        IF (xi(1)==1.OR.xi(2)==1) THEN
          D_1n(i,j,k)=e_n(k)+delta_DELTA_PH(i,k)+delta_DELTA_PC(i,k)- &
          delta_Tnw_TAH(i,j,k)+delta_Tnw_TAC(i,j,k)-delta_nwH(i,j,k)- &
          delta_nwC(i,j,k)-delta_MnwH(i,j,k)+delta_MnwC(i,j,k)-delta_onH(i,j,k)- &
          delta_onC(i,j,k)-(2*delta_VR) !mil
          ELSE IF (xi(3)==1) THEN
            D_1n(i,j,k)=e_n(k)+delta_DELTA_PC(i,k)-delta_Tnn_TAH(i,j,k)+ &
            delta_Tnw_TAC(i,j,k)-delta_nnH(i,j,k)-delta_nwC(i,j,k)- &
            delta_MnnH(i,j,k)+delta_MnwC(i,j,k)-delta_onH(i,j,k)- &
            delta_onC(i,j,k)-(2*delta_VR) !mil
            ELSE IF (xi(4)==1) THEN
              D_1n(i,j,k)=e_n(k)-delta_Tnn_TAH(i,j,k)+delta_Tnn_TAC(i,j,k)- &
              delta_nnH(i,j,k)-delta_nnC(i,j,k)-delta_MnnH(i,j,k)+ &
              delta_MnnC(i,j,k)-delta_onH(i,j,k)-delta_onC(i,j,k)-(2*delta_VR) !mil
        END IF
        D_n(i,j,k)=D_1n(i,j,k)-(DELTAD(k)*sin((2*pi*z(j,k))/lambda_g(k))) !mil
        !For the wide channel
        IF (xi(1)==1) THEN
          D_1w(i,j,k)=e_w(k)-delta_DELTA_PH(i,k)-delta_DELTA_PC(i,k)+ &
          delta_Tnw_TAH(i,j,k)-delta_Tnw_TAC(i,j,k)-delta_nwH(i,j,k)- &
          delta_nwC(i,j,k)+delta_MnwH(i,j,k)-delta_MnwC(i,j,k)-delta_owH(i,j,k)- &
          delta_owC(i,j,k)-(2*delta_VR) !mil
          ELSE IF (xi(2)==1) THEN
            D_1w(i,j,k)=e_w(k)-delta_DELTA_PH(i,k)+delta_Tnw_TAH(i,j,k)- &
            delta_Tww_TA(i,j,k)-delta_nwH(i,j,k)-delta_ww(i,j,k)+ &
            delta_MnwH(i,j,k)-delta_Mww(i,j,k)-delta_owH(i,j,k)- &
            delta_owC(i,j,k)-(2*delta_VR) !mil
            ELSE IF (xi(3)==1) THEN
              D_1w(i,j,k)=e_w(k)-delta_DELTA_PH(i,k)-delta_DELTA_PC(i,k)- &
              delta_Tnw_TAH(i,j,k)+delta_Tnw_TAC(i,j,k)-delta_nwH(i,j,k)- &
              delta_nwC(i,j,k)-delta_MnwH(i,j,k)+delta_MnwC(i,j,k)-delta_owH(i,j,k)- &
              delta_owC(i,j,k)-(2*delta_VR) !mil
        END IF
        D_w(i,j,k)=D_1w(i,j,k)-(DELTAD(k)*sin((2*pi*z(j,k))/lambda_g(k))) !mil
      END DO
    END DO
    !Product of channel thicknesses and axial increments; gets used in mass flow rate    
    DO i=1,m
      DO j=2,n
        D_nz(i,j,k)=(DELTAz(j,k)/1000)*((D_n(i,j-1,k)+D_n(i,j,k))/2) !in^2
        D_wz(i,j,k)=(DELTAz(j,k)/1000)*((D_w(i,j-1,k)+D_w(i,j,k))/2) !in^2
        D3_nz(i,j,k)=(DELTAz(j,k)*1000*12000**2)/((D_n(i,j-1,k)+ &
        D_n(i,j,k))/2)**3 !1/ft^2
        D3_wz(i,j,k)=(DELTAz(j,k)*1000*12000**2)/((D_w(i,j-1,k)+ &
        D_w(i,j,k))/2)**3 !1/ft^2
      END DO
    END DO
    !Mass flow rates down each strip.  Using Alternate Method II
    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_ne,U,oldw_n,phi_z,D_nz,D3_nz,w_no)
    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_we,U,oldw_w,phi_z,D_wz,D3_wz,w_wo)
    DO i=1,m
      w_n(i,k)=w_no(i) !lb_m/in-s
      w_w(i,k)=w_wo(i) !lb_m/in-s
    END DO
    !Bulk water temperature distribution for narrow coolant channel
    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_ne,U,w_n,z,D_n,phi_z,D_nz,D3_nz,T_no,rho_no,mu_no,Re_Dno)
    DO i=1,m
      DO j=1,n
        T_n(i,j,k)=T_no(i,j) !F
        mu_n(i,j,k)=mu_no(i,j) !lb_m/ft-hr
      END DO
    END DO
    !Bulk water temperature distribution for wide coolant channel
    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_we,U,w_w,z,D_w,phi_z,D_wz,D3_wz,T_wo,rho_wo,mu_wo,Re_Dwo)
    DO i=1,m
      DO j=1,n
        T_w(i,j,k)=T_wo(i,j) !F
        mu_w(i,j,k)=mu_wo(i,j) !lb_m/ft-hr
      END DO
    END DO
    !Fuel plate surface temperature at each position, in narrow and wide, for hot and cold
    CALL SurfaceTemp(fluid,k,m,n,U(8),P,w_n,z,D_n,T_n,mu_n,q_flux_H,T_SnHo)
    CALL SurfaceTemp(fluid,k,m,n,U(8),P,w_n,z,D_n,T_n,mu_n,q_flux_C,T_SnCo)
    CALL SurfaceTemp(fluid,k,m,n,U(8),P,w_w,z,D_w,T_w,mu_w,q_flux_H,T_SwHo)
    CALL SurfaceTemp(fluid,k,m,n,U(8),P,w_w,z,D_w,T_w,mu_w,q_flux_C,T_SwCo)
    DO i=1,m
      DO j=1,n
        T_SnH(i,j,k)=T_SnHo(i,j) !F
        T_SnC(i,j,k)=T_SnCo(i,j) !F
        T_SwH(i,j,k)=T_SwHo(i,j) !F
        T_SwC(i,j,k)=T_SwCo(i,j) !F
      END DO
    END DO
    !Oxide film thickness at each position
    DO i=1,m
      DO j=1,n
        IF(l==1) THEN
          X_nH(i,j,k)=443*U(9)*theta(l)**.778*exp(-8280/(T_SnH(i,j,k)+ &
          459.67)) !mil
          X_nC(i,j,k)=443*U(9)*theta(l)**.778*exp(-8280/(T_SnC(i,j,k)+ &
          459.67)) !mil
          X_wH(i,j,k)=443*U(9)*theta(l)**.778*exp(-8280/(T_SwH(i,j,k)+ &
          459.67)) !mil
          X_wC(i,j,k)=443*U(9)*theta(l)**.778*exp(-8280/(T_SwC(i,j,k)+ &
          459.67)) !mil
          ELSE
            psi_nH1(i,j,k)=(psi_nH2(i,j,k)+theta(l-1))*exp((10643* &
            (T_SnH1(i,j,k)-T_SnH(i,j,k)))/((T_SnH1(i,j,k)+459.67)* &
            (T_SnH(i,j,k)+459.67))) !hr
            psi_nC1(i,j,k)=(psi_nC2(i,j,k)+theta(l-1))*exp((10643* &
            (T_SnC1(i,j,k)-T_SnC(i,j,k)))/((T_SnC1(i,j,k)+459.67)* &
            (T_SnC(i,j,k)+459.67))) !hr
            psi_wH1(i,j,k)=(psi_wH2(i,j,k)+theta(l-1))*exp((10643* &
            (T_SwH1(i,j,k)-T_SwH(i,j,k)))/((T_SwH1(i,j,k)+459.67)* &
            (T_SwH(i,j,k)+459.67))) !hr
            psi_wC1(i,j,k)=(psi_wC2(i,j,k)+theta(l-1))*exp((10643* &
            (T_SwC1(i,j,k)-T_SwC(i,j,k)))/((T_SwC1(i,j,k)+459.67)* &
            (T_SwC(i,j,k)+459.67))) !hr
            X_nH(i,j,k)=443*U(9)*(theta(l)+psi_nH1(i,j,k))**.778*exp(-8280/ &
            (T_SnH(i,j,k)+459.67)) !mil
            X_nC(i,j,k)=443*U(9)*(theta(l)+psi_nC1(i,j,k))**.778*exp(-8280/ &
            (T_SnC(i,j,k)+459.67)) !mil
            X_wH(i,j,k)=443*U(9)*(theta(l)+psi_wH1(i,j,k))**.778*exp(-8280/ &
            (T_SwH(i,j,k)+459.67)) !mil
            X_wC(i,j,k)=443*U(9)*(theta(l)+psi_wC1(i,j,k))**.778*exp(-8280/ &
            (T_SwC(i,j,k)+459.67)) !mil
        END IF
      END DO
    END DO
    !Metal surface (metal-oxide interface) temperature at each position
    DO i=1,m
      DO j=1,n
        IF (X_nH(i,j,k)<=3) THEN
          T_MnH(i,j,k)=T_SnH(i,j,k)+((q_flux_H(i,j,k)*X_nH(i,j,k))/(12000* &
          k_oxide)) !F
          ELSE
            T_MnH(i,j,k)=T_SnH(i,j,k)+((3*q_flux_H(i,j,k))/(12000*k_oxide)) !F
        END IF
        IF (X_nC(i,j,k)<=3) THEN
          T_MnC(i,j,k)=T_SnC(i,j,k)+((q_flux_C(i,j,k)*X_nC(i,j,k))/(12000* &
          k_oxide)) !F
          ELSE
            T_MnC(i,j,k)=T_SnC(i,j,k)+((3*q_flux_C(i,j,k))/(12000*k_oxide)) !F
        END IF
        IF (X_wH(i,j,k)<=3) THEN
          T_MwH(i,j,k)=T_SwH(i,j,k)+((q_flux_H(i,j,k)*X_wH(i,j,k))/(12000* &
          k_oxide)) !F
          ELSE
            T_MwH(i,j,k)=T_SwH(i,j,k)+((3*q_flux_H(i,j,k))/(12000*k_oxide)) !F
        END IF
        IF (X_wC(i,j,k)<=3) THEN
          T_MwC(i,j,k)=T_SwC(i,j,k)+((q_flux_C(i,j,k)*X_wC(i,j,k))/(12000* &
          k_oxide)) !F
          ELSE
            T_MwC(i,j,k)=T_SwC(i,j,k)+((3*q_flux_C(i,j,k))/(12000*k_oxide)) !F
        END IF
      END DO
    END DO
    !Average metal temperature at each position
    DO i=1,m
      DO j=1,n
        T_MnnHi(i,j,k)=T_MnH(i,j,k) !F
        T_MnnCi(i,j,k)=T_MnC(i,j,k) !F
        T_MnwHi(i,j,k)=(T_MnH(i,j,k)+T_MwH(i,j,k))/2 !F
        T_MnwCi(i,j,k)=(T_MnC(i,j,k)+T_MwC(i,j,k))/2 !F
        T_Mwwi(i,j,k)=T_MwC(i,j,k) !F
      END DO 
    END DO
    !Average metal temperature down the length of the fuel plates
    DO i=4,m-2
      DO j=1,n
        T_MnnSH(i,j,k)=((T_MnnHi(i-1,j,k)+T_MnnHi(i,j,k))/2)*DELTAs(i,k) !F-in
        T_MnnSC(i,j,k)=((T_MnnCi(i-1,j,k)+T_MnnCi(i,j,k))/2)*DELTAs(i,k) !F-in
        T_MnwSH(i,j,k)=((T_MnwHi(i-1,j,k)+T_MnwHi(i,j,k))/2)*DELTAs(i,k) !F-in
        T_MnwSC(i,j,k)=((T_MnwCi(i-1,j,k)+T_MnwCi(i,j,k))/2)*DELTAs(i,k) !F-in
        T_MwwS(i,j,k)=((T_Mwwi(i-1,j,k)+T_Mwwi(i,j,k))/2)*DELTAs(i,k) !F-in
      END DO
    END DO
    DO j=1,n
      T_MnnH(j,k)=SUMMATION123(4,m-2,1,n,1,2,k,j,4,m-2,T_MnnSH)/ &
      SUMMATION12(1,m+1,1,2,k,4,m-2,DELTAs) !F
      T_MnnC(j,k)=SUMMATION123(4,m-2,1,n,1,2,k,j,4,m-2,T_MnnSC)/ &
      SUMMATION12(1,m+1,1,2,k,4,m-2,DELTAs) !F
      T_MnwH(j,k)=SUMMATION123(4,m-2,1,n,1,2,k,j,4,m-2,T_MnwSH)/ &
      SUMMATION12(1,m+1,1,2,k,4,m-2,DELTAs) !F
      T_MnwC(j,k)=SUMMATION123(4,m-2,1,n,1,2,k,j,4,m-2,T_MnwSC)/ &
      SUMMATION12(1,m+1,1,2,k,4,m-2,DELTAs) !F
      T_Mww(j,k)=SUMMATION123(4,m-2,1,n,1,2,k,j,4,m-2,T_MwwS)/ &
      SUMMATION12(1,m+1,1,2,k,4,m-2,DELTAs) !F
    END DO
  
    !Exit coolant density in the narrow & wide coolant channels
    DO i=1,m
      rho_ex_n(i,k)=rho_no(i,n) !lb_m/ft^3
      rho_ex_w(i,k)=rho_wo(i,n) !lb_m/ft^3
    END DO
    !Average differential pressure across a strip of unit width at s(i,k)
    DO i=1,m
      DELTAP_inlet(i,k)=((1.04*rho_in)/(2*g_c*144))*(((w_w(i,k)*12*12000)/ &
      (rho_in*e_we(k)))**2-((w_n(i,k)*12*12000)/(rho_in*e_ne(k)))**2) !psia
      DELTAP_exit(i,k)=((1-(1-(e_we(k)/(e_we(k)+t)))**2)*(rho_ex_w(i,k)/ &
      (2*g_c*144))*((w_w(i,k)*12*12000)/(rho_ex_w(i,k)*e_we(k)))**2) &
      -((1-(1-(e_ne(k)/(e_ne(k)+t)))**2)*(rho_ex_n(i,k)/(2*g_c*144))* &
      ((w_n(i,k)*12*12000)/(rho_ex_n(i,k)*e_ne(k)))**2) !psia
      DELTAP_nwi(i,k)=(DELTAP_inlet(i,k)+DELTAP_exit(i,k))/2 !psia
    END DO
    !Average differential pressure across a the entire width of a fuel plate between a narrow and wide coolant channel
    DO i=4,m-2
      DPs(i,k)=((DELTAP_nwi(i-1,k)+DELTAP_nwi(i,k))/2)*DELTAs(i,k) !psia-in
    END DO
    DELTAP_nw(k)=SUMMATION12(4,m-2,1,2,k,4,m-2,DPS)/ &
    SUMMATION12(1,m+1,1,2,k,4,m-2,DELTAs) !psia
    !Differences between previous & current values of average metal temperatures
    DO j=1,n
      DIFF_TMnnH(j)=ABS(T_MnnH(j,k)-oldT_MnnH(j,k)) !F
      DIFF_TMnnC(j)=ABS(T_MnnC(j,k)-oldT_MnnC(j,k)) !F
      DIFF_TMnwH(j)=ABS(T_MnwH(j,k)-oldT_MnwH(j,k)) !F
      DIFF_TMnwC(j)=ABS(T_MnwC(j,k)-oldT_MnwC(j,k)) !F
      DIFF_TMww(j)=ABS(T_Mww(j,k)-oldT_Mww(j,k)) !F
    END DO
    !Maximum differences from each type of plate
    DIFF_T(1)=MAXVAL(DIFF_TMnnH); DIFF_T(2)=MAXVAL(DIFF_TMnnC) !F
    DIFF_T(3)=MAXVAL(DIFF_TMnwH); DIFF_T(4)=MAXVAL(DIFF_TMnwC) !F
    DIFF_T(5)=MAXVAL(DIFF_TMww) !F
    !Maximum overall difference
    MAX_T=MAXVAL(DIFF_T)
  END DO
END DO

END SUBROUTINE