SUBROUTINE TYPE222 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)

!************************************************************************
! Object: External Receiver/Tower
! Simulation Studio Model: Type222
! 
! Author: Michael J. Wagner
! Date:	 March 14, 2008 
! Last modified: February 5th, 2009
! COPYRIGHT 2010 NATIONAL RENEWABLE ENERGY LABORATORY
!
! Parameters
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Input Units   |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  N_panels                      |  Number of individual panels on the receiver                    |  none          |  none
!  2|  D_rec                         |  The overall outer diameter of the receiver                     |  m             |  m
!  3|  H_rec                         |  The height of the receiver                                     |  m             |  m
!  4|  THT                           |  The height of the tower (hel. pivot pt to rec equator)         |  m             |  m
!  5|  D_out                         |  The outer diameter of an individual receiver tube              |  mm            |  m
!  6|  th_tu                         |  The wall thickness of a single receiver tube                   |  mm            |  m
!  7|  Material                      |  The material name of the receiver tubes                        |  none          |  none
!  8|  Coolant                       |  The name of the HTF used in the receiver                       |  none          |  none
!  9|  Flow_type                     |  A flag indicating which flow pattern is used                   |  none          |  none
! 10|  PLAT                          |  The lattitude of the power plant                               |  deg           |  deg
! 11|  LU_flux                       |  The logical unit of the flux map file                          |  none          |  none
! 12|  epsilon                       |  The emissivity of the receiver surface coating                 |  none          |  none
! 13|  hl_ffact                      |  The heat loss factor (thermal loss fudge factor)               |  none          |  none
!xxxxxxxxxxxx 14|  nuel = 8                      |  The number of solar elevation pts in the efficiency matrix     |  none          |  none
!xxxxxxxxxxxx 15|  nazm = 12                     |  The number of solar azimuth pts in the efficiency matrix       |  none          |  none

! Inputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Input Units   |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  azimuth                       |  Solar azimuth angle                                            |  deg           |  deg
!  2|  zenith                        |  Solar zenith angle                                             |  deg           |  deg
!  3|  T_salt_hot                    |  Desired HTF outlet temperature                                 |  C             |  K
!  4|  T_salt_cold                   |  Desired HTF inlet temperature                                  |  C             |  K
!  5|  V_wind                        |  Ambient wind velocity, ground level                            |  m/s           |  m/s
!  6|  P_amb                         |  Ambient atmospheric pressure                                   |  atm           |  Pa
!  7|  eta_pump                      |  Receiver HTF pump efficiency                                   |  none          |  none
!  8|  hour                          |  Hour of the day                                                |  hr            |  hr
!  9|  T_dp                          |  Ambient dew point temperature                                  |  C             |  C
! 10|  I_bn                          |  Direct (beam) normal radiation                                 |  kJ/m^2-hr     |  kW/m^2
! 11|  field_eff                     |  Heliostat field efficiency                                     |  none          |  none
! 12|  T_db                          |  Ambient dry bulb temperature                                   |  C             |  K
! 13|  night_recirc                  |  Flag to indicate night recirculation through the rec.          |  none          |  none
! 14|  hel_stow_deploy               |  Heliostat field stow/deploy solar angle                        |  deg           |  deg

! Outputs
!-----------------------------------------------------------------------------------------------------------------------------------------
!Nb |  Variable                      |  Description                                                    |  Output Units  |  Internal Units
!---|--------------------------------|-----------------------------------------------------------------|----------------|-----------------
!  1|  m_dot_salt_tot                |  Total HTF flow rate through the receiver                       |  kg/hr         |  kg/s
!  2|  eta_therm                     |  Receiver thermal efficiency                                    |  none          |  none
!  3|  W_dot_pump                    |  Receiver pump power                                            |  MWe           |  We
!  4|  q_conv_sum                    |  Receiver convective losses                                     |  MWt           |  Wt
!  5|  q_rad_sum                     |  Receiver radiative losses                                      |  MWt           |  Wt
!  6|  Q_thermal                     |  Receiver total thermal losses for timestep                     |  MWt           |  Wt
!  7|  T_salt_hot                    |  HTF outlet temperature                                         |  C             |  K
!  8|  -not named-                   |  Receiver power prior to thermal losses                         |  MWt           |  NA


!TRNSYS acess functions (allow to acess TIME etc.) 
USE TrnsysConstants
USE TrnsysFunctions

!DEC$ATTRIBUTES DLLEXPORT :: TYPE222				!SET THE CORRECT TYPE NUMBER HERE

!-----------------------------------------------------------------------------------------------------------------------
!    TRNSYS DECLARATIONS
IMPLICIT NONE			!REQUIRES THE USER TO DEFINE ALL VARIABLES BEFORE USING THEM

real(8):: XIN, OUT, TIME, PAR, STORED, T, DTDT 	
integer(4):: INFO(15),IUNIT,ITYPE,ICNTRL,NITEMS
integer,parameter:: NP=15,NI=14,NOUT=10,ND=0,NSTORED=0

DIMENSION:: XIN(NI),OUT(NOUT),PAR(NP),STORED(NSTORED),T(ND),DTDT(ND)
!-----------------------------------------------------------------------------------------------------------------------

!ADD DECLARATIONS AND DEFINITIONS FOR THE USER-VARIABLES HERE
real(8)::m_dot_salt_tot,eta_therm,W_dot_pump,q_conv_sum,q_rad_sum,Q_thermal, I_bn_inst,P_amb_kPa, q_inc_sum, q_abs_sum
real(8),allocatable:: T_s(:)

