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

C************************************************************************
C Object: External Receiver/Tower
C Simulation Studio Model: Type222
C 
C Author: Michael J. Wagner
C Date:	 March 14, 2008 
C Last modified: May 22, 2009
! COPYRIGHT 2009 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
! 14|  nuel                          |  The number of elevation angle points in the flux map           |  none          |  none
! 15|  nazm                          |  The number of azimuth angle points in the flux map             |  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


C (Comments and routine interface generated by TRNSYS Studio)
C************************************************************************

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

C-----------------------------------------------------------------------------------------------------------------------
C    REQUIRED BY THE MULTI-DLL VERSION OF TRNSYS
      !DEC$ATTRIBUTES DLLEXPORT :: TYPE222				!SET THE CORRECT TYPE NUMBER HERE
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C    TRNSYS DECLARATIONS
      IMPLICIT NONE			!REQUIRES THE USER TO DEFINE ALL VARIABLES BEFORE USING THEM

	DOUBLE PRECISION XIN	!THE ARRAY FROM WHICH THE INPUTS TO THIS TYPE WILL BE RETRIEVED
	DOUBLE PRECISION OUT	!THE ARRAY WHICH WILL BE USED TO STORE THE OUTPUTS FROM THIS TYPE
	DOUBLE PRECISION TIME	!THE CURRENT SIMULATION TIME - YOU MAY USE THIS VARIABLE BUT DO NOT SET IT!
	DOUBLE PRECISION PAR	!THE ARRAY FROM WHICH THE PARAMETERS FOR THIS TYPE WILL BE RETRIEVED
	DOUBLE PRECISION STORED !THE STORAGE ARRAY FOR HOLDING VARIABLES FROM TIMESTEP TO TIMESTEP
	DOUBLE PRECISION T		!AN ARRAY CONTAINING THE RESULTS FROM THE DIFFERENTIAL EQUATION SOLVER
	DOUBLE PRECISION DTDT	!AN ARRAY CONTAINING THE DERIVATIVES TO BE PASSED TO THE DIFF.EQ. SOLVER
	INTEGER*4 INFO(15)		!THE INFO ARRAY STORES AND PASSES VALUABLE INFORMATION TO AND FROM THIS TYPE
	INTEGER*4 NP,NI,NOUT,ND	!VARIABLES FOR THE MAXIMUM NUMBER OF PARAMETERS,INPUTS,OUTPUTS AND DERIVATIVES
	INTEGER*4 NPAR,NIN,NDER	!VARIABLES FOR THE CORRECT NUMBER OF PARAMETERS,INPUTS,OUTPUTS AND DERIVATIVES
	INTEGER*4 IUNIT,ITYPE	!THE UNIT NUMBER AND TYPE NUMBER FOR THIS COMPONENT
	INTEGER*4 ICNTRL		!AN ARRAY FOR HOLDING VALUES OF CONTROL FUNCTIONS WITH THE NEW SOLVER
	INTEGER*4 NSTORED		!THE NUMBER OF VARIABLES THAT WILL BE PASSED INTO AND OUT OF STORAGE
	CHARACTER*3 OCHECK		!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE OUTPUTS
	CHARACTER*3 YCHECK		!AN ARRAY TO BE FILLED WITH THE CORRECT VARIABLE TYPES FOR THE INPUTS
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    USER DECLARATIONS - SET THE MAXIMUM NUMBER OF PARAMETERS (NP), INPUTS (NI),
C    OUTPUTS (NOUT), AND DERIVATIVES (ND) THAT MAY BE SUPPLIED FOR THIS TYPE
      PARAMETER (NP=15,NI=14,NOUT=8,ND=0,NSTORED=0)
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    REQUIRED TRNSYS DIMENSIONS
      DIMENSION XIN(NI),OUT(NOUT),PAR(NP),YCHECK(NI),OCHECK(NOUT),
	1   STORED(NSTORED),T(ND),DTDT(ND)
      INTEGER NITEMS
