module global
    INTEGER, PARAMETER :: iknd=SELECTED_INT_KIND(8)  ! nobody does this
    INTEGER, PARAMETER :: DP=SELECTED_REAL_KIND(14)  ! this is double precision
    INTEGER, PARAMETER :: SP=SELECTED_REAL_KIND(6)  ! this is single precision   
    integer :: strat,downstream,furnace,ndays,nyears 
    integer :: opt,opt_method,iset,sens_map,sens_var
    integer :: yr,yr_opt !number of years of simulation
    integer :: site,var_vect(11),conv,conv_known
    integer :: opt_starter,opt_main,PC_index,print_stuff,guess,nyears_gen,select_file
    real(DP) :: PCon_sum(5000,1),xpass(10),Qset,Qset_ratio
    integer,parameter :: ilen = 8760 !length of input file, i.e. hours
    integer,parameter :: ilen2 = 175200
    !175200 !length of input file, i.e. hours
    integer,parameter :: day_len = 24 !length of day in hours
    real(DP) :: tground,comparison_results(2,20),energy_reduction
    !real(DP) :: qcoola(ilen),qheata(ilen),tdba(ilen),twba(ilen),ti(ilen)
    real(DP) :: qcoola(ilen2),qheata(ilen2),tdba(ilen2),twba(ilen2),ti(ilen2)
    real(DP) :: delt,delt_orig,L_hyb_m,C_CCCT,q_lh,q_lc,Qset_ash
    real(DP),allocatable :: qcool(:),qheat(:),tdb(:),twb(:),qcool_opt(:), &
                            qheat_opt(:),qcool_orig(:),qheat_orig(:),tdb_opt(:), &
                            twb_opt(:),tdb_orig(:),twb_orig(:)    
    logical :: first_thru
    character(len = 300) :: dir_loc
    real(DP) :: month_hours(12,2),day(365,2)
    integer :: print_vector(28)
    character*15 var_vect_text(8)
    character*15 var_vect_pl(8)
    real(DP) :: first_cost,e_source,g_source
    real(DP) :: total_fuel_cost
    
    
    !DP   THIS INCLUDE FILE MAKES IT EASY TO ADJUST THE SIZE OF THE ARRAYS
!DP   USED IN TRNVDST.FOR.
!DP   WRITTEN BY DANIEL PAHUD, OCT. 1996.

!DP---IDIM    = Dimension of the radial array for the global heat process.
	  integer, PARAMETER :: IDIM=26

!DP---JDIM    = Dimension of the vertical array for the global heat process.
      integer, PARAMETER :: JDIM=38

!DP---IDIMR   = Dimension of the radial array for the local heat process.
      integer, PARAMETER :: IDIMR=20

!DP---IDIMLO  = Dimension of the subregion array for the local solutions.
!DP---The maximum number of local solutions must not exceed (IDIMLO-1).
      integer, PARAMETER :: IDIMLO=121

!DP---MAXRLO  = Dimension of the radial subregion array.
      integer, PARAMETER :: MAXRLO=21

!DP---MRLOM1  = Maximum number of radial subregions.
!DP---MRLOM1 is also the maximum number of boreholes connected in series NSERIE.
!DP---Modify MAXRLO in case of necessity.
      integer, PARAMETER :: MRLOM1=MAXRLO-1

!DP---MAXZLO  = Dimension of the vertical subregion array.
      integer, PARAMETER :: MAXZLO=21

!DP---MZLOM1  = Maximum number of vertical subregions.
!DP---Modify MAXZLO in case of necessity.
      integer, PARAMETER :: MZLOM1=MAXZLO-1

!DP---NLAY    = Maximum number of ground layers.
      integer, PARAMETER :: NLAY=10

!DP---NSTLAY  = Maximum number of ground layers within the storage volume.
      integer, PARAMETER :: NSTLAY=5

!DP---MAXTR1  = Maximum number of global timesteps DT in TR1 (local transient).
      integer, PARAMETER :: MAXTR1=1500