!PARAMETERS
real(8):: N_panels, D_rec, H_rec, THT, D_out, th_tu,Material, Coolant, Flow_type, PLAT, LU_flux, epsilon,hl_ffact, nuel, nuaz

!INPUTS
real(8)::azimuth, zenith, T_salt_hot, T_salt_cold,V_wind, P_amb, eta_pump, hour, T_dp, I_bn, field_eff, T_db,&
         night_recirc, hel_stow_deploy
Logical::is_there
character::checkname*300, test*300
integer(4):: ios, i, qq


!-----------------------------------------------------------------------------------------------------------------------
!    SET THE VERSION INFORMATION FOR TRNSYS
IF(INFO(7).EQ.-2) THEN
    INFO(12)=16
    RETURN 1
ENDIF

!    DO ALL THE VERY LAST CALL OF THE SIMULATION MANIPULATIONS HERE
IF (INFO(8).EQ.-1) THEN

    !Deallocate the surface temperature array
    if(allocated(T_s)) deallocate(T_s)
    
    !Close the flux map file
    inquire(unit=LU_flux,opened=is_there,name=checkname)
    if (is_there) close(LU_flux)

    RETURN 1
ENDIF

!    PERFORM ANY 'AFTER-ITERATION' MANIPULATIONS THAT ARE REQUIRED HERE
!    e.g. save variables to storage array for the next timestep
IF (INFO(13).GT.0) THEN
    RETURN 1
ENDIF
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!    DO ALL THE VERY FIRST CALL OF THE SIMULATION MANIPULATIONS HERE
IF (INFO(7).EQ.-1) THEN
!-----------------------------------------------------------------------------------------------------------------------
    !Parameters
    N_panels=PAR(1)
    D_rec=PAR(2)
    H_rec=PAR(3)
    THT=PAR(4)
    D_out=PAR(5)
    th_tu=PAR(6)
    Material=PAR(7)
    Coolant=PAR(8)
    Flow_type=PAR(9)
    PLAT=PAR(10)
    LU_flux=PAR(11)
    epsilon=PAR(12)
    hl_ffact=PAR(13)

    !allocate space for the T_s array
    if(.not.allocated(T_s)) allocate(T_s(nint(N_panels)))

    !SET SOME INFO ARRAY VARIABLES TO TELL THE TRNSYS ENGINE HOW THIS TYPE IS TO WORK
    INFO(6)=NOUT				
    INFO(9)=1				
    INFO(10)=0	

    !Open the fluxmap file used for generating the flux distribution on the receiver
    inquire(unit=LU_flux,opened=is_there,name=checkname)
    if (.not.is_there) then
        open(unit=LU_flux,file=trim(checkname),iostat=ios,status="OLD",position="REWIND")
        !Check to see if the file was opened successfully
        if(ios.ne.0) then
            call messages(-1,"TRNSYS did not find the fluxmap file",'FATAL',INFO(1),INFO(2))
            stop
        endif
    endif


    !We need to check to make sure that the fluxmap file has been successfully
    !created by the PTGEN program.  Do this by counting the number of lines 
    !in the file to make sure that it is not zero or close to zero.
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    i=0 !initialize counting variable
    do 
        i=i+1
        if (i.gt.10) goto 400
        read(unit=LU_flux,fmt="(A)",err=400,eor=400,end=400, advance="YES") test
    enddo

    400 continue

    if(i.le.2) then
        call messages(-1,"TRNSYS did not find the fluxmap file",'FATAL',INFO(1),INFO(2))
        stop  !Quit the program here, returning through the hierarchy is cumbersome
    endif
    rewind(LU_flux)
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  
    
    !RETURN TO THE CALLING PROGRAM
    RETURN 1
ENDIF
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
!    DO ALL OF THE INITIAL TIMESTEP MANIPULATIONS HERE - THERE ARE NO ITERATIONS AT THE INTIAL TIME
IF (TIME .LT. (getSimulationStartTime() + getSimulationTimeStep()/2.D0)) THEN

    !SET THE UNIT NUMBER FOR FUTURE CALLS
    IUNIT=INFO(1)
    ITYPE=INFO(2)

    !PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE
    OUT = 0.d0
    OUT(7)=565.0
    
    !RETURN TO THE CALLING PROGRAM
    RETURN 1
ENDIF
!-----------------------------------------------------------------------------------------------------------------------

!-----------------------------------------------------------------------------------------------------------------------
! * * * Calculations for each timestep * * * 
!-----------------------------------------------------------------------------------------------------------------------

azimuth=XIN(1)    !Solar azimuth
zenith=XIN(2)     !Solar zenith
T_salt_hot=XIN(3)+273.15 !Desired hot temp, Convert from deg C to deg K
T_salt_cold=XIN(4)+273.15 !Specified inlet temp, Convert from C to K
V_wind=XIN(5)     !Wind velocity
P_amb=XIN(6)      !ambient pressure
eta_pump=XIN(7)   !Receiver HTF pump efficiency
hour=XIN(8)       !Hour of the day
T_dp=XIN(9)       !Dewpoint temp
I_bn=XIN(10)      !Beam normal radiation
field_eff=XIN(11) !Field efficiency value
T_db=XIN(12)      !Dry bulb temp
night_recirc=XIN(13) !Night recirculation control 0=empty receiver, 1=recirculate
hel_stow_deploy=XIN(14) !Solar elevation angle at which heliostats are stowed
    IUNIT=INFO(1)
    ITYPE=INFO(2)