C-----------------------------------------------------------------------------------------------------------------------
C-----------------------------------------------------------------------------------------------------------------------
C    ADD DECLARATIONS AND DEFINITIONS FOR THE USER-VARIABLES HERE
!	CHARACTER*40 matlist,coolist
!	DIMENSION matlist(1),coolist(2)
	double precision m_dot_salt_tot,eta_therm,W_dot_pump,
     & q_conv_sum,q_rad_sum,T_s,Q_thermal, I_bn_inst,P_amb_kPa
	DIMENSION T_s(50)

C    PARAMETERS
      DOUBLE PRECISION N_panels, D_rec, H_rec, THT, D_out, th_tu,
     &  Material, Coolant, Flow_type, PLAT, LU_flux, epsilon,
     &  hl_ffact, nuel, nuaz
C    INPUTS
      DOUBLE PRECISION 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*200, test*200
      integer ios, i

C-----------------------------------------------------------------------------------------------------------------------
C       READ IN THE VALUES OF THE PARAMETERS IN SEQUENTIAL ORDER
      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)
      nuel=PAR(14)
      nuaz=PAR(15)

C-----------------------------------------------------------------------------------------------------------------------
C    RETRIEVE THE CURRENT VALUES OF THE INPUTS TO THIS MODEL FROM THE XIN ARRAY IN SEQUENTIAL ORDER

      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)

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

C-----------------------------------------------------------------------------------------------------------------------
C    DO ALL THE VERY LAST CALL OF THE SIMULATION MANIPULATIONS HERE
      IF (INFO(8).EQ.-1) THEN
         !Close the flux map file
         inquire(unit=LU_flux,opened=is_there,name=checkname)
         if (is_there) close(LU_flux)
         
	   RETURN 1
	ENDIF
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    PERFORM ANY 'AFTER-ITERATION' MANIPULATIONS THAT ARE REQUIRED HERE
C    e.g. save variables to storage array for the next timestep
      IF (INFO(13).GT.0) THEN
	   NITEMS=0
C	   STORED(1)=... (if NITEMS > 0)
C        CALL setStorageVars(STORED,NITEMS,INFO)
	   RETURN 1
	ENDIF
C
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    DO ALL THE VERY FIRST CALL OF THE SIMULATION MANIPULATIONS HERE
      IF (INFO(7).EQ.-1) THEN

C       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	!STORAGE FOR VERSION 16 HAS BEEN CHANGED				

C       SET THE REQUIRED NUMBER OF INPUTS, PARAMETERS AND DERIVATIVES THAT THE USER SHOULD SUPPLY IN THE INPUT FILE
C       IN SOME CASES, THE NUMBER OF VARIABLES MAY DEPEND ON THE VALUE OF PARAMETERS TO THIS MODEL....
         NIN=NI
	   NPAR=NP
	   NDER=ND


C       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)
      !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~  

     
C       RETURN TO THE CALLING PROGRAM
         RETURN 1

      ENDIF
C-----------------------------------------------------------------------------------------------------------------------

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

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

C       CHECK THE PARAMETERS FOR PROBLEMS AND RETURN FROM THE SUBROUTINE IF AN ERROR IS FOUND
C         IF(...) CALL TYPECK(-4,INFO,0,"BAD PARAMETER #",0)

C       PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL VALUES OF THE OUTPUTS HERE
C		 m_dot_salt_tot
			OUT(1)=0.
C		 eta_therm
			OUT(2)=0.
C		 W_dot_pump
			OUT(3)=0.
C		 q_conv_sum
			OUT(4)=0.
C		 q_rad_sum
			OUT(5)=0.
C		 Q_thermal
			OUT(6)=0.
C		 T_salt_hot
			OUT(7)=565.0
C          Power before thermal losses
              OUT(8)=0.

C       PERFORM ANY REQUIRED CALCULATIONS TO SET THE INITIAL STORAGE VARIABLES HERE
         NITEMS=0
C	   STORED(1)=...

C       PUT THE STORED ARRAY IN THE GLOBAL STORED ARRAY
C         CALL setStorageVars(STORED,NITEMS,INFO)