!DP---MAXSPE  = Maximum number of specified output ground temperatures.
      integer, PARAMETER :: MAXSPE=10
      
    !chs1
    real(DP) :: AREAR(IDIM,JDIM)
    real(DP) :: AREAZ(IDIM,JDIM),CIN(IDIM,JDIM),FR(IDIM,JDIM),FZ(IDIM,JDIM)
    real(DP) :: GFR(IDIM,JDIM),GFZ(IDIM,JDIM),R(IDIM),RISOR(IDIM,JDIM)
    real(DP) :: VOL(IDIM,JDIM),RISOZ(IDIM,JDIM),RLAM(IDIM,JDIM),RM(IDIM)
    real(DP) :: T(IDIM,JDIM),Z(JDIM),ZM(JDIM),TSF(IDIM,JDIM)      
    integer :: IMAX,IMAX1,JMAX,JMAX1
    !chs2
    real(DP) :: TCMIN,TCMAX,TCPH,DRBREF,DZBREF,DTBREF,JREF
    !chs3
    real(DP) :: TAIRH,TAIR,TAIRM,ATAIR,TAPH,TSTMAX,ROW,TAHT,TAT
    !chs4
    real(DP) :: TVOHS,EVOHST,QVOHS
    integer :: IVOHSN,IVOHSX,JVOHSN,JVOHSX
    real(DP) :: GVOHST,TDEM,QDEM,LOMODE,QDEM1,VOLST,PLST,CAPSTO,EPLUS,ENEG
    real(DP) :: THPLUS,THNEG,T0,CST,DT,DTOKAY,TIME,LCOORD,TSF1(NSTLAY)
    real(DP) :: GVOHS(IDIM,JDIM)
    !chs5
    real(DP) :: IND,IMAXR1,ITR1,ITR1M1,IRAFLG,IRAMET,LPRT
    real(DP) :: LFPRT,R0,R1,TR1,DTRMAX,TIMGLO,DTDEM1,TSTAVE,RMPIPE,RAPIPE
    real(DP) :: RMPOLD,RMCONS,CW,TMSTOR,ETOP,ESIDE,EBOT,DETOP
    real(DP) :: DESIDE,DEBOT,DEPLUS,DENEG,TR(IDIMLO,IDIMR),RR(IDIMR)
    real(DP) :: RMR(IDIMR),GFRR(IDIMR,NSTLAY),CINR(IDIMR,NSTLAY)
    real(DP) :: ESUM(IDIMLO,MAXTR1),ESUMR(IDIMLO),DTEMP(IDIMLO),RLAMR(NSTLAY)
    real(DP) :: CR(NSTLAY),CRTOT(NSTLAY),TSFH(IDIMR),FRR(IDIMLO,IDIMR),TUTMED(0:MRLOM1)
    !chs6
    real(DP) :: NLOC,NLOC1,IFLAGG,CWQWGM,QWDEM1,QWDEM,CWQW
    real(DP) :: TIN,TOUT,QWGMS,QWGM,INVOL(IDIM,JDIM),VOLOC(IDIMLO)
    real(DP) :: PLSTL(IDIMLO),BETA(IDIMLO),SFBETA(IDIMLO),TGM(IDIMLO)
    real(DP) :: D1(NSTLAY)
    !chs7
    real(DP) :: IRR,JTZ,JBZ,ISO,ILAY,IRR1,JBZ1,IBAL
    integer :: NLAYST
    real(DP) :: RADIUS,HEIGHT,DEPTH,FAC,FACL,DRMIN,DRMAX,DRRMIN,ZDMIN,RIST
    real(DP) :: RISS,RISB,FRISO,DISPIP,VOLUME,DEEP,RBORE,THISO,RISLAM,TIMO3
    real(DP) :: TSTARG,TGRAD,NZLOCS(NSTLAY),JTZS(NSTLAY),JBZS(NSTLAY)
    real(DP) :: THLST(NSTLAY),RLAML(NLAY),CL(NLAY),D(NLAY),THL(NLAY)
    !chs8
    real(DP) :: TOLD(IDIM,JDIM),TROLD(IDIMLO,IDIMR),OLDE(IDIMLO,MAXTR1),OLDER(IDIMLO)
    !chs9
    real(DP) :: NSPEC,IBLOSS,IGSPEC(MAXSPE),JGSPEC(MAXSPE),ILSPEC(MAXSPE)
    !chs10
    real(DP) :: ITOPFL,JZLOC(MAXZLO),ILAYST(MAXZLO)
    integer :: nrloc,nzloc,nserie,IRLOC(MAXRLO)
    !chs11
    real(DP) :: NPIPES,DPIPES,Rb,Rpi,Rpo,L,Lb,Lpipe,Lgap,Gap,Rful0,Vref,Tref,Reynol
    integer :: Tflg,Pflg
    !chs12
    real(DP) :: LURRB,NXRB(2),XRB(2),YRB(5)
    !chsiou
    real(DP) :: LUR_DST,LUR1,LUW1
    real(DP) :: TUTME1(0:MRLOM1),TUTME2(0:MRLOM1)
    integer :: imut,ipre,IMAXR,LUW_DST
    parameter pi = 3.1415926535
    integer :: double
    logical :: first_run_opt