!Do necessary conversions
I_bn_inst = I_bn/3.6   !This converts from kJ/m^2-hr to kW/m^2

call CRS((hour-12.),azimuth,zenith,D_rec,H_rec,THT,D_out,th_tu,T_salt_hot,T_salt_cold,T_db,T_dp,V_wind,P_amb,I_bn_inst, &    !Inputs
        Material,Coolant,eta_pump,int(N_panels),epsilon,int(Flow_type),field_eff,hl_ffact,night_recirc, hel_stow_deploy,int(LU_flux), &
        m_dot_salt_tot,eta_therm,W_dot_pump,q_conv_sum,q_rad_sum,T_s,Q_thermal,time,info,qq,q_inc_sum,q_abs_sum)  !Outputs


! SET THE OUTPUTS FROM THIS MODEL IN SEQUENTIAL ORDER AND GET OUT

OUT(1)=(m_dot_salt_tot)*3600 !Convert to kg/hr
OUT(2)=eta_therm
OUT(3)=W_dot_pump/1.e6 !Convert We to MWe
OUT(4)=q_conv_sum/1.e6 !Convert Wt to MWt
OUT(5)=q_rad_sum/1.e6  !Convert Wt to MWt
OUT(6)=Q_thermal/1.e6  !Convert Wt to MWt
OUT(7)=	T_salt_hot - 273.15 !Convert from K back to C
OUT(8)=OUT(6)+abs(OUT(5)+OUT(4)) !Power before thermal losses
out(9)= 0. !
OUT(10) = q_inc_sum/1.e6 !MWt  !Total incident power on the receiver

RETURN 1
END
!-----------------------------------------------------------------------------------------------------------------------

!#######################################################################################################################
!	THE FOLLOWING SECTION CONTAINS THE SUBROUTINES CALLED BY TYPE222
!#######################################################################################################################


!#######################################################################################################################
subroutine CRS(hour,azimuth,zenith,D_rec,H_rec,THT,D_out,th_tu,T_salt_hot,T_salt_cold,T_db,T_dp,V_wind, P_amb,I_bn,Material,&
               Coolant,eta_pump,N_panels,epsilon,flowtype,field_eff,hl_ffact,night_recirc, hel_stow_deploy,LU_flux,&
               m_dot_salt_tot,eta_therm,W_dot_pump,q_conv_sum,q_rad_sum,T_s,Q_thermal,time,info,qq, q_inc_sum, q_abs_sum)

implicit none

!Define input variables ----------------------------------------------------
integer(4),intent(in)::N_panels,flowtype,LU_flux,INFO(15)
real(8)::azimuth,zenith,D_rec,H_rec,THT, T_salt_hot, T_salt_cold,T_db,T_dp,P_amb,eta_pump,D_out,th_tu,I_bn,&
         hour,field_eff,epsilon,hl_ffact, night_recirc, hel_stow_deploy
real(8)::V_wind,time,material,coolant

!Define output variables
real(8),intent(out)::m_dot_salt_tot,eta_therm,q_conv_sum, q_rad_sum, Q_thermal
real(8),dimension(N_panels),intent(out)::T_s

!Define local variables
integer(4),dimension(N_panels)::Panel,ceil,flo
integer(4),allocatable::Flow_pattern(:)
integer,dimension(2)::salt_out
integer::i,j,qq,nlines,daycall,hourcall,qq_max
real(8),parameter::pi=3.14159265,grav=9.81
real(8)::D_tube,th_tube,D_inner,D_in,A_tube,n_t,N_tube,A_receiver,q_guess,c_guess,A_rec_proj,A_node,P_atm,sigma,&
         flux_in(12),m_dot_salt,array(12),err,T_salt_hotX, m_dot_saltX,T_coolant_prop,T_s_ave,T_film_ave,&
         k_film,mu_film,rho_film,c_p_film,Re_for,ksD,Nusselt_for,h_for,conductivity,Viscosity,Density,specheat,Nusselt_FC,beta,nu_amb,&
         volexpcoef,m,T_sky,q_inc_sum,q_abs_sum,C_p_coolant,RelRough,mu_coolant,k_coolant,rho_coolant,u_coolant,&
         Re_inner,Pr_inner,h_inner,R_conv_inner,Nusselt_t,f,LoverD,m_dot_tube,L_e_45,L_e_90,DELTAP_tube,DELTAP_45,DELTAP_90,DELTAP,&
         DELTAP_THT,DELTAP_net,Pres_D,W_dot_pump,W_dot_pump_hp,T_amb, azi_adj,skytemp,tol
         
real(8),dimension(N_panels)::ind,Psp_field,P_field,T_panel_out,T_panel_in,T_sX,T_panel_outX,T_panel_inX,T_film,&
                             h_rad_amb,h_rad_sky,q_dot_amb,q_dot_sky,q_dot_rad,q_dot_conv,q_dot_loss,q_loss_sum,&
                             q_dot_inc,q_dot_abs,T_wall,k_tube,R_tube_wall,ppos,T_panel_ave,Gr_nat,Nusselt_nat,h_nat,h_mixed
save::flow_pattern,salt_out,nlines
!Initialize sensitive variables
daycall=0 ; hourcall=0 ; ppos=0
!---------------------------------------------------------------------------

!Receiver dimensions, parameters
D_tube = D_out/1000.     !Value in meters
th_tube = th_tu/1000.     !Thickness of the tube
D_inner = D_tube - 2.*th_tube	      !Diameter of each receiver tube
D_in = D_inner*1000.      !Inner diameter of the tube, D_in set from D_out and th_tube
A_tube = pi*D_tube/2.*H_rec           !Outer area of each tube
n_t = Nint((pi*D_rec)/(D_tube*N_panels))  !The number of tubes per panel, as a function of the number
									!  of panels and the desired diameter of the receiver
