module perturb
contains

    subroutine perturb_TMY(mu,sigma,Tdb,Tdbmean,rnumsin,out)
    use global, only: DP,month_hours
    use fit_interp_etc
    
    implicit none
    real(DP) :: out(8760)
    integer :: ii,jj
    real(DP) :: Tdbmean(365)
    real(DP) :: mu,sigma,Z1(12)
    real(DP) :: rnumsin(12)
    real(DP) :: Tdb(8760),std_perturb_day(8760),Tdb2(8760),Tdb3(8760)
    integer :: mwt
    integer,parameter :: ndata = 8760
    real(DP) :: a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
    real(DP) :: x1(365),y1(365)
    real(DP),allocatable :: coeff_month(:)
       
    allocate(coeff_month(3))
    x1 = 0
    y1 = 0
       
    call perturb_data(mu,sigma,rnumsin,Z1)    
   
    do ii = 1,12
        std_perturb_day(month_hours(ii,1):month_hours(ii,2)) = Z1(ii)
    enddo
    Tdb2 = Tdb + std_perturb_day
    Tdb3 = Tdb2
    !smooth data at each month. Will need to do the same at the change of each year
    do ii = 1,11
        do jj = 1,12
            x1(jj) = month_hours(ii,2) - 6 + jj - 1
        enddo
        y1(1:12) = Tdb2(x1(1:12))
        call LEASTSQR(x1,y1,12,2,coeff_month)
        Tdb3(x1(3:10)) = coeff_month(3)*x1(3:10)*x1(3:10)+coeff_month(2)*x1(3:10)+ &
                         coeff_month(1)             
    enddo

    out = Tdb3
    end subroutine perturb_TMY
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine perturb_data(mu,sigma,rnumsin,Z1)
    use global
    
    implicit none
    integer :: nseed,ii
    real(DP) :: rnumsin(12)
    real(DP) :: mu,sigma,Z1(12)
    
    !open(unit = 21,file = 'E:\Research\GUI\Weather_histogram\output_files\ran1.txt',ACCESS = 'APPEND',recl = 2140000)  
    do ii = 1,12
        Z1(ii) = mu + sigma*log((1-rnumsin(ii))/rnumsin(ii))
    enddo   
    !close(21)
    end subroutine perturb_data 
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine test_random_number(nseed,out)
    use global
    
    implicit none
    integer :: nseed
    REAL(DP) :: out(nseed)  
        CALL init_random_seed_orig(nseed)         ! see example of RANDOM_SEED
        CALL RANDOM_NUMBER(out)
    end subroutine test_random_number
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    SUBROUTINE init_random_seed_orig(n)
    INTEGER :: i, n, clock,clockrate,clockmax
    INTEGER, DIMENSION(:), ALLOCATABLE :: seed
    logical,save :: first_call
    integer,save :: oldseed

        ALLOCATE(seed(n))

        CALL SYSTEM_CLOCK(COUNT=clock,count_rate = clockrate,count_max = clockmax)

        seed = clock + 37 * (/ (i - 1, i = 1, n) /)
        CALL RANDOM_SEED(PUT = seed)

        oldseed = 2*seed(n)
        DEALLOCATE(seed)
        
    END SUBROUTINE init_random_seed_orig
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
    subroutine DD(Tdb,HDD,CDD)
    use global
    
    implicit none
    real(DP) :: Tmax(365),Tmin(365),Tmean(365),Tdb(8760)
    real(DP),parameter :: Tbase = 18.3
    real(DP) :: HDD,CDD,temph(365),tempc(365)
    integer :: ii
    
    do ii = 1,365
        Tmax(ii) = maxval(Tdb(day(ii,1):day(ii,2)))
        Tmin(ii) = minval(Tdb(day(ii,1):day(ii,2)))
        Tmean(ii) = (Tmax(ii) + Tmin(ii))/2;
    enddo
    do ii = 1,365
        temph(ii) = max(0.,Tbase - Tmean(ii))
        tempc(ii) = max(0.,Tmean(ii) - Tbase)
    enddo
    
    HDD = sum(temph)
    CDD = sum(tempc)
    
    end subroutine DD
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine std_mean(ydata,std,mean)
    use global
    
    implicit none
    real(DP) :: ydata(nyears),std,mean,delta(nyears)
    
    mean = sum(ydata)/nyears
    delta = (ydata - mean)**2
    std = sqrt( sum(delta)/(nyears - 1) )/mean*100    
    
    end subroutine std_mean
    
end module perturb