end module
module subplex_global
    use global, only : DP,double,first_cost,total_fuel_cost
    real(DP) :: xkeep(5000,10),fmin_subplex(5000,1),fmin(2,1)
    integer, parameter :: subplex_log = 2
    integer :: fmin_index,count    
end module
module global_print
    use global, only : dir_loc
    character(len = 700) :: file1,file2,file3,file4,file5,file6,file7,file8, &
                            file9,file10,file11,file12,file13,file14,file15, &
                            file16,file17,file17a,file15a, &
                            file50,file51,file52,file53,file54   
    contains
    subroutine file_names(opt2)
    integer,intent(in),optional :: opt2
    character(len = 1) :: opt2c
    logical :: file_exists  
    character(len = 500) :: makedirectory
    
    !Check if the directory exists
    inquire(directory = trim(dir_loc),exist = file_exists)
    !if the directory doesn't exist, create it
    if(.not.file_exists)then   
        makedirectory = 'mkdir ' // dir_loc
        print*, makedirectory
        call system(trim(makedirectory))
    endif

    if(present(opt2))then  
        write(opt2c,'(i1)') opt2
        file1 = trim(dir_loc) // "\optimization_summary" // "_" // opt2c // ".txt"
        file2 = trim(dir_loc) // "\summary" // "_" // opt2c // ".txt"
        file3 = trim(dir_loc) // "\temps_day" // "_" // opt2c // ".txt"
        file4 = trim(dir_loc) // "\loads_day" // "_" // opt2c // ".txt"
        file5 = trim(dir_loc) // "\time_summary_day" // "_" // opt2c // ".txt"
        file6 = trim(dir_loc) // "\results_day" // "_" // opt2c // ".txt"
        file7 = trim(dir_loc) // "\inputs" // "_" // opt2c // ".txt"
        file8 = trim(dir_loc) // "\solver_stats" // "_" // opt2c // ".txt"
        file9 = trim(dir_loc) // "\power_day" // "_" // opt2c // ".txt"
        file10 = trim(dir_loc) // "\economics_day" // "_" // opt2c // ".txt"
        file11 = trim(dir_loc) // "\tempfield_day" // "_" // opt2c // ".txt"
        file12 = trim(dir_loc) // "\troubleshooting" // "_" // opt2c // ".txt"
        file13 = trim(dir_loc) // "\fluxfield_day" // "_" // opt2c // ".txt"
        file14 = trim(dir_loc) // "\time_summary" // "_" // opt2c // ".txt"
        file15 = trim(dir_loc) // "\direct" // "_" // opt2c // "_starter.txt"
        file15a = trim(dir_loc) // "\direct" // "_" // opt2c // "_main.txt"
        file16 = trim(dir_loc) // "\results" // "_" // opt2c // ".txt"
        file17 = trim(dir_loc) // "\subplex" // "_" // opt2c // "_starter.txt"
        file17a = trim(dir_loc) // "\subplex" // "_" // opt2c // "_main.txt"
        file50 = trim(dir_loc) // "\temperature_summary" // "_" // opt2c // ".txt"
        file51 = trim(dir_loc) // "\flow_rate_summary" // "_" // opt2c // ".txt"
        file52 = trim(dir_loc) // "\power_summary" // "_" // opt2c // ".txt"
        file53 = trim(dir_loc) // "\heat_transfer_summary" // "_" // opt2c // ".txt"
        file54 = trim(dir_loc) // "\misc_summary" // "_" // opt2c // ".txt"
    else
        file1 = trim(dir_loc) // "\optimization_summary.txt"
        file2 = trim(dir_loc) // "\summary.txt"
        file3 = trim(dir_loc) // "\temps_day.txt"
        file4 = trim(dir_loc) // "\loads_day.txt"
        file5 = trim(dir_loc) // "\time_summary_day.txt"
        file6 = trim(dir_loc) // "\results_day.txt"
        file7 = trim(dir_loc) // "\inputs.txt"
        file8 = trim(dir_loc) // "\solver_stats.txt"
        file9 = trim(dir_loc) // "\power_day.txt"
        file10 = trim(dir_loc) // "\economics_day.txt"
        file11 = trim(dir_loc) // "\tempfield_day.txt"
        file12 = trim(dir_loc) // "\troubleshooting.txt"
        file13 = trim(dir_loc) // "\fluxfield_day.txt"
        file14 = trim(dir_loc) // "\time_summary.txt"
        file15 = trim(dir_loc) // "\direct_starter.txt"
        file15a = trim(dir_loc) // "\direct_main.txt"
        file16 = trim(dir_loc) // "\results.txt"
        file17 = trim(dir_loc) // "\subplex_starter.txt"
        file17a = trim(dir_loc) // "\subplex_main.txt"
        file50 = trim(dir_loc) // "\temperature_summary.txt"
        file51 = trim(dir_loc) // "\flow_rate_summary.txt"
        file52 = trim(dir_loc) // "\power_summary.txt"
        file53 = trim(dir_loc) // "\heat_transfer_summary.txt"
        file54 = trim(dir_loc) // "\misc_summary.txt"
    endif    
    end subroutine file_names                                           
