subroutine weather_gen(qchour,qhhour,tdb,twb,rnums,mu,sigma,out_syn)
use global, only: DP,nyears,month_hours,day,yr
use perturb
use fit_interp_etc

implicit none
real(DP) :: Tdb(8760),Twb(8760),month(12,2),vector(40,2)
real(DP) :: mu,sigma,rnums(12*40)
integer :: ii,jj,kk,nseed
real(DP) :: test,check(365)
real(DP) :: time(8760),Qchour(8760),Qhhour(8760)
real(DP) :: Tdbmax(365),Tdbmin(365),Tdbmean(365),Qcday(365),Qhday(365)
real(DP) :: Twbmax(365),Twbmin(365),Twbmean(365),coeff_wb(2)
real(DP) :: coeff_cool(3),coeff_heat(3)
real(DP) :: x1(365),x1a(365),c(2),x1b(365),y1b(365)
real(DP) :: DTwb(365),DTwbhour(8760)
real(DP) :: Qcfit(365),Qhfit(365),Qcfit_day(8760),Qhfit_day(8760)
integer :: kkc,kkh,nyears_gen
real(DP) :: vtempc(24),vtemph(24),rcday(365),rhday(365)
real(DP),allocatable :: Qc_syn(:,:)
real(DP),allocatable :: Qh_syn(:,:),Tdb_syn(:,:),Tdb_syn2(:,:)
real(DP),allocatable :: Twb_syn(:,:),Tdbmean2(:,:),Twbmean2(:,:)
real(DP),allocatable :: yr_index(:,:),qc_total(:),qh_total(:)
real(DP),allocatable :: tdb_total(:),twb_total(:)
real(DP) :: out_syn(40*8760,4)
real(DP) :: out_TMY(8760,1)

nyears_gen = yr - nyears

allocate(Qc_syn(8760,nyears_gen))
allocate(Qh_syn(8760,nyears_gen))
allocate(Tdb_syn(8760,nyears_gen))
allocate(Tdb_syn2(8760,nyears_gen))
allocate(Twb_syn(8760,nyears_gen))
allocate(Tdbmean2(365,nyears_gen))
allocate(Twbmean2(365,nyears_gen))
allocate(yr_index(nyears_gen,2))
allocate(qc_Total(nyears_gen*8760))
allocate(qh_Total(nyears_gen*8760))
allocate(tdb_Total(nyears_gen*8760))
allocate(twb_Total(nyears_gen*8760))

out_syn = 0
Qc_syn = 0
Qh_syn = 0
Tdb_syn = 0
Tdb_syn2 = 0
Twb_syn = 0
Tdbmean2 = 0
Twbmean2 = 0

vector(1,1) = 1
vector(1,2) = 12
do ii = 2,nyears_gen
    vector(ii,:) = vector(ii-1,:) + 12
enddo

month(:,1) = (/1,32,60,91,121,152,182,213,244,274,305,335/)
month(:,2) = (/31,59,90,120,151,181,212,243,273,304,334,365/)
month_hours = month*24
month_hours(:,1) = month_hours(:,1) - 23

day(1,:) = (/1,24/)
do ii = 2,365
    day(ii,:) = day(ii-1,:) + 24;
enddo

!read in load file
!Get data into the form needed for sub-programs
!!!!!!!get rid of this stuff in final app!!!!!!!
!open(unit = 11, file = 'E:\Research\GUI\program_files\Building Loads\ATL\ATLTMY3.txt')
!read(11,*)
!read(11,*)
!do ii = 1,8760
!    read(11,*) time(ii),Qchour(ii),Qhhour(ii),Tdb(ii),Twb(ii)
!enddo

do ii = 1,365
    Tdbmax(ii) = maxval(Tdb(day(ii,1):day(ii,2)))
    Tdbmin(ii) = minval(Tdb(day(ii,1):day(ii,2)))
    Tdbmean(ii) = (Tdbmax(ii) + Tdbmin(ii))/2
    
    Twbmax(ii) = maxval(Twb(day(ii,1):day(ii,2)))
    Twbmin(ii) = minval(Twb(day(ii,1):day(ii,2)))
    Twbmean(ii) = (Twbmax(ii) + Twbmin(ii))/2
    
    Qcday(ii) = sum(Qchour(day(ii,1):day(ii,2)))
    Qhday(ii) = sum(Qhhour(day(ii,1):day(ii,2)))
enddo

!Get the curve fit for wet bulb to dry bulb
!coeff_wb is the coefficient vector
!cw(2)*x + cw(1)
call fit(Tdbmean,Twbmean,365,coeff_wb)
!Call leastsqr method to find curve fit for Qday vs. Tmean
!coeff_cool/coeff_heat returns the coefficients as:
!cc(3)*x^2 + cc(2)*x + cc(1)
!allocate(coeff_cool(3))
!allocate(coeff_heat(3))
call LEASTSQR(Tdbmean,Qcday,365,2,coeff_cool)
call LEASTSQR(Tdbmean,Qhday,365,2,coeff_heat)

!Get synthetic dry bulb temperatures
!this is where you would change the code to accomodate user input
!something like:
!if(specific_fit == 1)then
!   get mu and sigma from namelist file, which is read in calling code because
!   this file will also contain the info about whether or not weather needs
!   to be generated
!else
if((mu > 0.385*1.001).or.(mu < 0.385*0.999))then
    mu = 0.9206*mu + 0.0555
endif
if((sigma > 1.186*1.001).or.(sigma < 1.186*0.999))then
    sigma = 0.5055*sigma + 0.1235