N_tube =n_t*N_panels                  !Number of tubes in the system
A_receiver = pi**2*D_rec*H_rec        !The area of the receiver exposed to the ambient air, [S&K]
A_rec_proj = D_tube*H_rec*N_tube      !The projected area of the tubes on a plane parallel to the center lines of the tubes
T_amb = T_db+273.15           !Ambient temperature
T_sky = skytemp(T_amb,(T_dp+273.15),hour)     !The effective sky temperature [K]
A_node = pi*D_rec/N_panels*H_rec      !The area associated with each node
P_atm=P_amb*101325.     !Ambient pressure, in [Pa]
azi_adj = azimuth + 180.              !By TRNSYS convention, the azimuth angle is 0 at due south, negative to the east,
								! and positive to the west. The range is then -180 to 180. By the convention used
								! here, the azimuth is 0 at due north, and ranges clockwise from 0 to 360. This adjusts.
!Other parameters which are useful for calculating heat flow
sigma = 5.670e-8                       ![W/m^2-K^4] Stefan-Boltzmann constant
!epsilon = set as parameter                         !Taumoefolau, T. et al., 2004

!Correct the windspeed using the logarithmic profile law (HOMER documentation, Wind Shear Inputs section)
V_wind = log((THT+H_rec/2.)/.003)/log(10./.003)*V_wind

!Get flow pattern
if(.not.allocated(Flow_pattern)) then
    allocate(Flow_pattern(N_panels))
    call flowPatterns(N_panels,flowtype,Flow_pattern,salt_out,nlines)
else
    if(time>=8760.) deallocate(Flow_pattern)
endif


!Do an initial check to make sure the solar position called is valid.
!  If its not, return the output equal to zeros. Also check to make sure 
!  the solar flux is at a certain level, otherwise the correlations
!  aren't valid
if(((zenith.gt.(90.-hel_stow_deploy)).or.(I_bn.le.1.e-6)).or.((zenith.eq.0.).and.(azimuth.eq.0))) then
    if(night_recirc.eq.1.)then
        I_bn=0.
    else
        m_dot_salt_tot=0 ; eta_therm=0 ; W_dot_pump=0  
        q_conv_sum=0 ; q_rad_sum=0 ; T_s=0 ;  Q_thermal=0
        q_inc_sum=0; q_abs_sum=0
        return
    endif
endif

!Get the values of the flux from the fluxmap and store them as flux_in(col,row)
if(I_bn.gt.1.) then
    call fluxinterp(LU_flux,zenith,azi_adj,array,daycall,hourcall,info)
    flux_in(:)=array(:)/(950.)*I_bn*field_eff  !The weather-adjusted, efficiency-adjusted flux values
else
    flux_in=0.
endif

!Translate to the number of panels, so each panel has its own linearly interpolated flux value
do j=1,N_panels
    Panel(j) = j	                       !The position of each panel
    ppos(j)=(12./dble(N_panels)*(j-1.)+6./dble(N_panels))+1.
    flo(j)=floor(ppos(j))
    ceil(j)=ceiling(ppos(j))
    ind(j)= (ppos(j)-flo(j))/max(dble(ceil(j)-flo(j)),1.e-6)
    if(ceil(j).gt.12.) ceil(j)=1
    Psp_field(j)=ind(j)*(flux_in(ceil(j))-flux_in(flo(j)))+flux_in(flo(j))  !Average area-specific power for each node
    P_field(j)=A_node*Psp_field(j)	        !The power incident on each node
enddo

!Guess values -------------------------------------------------------------
if(night_recirc.eq.1.) then
    T_sX = T_salt_hot                !Guess temperature for the surface nodes
    T_panel_outX = (T_salt_hot+T_salt_cold)/2.       !Guess values for the fluid temp coming out of the control volume
    T_panel_inX = (T_salt_hot+T_salt_cold)/2.    !Guess values for the fluid temperature coming into the control volume
else
    T_sX = T_salt_hot                !Guess temperature for the surface nodes
    T_panel_outX = T_salt_cold       !Guess values for the fluid temp coming out of the control volume
    T_panel_inX = T_salt_cold        !Guess values for the fluid temperature coming into the control volume
endif
      
if(I_bn.gt.1.e-6) then
    q_guess = 0.5*sum(P_field(1:N_panels)) !Estimate the thermal power produced by the receiver [kWt].
                                         !0.5 coefficient is optimized based on run-time 
                                         ![0.2|27.77s, 0.4|25.05s, 0.6|25.87s, 0.8|26.27s, 0.9|26.3s, 1.0|26.6s]
    c_guess = specheat(Coolant, (T_salt_hot+T_salt_cold)/2.) !Estimate the specific heat of the fluid in the receiver [kJ/kg-K]
    m_dot_saltX = q_guess/(c_guess*(T_salt_hot-T_salt_cold)*dble(nlines))     !coolant mass flow rate guess value
else !The tower recirculates at night (based on earlier condition) 
    !Enter recirculation mode, where inlet/outlet temps switch
    T_salt_hot = T_salt_cold
    T_salt_cold = T_sX(1) !T_sX is set to T_salt_hot before.. just a way of switching
    c_guess = specheat(Coolant, (T_salt_hot+T_salt_cold)/2.) !Estimate the specific heat of the fluid in the receiver [kJ/kg-K]
    m_dot_saltX = -3500./(c_guess*(T_salt_hot-T_salt_cold)/2.) !Estimate the flow rate to be some low number