C       RETURN TO THE CALLING PROGRAM
         RETURN 1

      ENDIF
C-----------------------------------------------------------------------------------------------------------------------

C-----------------------------------------------------------------------------------------------------------------------
C    *** ITS AN ITERATIVE CALL TO THIS COMPONENT ***
C-----------------------------------------------------------------------------------------------------------------------

	    
C-----------------------------------------------------------------------------------------------------------------------
C    RETRIEVE THE VALUES IN THE STORAGE ARRAY FOR THIS ITERATION
C      NITEMS=
C	CALL getStorageVars(STORED,NITEMS,INFO)
C      STORED(1)=
C-----------------------------------------------------------------------------------------------------------------------

C	Do necessary conversions
	I_bn_inst = I_bn/3.6   !This converts from kJ/m^2-hr to kW/m^2
      continue
	call CRS((hour-12.),azimuth,zenith,D_rec,H_rec,THT,D_out,th_tu, !Inputs
     & 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,nuaz,nuel,
     & 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)  !Outputs

      continue


C-----------------------------------------------------------------------------------------------------------------------
C    SET THE STORAGE ARRAY AT THE END OF THIS ITERATION IF NECESSARY
C      NITEMS=
C      STORED(1)=
C	CALL setStorageVars(STORED,NITEMS,INFO)
C-----------------------------------------------------------------------------------------------------------------------


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

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

C-----------------------------------------------------------------------------------------------------------------------
C    EVERYTHING IS DONE - RETURN FROM THIS SUBROUTINE AND MOVE ON
      RETURN 1
      END
C-----------------------------------------------------------------------------------------------------------------------

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


C#######################################################################################################################
	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,
     & nuel,nuaz,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)

	implicit none

	!Define input variables ----------------------------------------------------
	integer,intent(in)::N_panels,flowtype,LU_flux
	INTEGER*4,intent(in):: INFO(15)
	double precision::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
      double precision::V_wind,time,material,coolant