end module global_print

program main_call
    !use consts
    use global
    use subplex_global
    use opt_v_m
    use HyGSHP_calc
    use HyGSHP_calc_down
    !use weather_gen
    implicit none    
       
    integer :: i,runs,k,ii,jj,qq,hh,kk,nspeeds,wgen
    character(len = 50) :: a 
    character(len = 50) :: b 
    character(len = 50) :: c,gr
    real(DP) :: in(10),LCC,totallength,size_tower,dt1,tc1,tc2,tpc,dtpc,testvar(4),jindex(7)
    
    integer :: niwmax,nwmax,nxmax
    parameter (niwmax=300,nwmax=300,nxmax=300)
    integer :: iflag,iwork(niwmax),maxnfe,mdcont,mdsing,mduser,mode, &
               nf1,nf2,nfe,nfinc,n         
    real(DP) :: fx,scale(nxmax),scl, tol,tol1,tol2,tolfac,work(nwmax), &
                x(nwmax),junk,lower(nxmax),upper(nxmax)
    real(DP) :: tol1_main,tol1_starter,dfmin_lim_main,dfmin_lim_starter,dfmin_lim
    integer :: count_lim_main,count_lim_starter,count_lim
    
    character*26 flgdsc(-2:2)
    character*40 fname(5)
    character(len = 400) ::copy_file
    external subplx
    Integer :: findmin,mindex 
    data flgdsc /'invalid input', &
                  'maxnfe exceeded', &
                  'tol satisfied', &
                  'limit of machine precision', &
                  'fstop reached'/
    Integer :: DIRmaxf
    real(DP) :: DIRf,DIRx(128),opt_time,DIRx_print(15),sens_keep(5,9),sens_scale(5,2) 
    integer :: t1,t2,clock_rate,clock_max,sens_min      
    integer :: yr_orig,m_or_s,s
    real(DP) :: GHX_next, CCCT_next      
    real(DP) :: opt1,mean_input,std_input
    integer :: opt2
    integer :: temp_min,temp_max,temp_size,var_vect2(11)
    integer :: compare_conv,doloop,wgen_hot_cold,wgen_cold_hot
    real(DP) :: MULT_CTFANP,MULT_CTFANFL,MULT_CTFLUIDFL,Lo_SPD,lo_spd_tset,Lo_SPD_PC,spd1,spd2,spd3,spd4
    real(DP), allocatable :: n_ct_cont(:)    
    integer :: nseed
    integer :: Tghx_in_p,Thgx_out_p,Tghx_store,Tct_in_p,Tct_out_p,Thp_in_p,Thp_out_p1,Thp_out_p2,Thp_down,Twb_p,Tdb_p, &
               Flow_pump,flow_hpc,Flow_hph,Flow_CT,Flow_GHX,Pct_p,Php_p,Ppump_p,Qghx_h_p,Qghx_c_p, &
               Qct_p,Qhp_h_p,Qhp_c_p,Qfurnace_p,Qheat_p,Qcool_p,CT_speed_p
                  
    namelist /operation_settings/ strat,opt,yr,site,delt,opt_starter,opt_main,opt_method, &
                                  guess,print_stuff,sens_map,sens_var,downstream, &
                                  yr_opt,furnace,conv,compare_conv,energy_reduction, &
                                  conv_known,wgen,wgen_hot_cold,wgen_cold_hot,mean_input,std_input
    namelist /CT/ MULT_CTFANP,MULT_CTFANFL,MULT_CTFLUIDFL,Lo_SPD,lo_spd_tset,Lo_SPD_PC,spd1,spd2,nspeeds
    namelist /print_options/ Tghx_in_p,Thgx_out_p,Tghx_store,Tct_in_p,Tct_out_p, &
                             Thp_in_p,Thp_out_p1,Thp_out_p2,Thp_down,Twb_p,Tdb_p, &
                             Flow_pump,flow_hpc,Flow_hph,Flow_CT,Flow_GHX, &
                             Pct_p,Php_p,Ppump_p,Qghx_h_p,Qghx_c_p, &
                             Qct_p,Qhp_h_p,Qhp_c_p,Qfurnace_p,Qheat_p, &
                             Qcool_p,CT_speed_p
    
    call system_clock(t1,clock_rate,clock_max)
    
    jindex = 0
    comparison_results = 0
    
    a = "..\..\namelist_files\op_set.txt      "
    open(unit = 1,file = trim(a))
    read(1,nml = operation_settings)
    close(1)

    open(unit = 1,file = "..\..\namelist_files\component_inputs.txt")
    read(1,nml = CT)
    close(1)
    
    open(unit = 1,file = "..\..\namelist_files\print_options.txt")
    read(1,nml = print_options)
    close(1)
    
    print_vector = (/Tghx_in_p,Thgx_out_p,Tghx_store,Tct_in_p,Tct_out_p, &
                             Thp_in_p,Thp_out_p1,Thp_out_p2,Thp_down,Twb_p,Tdb_p, &
                             Flow_pump,flow_hpc,Flow_hph,Flow_CT,Flow_GHX, &
                             Pct_p,Php_p,Ppump_p,Qghx_h_p,Qghx_c_p, &
                             Qct_p,Qhp_h_p,Qhp_c_p,Qfurnace_p,Qheat_p, &
                             Qcool_p,CT_speed_p/)
                             
    if(conv == 1)then
        Tghx_in_p = 0
        Thgx_out_p = 0
        Tghx_store = 0
        Twb_p = 0
        Tdb_p = 0
        Flow_pump = 0
        Flow_hph = 0
        Flow_GHX = 0
        Qghx_h_p = 0
        Qghx_c_p = 0
        Qhp_h_p = 0
    endif
    
    var_vect_text = (/"GHX Length (m)","CT Size (kW)","CT1 (C)","CT2 (C)","GHX1 (C)","TPC (C)","DTPC (C)","Furnace (ratio)"/)
    
    if((compare_conv > 0).or.((energy_reduction > 0).and.(energy_reduction < 1000)))then
        doloop = 2
    else
        doloop = 1
    endif
    !wgen = 1
    !yr = 5
    call opening_weather(wgen,wgen_hot_cold,wgen_cold_hot,mean_input,std_input) 
    call ASHRAE_Sizing    
    
    do kk = 1,doloop
        var_vect = 0
        xpass = 0
        xkeep = 0
        first_thru = .true.
        
        if((kk == 1).and.(doloop == 2))then
            conv = 1
        elseif(kk == 2)then
            conv = 0
            first_run_opt = .false.
        endif
        
        if(nspeeds == 1) goto 103
        if(nspeeds == 2) goto 104
        if(nspeeds > 2) goto 105        
        
        !!!1 speed 
 103    if(conv)then
            var_vect(1:3) = (/2,2,3/)
            goto 100
        endif            
        
        if(strat == 0)then
            var_vect(1:3) = (/2,1,5/)    !CS0, GSHP only  
            if(furnace) var_vect(1:4) = (/3,1,5,8/)      
        elseif(strat == 3)then 
            var_vect(1:5) = (/4,1,2,3,5/) !CS2, tdiff
            if(furnace) var_vect(1:6) = (/5,1,2,3,5,8/) !CS2, tdiff
        elseif(strat == 5)then 
            var_vect(1:6) = (/5,1,2,5,6,7/)  !CS3, PC only
            if(furnace) var_vect(1:7) = (/6,1,2,5,6,7,8/)
        elseif(strat == 7)then 
            var_vect(1:7) = (/6,1,2,3,5,6,7/)  !CS5, Tdiff + PC
            if(furnace) var_vect(1:8) = (/7,1,2,3,5,6,7,8/)
        elseif(strat == 8)then 
            var_vect(1:7) = (/6,1,2,3,5,6,7/)
            if(furnace) var_vect(1:8) = (/7,1,2,3,5,6,7,8/)
        elseif(strat == 10)then 
            var_vect(1:3) = (/2,3,5/)
            !if(furnace) var_vect(1:4) = (/3,3,5,8/)
        elseif(strat == 13)then
            var_vect(1:5) = (/4,3,5,6,7/)  
            !if(furnace) var_vect(1:6) = (/5,3,5,6,7,8/) 
        endif    
        goto 100
        
        !!!2 speed