endif
T_salt_hotX = 9999.                   !Initial value for error calculation
!--------------------------------------------------------------------------

!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
!                            ITERATION STARTS HERE
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
qq_max=50
do qq=1,qq_max+1
    !if the problem fails to converge after 50 iterations, then the power is likely negligible and
    !..the zero set can be returned
    if(qq.gt.qq_max)then
        m_dot_salt_tot=0 ; eta_therm=0 ; W_dot_pump=0
        q_conv_sum=0 ; q_rad_sum=0 ; T_s=0 ;  Q_thermal=0
        return
    endif

    err=(T_salt_hotX - T_salt_hot)/T_salt_hot  !The basis on which convergence is determined
    if(night_recirc.eq.1.) then
        tol= 0.005
    else
        tol= 0.001 
    endif

    if(abs(err).lt.tol) exit        !Check to see if the problem has converged
    
    !Set the variables equal to their newly calculated guess values----------
    T_s(:) = T_sX(:)
    T_panel_out(:) = T_panel_outX(:)
    T_panel_in(:) = T_panel_inX(:)

    m_dot_salt = m_dot_saltX
    !------------------------------------------------------------------------

    !Now do the actual calculations
    T_panel_ave(:) = (T_panel_in(:)+T_panel_out(:))/2.  !The average coolant temperature in each control volume
    T_film(:) = (T_s(:)+T_amb)/2.
    T_coolant_prop = (T_salt_hot + T_salt_cold)/2.             !The temperature at which the coolant properties are evaluated. Validated as constant
    T_salt_hotX = sum(T_panel_out(salt_out(1:nlines)))/dble(nlines)!Calculates the mixed outlet temperature of the salt

    !Calculate the average surface temperature so that the forced convection coefficient can be determined
    T_s_ave = sum(T_s(1:N_panels))/(dble(N_panels))
    T_film_ave = (T_amb+T_salt_hot)/2.

    !Convection coefficient for external forced convection using Siebers & Kraabel
    k_film = Conductivity(1.d0,T_film_ave)                    !The conductivity of the air
    mu_film = Viscosity(1.d0,T_film_ave)                      !Dynamic viscosity of the air
    rho_film = Density(1.d0, T_film_ave, P_atm)               !Density of the air
    c_p_film = specheat(1.d0,T_film_ave)*1000.    !Specific heat of the air
    Re_for = rho_film*V_wind*D_rec/mu_film                     !Reynolds number
    ksD = (D_tube/2.)/D_rec                                    !The effective roughness of the cylinder [Siebers, Kraabel 1984]
    Nusselt_for=Nusselt_FC(ksD,Re_for)                         ![S&K]
    h_for = Nusselt_for*k_film/D_rec*hl_ffact                  !The forced convection heat transfer coefficient

    !Convection coefficient for external natural convection using Siebers & Kraabel"
    !Note: This relationship applies when the surrounding properties are evaluated at ambient conditions [S&K]"
    beta = 1./T_amb                                          !Volumetric expansion coefficient
    nu_amb = Viscosity(1.d0,T_amb)/Density(1.d0,T_amb,P_atm) !Kinimatic viscosity
    Gr_nat(:) = grav*beta*(T_s(:)-T_amb)*H_rec**3/nu_amb**2      !Grashof Number at ambient conditions
    Nusselt_nat(:) = .098*Gr_nat(:)**(1./3.)*(T_s(:)/T_amb)**(-.14) !Nusselt number
    h_nat(:) = Nusselt_nat(:)*k_film/H_rec*hl_ffact                  !The natural convection cofficient

    !Mixed convection
    m=3.2
    h_mixed(:) = (h_for**m+h_nat(:)**m)**(1./m)              !The mixed convection relationship
    q_dot_conv(:) = h_mixed(:)*A_node*(T_s(:) - T_film(:))   !Convection losses per node

    !Radiation from the receiver
    !Calculate the radiation node by node
    h_rad_amb(:) = sigma*epsilon*(T_s(:)**2+T_amb**2)*(T_s(:)+T_amb)*hl_ffact  !The radiation coefficient for amb
    h_rad_sky(:) = sigma*epsilon*(T_s(:)**2+T_sky**2)*(T_s(:)+T_sky)*hl_ffact  !The radiation coef. for sky
    q_dot_amb(:) = .5*h_rad_amb(:)*A_node*(T_s(:)-T_amb)  !amb losses per node
    q_dot_sky(:) = .5*h_rad_sky(:)*A_node*(T_s(:) - T_sky)  !sky losses per node

    !Calculate the losses from the surface"
    q_dot_rad(:) = q_dot_amb(:)+q_dot_sky(:)		     !Total radiation losses per node
    q_dot_loss(:) = q_dot_rad(:)+q_dot_conv(:)		     !Total overall losses per node

    q_conv_sum = sum(q_dot_conv(1:N_panels))                   !Receiver convective losses
    q_rad_sum = sum(q_dot_rad(1:N_panels))                     !Receiver radiation losses
  
    !Calculate the flux incident on the surface
    q_dot_inc(:)=P_field(:)*1000.     ![kW]
    q_dot_abs(:) = q_dot_inc(:) - q_dot_loss(:)     !The absorbed flux at each node
    q_inc_sum = sum(q_dot_inc(1:N_panels))	             !The total power incident on the surface
    q_abs_sum = sum(q_dot_abs(1:N_panels))	             !The total power absorbed by the receiver"

    !Calculate the temperature drop across the receiver tube wall.. assume a cylindrical thermal resistance
    T_wall(:) = (T_s(:) + T_panel_ave(:))/2.	             !The temperature at which the conductivity of the wall is evaluated
    k_tube(:) = Conductivity(Material,T_wall(:))	     !The conductivity of the wall
    R_tube_wall(:) = th_tube/(k_tube(:)*H_rec*D_rec*pi**2/2.) !The thermal resistance of the wall
 
    !Calculations for the inside of the tube
    C_p_coolant = specheat(Coolant, T_coolant_prop)*1000.    !Specific heat of the coolant
    LoverD = H_rec/D_inner
    RelRough = (1.5e-6)/D_inner                                !Relative roughness of the tubes. http:www.efunda.com/formulae/fluids/roughness.cfm
    mu_coolant = viscosity(Coolant, T_coolant_prop)            !Absolute viscosity of the coolant
    k_coolant = conductivity(Coolant,  T_coolant_prop)         !Conductivity of the coolant fluid
    rho_coolant = density(Coolant, T_coolant_prop, dble(0))         !Density of the coolant

    u_coolant = m_dot_salt/(n_t*rho_coolant*(D_inner/2.)**2*pi)!Average velocity of the coolant through the receiver tubes.
    Re_inner = rho_coolant*u_coolant*D_inner/mu_coolant        !Reynolds number for internal flow
    Pr_inner = C_p_coolant*mu_coolant/k_coolant                !Prandtl number for internal flow
    call PipeFlow(Re_inner,Pr_inner,LoverD,relRough,Nusselt_t,f)  !The internal convection correlation. Petukhov, Gneilinski
    if(Nusselt_t.le.0) then
        m_dot_salt_tot=0 ; eta_therm=0 ; W_dot_pump=0 
        q_conv_sum=0 ; q_rad_sum=0 ; T_s=0 ;  Q_thermal=0
        return
    endif
    h_inner = Nusselt_t*k_coolant/D_inner                      !Convection coefficient between the inner tube wall and the coolant
    R_conv_inner = 1./(h_inner*pi*D_inner/2.*H_rec*n_t)        !The thermal resistance associated with this value

    !Set up numerical flow grid
    do j=1,N_panels
        if(Flow_pattern(j).lt.1) then    !Declare equation duplicities
            T_panel_in(j) = T_salt_cold
        else
            T_panel_in(j) = T_panel_out(Flow_pattern(j))       !The panel inlet temp is equal to the panel outlet temp from the previous panel, according to the flow diagram
        endif
    enddo
    T_panel_outX(:) = T_panel_in(:) + q_dot_abs(:)/(m_dot_salt*c_p_coolant) !The energy balance for each node
    T_sX(:)=T_panel_ave(:)+q_dot_abs(:)*(R_conv_inner+R_tube_wall(:)) !Calculate the surface temperature based on the absorbed heat
    if(q_inc_sum.gt.0.) then
        eta_therm = q_abs_sum/q_inc_sum
    else
        eta_therm = 0.
    endif
        
    !Final calculations
    m_dot_saltX = q_abs_sum/(nlines*C_p_coolant*(T_salt_hot - T_salt_cold))
    !Do a check to make sure the mass flow rate is reasonable
    if(m_dot_saltX.lt.1.e-5) then  !Modified 3-30-09, MJW
        m_dot_salt_tot=0 ; eta_therm=0 ; W_dot_pump=0  
        q_conv_sum=0 ; q_rad_sum=0 ; T_s=0 ;  Q_thermal=0
        return
    endif
