subroutine opening_weather(wgen,wgen_hot_cold,wgen_cold_hot,mu,sigma)
!use consts
use global
use fit_interp_etc
use perturb

implicit none
character(len = 500) :: loadfilename,makedirectory
integer,parameter :: headerlines = 2
integer :: i,ii,jj,dummy,yra,nyears2,j,tempint,wgen,nseed,wgen_hot_cold,wgen_cold_hot
real(DP) :: qcoola_temp(40*8760),qheata_temp(40*8760),tdba_temp(40*8760), &
            twba_temp(40*8760),ti_temp(40*8760),time_gr,time_step,timeold
real(DP) :: out_syn(40*8760,4),rnums(12*40),a,mu,sigma
real(DP), allocatable :: qcool2(:),qheat2(:),tdb2(:),twb2(:), &
                         qcool3(:),qheat3(:),tdb3(:),twb3(:)
real(DP) :: HDD(yr),CDD(yr),yr_index(yr,2),rat(yr),out_sort(yr,2)
real(DP) :: qcool_temp(40*8760),qheat_temp(40*8760),tdb_temp(40*8760),twb_temp(40*8760)
real(DP) :: temp
logical :: file_exists   
character :: delimiter
character*150 :: command                                            
                         
open (unit = 17, status = "old", file = "..\..\Building Loads\loadfilename.txt")  
read(17,*)
read(17,*)
read(17,*) loadfilename
read(17,*) dir_loc
close(17)


!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

open (unit = 17, status = "old", file = trim(loadfilename))  

do i = 1,headerlines
    read(17,*)
enddo 

ti = 0
qcoola = 0
qheata = 0
tdba = 0
twba = 0

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Run this portion if user has input a multi-yr file
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!read in all of the data in the file and store it in inputs
!if multi_yr == 1
    dummy = 0
    ii = 1
    jj = 1
    do while (jj == 1)
        read(17,*,Iostat = dummy) ti_temp(ii),qcoola_temp(ii),qheata_temp(ii),tdba_temp(ii),twba_temp(ii)     
        ii = ii + 1
        if(ti_temp(ii-1) < 1)then 
            jj = ii - 1
        endif
    enddo
    10  ndays = ii - 1
    nyears = ndays/8760

    ti = ti_temp(1:nyears*8760)
    qcoola = qcoola_temp(1:nyears*8760)
    qheata = qheata_temp(1:nyears*8760)
    tdba = tdba_temp(1:nyears*8760)
    twba = twba_temp(1:nyears*8760)
    
    !write(*,*) "qcoola = ", size(qcoola)/8760
    
    allocate(qcool2(ilen*(yr)+1))
    allocate(qheat2(ilen*(yr)+1))
    allocate(tdb2(ilen*(yr)+1))
    allocate(twb2(ilen*(yr)+1))  
    
    allocate(qcool3(ilen*(yr)+1))
    allocate(qheat3(ilen*(yr)+1))
    allocate(tdb3(ilen*(yr)+1))
    allocate(twb3(ilen*(yr)+1))  
       
    qcool2 = 0
    qheat2 = 0
    tdb2 = 0
    twb2 = 0
    
    if(nyears > yr)then
        qcool2(1:yr*8760) = qcoola(1:yr*8760)
        qheat2(1:yr*8760) = qheata(1:yr*8760)
        tdb2(1:yr*8760) = tdba(1:yr*8760)
        twb2(1:yr*8760) = twba(1:yr*8760)      
    !    write(*,*) size(qcool2)  
    elseif(nyears == yr)then
        qcool2 = qcoola
        qheat2 = qheata
        tdb2 = tdba
        twb2 = twba        
    elseif((nyears < yr).and.(wgen /= 1))then      
   ! elseif((nyears < yr).and.(wgen == 0))then       
        yra = nyears
        do while(yra <= yr)
            qcool2((yra - nyears)*8760 + 1:yra*8760) = qcoola
            qheat2((yra - nyears)*8760 + 1:yra*8760) = qheata
            tdb2((yra - nyears)*8760 + 1:yra*8760) = tdba
            twb2((yra - nyears)*8760 + 1:yra*8760) = twba
            yra = yra + nyears
        enddo
    elseif((nyears < yr).and.(wgen == 1))then    
       nyears_gen = yr - nyears
       
        nseed = 12*nyears_gen
        call test_random_number(nseed,rnums)
       !allocate(out_syn(nyears_gen*8760,4))

       call weather_gen(qcoola(1:8760),qheata(1:8760),tdba(1:8760),twba(1:8760),rnums,mu,sigma,out_syn)
       qcool2(1:nyears*8760) = qcoola(1:nyears*8760)
       qcool2(nyears*8760+1:yr*8760) = out_syn(1:nyears_gen*8760,1)
       qheat2(1:nyears*8760) = qheata(1:nyears*8760)
       qheat2(nyears*8760+1:yr*8760) = out_syn(1:nyears_gen*8760,2)
       tdb2(1:nyears*8760) = tdba(1:nyears*8760)
       tdb2(nyears*8760+1:yr*8760) = out_syn(1:nyears_gen*8760,3)
       twb2(1:nyears*8760) = twba(1:nyears*8760)
       twb2(nyears*8760+1:yr*8760) = out_syn(1:nyears_gen*8760,4)
    endif
    
    if(yr_opt > yr) yr_opt = yr
    
    allocate(qcool_opt(int(ilen*yr_opt)+1))
    allocate(qheat_opt(int(ilen*yr_opt)+1))
    allocate(tdb_opt(int(ilen*yr_opt)+1))
    allocate(twb_opt(int(ilen*yr_opt)+1))
    
    qcool_opt = qcool2(1:(yr_opt*8760+1))
    qheat_opt = qheat2(1:(yr_opt*8760+1))
    tdb_opt = tdb2(1:(yr_opt*8760+1))
    twb_opt = twb2(1:(yr_opt*8760+1))
    
    a = size(Qcool2)/8760