104     if(conv)then
            var_vect(1:4) = (/3,2,3,4/)
            goto 100
        endif            
        
        if(strat == 0)then
            var_vect(1:3) = (/2,1,5/)    !CS0, GSHP only  
            if(furnace) var_vect(1:4) = (/3,1,5,8/)      
        elseif(strat == 3)then 
            var_vect(1:6) = (/5,1,2,3,4,5/) !CS2, tdiff
            if(furnace) var_vect(1:7) = (/6,1,2,3,4,5,8/) !CS2, tdiff
        elseif(strat == 5)then 
            var_vect(1:6) = (/5,1,2,5,6,7/)  !CS3, PC only
            if(furnace) var_vect(1:7) = (/6,1,2,5,6,7,8/)
        elseif(strat == 7)then 
            var_vect(1:8) = (/7,1,2,3,4,5,6,7/)  !CS5, Tdiff + PC
            if(furnace) var_vect(1:9) = (/8,1,2,3,4,5,6,7,8/)
        elseif(strat == 8)then 
            var_vect(1:8) = (/7,1,2,3,4,5,6,7/)
            if(furnace) var_vect(1:9) = (/8,1,2,3,4,5,6,7,8/)
        elseif(strat == 10)then 
            var_vect(1:4) = (/3,3,4,5/)
            !if(furnace) var_vect(1:5) = (/4,3,4,5,8/)
        elseif(strat == 13)then
            var_vect(1:6) = (/5,3,4,5,6,7/)  
            !if(furnace) var_vect(1:7) = (/6,3,4,5,6,7,8/) 
        endif    
        goto 100
        
         !!!Continuous