!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    subroutine ran1(nseed,out)
!    use global
!    !adapted from numerical recipes in FORTRAN 90
!    implicit none
!    integer :: idum,IA,IM,IQ,IR,NTAB,NDIV
!    real(DP) :: AM,EPS,RNMX
!    Parameter (IA = 16807,IM = 2147483647,AM = 1./IM,IQ = 127773,IR = 2836)
!    parameter (NTAB = 32,NDIV = 1+(IM-1)/NTAB,EPS = 1.2e-7,RNMX = 1.-EPS)
!    !"Minimal" random number generator of Park and Miller with Bays-Durham shuffle
!    !and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 
!    !(exclusive of the endpoint values). Call with idum a negative integer to initialize;
!    !thereafter, do not alter idum between successive deviates in a sequence. 
!    !RNMX should approximate the largest floating value that is less than 1.
!    integer :: j,k,iv(NTAB),iy,ii
!    save iv,iy
!    data iv /Ntab*0/, iy /0/
!    integer :: nseed
!    integer :: seed(nseed)
!    real(DP) :: out(nseed),prob
!   
!    call init_random_seed(nseed,seed)
!
!    do ii = 1,nseed
!        idum = seed(ii)
!        if((idum <= 0.).or.(iy == 0))then   !Prevent idum = 0
!            idum = max(-idum,1)             !Don't allow idum to be less than 0
!            do j = ntab+8,1,-1              !load the shuffle table after 8 warm-ups
!                k = idum/IQ
!                idum = IA*(idum-k*IQ) - IR*k
!                if(idum < 0) idum = idum + IM
!                if(j <= ntab) iv(j) = idum
!            enddo
!            iy = iv(1)
!        endif
!        k = idum/IQ     !start here when not initializing
!        idum = IA*(idum-k*IQ)-IR*k  !compute idum = mod(ia*idum,im) without overflows by Schrage's method
!        if(idum < 0) idum = idum + IM
!        j = 1+ iy/NDIV  !will be in the range 1:ntab
!        iy = iv(j)      !output previously stored value and refill the shuffle table
!        iv(j) = idum    
!        out(ii) = min(am*iy,rnmx)   !users don't expect endpoint values
!    enddo
!    
!    return
!    end subroutine ran1
!
!    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    subroutine ran0(nseed,out)
!    use global
!    
!    implicit none
!    integer :: nseed
!    real(DP) :: out(nseed)
!    integer :: IA,IM,IQ,IR,MASK,seed(nseed),idum
!    real(DP) :: AM
!    Parameter (IA = 16807,IM = 2147483647,AM=1./IM,IQ=127773,IR=2836,MASK=123459876)
!    integer :: k,ii
!    
!    call init_random_seed(nseed,seed)
!        
!    do ii = 1,nseed
!        idum = seed(ii)
!        idum = ieor(idum,MASK)
!        k = idum/IQ
!        idum = IA*(idum-k*IQ)-IR*k
!        if (idum < 0.) idum = idum + IM
!        out(ii) = AM*idum
!        idum = ieor(idum,MASK)
!    enddo    
!        
!    end subroutine ran0