!    open(unit = 21,file = 'E:\Research\GUI\program_files\results\weather_Test.txt',recl = 210000)
!    do ii = 1,yr*8760
!        write(21,*) ii,qcool2(ii),qheat2(ii),tdb2(ii),twb2(ii)
!    enddo
!    close(unit = 21)
    nyears2 = size(qcool2)/8760
     
    if(nyears2 > yr)then
        qcool3 = qcool2(1:yr*8760)
        qheat3 = qheat2(1:yr*8760)
        tdb3 = tdb2(1:yr*8760)
        twb3 = twb2(1:yr*8760) 
    else
        qcool3 = qcool2
        qheat3 = qheat2
        tdb3 = tdb2
        twb3 = twb2
    endif
    !write(*,*) tdb3(1),tdb3(8761),tdb3(17521),tdb3(26281),tdb3(35041)
    !Calculate CDD and HD and order the years from greatest HDD/CDD to least and 
    !vice versa
    !wgen_hot_cold = 1
    if((wgen_hot_cold + wgen_cold_hot) >= 1)then
        yr_index(1,:) = (/1,8760/)
        do ii = 2,yr
            yr_index(ii,:) = yr_index(ii-1,:) + 8760
        enddo
        do ii = 1,yr
            call DD(tdb3(yr_index(ii,1):yr_index(ii,2)),HDD(ii),CDD(ii))
        enddo
        !calculate ratio
        rat = HDD/CDD
        !sort data in ascending order
        if(wgen_hot_cold == 1)then
            call PIKSRT(yr,rat,out_sort)
        elseif(wgen_cold_hot == 1)then
            rat = 1/rat
            call PIKSRT(yr,rat,out_sort)
        endif
        do ii = 1,yr
            tdb_temp(yr_index(ii,1):yr_index(ii,2)) = tdb3(yr_index(out_sort(ii,2),1):yr_index(out_sort(ii,2),2))
            twb_temp(yr_index(ii,1):yr_index(ii,2)) = twb3(yr_index(out_sort(ii,2),1):yr_index(out_sort(ii,2),2))
            qcool_temp(yr_index(ii,1):yr_index(ii,2)) = qcool3(yr_index(out_sort(ii,2),1):yr_index(out_sort(ii,2),2))
            qheat_temp(yr_index(ii,1):yr_index(ii,2)) = qheat3(yr_index(out_sort(ii,2),1):yr_index(out_sort(ii,2),2))
        enddo   
        tdb3 = tdb_temp
        twb3 = twb_temp
        qcool3 = qcool_temp
        qheat3 = qheat_temp     
    endif
    
    !write(*,*) tdb3(1),tdb3(8761),tdb3(17521),tdb3(26281),tdb3(35041)
              
    allocate(qcool(int(ilen*yr/delt)+1))
    allocate(qheat(int(ilen*yr/delt)+1))
    allocate(qcool_orig(int(ilen*yr/delt)+1))
    allocate(qheat_orig(int(ilen*yr/delt)+1))
    allocate(tdb_orig(int(ilen*yr/delt)+1))
    allocate(twb_orig(int(ilen*yr/delt)+1))
    allocate(tdb(int(ilen*yr/delt)+1))
    allocate(twb(int(ilen*yr/delt)+1)) 
    
    qcool = 0
    qheat = 0
    tdb = 0
    twb = 0   
    qcool_orig = 0
    qheat_orig = 0
    tdb_orig = 0
    twb_orig = 0 
         
    time_step = ilen*yr/delt !total simulation time considering a time step         
    if (delt /= 1.) then
        time_gr = 1
        qcool(1) = qcool3(1)
        qheat(1) = qheat3(1)
        tdb(1) = tdb3(1)
        twb(1) = twb3(1)
        do j = 2,time_step+1  
           timeold = time_gr
           tempint = int(time_gr)
           qcool(j) = qcool3(tempint)
           qheat(j) = qheat3(tempint)
           tdb(j) = tdb3(tempint)
           twb(j) = twb3(tempint)
           time_gr = time_gr + delt  
        enddo
    else
        qcool = qcool2
        qheat = qheat2
        tdb = tdb2
        twb = twb2
    endif
    close(unit = 17)    
    
    qcool_orig = qcool
    qheat_orig = qheat
    tdb_orig = tdb
    twb_orig = twb
    
    deallocate(qcool2)
    deallocate(qheat2)
    deallocate(tdb2)
    deallocate(twb2)  
    
    q_lc = -maxval(qcoola)*0.9478    !peak cooling load, kW converted to Btu/h
!    if(furnace == 0)then
    q_lh = maxval(qheata)*0.9478     !peak heating load, kW converted to Btu/h
!    else
!        q_lh = 0.85*maxval(qheata)*0.9478     !peak heating load assuming the furnace picks up 15% of the peak, kW converted to Btu/h
!    endif
    !write(*,*) maxval(Qcoola),maxval(qcool_opt),maxval(qcool),maxval(Qcool_orig),q_lc/.9478,maxval(qcoola_temp)
    !write(*,*) maxval(Qheata),maxval(qheat_opt),maxval(qheat),maxval(Qheat_orig),q_lh/.9478,maxval(qheata_temp)
end subroutine opening_weather