enddo

!Now we can calculate some of the parasitics associated with pumping the coolant fluid
!Calculating the pressure drop across the receiver panels
m_dot_salt_tot = m_dot_salt*nlines
m_dot_tube = m_dot_salt/dble(n_t)                  !The mass flow through each individual tube
L_e_45 = 16.                                       !The equivalent length produced by the bends in the tubes.
											 !'Intro to Fluid Mechanics, Fox, et al.'
L_e_90 = 30.

DELTAP_tube = rho_coolant*(f*H_rec/D_inner*u_coolant**2/2.)  !Pressure drop across the tube, straight length
DELTAP_45 = rho_coolant*(f*L_e_45*u_coolant**2/2.) !Pressure drop across 45 degree bends
DELTAP_90 = rho_coolant*(f*L_e_90*u_coolant**2/2.) !Pressure drop across 90 degree bends
DELTAP = DELTAP_tube + 2*DELTAP_45 + 4*DELTAP_90   !Total pressure drop across the tube with (4) 90-deg bends, (2) 45-deg bends
DELTAP_THT = rho_coolant*THT*grav                  !The pressure drop from pumping up to the receiver
DELTAP_net = DELTAP*N_panels/dble(nlines)+DELTAP_THT !The net pressure drop across the receiver panels
Pres_D = DELTAP_net*1.e-6 ![MPa]
W_dot_pump = DELTAP_net*u_coolant*pi*D_inner**2/4.*n_t/eta_pump !The work done by the pump to move the coolant through the receiver
!Calculate the thermal output of the tower. This includes all thermal losses.
Q_thermal = m_dot_salt_tot*C_p_coolant*(T_salt_hot - T_salt_cold)
!Allow a thermal loss factor, supplied by the user, to tweak the thermal losses
Q_thermal = q_inc_sum-(q_inc_sum-Q_thermal) 
q_conv_sum = q_conv_sum 
q_rad_sum = q_rad_sum 