!	character(*),intent(in)::Material,Coolant

	!Define output variables
	double precision,intent(out)::m_dot_salt_tot,eta_therm,q_conv_sum,
     & q_rad_sum, Q_thermal
	double precision,dimension(N_panels),intent(out)::T_s

	!Define local variables
	integer,dimension(N_panels)::Flow_pattern
	integer,dimension(2)::salt_out
	integer::i,j,qq,nlines,daycall,hourcall,qq_max
	integer,dimension(N_panels)::Panel,ceil,flo
	double precision,parameter::pi=3.14159265,grav=9.81
	double precision::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,myconvert,converttemp,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,nuaz,nuel
	double precision,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

	!Initialize sensitive variables
	daycall=0 ; hourcall=0 ; ppos=0
	!---------------------------------------------------------------------------

	!Receiver dimensions, parameters
	D_tube = D_out*myconvert('mm','m')      !Value in meters
	th_tube = th_tu*myconvert('mm','m')     !Thickness of the tube
	D_inner = D_tube - 2.*th_tube	      !Diameter of each receiver tube
	D_in = D_inner*myconvert('m','mm')      !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 = converttemp('C','K',T_db)    !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*myconvert('atm','Pa')     !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

	!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.
     1   ((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
	      return
	    endif
	endif

	!Get flow pattern
	call flowPatterns(N_panels,flowtype,Flow_pattern,salt_out,nlines)

	!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,nuel,array,daycall,hourcall
     &                    ,info)
	    do j=1,int(nuel)
	      flux_in(j)=array(j)/(950.)*I_bn*field_eff  !The weather-adjusted, efficiency-adjusted flux values
	    enddo
	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)=(nuel/dble(N_panels)*(j-1.)+nuel/2./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.nuel) then
	    ceil(j)=1
	  endif
	  Psp_field(j)=ind(j)*(flux_in(ceil(j))-flux_in(flo(j)))+
	1	  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=5.0e-3
	  else
	      tol=1.0e-4
	  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----------
	  do j=1,N_panels
	    T_s(j) = T_sX(j)
	    T_panel_out(j) = T_panel_outX(j)
	    T_panel_in(j) = T_panel_inX(j)
	  enddo
	  m_dot_salt = m_dot_saltX
	  !------------------------------------------------------------------------

	  !Now do the actual calculations
	  do j=1,N_panels
	    T_panel_ave(j) = (T_panel_in(j)+T_panel_out(j))/2.  !The average coolant temperature in each control volume
	    T_film(j) = (T_s(j)+T_amb)/2.
	  enddo
	  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)*myconvert('kJ','J')    !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                           !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 = volexpcoef('Air',T_amb)                             !Volumetric expansion coefficient
	  nu_amb = Viscosity(1.d0,T_amb)/Density(1.d0,T_amb,P_atm) !Kinimatic viscosity
	  do j=1,N_panels
	    Gr_nat(j) = grav*beta*(T_s(j)-T_amb)*H_rec**3/nu_amb**2      !Grashof Number at ambient conditions
	    Nusselt_nat(j) = 
	1	.098*Gr_nat(j)**(1./3.)*(T_s(j)/T_amb)**(-.14) !Nusselt number
	    h_nat(j) = Nusselt_nat(j)*k_film/H_rec                           !The natural convection cofficient
	  enddo

	  !Mixed convection
	  m=3.2
	  do j=1,N_panels
	    h_mixed(j) = (h_for**m+h_nat(j)**m)**(1./m)              !The mixed convection relationship
	    q_dot_conv(j) = h_mixed(j)*A_node*(T_s(j) - T_film(j))   !Convection losses per node
	  enddo
	  !Radiation from the receiver
	    !Calculate the radiation node by node
	  do j=1,N_panels
	    h_rad_amb(j) = 
	1	sigma*epsilon*(T_s(j)**2+T_amb**2)*(T_s(j)+T_amb)  !The radiation coefficient for amb
	    h_rad_sky(j) = 
	1	sigma*epsilon*(T_s(j)**2+T_sky**2)*(T_s(j)+T_sky)  !The radiation coef. for sky
	    q_dot_amb(j) = .5*h_rad_amb(j)*A_node*(T_s(j)-T_amb)  !amb losses per node
	    q_dot_sky(j) = .5*h_rad_sky(j)*A_node*(T_s(j) - T_sky)  !sky losses per node
	  enddo

	  !Calculate the losses from the surface"
	  do j=1,N_panels
	    q_dot_rad(j) = q_dot_amb(j)+q_dot_sky(j)		     !Total radiation losses per node
	    q_dot_loss(j) = q_dot_rad(j)+q_dot_conv(j)		     !Total overall losses per node
	  enddo
	  !q_loss_sum = sum(q_dot_loss(1:N_panels))                   !Receiver total losses
	  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
	  do j=1,N_panels
	    q_dot_inc(j)=P_field(j)*myconvert('kW','W')
	    q_dot_abs(j) = q_dot_inc(j) - q_dot_loss(j)     !The absorbed flux at each node
	  enddo
	  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
	  do j=1,N_panels
	    T_wall(j) = (T_s(j) + T_panel_ave(j))/2.	             !The temperature at which the conductivity of the wall is evaluated
	    k_tube(j) = Conductivity(Material,T_wall(j))	     !The conductivity of the wall
	    R_tube_wall(j) = th_tube/(k_tube(j)*H_rec*D_rec*pi**2/2.) !The thermal resistance of the wall
	  enddo

	  !Calculations for the inside of the tube
	  C_p_coolant = specheat(Coolant, T_coolant_prop)
	1  *myconvert('kJ','J') !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
	    T_panel_outX(j) = T_panel_in(j) + q_dot_abs(j)/
	1	(m_dot_salt*c_p_coolant) !The energy balance for each node
	    T_sX(j)=T_panel_ave(j)+q_dot_abs(j)*
	1	(R_conv_inner+R_tube_wall(j)) !Calculate the surface temperature based on the absorbed heat
	  enddo
	  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*
	1  (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*myconvert('Pa','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)*hl_ffact 
      q_conv_sum = q_conv_sum*hl_ffact
      q_rad_sum = q_rad_sum*hl_ffact
	end subroutine






C#######################################################################################################################
	double precision function myconvert(units1,units2)
	implicit none
	!This function takes as an input two text strings and outputs the 
	!factor by which the two units differ (a real number)
	character(*),intent(in)::units1,units2
	character(len=10)::units(2)
	integer::i
	double precision::factor(2),diff(2)
	factor=1 ; diff=0  !initialize variables to default values
	units(1)=units1
	units(2)=units2

	!get unit in terms of its base units
	do i=1,2
	!Length --> m
	  if(trim(units(i)).eq.'m') then
	    factor(i)=1 ; diff(i)=0
	  elseif(trim(units(i)).eq.'cm') then
	    factor(i)=1./100. ; diff(i)=0
	  elseif(trim(units(i)).eq.'mm') then
	    factor(i)=1./1000. ;diff(i)=0
	  elseif(trim(units(i)).eq.'in') then
	    factor(i)=.0254 ; diff(i)=0
	  endif
	!velocity --> m/s
	  if(trim(units(i)).eq.'mph') then
	    factor(i)=.44704 ; diff(i)=0
	  elseif(trim(units(i)).eq.'m/s') then
	    factor(i)=1 ; diff(i)=0
	  elseif(trim(units(i)).eq.'km/s') then
	    factor(i)=1./1000. ; diff(i)=0
	  elseif(trim(units(i)).eq.'km/hr') then
	    factor(i)=.27777777778 ; diff(i)=0
	  endif
	!time --> s
	  if(trim(units(i)).eq.'hr') then
	    factor(i)=3600. ; diff(i)=0
	  elseif(trim(units(i)).eq.'min') then
	    factor(i)=60. ; diff(i)=0
	  elseif(trim(units(i)).eq.'day') then
	    factor(i)= 86400. ; diff(i)=0
	  elseif(trim(units(i)).eq.'year') then
	    factor(i)=31556926. ; diff(i)=0
	  endif
	!energy --> J
	  if(trim(units(i)).eq.'kJ') then
	    factor(i)=1000. ; diff(i)=0
	  elseif(trim(units(i)).eq.'J') then
	    factor(i)=1. ; diff(i)=0
	  elseif(trim(units(i)).eq.'MJ') then
	    factor(i)=1.0e6 ; diff(i)=0
	  elseif(trim(units(i)).eq.'GJ') then
	    factor(i)=1.0e9 ; diff(i)=0
	  endif
	!energy rate --> W (J/s)
	  if(trim(units(i)).eq.'kW') then
	    factor(i)=1000. ; diff(i)=0
	  elseif(trim(units(i)).eq.'W') then
	    factor(i)=1. ; diff(i)=0
	  elseif(trim(units(i)).eq.'MW') then
	    factor(i)=1.0e6 ; diff(i)=0
	  elseif(trim(units(i)).eq.'GW') then
	    factor(i)=1.0e9 ; diff(i)=0
	  elseif(trim(units(i)).eq.'hp') then     !horsepower
	    factor(i)=745.699872 ; diff(i)=0
	  endif
	!pressure --> Pa
	  if(trim(units(i)).eq.'kPa') then
	    factor(i)=1000. ; diff(i)=0
	  elseif(trim(units(i)).eq.'Pa') then
	    factor(i)=1. ; diff(i)=0
	  elseif(trim(units(i)).eq.'MPa') then
	    factor(i)=1.0e6 ; diff(i)=0
	  elseif(trim(units(i)).eq.'GPa') then
	    factor(i)=1.0e9 ; diff(i)=0
	  elseif(trim(units(i)).eq.'psi') then
	    factor(i)=6894.75729 ; diff(i)=0
	  elseif(trim(units(i)).eq.'millibar') then
	    factor(i)=100. ; diff(i)=0
	  elseif(trim(units(i)).eq.'mbar') then     !also millibar
	    factor(i)=100. ; diff(i)=0
	  elseif(trim(units(i)).eq.'atm') then
	    factor(i)=101325. ; diff(i)=0
	  endif
	!viscosity --> Pa-s (kg/m-s)
	  if(trim(units(i)).eq.'kPa-s') then
	    factor(i)=1000. ; diff(i)=0
	  elseif(trim(units(i)).eq.'Pa-s') then
	    factor(i)=1. ; diff(i)=0
	  elseif(trim(units(i)).eq.'cP') then
	    factor(i)=1./1000. ; diff(i)=0
	  endif

	enddo

	myconvert=factor(1)/factor(2)
	end function
 
	double precision function converttemp(units1,units2,value)
	character(len=*),intent(in)::units1,units2
	character(len=10)::units(2)
	double precision,intent(in)::value
	double precision::factor(2),diff(2)
	integer::i
	factor=1 ; diff=0
	units(1)=units1
	units(2)=units2

 
	!temperature --> C
	do i=1,2
	  if(trim(units(i)).eq.'K') then
	    factor(i)=1 ; diff(i)=-273.15
	  elseif(trim(units(i)).eq.'C') then
	    factor(i)=1 ; diff(i)=0
	  elseif(trim(units(i)).eq.'F') then
	    factor(i)=5./9. ; diff(i)=-32.
	  elseif(trim(units(i)).eq.'R') then
	    factor(i)=5./9. ; diff(i)=-491.69
	  endif
	enddo
	converttemp= (value + diff(1))*factor(1)/factor(2) - diff(2)
 
	end function
C#######################################################################################################################


	subroutine flowPatterns(N_panels,flowtype,Flow_pattern,salt_out,
     1 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


	 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
	   if(dble(floor(N_panels/2.)).lt.(dble(N_panels)/2.)) then
	     print *,'The number of panels for this flow configuration'//
	1	 ' must be divisible by 2'
	     print *,'Hit [Enter]'
	     read (5,*)
	   else
	     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)/)
	   endif
	 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.
	   if(dble(floor(N_panels/2.)).lt.(dble(N_panels)/2.)) then
	     print *,'The number of panels for this flow configuration'//
	1	 ' must be divisible by 2'
	     print *,'Hit [Enter]'
	     read (5,*)
	   else
	   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)/)
	   endif
	 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
	   if(dble(floor(N_panels/2.)).lt.(dble(N_panels)/2.)) then
	     print *,'The number of panels for this flow configuration'//
	1	 ' must be divisible by 2'
	     print *,'Hit [Enter]'
	     read (5,*)
	   else
	   nlines=2
	   salt_out(1)=1
	   salt_out(2)=N_panels
	   Flow_Pattern = (/(i+1,i=1,N_panels/2-1),0,0,
     1   (i-1,i=(N_panels/2+2),N_panels),(0,i=1,100-N_panels)/)
	   endif
	 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
	   if(dble(floor(N_panels/2.)).lt.(dble(N_panels)/2.)) then
	     print *,'The number of panels for this flow configuration'//
	1	 ' must be divisible by 2'
	     print *,'Hit [Enter]'
	     read (5,*)
	   else
	   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),
	1   (i+1,i=N_panels/2+1,N_panels-1),0,(0,i=1,100-N_panels)/)
	   endif
	 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.)),
	1   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),
	1   0,(i+1,i=floor(N_panels/2.)+1,N_panels-1),1,
	2   (0,i=1,100-N_panels)/)
	 end select

	end subroutine