!subroutine fit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q)
!    use global
!    
!    implicit none
!    integer :: mwt,ndata
!    real(DP) :: a,b,chi2,q,siga,sigb,sig(ndata),x(ndata),y(ndata)
!    !uses gammq; from Numerical Recipes in FORTRAN
!    !Given a set of data points x(1:ndata),y(1:ndata) with individual standard deviations
!    !sig(1:ndata), fit them to a straight line y = a + bx by minimizing chi-square. Returned are
!    !a,b and their respective probabl uncertainties siga and sigb, te chi-square chi2,j and
!    !the goodness-of-fit probability q (that the fit would have chi-square this large or larger).
!    !If mwt = 0 on input, then the standard deviations are assumed to be unavailable: q is returned
!    !as 1.0 and the normalization of chi2 is to unit standard deviation on all points
!    integer :: i
!    real(DP) :: sigdat,ss,st2,sx,sxoss,sy,t,wt,gammqj,temp1,temp2
!        sig = max(0.1*x,0.0001)
!        sx = 0. !initialize sums to zero
!        sy = 0.
!        st2 = 0.
!        b = 0.
!        if(mwt /= 0)then    !Accumulate sums...
!            ss = 0.
!            do 11 i = 1,ndata       !...with weights
!                wt = 1./(sig(i)**2)
!                ss = ss + wt
!                sx = sx + x(i)*wt
!                sy = sy + y(i)*wt
!            11 enddo 
!        else
!            do 12 i = 1,ndata   !...or without weights.
!                sx = sx + x(i)
!                sy = sy + y(i)
!            12 enddo 
!            ss = float(ndata)
!        endif
!        sxoss = sx/ss
!        if(mwt /= 0)then
!            do 13 i = 1,ndata
!                t = (x(i) - sxoss)/sig(i)
!                st2 = st2 + t*t
!                b = b + t*y(i)/sig(i)
!            13 enddo 
!        else
!            do 14 i = 1,ndata
!                t = x(i) - sxoss
!                st2 = st2 + t*t
!                b = b + t*y(i)
!            14 enddo 
!        endif
!        b = b/st2   !solve for a, b, sigma_a and sigma_b
!        a = (sy - sx*b)/ss
!        siga = sqrt((1.+sx*sx/(ss*st2))/ss)
!        sigb = sqrt(1./st2)
!        chi2 = 0.   !calculate chi-squared
!        q = 1.
!        if(mwt == 0)then
!            do 15 i = 1,ndata
!                chi2 = chi2 + (y(i) - a - b*x(i))**2
!            15 enddo 
!            sigdat = sqrt(chi2/(ndata-1))  !for unweighted data evaluate typical sig using
!            siga = siga*sigdat             !chi2, and adjust the standard deviations.
!            sigb = sigb*sigdat
!        else
!            do 16 i = 1,ndata
!                chi2 = chi2 + ((y(i) - a - b*x(i))/sig(i))**2
!            16 enddo 
!            temp1 = 0.5*(ndata - 2)
!            temp2 = 0.5*chi2
!            if(ndata > 2) q = gammq(temp1,temp2)        
!        endif
!        return
!    end subroutine fit
!    
!    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    function gammq(a,x)
!    use global
!    
!    implicit none
!    real(DP) :: a,gammq,x
!    !uses gcf,gser
!    !returns the incomplete gamma function Q(a,x) = 1-P(a,x)
!    real(DP) :: gammcf,gamser,gln
!    
!    if((x < 0.).or.(a <= 0.))pause 'bad arguments in gammq'
!    if(x < a+1.)then    !use the series representation
!        call gser(gamser,a,x,gln)
!        gammq = 1. - gamser !and take its complement.
!    else    !use the continued fraction representation
!        call gcf(gammcf,a,x,gln)
!        gammq = gammcf
!    endif
!    end function gammq
!    
!    subroutine gser(gamser,a,x,gln)
!    use global
!    
!    implicit none
!    integer :: itmax
!    real(DP) :: a,gamser,gln,x,EPS
!    parameter (itmax = 100,eps = 3.e-7)
!    !uses gamm1n
!    !returns the incomplete gamma function P(a,x) evaluated by its series represntation as
!    !gamser. Also returns ln gamma(a) as gln
!    integer :: n
!    real(DP) :: ap,del,sum,gammln
!    gln = gammln(a)
!    if(x <= 0.)then
!        if(x < 0.) pause 'x < 0 in gser'
!        gamser = 0.
!        return
!    endif
!    ap = a
!    sum = 1./a
!    del = sum
!    do 11 n = 1,itmax
!        ap = ap + 1.
!        del = del*x/ap
!        sum = sum + del
!        if(abs(del) < abs(sum)*eps) goto 1
!    11 enddo
!    pause 'a too large, itmax too small in gser'
!    1 gamser = sum*exp(-x+a*log(x) - gln)
!    return
!    end subroutine gser
!    
!    subroutine gcf(gammcf,a,x,gln)
!    use global
!    
!    implicit none
!    integer :: itmax
!    real(DP) :: a,gammcf,gln,x,eps,fpmin
!    parameter(itmax = 100,eps = 3.e-7,fpmin = 1.e-30)
!    !returns the incomplete gamma function Q(a,x) evaluated by its continued fraction
!    !representation as gammcf. also returns gln. itmax is the max allowed number of iterations
!    !eps is the relative accuracy, fpmin is the number near the smallest represntable floating point number
!    integer :: i
!    real(DP) :: an,b,c,d,del,h,gammln
!    gln = gammln(a)
!    b = x + 1. - a
!    c = 1./fpmin
!    d = 1./b
!    h = d
!    do 11 i = 1,itmax
!        an = -i*(i-a)
!        b = b+2
!        d = an*d + b
!        if(abs(d) < fpmin) d = fpmin
!        c = b + an/c
!        if(abs(c) < fpmin) c = fpmin
!        d = 1./d
!        del = d*c
!        h = h*del
!        if(abs(del - 1.) < eps) goto 1
!    11 enddo
!    pause 'a too large, itmax too small in gcf'
!    1 gammcf = exp(-x+a*log(x)-gln)*h
!    return     
!    end subroutine gcf
!    
!    function gammln(XX)
!    use global
!    
!    implicit none
!    real(DP) :: gammln,xx
!        !returns ln(gamma) for xx > 0
!    integer :: j
!    real(DP) :: ser,stp,tmp,x,y,cof(6)
!    save cof,stp
!    data cof,stp/76.18009172947146d0,-86.50532032941677d0, &
!        24.01409824083091d0,-1.231739572450155d0,0.1208650973866179d-2, &
!        -0.5395239384953d-5,2.5066282746310005d0/
!        
!        x = xx
!        y = x
!        tmp = x + 5.5d0
!        tmp = (x + 0.5d0)*log(tmp) - tmp
!        ser = 1.000000000190015d0
!        do 11 j = 1,6
!            y = 1 + 1.d0
!            ser = ser + cof(j)/y
!        11 enddo
!        gammln = tmp + log(stp*ser/x)
!        return    
!    end function gammln