end subroutine

!#######################################################################################################################



subroutine flowPatterns(N_panels,flowtype,Flow_pattern,salt_out,nlines)
implicit none

!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! This subroutine takes the number of panels, the requested flow type, and 
! returns the corresponding flow pattern and panel numbers where the salt 
! exits.
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
integer,intent(in)::N_panels
integer,dimension(N_panels)::Flow_pattern
integer,dimension(2)::salt_out
integer::nlines,flowtype,i

if((flowtype<=4).and.(dble(floor(N_panels/2.)).lt.(dble(N_panels)/2.))) then
    call messages(-1,'The number of panels for this flow configuration must be divisible by 2',"Fatal",0,222)
endif

select case(flowtype)
case(1)
    !This flow pattern begins at the northmost 2 panels, splits into 2 flows, and crosses over
    !  at the quarter position, exiting in 2 flows on the southmost 2 panels. This is the flow
    !  configuration that was used for SOLAR II
    nlines = 2
    salt_out(1)= 1
    salt_out(2)= N_panels
    Flow_Pattern = (/(i+1,i=1,(floor(N_panels/4.)-1)),(N_panels-floor(N_panels/4.)),(i+1,i=(floor(N_panels/4.)+1),&
                     (N_panels/2)-1),0,0,(i,i=N_panels/2+1,(N_panels-floor(N_panels/4.)-1)),(floor(N_panels/4.)+1),&
                     (i-1,i=(N_panels-floor(N_panels/4.)+2),N_panels),(0,i=1,100-N_panels)/)
case(2)
    !This flow pattern is the same as flow pattern #1, but in reverse. The salt enters
    !  on the 2 southmost panels, crosses over, and exits on the 2 northmost panels.
    nlines = 2
    salt_out(1)=N_panels/2
    salt_out(2)=N_panels/2+1
    Flow_Pattern = (/0,(i-1,i=2,floor(N_panels/4.)),(N_panels+1-floor(N_panels/4.)),(i-1,i=floor(N_panels/4.)+2,N_panels/2),&
                    (i+1,i=(N_panels/2+1),(N_panels-floor(N_panels/4.)-1)),(floor(N_panels/4.)),(i+1,i=(N_panels-floor(N_panels/4.)+1),&
                    (N_panels-1)),0,(0,i=1,100-N_panels)/)
case(3)
!This flow pattern has 2 separate flows that enter in 2 of the northmost panels
!  and flow around (without crossing over), exiting at the 2 southmost panels
    nlines=2
    salt_out(1)=1
    salt_out(2)=N_panels
    Flow_Pattern = (/(i+1,i=1,N_panels/2-1),0,0,(i-1,i=(N_panels/2+2),N_panels),(0,i=1,100-N_panels)/)
case(4)
!This flow pattern has 2 separate flows that enter in 2 of the southmost panels
!  and flow around (without crossing over), exiting at the 2 northmost panels
    nlines=2
    salt_out(1)=N_panels/2
    salt_out(2)=N_panels/2+1
    Flow_Pattern = (/0,(i-1,i=2,N_panels/2),(i+1,i=N_panels/2+1,N_panels-1),0,(0,i=1,100-N_panels)/)
case(5)
    !This flow type enters on a panel at the southmost side of the receiver,
    !  travels completely around the receiver in a clockwise direction,
    !  and exits again on the south side
    nlines=1
    salt_out(1)=N_panels
    Flow_Pattern = (/0,(i-1,i=2,N_panels),(0,i=1,100-N_panels)/)
case(6)
    !This flow type enters on a panel at the southmost side of the receiver,
    !  travels completely around the receiver in a counter-clockwise direction,
    !  and exits again on the south side
    nlines=1
    salt_out(1)=1
    Flow_Pattern = (/(i+1,i=1,N_panels-1),0,(0,i=1,100-N_panels)/)
case(7)
    !This flow type enters on a panel at the northmost side of the receiver,
    !  travels completely around the receiver in a clockwise direction,
    !  and exits again on the north side
    nlines=1
    salt_out(1)=floor(N_panels/2.)
    Flow_Pattern = (/N_panels,(i-1,i=2,floor(N_panels/2.)),0,(i-1,i=floor(N_panels/2.)+2,N_panels),(0,i=1,100-N_panels)/)
case(8)
    !This flow type enters on a panel at the northmost side of the receiver,
    !  travels completely around the receiver in a counter-clockwise direction,
    !  and exits again on the north side
    nlines=1
    salt_out(1)=floor(N_panels/2.)+1
    Flow_Pattern = (/(i+1,i=1,floor(N_panels/2.)-1),0,(i+1,i=floor(N_panels/2.)+1,N_panels-1),1,(0,i=1,100-N_panels)/)
end select

end subroutine

!#######################################################################################################################

subroutine fluxinterp(LU_flux,zen_in,azi_in,array,dayp1,hourp1,info)

!**********************************************************************
! This subroutine reads the 1-D fluxmap output by the PTGEN.exe       *
! program and locates the flux map closest to the requested azimuth/  *
! zenith angle, returning it in an array.  For justification of this  *
! selection method, see the PTGEN portion of the Plant Sizing and     *
! Optimization chapter.                                               *
! Inputs and outputs for this model include:                          *
! INPUTS:                                                             *
!   LU_flux :: The logical unit to be used for the fluxmap file       *
!   zen_in  :: The solar zenith angle                                 *
!   azi_in  :: The solar azimuth angle                                *
! OUTPUTS:                                                            *
!   array   :: a 12x1 array of the flux distrubution around the       *
!              circumference of the receiver                          *  
!   dayp1,hourp1 :: The day and hour associated with the chosen flux  *
!                   map                                               *
!**********************************************************************