C#######################################################################################################################


	subroutine fluxinterp(LU_flux,zen_in,azi_in,nuel,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,intent(in)::LU_flux
	double precision,intent(in)::zen_in,azi_in,nuel
	INTEGER*4,intent(in):: INFO(15)
	integer::mylen
	double precision,intent(out)::array(int(nuel)) 
	double precision::azimuth(100),zenith(100),xdist,save,
	1   day(100),time(100)
	   !All of the flux maps are 1x12
	integer::i,j,ios,numrec,p1,dayp1,hourp1
      logical::is_there
      character::checkname*100, test*200

	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(is_there.eq..FALSE.) 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 while(ios.eq.0)
	  i=i+1  !counting variable
	  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=i-1 !The number of flux arrays included


	!Find the closest neighbor to the inputted value
	save = 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.save) then
	    save=xdist
	    p1=i
	  endif
	enddo

	!Read in the flux map for the point p1
	!The actual flux map records begin at line 101
	!The flux array data begins on the 13th line of the record
	!Each record is 6 lines long, we want the p1'th record
	rewind(LU_flux)
	do i=1,(numrec+7+6*(p1-1))
	  read(LU_flux,*)
	enddo
	do j=1,int(nuel)
	  read(LU_flux,advance="NO",fmt='(F8.1,1X)') array(j)
      enddo
      if(is_there.eq..FALSE.) 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