endif
do ii = 1,nyears_gen        
    call perturb_TMY(mu,sigma,Tdb,Tdbmean,rnums(vector(ii,1):vector(ii,2)),out_TMY)  
    Tdb_syn(1:8760,ii) = out_TMY(1:8760,1)    
enddo

!smooth data at each year
x1(1:26) = (/8748,8749,8750,8751,8752,8753,8754,8755,8756,8757,8758,8759,8760,&
             1,2,3,4,5,6,7,8,9,10,11,12,13/) 
x1a(1:26) = (/8748,8749,8750,8751,8752,8753,8754,8755,8756,8757,8758,8759,8760,&
              8761,8762,8763,8764,8765,8766,8767,8768,8769,8770,8771,8772,8773/)          
Tdb_syn2 = Tdb_syn

do ii = 1,nyears_gen-1    
    x1b(1) = 8756
    x1b(2) = 8765
    y1b(1) = Tdb_syn(8756,ii)   
    y1b(2) = Tdb_syn(5,ii+1)   
    c(2) = (y1b(2) - y1b(1))/(x1b(2) - x1b(1))     
    c(1) = y1b(1) - c(2)*x1b(1)             
    Tdb_syn2(x1(9:13),ii) = c(2)*x1a(9:13) + c(1)  
    Tdb_syn2(x1(14:18),ii+1) = c(2)*x1a(14:18) + c(1)                                             
enddo

!Now generate the wet bulb temperatures   
do ii = 1,nyears_gen
    do jj = 1,365
        Tdbmax(jj) = maxval(Tdb_syn2(day(jj,1):day(jj,2),ii))
        Tdbmin(jj) = minval(Tdb_syn2(day(jj,1):day(jj,2),ii))
        Tdbmean2(jj,ii) = (Tdbmax(jj) + Tdbmin(jj))/2
    enddo
    Twbmean2(:,ii) = coeff_wb(2)*Tdbmean2(:,ii) + coeff_wb(1)
    DTwb = Twbmean - Twbmean2(:,ii)
    do jj = 1,365
        DTwbhour(day(jj,1):day(jj,2)) = DTwb(jj)
    enddo
    Twb_syn(:,ii) = Twb + DTwbhour
enddo

!open(unit = 21,file = 'E:\Research\GUI\Weather_histogram\output_files\wet_bulb.txt',recl = 2100000)
!do ii = 1,8760
!    write(21,*) ii,Twb(ii),Twb_syn(ii,:)
!enddo
!close (unit = 21)

do kk = 1,nyears_gen
    !Now generate the cooling and heating loads 
    Qcfit = max(0.,coeff_cool(3)*Tdbmean2(:,kk)*Tdbmean2(:,kk) + coeff_cool(2)*Tdbmean2(:,kk) + coeff_cool(1))
    Qhfit = max(0.,coeff_heat(3)*Tdbmean2(:,kk)*Tdbmean2(:,kk) + coeff_heat(2)*Tdbmean2(:,kk) + coeff_heat(1))
    !figure out how many hours a day there is a load    
    do ii = 1,365
        kkc = 0
        kkh = 0
        vtempc = Qchour(day(ii,1):day(ii,2))
        vtemph = Qhhour(day(ii,1):day(ii,2))
        do jj = 1,24
            if(vtempc(jj) > 0)then
                kkc = kkc + 1
            endif
            if(vtemph(jj) > 0)then
                kkh = kkh + 1
            endif
        enddo
        if(kkc > 0)then
            rcday(ii) = kkc - 1
        else
            rcday(ii) = 0
        endif
        if(kkh > 0)then
            rhday(ii) = kkh - 1
        else
            rhday(ii) = 0
        endif
    enddo

    do ii = 1,365
        if(rcday(ii) > 0)then
            Qcfit_day(day(ii,1):day(ii,2)) = (Qcfit(ii) - Qcday(ii))/rcday(ii)
        else
            Qcfit_day(day(ii,1):day(ii,2)) = 0.
        endif
        if(rhday(ii) > 0)then
            Qhfit_day(day(ii,1):day(ii,2)) = (Qhfit(ii) - Qhday(ii))/rhday(ii)
        else
            Qhfit_day(day(ii,1):day(ii,2)) = 0.
        endif
    enddo
    do ii = 1,8760
        if(Qchour(ii) > 0)then
            Qc_syn(ii,kk) = max(0.,Qcfit_day(ii) + Qchour(ii))
        else
            Qc_syn(ii,kk) = Qchour(ii)
        endif
        if(Qhhour(ii) > 0)then
            Qh_syn(ii,kk) = max(0.,Qhfit_day(ii) + Qhhour(ii))
        else
            Qh_syn(ii,kk) = Qhhour(ii)
        endif
    enddo
enddo
yr_index(1,:) = (/1,8760/)
do ii = 2,nyears_gen
    yr_index(ii,:) = yr_index(ii-1,:) + 8760
enddo
do ii = 1,nyears_gen
    Qc_total(yr_index(ii,1):yr_index(ii,2)) = Qc_syn(:,ii)
    Qh_total(yr_index(ii,1):yr_index(ii,2)) = Qh_syn(:,ii)
    Tdb_total(yr_index(ii,1):yr_index(ii,2)) = Tdb_syn2(:,ii)
    Twb_total(yr_index(ii,1):yr_index(ii,2)) = Twb_syn(:,ii)
enddo
out_syn(1:8760*nyears_gen,1) = Qc_total
out_syn(1:8760*nyears_gen,2) = Qh_total
out_syn(1:8760*nyears_gen,3) = Tdb_total
out_syn(1:8760*nyears_gen,4) = Twb_total
end subroutine weather_gen