105     if(conv)then
            var_vect(1:2) = (/1,2/)
            goto 100
        endif            
        
        if(strat == 0)then
            var_vect(1:3) = (/2,1,5/)    !CS0, GSHP only  
            if(furnace) var_vect(1:4) = (/3,1,5,8/)      
        elseif(strat == 3)then 
            var_vect(1:6) = (/5,1,2,3,4,5/) !CS2, tdiff
            if(furnace) var_vect(1:7) = (/6,1,2,3,4,5,8/) !CS2, tdiff
        elseif(strat == 5)then 
            var_vect(1:6) = (/5,1,2,5,6,7/)  !CS3, PC only
            if(furnace) var_vect(1:7) = (/6,1,2,5,6,7,8/)
        elseif(strat == 7)then 
            var_vect(1:8) = (/7,1,2,3,4,5,6,7/)  !CS5, Tdiff + PC
            if(furnace) var_vect(1:9) = (/8,1,2,3,4,5,6,7,8/)
        elseif(strat == 8)then 
            var_vect(1:8) = (/7,1,2,3,4,5,6,7/)
            if(furnace) var_vect(1:9) = (/8,1,2,3,4,5,6,7,8/)
        elseif(strat == 10)then 
            var_vect(1:4) = (/3,3,4,5/)
            !if(furnace) var_vect(1:5) = (/4,3,4,5,8/)
        elseif(strat == 13)then
            var_vect(1:6) = (/5,3,4,5,6,7/)  
            !if(furnace) var_vect(1:7) = (/6,3,4,5,6,7,8/) 
        endif  
        goto 100
                                
    100 if((kk == 1).and.(doloop == 2))then
            open(unit = 220,file = "..\..\results\opt_design_conv.txt",recl = 2140000) 
        else
           open(unit = 220,file = "..\..\results\opt_design.txt",recl = 2140000) 
        endif
        if(conv == 0)then 
            write(220,15) "totallength","size_towerkw","Furnace","ct1","ct2","ghx1","tpc","dtpc","LCC"
            write(220,15) "[m]","[kW]","[kBtu/hr]","[C]","[C]","[C]","[C]","[C]","[$]"
            15 format(a,',',a,',',a,',',a,',',a,',',a,',',a,',',a,',',a)  
        elseif(conv == 1)then
            write(220,25) "size_towerkw","Furnace","ct1","ct2","LCC"
            write(220,25) "[kW]","[kBtu/hr]","[C]","[C]","[$]"
            25 format(a,',',a,',',a,',',a,',',a)  
        endif
        !!!!!!!!!!!!!!!Print out a header saying what's going to occur!!!!!!!!!!!!!!!!!!!!
        if(opt == 0)then
            write(*,*) '************Single simulation************'
            write(*,*)
        else
            write(*,*) '************Optimization will be performed************'
            write(*,*)
        endif
        
        if((kk == 1).and.(conv_known == 1).and.(opt == 1))then  
            !Going to optimize hybrid design, but already know conventional design and want to
            !compare hybrid design to conventional, so start by calculating stats for conventional
            !design 
            call no_opt(kk,compare_conv) 
            goto 101
        endif
         
        if(opt == 0)then
            call no_opt(kk,compare_conv) 
        elseif(opt == 1)then
            if((kk == 2).or.(doloop == 1))then
                call opt_v(DIRx_print,DIRf,nspeeds)
            elseif(kk == 1)then
                call opt_v(DIRx_print,DIRf,nspeeds,opt2 = 1)
            endif
            
            if(conv == 0)then
                write(220,20) DIRx_print(1),DIRx_print(2),q_lh/1000*(1-DIRx_print(8)),DIRx_print(3),DIRx_print(4), &
                              DIRx_print(5),DIRx_print(6),DIRx_print(7),DIRf
            elseif(conv == 1)then
                write(220,30) DIRx_print(2),q_lh/1000*(1-DIRx_print(8)),DIRx_print(3),DIRx_print(4),DIRf 
            endif
                    
            if((sens_map == 1).and.(conv == 0))then
                if((sens_var <= 2).or.(sens_var == 8))then
                    testvar(1) = .8*DIRx_print(sens_var)
                    testvar(2) = .9*DIRx_print(sens_var)      
                    testvar(3) = 1.1*DIRx_print(sens_var) 
                    testvar(4) = 1.2*DIRx_print(sens_var) 
                else
                    testvar(1) = DIRx_print(sens_var) - 2
                    testvar(2) = DIRx_print(sens_var) - 1
                    testvar(3) = DIRx_print(sens_var) + 1
                    testvar(4) = DIRx_print(sens_var) + 2
                endif
                jj = 2;
                do ii = 2,var_vect(1)+1
                    if(var_vect(ii) /= sens_var)then
                        var_vect(jj) = var_vect(ii);
                        jj = jj+1;
                    endif
                enddo
                var_vect(var_vect(1)+1) = 0
                var_vect(1) = var_vect(1) - 1
                sens_keep(1,1:9) = (/DIRx_print(1),DIRx_print(2),DIRx_print(3),DIRx_print(4), &
                                     DIRx_print(5),DIRx_print(6),DIRx_print(7),DIRx_print(8),DIRf/)
                do ii = 1,4
                    call opt_v(DIRx_print,DIRf,nspeeds,opt1 = testvar(ii),opt2 = ii)
                    write(220,20) DIRx_print(1),DIRx_print(2),DIRx_print(3),DIRx_print(4), &
                                  DIRx_print(5),DIRx_print(6),DIRx_print(7),DIRx_print(8),DIRf
                    sens_keep(ii+1,1:9) = (/DIRx_print(1),DIRx_print(2),DIRx_print(3),DIRx_print(4), &
                                            DIRx_print(5),DIRx_print(6),DIRx_print(7),DIRx_print(8),DIRf/)             
                enddo
            endif
        endif
        16 format(a,',',a)
        20 format(f25.6,',',f12.7,',',f12.7,',',f12.7,',',f12.7,',',f12.7,',',f12.7,',',f12.7,',',f25.2)
        30 format(f12.7,',',f12.7,',',f12.7,',',f12.7,',',f25.2)
        21 format(f25.6,',',f12.2)
        close(unit = 220)
        close(unit = 221)
        
 101    conv = 0
    enddo
    
    !copy namelist files to the output directory:
    !component_inputs.txt,direct_optim.txt,op_set.txt,optim_settings.txt,print_options.txt
    fname(1) = "component_inputs.txt"
    fname(2) = "direct_optim.txt"
    fname(3) = "op_set.txt"
    fname(4) = "optim_settings.txt"
    fname(5) = "print_options.txt"
    do ii = 1,5
        copy_file = 'copy ' // '..\..\namelist_files\' // trim(fname(ii)) // ' ' // trim(dir_loc) // '\' // trim(fname(ii))   
        !write(*,*) copy_file
        call system(copy_file)
    enddo
    
    call system_clock(t2,clock_rate,clock_max)
    opt_time = real(t2 - t1)/real(clock_rate) 
    
    open(unit = 222,file = "complete.txt")
    write(222,*) "done"
    close(unit = 222)
    
end program main_call
















   