!        integer function mylen(text)
!        implicit none
!
!        !**********************************************************************
!        ! This function determines the length of a text string that includes  *
!        ! spaces. The return value for this function is an integer value.     *
!        !**********************************************************************
!
!        character(len=*),intent(in)::text
!        integer::length,i
!
!        do i = 1,len(text)
!         IF(.not.((text(i:i).eq.'').or.(text(i:i).eq.' '))) THEN
!            length = i
!         ENDIF
!        enddo
!        mylen=length
!        end function
C#######################################################################################################################


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

C#######################################################################################################################

	subroutine PipeFlow(Re, Pr, LoverD,relRough, Nusselt, f)
	implicit none

	!*********************************************************************
	!* PipeFlow_turbulent:                                               *
	!* This procedure calculates the average Nusselt number and friction *
	!* factor for turbulent flow in a pipe given Reynolds number (Re),   *
	!* Prandtl number (Pr), the pipe length diameter ratio (LoverD) and  *
	!* the relative roughness}                                           *
	!*********************************************************************

	double precision,intent(in)::Re,Pr,LoverD,relRough
	double precision,intent(out)::Nusselt,f
	double precision::f_fd,Nusselt_L, Gz, Gm, Nusselt_T, Nusselt_H,fR,X

      !Correlation for laminar flow.. Note that no transitional effects are considered
	if (Re.lt.2300.) then 
	    !This procedure calculates the average Nusselt number and friction factor for laminar flow in a pipe 
	    !..given Reynolds number (Re), Prandtl number (Pr), the pipe length diameter ratio (LoverD) 
	    !..and the relative roughness}
	    Gz=Re*Pr/LoverD	
	    x=LoverD/Re
	    fR=3.44/sqrt(x)+(1.25/(4*x)+16-3.44/sqrt(x))/(1+0.00021*x**(-2))
	    f=4.*fR/Re
	    !{f$='Shah' {Shah, R.K.  and London, A.L. "Laminar Flow Forced Convection in Ducts", 
	    !..Academic Press, 1978 ,Eqn 192, p98}}
	    Gm=Gz**(1./3.)
   	    Nusselt_T=3.66+((0.049+0.02/Pr)*Gz**1.12)/(1+0.065*Gz**0.7)
	    Nusselt_H=4.36+((0.1156 +0.08569 /Pr**0.4)*Gz)/(1+0.1158*Gz**0.6) 
	    !{Nusselt$='Nellis and Klein fit to Hornbeck'  {Shah, R.K.  and London, A.L. "Laminar Flow Forced Convection in Ducts",
	    !..Academic Press, 1978 ,Tables  20 and 22}}
	    Nusselt = Nusselt_T  !Constant temperature Nu is better approximation
	    
	else  !Correlation for turbulent flow
     	    f_fd = (0.79*log(Re)-1.64)**(-2) !Petukhov, B.S., in Advances in Heat Transfer, Vol. 6, Irvine and Hartnett, Academic Press, 1970
	    Nusselt_L= ((f_fd/8.)*(Re-1000)*Pr)/(1.+12.7*sqrt(f_fd/8.)*
     1                (Pr **(2/3.)-1.)) !Gnielinski, V.,, Int. Chem. Eng., 16, 359, 1976

	    if (relRough.gt.1e-5) then

	      f=8.*((8./Re)**12+((2.457*log(1./((7./Re)**0.9+0.27*
	1           (RelRough))))**16+(37530./Re)**16)**(-1.5))**(1./12.)

	      f_fd=(-2.*log10(2*relrough/7.4-5.02*
	1           log10(2*relrough/7.4+13/Re)/Re))**(-2)
    	
	      Nusselt_L= ((f_fd/8.)*(Re-1000.)*Pr)/(1.+12.7*sqrt(f_fd/8.)*
     1                 (Pr **(2./3.)-1.)) !Gnielinski, V.,, Int. Chem. Eng., 16, 359, 1976}
	    endif
	    f=f_fd*(1.+(1./LoverD)**0.7) !account for developing flow
	    Nusselt= Nusselt_L*(1.+(1./LoverD)**0.7)  !account for developing flow
      endif
      