implicit none
integer(4),intent(in)::LU_flux,INFO(15)
real(8),intent(in)::zen_in,azi_in

real(8),intent(out)::array(12) 
real(8)::xdist,hold,array_all(96,12)
real(8),allocatable::azimuth(:),zenith(:), day(:), time(:)
integer(4)::mylen,i,j,ios,numrec,p1,dayp1,hourp1
logical::is_there
character::checkname*300, test*300

!Save the variables that are used in future runs
save::numrec, azimuth, zenith, day, time, array_all

if(.not.allocated(azimuth)) then
    allocate(azimuth(96),zenith(96),day(96),time(96))
    azimuth=0.d0; zenith=0.d0; day=0.d0; time=0.d0
    
    ios=0 !Initialize
    !Load in the array of azimuth and zenith values from fluxmap.csv
    inquire(unit=LU_flux,opened=is_there,name=checkname)
    if(.not.is_there) then
        open (unit=LU_flux,file=trim(checkname),status="OLD",iostat=ios,position="rewind")
    endif

    !The data begins on line #5 of the fluxmap.csv file
    rewind(LU_flux)
    do i=1,4
        read (LU_flux,*)
    enddo
    i=0
    do i=1,96 !while(ios.eq.0)
        read (LU_flux,fmt=500,iostat=ios) day(i),time(i),azimuth(i),zenith(i)
        500     format(3(F5.1,1X),F5.1)
    enddo
    numrec=96 !The number of flux arrays included
    
    !Read in all of the data into the array_all variable
    do i=1,3
        read(LU_flux,*)
    enddo
    do i=1,96
        read(LU_flux,fmt='(12(F8.1,1X))') (array_all(i,j),j=1,12)  !array is always 12 wide
        if(i>=96) exit
        do j=1,5
            read(LU_flux,*)
        enddo
    enddo
endif

!Find the closest neighbor to the inputted value
hold = 1000  !set it to a high value so its sure to be reset in the loop
do i=1,numrec
    xdist = sqrt((azimuth(i)-azi_in)**2 + (zenith(i) - zen_in)**2)
    if(xdist.le.hold) then
        hold=xdist
        p1=i
    endif
enddo

!Read in the flux map for the point p1
array(:)=array_all(p1,:)


if(.not.is_there) close(LU_flux)  !If the flux file was closed when the subroutine was called, return to this condition
dayp1=day(p1)
hourp1=time(p1)


end subroutine


!#######################################################################################################################


double precision function Nusselt_FC(ksDin,Re)
implicit none

double precision,intent(in)::ksDin,Re
double precision::Nomval,ValHi,ValLo,ValHi2,ValLo2, Nu_Hi,Nu_Lo,chi,ksD
integer::rerun

!This is the forced convection correlation presented in [Siebers & Kraabel, 1984]
!The value of ks\D determines the interpolation that takes place between these 4 data points.
ksD = ksDin
Nomval = ksD
rerun=0
!Select the bounding conditions

10	continue
!Point 1: ksD=0.0
if(ksD.lt.75e-5) then
    Nusselt_FC=.3+.488*Re**.5*(1.+(Re/282000.)**.625)**.8
    ValHi=75.e-5
    ValLo=0.0
else
    !Point 2: ksD= 75 x 10**(-5)
    if((ksD>=75.e-5).and.(ksD<300.e-5)) then
        ValHi=300.e-5
        ValLo=75.e-5
        if(Re<=7e5) then
            Nusselt_FC=.3+.488*Re**.5*(1.+(Re/282000.)**.625)**.8
        else
            if((Re>7.0e5).and.(Re<2.2e7)) then
                Nusselt_FC=2.57e-3*Re**.98
            else
                Nusselt_FC=.0455*Re**.81
            endif
        endif
    else
        !Point 3: ksD= 300 x 10**(-5)
        if((ksD>=300.0e-5).and.(ksD<900.0e-5)) then
            ValHi=900.0e-5
            ValLo=300.0e-5
            if(Re<=1.8e5) then
                Nusselt_FC=.3+.488*Re**.5*(1.+(Re/282000.)**.625)**.8
            else
                if((Re>1.8e5).and.(Re<4.0e6)) then
                    Nusselt_FC=.0135*Re**.89
                else
                    Nusselt_FC=.0455*Re**.81
                endif
            endif
        else
            !Point 4: ksD = 900 x 10**(-5)
            if(ksD>=900.0e-5) then
                ValHi=900.0e-5
                ValLo=900.0e-5
                if(Re<=1e5) then
                    Nusselt_FC=.3+.488*Re**.5*(1.+(Re/282000.0)**.625)**.8
                else
                    Nusselt_FC=.0455*Re**.81
                endif
            endif
        endif
    endif
endif

if (rerun.eq.1) then
    goto 20
else
    rerun=1
    Nu_Lo=Nusselt_FC
    ksD=ValHi
    ValLo2=ValLo
    ValHi2=ValHi
    goto 10
endif

20 continue
Nu_hi=Nusselt_FC

if (Nomval.ge.900.0e-5) then
    chi=0
else
    chi=(Nomval-ValLo2)/(ValHi2-ValLo2)
endif

Nusselt_FC=Nu_Lo+(Nu_hi-Nu_Lo)*chi

end function

!#######################################################################################################################