99	continue

	end

C#######################################################################################################################



	!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

	double precision function volexpcoef(substance,T)
	implicit none
	!This function outputs in units of [J/kg-K]
	double precision,intent(in)::T
	character(*),intent(in)::substance

	if((trim(substance).eq.'air').or.(trim(substance).eq.'Air')) then
	  volexpcoef = 1./T
	else
	  print *,"The substance that was called is not available."
	  print *,
	1  "Please make sure the substance is supported.  Press [Enter]..."
	  read (5,*)
	  stop
	endif

	end function

	!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

	subroutine substances
	implicit none
	integer,parameter::N=14
	integer::i
	character(len=40)::subst(N)

	data subst/
     &'air',
     &'Air',
     &'Salt (68% KCl, 32% MgCl2)',
     &'Salt (8% NaF, 92% NaBF4)',
     &'Salt (25% KF, 75% KBF4)',
     &'Salt (31% RbF, 69% RbBF4)',
     &'Salt (46.5% LiF, 11.5%NaF, 42%KF)',
     &'Salt (49% LiF, 29% NaF, 29% ZrF4)',
     &'Salt (58% KF, 42% ZrF4)',
     &'Salt (58% LiCl, 42% RbCl)',
     &'Salt (58% NaCl, 42% MgCl2)',
     &'Salt (59.5% LiCl, 40.5% KCl)',
     &'Salt (59.5% NaF, 40.5% ZrF4)',
     &'Salt (60% NaNO3, 40% KNO3)'/

	do i=1,N
	  print *,subst(i)
	enddo

	end subroutine



	double precision function skytemp(T_amb,T_dp,hour)
	implicit none

	!**********************************************************************
	!  This function uses the correlation for Sky Temperature             *
	!  that was provided in Duffie & Beckman (2006), and was              *
	!  also implemented in EES.                                           *
	!                                                                     *
	!  This function takes as inputs:                                     *
	!    - T_amb -> ambient air temperature, dry bulb [K]                 *
	!    - T_dp  -> the ambient dewpoint temperature [K]                  *
	!    - hour  -> the hour, in solar time, with noon as 0               *
	!  The function outputs:                                              *
	!    - skytemp -> the effective temperature of the sky, in degrees [K]*
	!                                                                     *
	!  Note that the "converttemp" subroutine is required in this function*
	!**********************************************************************

	double precision,intent(in)::T_amb,T_dp,hour
	double precision::T_sky,converttemp,T_ambC,T_dpC,T_skyC,time
	double precision,parameter::pi=3.1415

	!express "time" in terms of an angle  [rad]
	time = (180.-hour*15.)*pi/180.

	!The inputs are in terms of degrees K, but the equation is in terms of degrees C
	T_ambC = converttemp('K','C',T_amb)
	T_dpC = converttemp('K','C',T_dp)

	!The sky temperature relationship
	T_skyC = (T_ambC*(.711+.0056*T_dpC+.000073*T_dpC**2+.013*
	1  cos(time))**.25)

	!Convert the temperature back to degrees K for output
	skytemp = converttemp('C','K',T_skyC)


	end function