!
! Test for trend-variability relation in high-resolution EDC time series from MIS 1, MIS 5 and MIS 7
!
! f:/EDC-VAR/test/
!
!=======================================================================================================================
!
! Author
! ======
!
! Manfred Mudelsee
! Climate Risk Analysis
! Schneiderberg 26
! 30167 Hannover
! Germany
! Email: mudelsee@mudelsee.com
! URL: http://www.mudelsee.com
!
!=======================================================================================================================
!
! Change log
! ==========
!
! Version       Date                                    Comments
!
! 1.1            3 May   2011                           MIS1 added
! 1.0           18 April 2011 to 25 April 2011          original version
!
!=======================================================================================================================
!
! NR files
! ========
!
include 'd:\nr7790\double\nrtype.f90' ! uncomment this line if using gfortran instead of AbSoft
include 'd:\nr7790\double\nrutil.f90' ! uncomment this line if using gfortran instead of AbSoft
include 'd:\nr7790\double\nr.f90' ! uncomment this line if using gfortran instead of AbSoft
!
!=======================================================================================================================
!
module datatx
        use nrtype
        implicit none
!       Data: t, x, t_rm_sim, x_rm_sim.
        save
!
        real(sp), allocatable, dimension(:) :: t ! time
        real(sp), allocatable, dimension(:) :: x ! x-value
        real(sp), allocatable, dimension(:) :: t_rm_sim ! time
        real(sp), allocatable, dimension(:) :: x_rm_sim ! x-value
end module datatx
!
!=======================================================================================================================
!
module parameters
        use nrtype
        implicit none
!       Parameters.
        save
!
! 1.    Data
!       ====
!
        character (len=4), parameter :: mis1 = 'MIS1' ! attribute
        character (len=4), parameter :: mis5 = 'MIS5' ! attribute
        character (len=4), parameter :: mis7 = 'MIS7' ! attribute
        character (len=16), parameter :: infile_mis1 = '../data/mis1.txt' ! file name input
        character (len=16), parameter :: infile_mis5 = '../data/mis5.txt' ! file name input
        character (len=16), parameter :: infile_mis7 = '../data/mis7.txt' ! file name input
        character (len=12), parameter :: outfile_mis1 =    'testmis1.txt' ! file name output
        character (len=12), parameter :: outfile_mis5 =    'testmis5.txt' ! file name output
        character (len=12), parameter :: outfile_mis7 =    'testmis7.txt' ! file name output
!
        integer, parameter :: n_mis1=639 ! data size
        integer, parameter :: n_mis5=2416 ! data size
        integer, parameter :: n_mis7=505 ! data size
!
! 2.    Trend estimation 1: Running mean, number of window points
!       =========================================================
!
        integer, parameter :: k_rm_mis1=163
        integer, parameter :: k_rm_mis5=335
        integer, parameter :: k_rm_mis7=153
!
! 3.    Trend estimation 2: CHECK
!       =========================
!
!
! 4.    Segmentation
!       ============
!
!       Note: n points are put into m-1 segments of k_seg points and 1 last segment with the remaining points.
!             Average trend and standard deviation are then calculated over segments points.
!             => Adjust k_seg start, end  and step values such to avoid the last segment having only 1 point.
!
!
        integer, parameter :: k_seg_start_mis1=10 ! number of points per segment, start
        integer, parameter :: k_seg_start_mis5=10 ! number of points per segment, start
        integer, parameter :: k_seg_start_mis7=10 ! number of points per segment, start
!
        integer, parameter :: k_seg_end_mis1=60 ! number of points per segment, end
        integer, parameter :: k_seg_end_mis5=250 ! number of points per segment, end
        integer, parameter :: k_seg_end_mis7=50 ! number of points per segment, end
!
        integer, parameter :: k_seg_delta_mis1=10 ! number of points per segment, step
        integer, parameter :: k_seg_delta_mis5=10 ! number of points per segment, step
        integer, parameter :: k_seg_delta_mis7=10 ! number of points per segment, step
!
! 5.    Bootstrap simulations
!       =====================
!
        integer, parameter :: n_sim=100000
!
!
! x.    Other
!       =====
!
        integer :: ibig=10000
end module parameters
!
!=======================================================================================================================
!
module residuals
        use nrtype
        implicit none
!       Residuals.
        save
        real(sp), allocatable, dimension(:) :: e_rm ! running-mean residuals
        real(sp), allocatable, dimension(:) :: e_rm_sim ! running-mean residuals, simulation
end module residuals
!
!=======================================================================================================================
!
module segmented
        use nrtype
        implicit none
!       Segmented trends and variabilities.
        save
        real(sp), allocatable, dimension(:) :: x_trend_rm_seg ! trend, running mean
        real(sp), allocatable, dimension(:) :: x_trend_rm_sim_seg ! trend, running mean, simulation
        real(sp), allocatable, dimension(:) :: s_rm_seg ! variability around trend, running mean
        real(sp), allocatable, dimension(:) :: s_rm_sim_seg ! variability around trend, running mean, simulation
end module segmented
!
!=======================================================================================================================
!
module setting
        use nrtype
        implicit none
!       Setting.
        save
        character (len=4) :: mis ! input: MIS1, MIS5 or MIS7
        character (len=16) :: infile ! input file
        character (len=12) :: outfile ! output file
        integer :: k_rm=-999 ! x_trend_rm estimation: number of window points for running mean
        integer, allocatable, dimension(:) :: k_seg ! segmentation, number of points per segment
        integer :: n_segmentation=-999 ! segmentation, number of segmentations
        integer, allocatable, dimension(:) :: m ! segmentation, number of segments
        integer :: n=-999 ! data size
!
        real(sp) :: a_rm=-999.0_sp ! equivalent autocorrelation coefficient, residuals from trend 1: running mean
        real(sp) :: tau_rm=-999.0_sp ! persistence time, residuals from trend 1: running mean
        integer :: lopt_rm=-999 ! optimal block length, residuals from trend 1: running mean
end module setting
!
!=======================================================================================================================
!
module statistics
        use nrtype
        implicit none
!       Test statistics.
        save
        real(sp), allocatable, dimension(:) :: r_rm ! correlation, trend-variability, running mean
        real(sp), allocatable, dimension(:) :: r_rm_900 ! correlation, false-alarm level, 90.0%
        real(sp), allocatable, dimension(:) :: r_rm_910 ! correlation, false-alarm level, 91.0%
        real(sp), allocatable, dimension(:) :: r_rm_920 ! correlation, false-alarm level, 92.0%
        real(sp), allocatable, dimension(:) :: r_rm_930 ! correlation, false-alarm level, 93.0%
        real(sp), allocatable, dimension(:) :: r_rm_940 ! correlation, false-alarm level, 94.0%
        real(sp), allocatable, dimension(:) :: r_rm_950 ! correlation, false-alarm level, 95.0%
        real(sp), allocatable, dimension(:) :: r_rm_960 ! correlation, false-alarm level, 96.0%
        real(sp), allocatable, dimension(:) :: r_rm_970 ! correlation, false-alarm level, 97.0%
        real(sp), allocatable, dimension(:) :: r_rm_980 ! correlation, false-alarm level, 98.0%
        real(sp), allocatable, dimension(:) :: r_rm_990 ! correlation, false-alarm level, 99.0%
        real(sp), allocatable, dimension(:) :: r_rm_991 ! correlation, false-alarm level, 99.1%
        real(sp), allocatable, dimension(:) :: r_rm_992 ! correlation, false-alarm level, 99.2%
        real(sp), allocatable, dimension(:) :: r_rm_993 ! correlation, false-alarm level, 99.3%
        real(sp), allocatable, dimension(:) :: r_rm_994 ! correlation, false-alarm level, 99.4%
        real(sp), allocatable, dimension(:) :: r_rm_995 ! correlation, false-alarm level, 99.5%
        real(sp), allocatable, dimension(:) :: r_rm_996 ! correlation, false-alarm level, 99.6%
        real(sp), allocatable, dimension(:) :: r_rm_997 ! correlation, false-alarm level, 99.7%
        real(sp), allocatable, dimension(:) :: r_rm_998 ! correlation, false-alarm level, 99.8%
        real(sp), allocatable, dimension(:) :: r_rm_999 ! correlation, false-alarm level, 99.9%
        real(sp), allocatable, dimension(:,:) :: r_rm_sim ! correlation, trend-variability, running mean, simulation
end module statistics
!
!=======================================================================================================================
!
module trends
        use nrtype
        implicit none
!       Trends.
        save
        real(sp), allocatable, dimension(:) :: x_trend_rm ! running mean
        real(sp), allocatable, dimension(:) :: x_trend_rm_sim ! running mean, simulation
end module trends
!
!=======================================================================================================================
!
program t1
        use datatx
        use parameters, only: n_sim
        use residuals
        use segmented
        use setting
        use statistics
        use trends
        implicit none
!       Main program.
        integer :: i
        integer :: i_segmentation
        integer :: i_sim
        integer :: j
!
! 1.    Initial
!       =======
!
        call info1 ! welcome
        call chsett1 ! changes setting: mis, infile, outfile, k_rm, n
        call allocate1 ! allocates and initializes t, x
        call allocate2 ! allocates and initializes x_trend_rm, e_rm
        call chsett2 ! changes setting: n_segmentation
        call allocate3 ! allocates and initializes k_seg, m, r_rm, false-alarm levels
        call chsett3 ! changes setting: k_seg, m
        call read1 ! reads t, x
!
! 2.    Calculation of test statistics
!       ==============================
!
!
! 2.1   Trend estimation 1: Running mean
!       ===================
!
        call trend_rm(x,n,k_rm,x_trend_rm)
!
!
! 2.2   Trend estimation 2: CHECK
!       ===================
!
!       call trend_rm(x,n,k_rm,x_trend_rm)
!
!
! 2.3   Residuals, trend estimation 1: Running mean
!       =============================
!
        call calc_residuals(n,x,x_trend_rm,e_rm)
!
!
! 2.4   Residuals, trend estimation 2: CHECK
!       =============================
!
!       call calc_residuals(n,x,x_trend_rm,e_rm)
!
!
! 2.5   Persistence time estimation 1: Running-mean residuals
!       =============================
!
        call tauest(t,e_rm,n,'rm',tau_rm)
!
!
! 2.6   Persistence time estimation 2: CHECK
!       =============================
!
!       call tauest(t,e_rm,n,'rm',tau_rm)
!
! test
!
!       print *, 'tau_rm: ',tau_rm
!       open (unit=7, file='e.dat', status='unknown', form='formatted', action='write')
!       do i=1,n
!          write (unit=7, fmt='(2(2x,f12.6))') t(i),e_rm(i)
!       end do
!       close (unit=7, status='keep')
!
!
! 2.7   Segmentation loop, start
!       ========================
!
seg:    do i_segmentation=1,n_segmentation
!
           call allocate4(m(i_segmentation)) ! allocates and initializes x_trend_rm_seg, s_rm_seg
!
!
! 2.8      Segmentation, trend 1: Running mean
!          =====================
!
           call trend_seg(n,x_trend_rm,k_seg(i_segmentation),m(i_segmentation),x_trend_rm_seg)
!
!
! 2.9      Segmentation, trend 2: CHECK
!          =====================
!
!          call trend_seg(n,x_trend_rm,k_seg(i_segmentation),m(i_segmentation),x_trend_rm_seg)
!
!
! 2.10     Segmentation, variability around trend 1: Running mean
!          ========================================
!
           call variability_seg(n,e_rm,k_seg(i_segmentation),m(i_segmentation),s_rm_seg)
!
!
! 2.11     Segmentation, variability around trend 2: CHECK
!          ========================================
!
!          call variability_seg(n,e_rm,k_seg(i_segmentation),m(i_segmentation),s_rm_seg)
!
!
! 2.12     Calculation of test statistic for measures from trend 1: Running mean
!          =======================================================
!
           call calc_statistic(m(i_segmentation),x_trend_rm_seg,s_rm_seg,r_rm(i_segmentation))
!
!
! 2.13     Calculation of test statistic for measures from trend 2: CHECK
!          =======================================================
!
!          call calc_statistic(m(i_segmentation),x_trend_rm_seg,s_rm_seg,r_rm(i_segmentation))
!
!
! 2.14     Final stuff
!          ===========
!
! test
! optional output of segmented trend and variability values
!
          open (unit=1, file='seg_rm.dat', status='unknown', form='formatted', action='write')
          write(unit=1, fmt='(a,2x,a)') 'x_trend_rm_seg','s_rm_seg'
          do j=1,m(i_segmentation)
             write(unit=1, fmt='(2(2x,f12.6))') x_trend_rm_seg(j),s_rm_seg(j)
          end do
          close (unit=1, status='keep')
          print *
          print *
          print *
          print *, 'Currently ',k_seg(i_segmentation),' values per segment'
          print *
          print *, 'xtrend_rm_seg vs s_rm_seg data written into seg_rm.dat'
          read (*,'()')
!
           call deallocate4
!
!
! 2.14  Segmentation loop, end
!       ======================
!
        end do seg
!
!
! 3.    Calculation of null distribution of test statistics
!       ===================================================
!
!
! 3.1   Initial stuff
!       =============
!
        call ranseed ! seed random number generator
        call allocate5 ! allocates and initializes t_rm_sim, x_rm_sim
        call allocate6 ! allocates and initializes x_trend_rm_sim
        call allocate7 ! allocates and initializes e_rm_sim
        call calc_lopt(n,t,tau_rm,a_rm,lopt_rm) ! calculates a_rm, lopt_rm
        call allocate8 ! allocates and initializes r_rm_sim
!
!
! 3.2   Bootstrap simulation loop, start
!       ================================
!
sim:    do i_sim=1,n_sim
!
           if (mod(i_sim,100) .eq. 0)                                                             &
           print '(2(a,i6))', 'Bootstrap simulation #    ',i_sim,'    from    ',n_sim
!
!
! 3.3      Simulation x data, using trend 1: Running mean
!          ================================
!
           call simulation_x(n,x_trend_rm,e_rm,lopt_rm,x_rm_sim)
!
!
! 3.4      Simulation x data, using trend 2: CHECK
!          ================================
!
!          call simulation_x(n,x_trend_rm,e_rm,lopt_rm,x_rm_sim)
!
!
! 3.5      Simulation t data, using trend 1: Running mean
!          ================================
!
           call simulation_t(n,t,t_rm_sim) ! currently: t_rm_sim = t
!
!
! 3.6      Trend estimation 1: Running mean
!          ===================
!
           call trend_rm(x_rm_sim,n,k_rm,x_trend_rm_sim)
!
!
! 3.7      Trend estimation 2: CHECK
!          ===================
!
!          call trend_rm(x_rm_sim,n,k_rm,x_trend_rm_sim)
!
!
! 3.8      Residuals, trend estimation 1: Running mean
!          =============================
!
           call calc_residuals(n,x_rm_sim,x_trend_rm_sim,e_rm_sim)
!
!
! 3.9      Residuals, trend estimation 2: CHECK
!          =============================
!
!          call calc_residuals(n,x_rm_sim,x_trend_rm_sim,e_rm_sim)
!
!
! 3.10     Segmentation (simulations) loop, start
!          ======================================
!
segsim:    do i_segmentation=1,n_segmentation
!
              call allocate9(m(i_segmentation)) ! allocates and initializes x_trend_rm_sim_seg, s_rm_sim_seg
!
!
! 3.11        Segmentation, trend 1: Running mean
!             =====================
!
              call trend_seg(n,x_trend_rm_sim,k_seg(i_segmentation),m(i_segmentation),x_trend_rm_sim_seg)
!
!
! 3.12        Segmentation, trend 2: CHECK
!             =====================
!
!             call trend_seg(n,x_trend_rm_sim,k_seg(i_segmentation),m(i_segmentation),x_trend_rm_sim_seg)
!
!
! 3.13        Segmentation, variability around trend 1: Running mean
!             ========================================
!
              call variability_seg(n,e_rm_sim,k_seg(i_segmentation),m(i_segmentation),s_rm_sim_seg)
!
!
! 3.14        Segmentation, variability around trend 2: CHECK
!             ========================================
!
!             call variability_seg(n,e_rm_sim,k_seg(i_segmentation),m(i_segmentation),s_rm_sim_seg)
!
!
! 3.15        Calculation of test statistic for measures from trend 1: Running mean
!             =======================================================
!
              call calc_statistic(m(i_segmentation),x_trend_rm_sim_seg,s_rm_sim_seg,r_rm_sim(i_segmentation,i_sim))
!
!
! 3.16        Calculation of test statistic for measures from trend 2: CHECK
!             =======================================================
!
!             call calc_statistic(m(i_segmentation),x_trend_rm_sim_seg,s_rm_sim_seg,r_rm_sim(i_segmentation,i_sim))
!
!
! 3.17        Final stuff
!             ===========
!
! test
! optional output of segmented trend and variability values, simulation
!
!             open (unit=1, file='seg_rm_sim.dat', status='unknown', form='formatted', action='write')
!             write(unit=1, fmt='(a,2x,a)') 'x_trend_rm_sim_seg','s_rm_sim_seg'
!             do j=1,m(i_segmentation)
!                write(unit=1, fmt='(2(2x,f12.6))') x_trend_rm_sim_seg(j),s_rm_sim_seg(j)
!             end do
!             close (unit=1, status='keep')
!             print *
!             print *
!             print *
!             print *
!             print *, 'Currently ',k_seg(i_segmentation),' values per segment'
!             print *
!             print *, 'xtrend_rm_sim_seg vs s_rm_sim_seg data written into seg_rm_sim.dat'
!             read (*,'()')
!
              call deallocate9
!
!
! 3.18     Segmentation (simulations) loop, end
!          ====================================
!
           end do segsim
!
!
! 3.19  Bootstrap simulation loop, end
!       ===============================
!
        end do sim
!
!
! 4.    Calculate false-alarm levels
!       ============================
!
!
! 4.1   For measures from trend 1: Running mean
!       =========================
!
        call calc_statistics_fal(n_segmentation,n_sim,r_rm_sim,r_rm_900,r_rm_910,r_rm_920,r_rm_930,r_rm_940,     &
                                                               r_rm_950,r_rm_960,r_rm_970,r_rm_980,r_rm_990,     &
                                                               r_rm_991,r_rm_992,r_rm_993,r_rm_994,r_rm_995,     &
                                                               r_rm_996,r_rm_997,r_rm_998,r_rm_999)
!
!
! 5.    Output
!       ======
!
        call output
!
!
! 6.    Exit
!       ====
!
        call deallocate8
        call deallocate7
        call deallocate6
        call deallocate5
        call deallocate3
        call deallocate2
        call deallocate1
!
end program t1
!
!=======================================================================================================================
!
subroutine allocate1
        use nrtype
        use datatx
        use setting, only: n
        implicit none
!       Allocates and initializes t, x.
        allocate(t(n),x(n))
        t=-999.0_sp
        x=-999.0_sp
end subroutine allocate1
!
!=======================================================================================================================
!
subroutine allocate2
        use nrtype
        use residuals
        use trends
        use setting, only: n
        implicit none
!       Allocates and initializes x_trend_rm, e_rm.
        allocate(x_trend_rm(n))
        x_trend_rm=-999.0_sp
        allocate(e_rm(n))
        e_rm=-999.0_sp
end subroutine allocate2
!
!=======================================================================================================================
!
subroutine allocate3
        use nrtype
        use setting, only: k_seg,m,n_segmentation
        use statistics
        implicit none
!       Allocates and initializes k_seg, m, r_rm, false-alarm levels.
        allocate(k_seg(n_segmentation))
        k_seg=-999
        allocate(m(n_segmentation))
        m=-999
        allocate(r_rm(n_segmentation))
        r_rm=-999.0_sp
        allocate(r_rm_900(n_segmentation))
        r_rm_900=-999.0_sp
        allocate(r_rm_910(n_segmentation))
        r_rm_910=-999.0_sp
        allocate(r_rm_920(n_segmentation))
        r_rm_920=-999.0_sp
        allocate(r_rm_930(n_segmentation))
        r_rm_930=-999.0_sp
        allocate(r_rm_940(n_segmentation))
        r_rm_940=-999.0_sp
        allocate(r_rm_950(n_segmentation))
        r_rm_950=-999.0_sp
        allocate(r_rm_960(n_segmentation))
        r_rm_960=-999.0_sp
        allocate(r_rm_970(n_segmentation))
        r_rm_970=-999.0_sp
        allocate(r_rm_980(n_segmentation))
        r_rm_980=-999.0_sp
        allocate(r_rm_990(n_segmentation))
        r_rm_990=-999.0_sp
        allocate(r_rm_991(n_segmentation))
        r_rm_991=-999.0_sp
        allocate(r_rm_992(n_segmentation))
        r_rm_992=-999.0_sp
        allocate(r_rm_993(n_segmentation))
        r_rm_993=-999.0_sp
        allocate(r_rm_994(n_segmentation))
        r_rm_994=-999.0_sp
        allocate(r_rm_995(n_segmentation))
        r_rm_995=-999.0_sp
        allocate(r_rm_996(n_segmentation))
        r_rm_996=-999.0_sp
        allocate(r_rm_997(n_segmentation))
        r_rm_997=-999.0_sp
        allocate(r_rm_998(n_segmentation))
        r_rm_998=-999.0_sp
        allocate(r_rm_999(n_segmentation))
        r_rm_999=-999.0_sp
end subroutine allocate3
!
!=======================================================================================================================
!
subroutine allocate4(m)
        use nrtype
        use segmented
        implicit none
        integer, intent(in) ::m
!       Allocates and initializes x_trend_rm_seg, s_rm_seg.
        allocate(x_trend_rm_seg(m))
        x_trend_rm_seg=-999.0_sp
        allocate(s_rm_seg(m))
        s_rm_seg=-999.0_sp
end subroutine allocate4
!
!=======================================================================================================================
!
subroutine allocate5
        use nrtype
        use datatx
        use setting, only: n
        implicit none
!       Allocates and initializes t_rm_sim, x_rm_sim.
        integer :: i
        allocate(t_rm_sim(n),x_rm_sim(n))
        t_rm_sim=-999.0_sp
        x_rm_sim=-999.0_sp
end subroutine allocate5
!
!=======================================================================================================================
!
subroutine allocate6
        use nrtype
        use trends
        use setting, only: n
        implicit none
!       Allocates and initializes x_trend_rm_sim.
        integer :: i
        allocate(x_trend_rm_sim(n))
        x_trend_rm_sim=-999.0_sp
end subroutine allocate6
!
!=======================================================================================================================
!
subroutine allocate7
        use nrtype
        use residuals
        use setting, only: n
        implicit none
!       Allocates and initializes e_rm_sim.
        integer :: i
        allocate(e_rm_sim(n))
        e_rm_sim=-999.0_sp
end subroutine allocate7
!
!=======================================================================================================================
!
subroutine allocate8
       use nrtype
       use parameters, only: n_sim
       use setting, only: n_segmentation
       use statistics
       implicit none
!      Allocates and initializes r_rm_sim.
       allocate(r_rm_sim(n_segmentation,n_sim))
       r_rm_sim=-999.0_sp
end subroutine allocate8
!
!=======================================================================================================================
!
subroutine allocate9(m)
        use nrtype
        use segmented
        implicit none
        integer, intent(in) ::m
!       Allocates and initializes x_trend_rm_sim_seg, s_rm_sim_seg.
        allocate(x_trend_rm_sim_seg(m))
        x_trend_rm_sim_seg=-999.0_sp
        allocate(s_rm_sim_seg(m))
        s_rm_sim_seg=-999.0_sp
end subroutine allocate9
!
!=======================================================================================================================
!
subroutine bootstrap1(n,x,l,x_resample)
        use nrtype
        use nr, only: ran
        implicit none
        integer, intent(in) :: n
        real(sp), dimension(n), intent(in)  :: x
        integer, intent(in) :: l
        real(sp), dimension(n), intent(out) :: x_resample
!       MBB.
        integer, dimension(n) :: indxx
        integer :: i=0
        integer :: j=0
        integer :: k=0
        integer :: idum=1
        integer :: n_block_start
        n_block_start=n-l+1
        if (l == 1) then
           do i=1,n
              indxx(i)=int(n*ran(idum))+1
           end do
        else
           k=1
outer:     do i=1,n
              indxx(k)=int(n_block_start*ran(idum))+1
              k=k+1
              if (k > n) exit outer
              do j=1,l-1
                 indxx(k)=indxx(k-1)+1
                 k=k+1
                 if (k > n) exit outer
              end do
           end do outer
        end if
        x_resample(:)=x(indxx(:))
end subroutine bootstrap1
!
!=======================================================================================================================
!
function brent(ax,bx,cx,lstau,tol,xmin,xfunc,yfunc,nfunc)
        use nrtype; use nrutil, only : nrerror
        implicit none
!       Brents minimization (Numerical Recipes, modified: last three arguments).
        integer, intent(in) :: nfunc
        real(sp), dimension(nfunc) :: xfunc
        real(sp), dimension(nfunc) :: yfunc
        real(sp), intent(in) :: ax,bx,cx,tol
        real(sp), intent(out) :: xmin
        real(sp) :: brent
        interface
        function lstau(a,t,x,n)
               use nrtype
               implicit none
               integer, intent(in) :: n
               real(sp), dimension(n), intent(in) :: t
               real(sp), dimension(n), intent(in) :: x
               real(sp), intent(in) :: a
               real(sp) :: lstau
        end function lstau
        end interface
        integer(i4b), parameter :: itmax=100
        real(sp), parameter :: cgold=0.3819660_sp,zeps=1.0e-3_sp*epsilon(ax)
        integer(i4b) :: iter
        real(sp) :: a,b,d,e,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm
        a=min(ax,cx)
        b=max(ax,cx)
        v=bx
        w=v
        x=v
        e=0.0
        fx=lstau(x,xfunc,yfunc,nfunc)
        fv=fx
        fw=fx
        do iter=1,itmax
                xm=0.5_sp*(a+b)
                tol1=tol*abs(x)+zeps
                tol2=2.0_sp*tol1
                if (abs(x-xm) <= (tol2-0.5_sp*(b-a))) then
                        xmin=x
                        brent=fx
                        return
                end if
                if (abs(e) > tol1) then
                        r=(x-w)*(fx-fv)
                        q=(x-v)*(fx-fw)
                        p=(x-v)*q-(x-w)*r
                        q=2.0_sp*(q-r)
                        if (q > 0.0) p=-p
                        q=abs(q)
                        etemp=e
                        e=d
                        if (abs(p) >= abs(0.5_sp*q*etemp) .or. &
                                p <= q*(a-x) .or. p >= q*(b-x)) then
                                e=merge(a-x,b-x, x >= xm )
                                d=cgold*e
                        else
                                d=p/q
                                u=x+d
                                if (u-a < tol2 .or. b-u < tol2) d=sign(tol1,xm-x)
                        end if
                else
                        e=merge(a-x,b-x, x >= xm )
                        d=cgold*e
                end if
                u=merge(x+d,x+sign(tol1,d), abs(d) >= tol1 )
                fu=lstau(u,xfunc,yfunc,nfunc)
                if (fu <= fx) then
                        if (u >= x) then
                                a=x
                        else
                                b=x
                        end if
                        call shft(v,w,x,u)
                        call shft(fv,fw,fx,fu)
                else
                        if (u < x) then
                                a=u
                        else
                                b=u
                        end if
                        if (fu <= fw .or. w == x) then
                                v=w
                                fv=fw
                                w=u
                                fw=fu
                        else if (fu <= fv .or. v == x .or. v == w) then
                                v=u
                                fv=fu
                        end if
                end if
        end do
        call nrerror('brent: exceed maximum iterations')
        contains
!bl
        subroutine shft(a,b,c,d)
        real(sp), intent(out) :: a
        real(sp), intent(inout) :: b,c
        real(sp), intent(in) :: d
        a=b
        b=c
        c=d
        end subroutine shft
end function brent
!
!=======================================================================================================================
!
subroutine calc_lopt(n,t,tau,a,lopt)
        use nrtype
        implicit none
        integer, intent(in) :: n ! number of data points
        real(sp), dimension(n), intent(in) :: t ! time
        real(sp), intent(in) :: tau ! persistence time
        real(sp), intent(out) :: a ! equivalent autocorrelation coefficient
        integer, intent(out) :: lopt ! optimal block length
!
!       Calculates optimal block length after Mudelsee (2010: Eq. 3.28).
!
        real(sp) :: d ! average spacing
        real(sp), parameter :: zero=0.0_sp
        real(sp), parameter :: onethird=0.3333333333333333333333333333333333333333_sp
        real(sp), parameter :: twothird=0.6666666666666666666666666666666666666666_sp
        real(sp), parameter ::  half=0.5_sp
        real(sp), parameter ::   one=1.0_sp
        real(sp), parameter ::   two=2.0_sp
        real(sp), parameter :: six=6.0_sp
        real(sp) :: dummy1
        real(sp) :: dummy2
        real(sp) :: dummy3
!
! 1.    Equivalent autocorrelation coefficient
!       ======================================
!
        d=(t(n)-t(1))/(n-one)
        if (tau/d .le. 0.001_sp) then
           a=zero
        else
           a=exp(-d/tau)
        end if
!
! 2.    Optimal block length
!       ====================
!
        dummy1=six**half*a
        dummy2=one-a**two
        dummy3=(dummy1/dummy2)**twothird*n**onethird
        lopt=anint(dummy3)
        lopt=max(1,lopt)
        lopt=min(n-1,lopt)
!
! test
!       print '(a,i4,a,f8.2,a,f8.2)', 'lopt: ',lopt,'    a:',a,'    tau: ',tau
end subroutine calc_lopt
!
!=======================================================================================================================
!
subroutine calc_residuals(n,x,x_trend,e)
        use nrtype
        implicit none
        integer, intent(in) :: n ! number of data points
        real(sp), dimension(n), intent(in) :: x ! x-value
        real(sp), dimension(n), intent(in) :: x_trend ! trend
        real(sp), dimension(n), intent(out) :: e ! residual
!       Calculates residuals.
        integer :: i
        e=x-x_trend
!
! test
!       do i=1,n
!          print '(a,i4,a,f8.2,a,f8.2,a,f8.2)',                                         &
!                'i = ',i,'    x(i) = ',x(i),'    x_trend(i) = ',x_trend(i),'    e(i) = ',e(i)
!       end do
end subroutine calc_residuals
!
!=======================================================================================================================
!
subroutine calc_statistic(m,x_trend_seg,s_seg,r)
        use nrtype
        implicit none
        integer, intent(in) :: m ! number of segments
        real(sp), dimension(m), intent(in) :: x_trend_seg ! trend, segmented
        real(sp), dimension(m), intent(in) :: s_seg ! variability, segmented
        real(sp), intent(out) :: r ! test statistic = correlation coefficient
!       Calculates test statistic.
        call pearsn(x_trend_seg,s_seg,m,r)
end subroutine calc_statistic
!
!=======================================================================================================================
!
subroutine calc_statistics_fal(n_segmentation,n_sim,r_sim,r_900,r_910,r_920,r_930,r_940,     &
                                                          r_950,r_960,r_970,r_980,r_990,     &
                                                          r_991,r_992,r_993,r_994,r_995,     &
                                                          r_996,r_997,r_998,r_999)
        use nrtype
        use nr, only: sort
        implicit none
        integer, intent(in) :: n_segmentation ! number of segmentations
        integer, intent(in) :: n_sim ! number of simulations
        real(sp), dimension(n_segmentation,n_sim), intent(in) :: r_sim ! simulated test statistic
        real(sp), dimension(n_segmentation), intent(out) :: r_900 ! correlation, false-alarm level, 90.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_910 ! correlation, false-alarm level, 91.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_920 ! correlation, false-alarm level, 92.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_930 ! correlation, false-alarm level, 93.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_940 ! correlation, false-alarm level, 94.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_950 ! correlation, false-alarm level, 95.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_960 ! correlation, false-alarm level, 96.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_970 ! correlation, false-alarm level, 97.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_980 ! correlation, false-alarm level, 98.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_990 ! correlation, false-alarm level, 99.0%
        real(sp), dimension(n_segmentation), intent(out) :: r_991 ! correlation, false-alarm level, 99.1%
        real(sp), dimension(n_segmentation), intent(out) :: r_992 ! correlation, false-alarm level, 99.2%
        real(sp), dimension(n_segmentation), intent(out) :: r_993 ! correlation, false-alarm level, 99.3%
        real(sp), dimension(n_segmentation), intent(out) :: r_994 ! correlation, false-alarm level, 99.4%
        real(sp), dimension(n_segmentation), intent(out) :: r_995 ! correlation, false-alarm level, 99.5%
        real(sp), dimension(n_segmentation), intent(out) :: r_996 ! correlation, false-alarm level, 99.6%
        real(sp), dimension(n_segmentation), intent(out) :: r_997 ! correlation, false-alarm level, 99.7%
        real(sp), dimension(n_segmentation), intent(out) :: r_998 ! correlation, false-alarm level, 99.8%
        real(sp), dimension(n_segmentation), intent(out) :: r_999 ! correlation, false-alarm level, 99.9%
!
!       Calculates false-alarm levels.
!
        integer :: i
        real(sp), dimension(n_sim) :: dumr ! dummy
        real(sp), parameter :: z900=0.900_sp
        real(sp), parameter :: z910=0.910_sp
        real(sp), parameter :: z920=0.920_sp
        real(sp), parameter :: z930=0.930_sp
        real(sp), parameter :: z940=0.940_sp
        real(sp), parameter :: z950=0.950_sp
        real(sp), parameter :: z960=0.960_sp
        real(sp), parameter :: z970=0.970_sp
        real(sp), parameter :: z980=0.980_sp
        real(sp), parameter :: z990=0.990_sp
        real(sp), parameter :: z991=0.991_sp
        real(sp), parameter :: z992=0.992_sp
        real(sp), parameter :: z993=0.993_sp
        real(sp), parameter :: z994=0.994_sp
        real(sp), parameter :: z995=0.995_sp
        real(sp), parameter :: z996=0.996_sp
        real(sp), parameter :: z997=0.997_sp
        real(sp), parameter :: z998=0.998_sp
        real(sp), parameter :: z999=0.999_sp
        real(sp), parameter :: one=1.0_sp
!
        if ((one-z999)*n_sim .le. 10) print *,'Too few simulations---increase n_sim and re-run.'
!
        do i=1,n_segmentation
           dumr(1:n_sim)=r_sim(i,1:n_sim)
           call sort(dumr)
           r_900(i)=dumr(nint(z900*n_sim))
           r_910(i)=dumr(nint(z910*n_sim))
           r_920(i)=dumr(nint(z920*n_sim))
           r_930(i)=dumr(nint(z930*n_sim))
           r_940(i)=dumr(nint(z940*n_sim))
           r_950(i)=dumr(nint(z950*n_sim))
           r_960(i)=dumr(nint(z960*n_sim))
           r_970(i)=dumr(nint(z970*n_sim))
           r_980(i)=dumr(nint(z980*n_sim))
           r_990(i)=dumr(nint(z990*n_sim))
           r_991(i)=dumr(nint(z991*n_sim))
           r_992(i)=dumr(nint(z992*n_sim))
           r_993(i)=dumr(nint(z993*n_sim))
           r_994(i)=dumr(nint(z994*n_sim))
           r_995(i)=dumr(nint(z995*n_sim))
           r_996(i)=dumr(nint(z996*n_sim))
           r_997(i)=dumr(nint(z997*n_sim))
           r_998(i)=dumr(nint(z998*n_sim))
           r_999(i)=dumr(nint(z999*n_sim))
        end do
end subroutine calc_statistics_fal
!
!=======================================================================================================================
!
subroutine chsett1
        use parameters
        use setting
        implicit none
!       Changes setting: mis, infile, outfile, k_rm, n.
        character (len=1) :: answ1
        integer :: i
outer:  do i=1,ibig
           print '(a)', 'Input data:                                          MIS 1    [1] '
           print '(a)', '                                                     MIS 5    [5] '
           print '(a)', '                                                     MIS 7    [7] '
           read (unit=5, fmt='(a1)') answ1
           if (answ1 .eq. '1') then
              mis=mis1
              infile=infile_mis1
              outfile=outfile_mis1
              k_rm=k_rm_mis1
              n=n_mis1
              exit outer
           else if (answ1 .eq. '5') then
              mis=mis5
              infile=infile_mis5
              outfile=outfile_mis5
              k_rm=k_rm_mis5
              n=n_mis5
              exit outer
           else if (answ1 .eq. '7') then
              mis=mis7
              infile=infile_mis7
              outfile=outfile_mis7
              k_rm=k_rm_mis7
              n=n_mis7
              exit outer
           end if
        end do outer
        if (mod(k_rm,2) .eq. 0) then
           print *,' k_rm: ',k_rm, ' running mean # window points is even --- please check.'
           stop
        end if
end subroutine chsett1
!
!=======================================================================================================================
!
subroutine chsett2
        use nrtype
        use parameters
        use setting, only: mis,n_segmentation
        implicit none
!       Changes setting: n_segmentation.
        integer :: k_seg_start
        integer :: k_seg_end
        integer :: k_seg_delta
        if (mis .eq. mis1) then
           k_seg_start=k_seg_start_mis1
           k_seg_end=k_seg_end_mis1
           k_seg_delta=k_seg_delta_mis1
        else if (mis .eq. mis5) then
           k_seg_start=k_seg_start_mis5
           k_seg_end=k_seg_end_mis5
           k_seg_delta=k_seg_delta_mis5
        else if (mis .eq. mis7) then
           k_seg_start=k_seg_start_mis7
           k_seg_end=k_seg_end_mis7
           k_seg_delta=k_seg_delta_mis7
        end if
        n_segmentation=int((k_seg_end-k_seg_start)/k_seg_delta)+1
end subroutine chsett2
!
!=======================================================================================================================
!
subroutine chsett3
        use parameters
        use setting, only: k_seg,m,mis,n,n_segmentation
        implicit none
!       Changes setting: k_seg.
        integer :: i
        if (mis .eq. mis1) then
           do i=1,n_segmentation
              k_seg(i)=k_seg_start_mis1+(i-1)*k_seg_delta_mis1
           end do
        else if (mis .eq. mis5) then
           do i=1,n_segmentation
              k_seg(i)=k_seg_start_mis5+(i-1)*k_seg_delta_mis5
           end do
        else if (mis .eq. mis7) then
           do i=1,n_segmentation
              k_seg(i)=k_seg_start_mis7+(i-1)*k_seg_delta_mis7
           end do
        end if
        do i=1,n_segmentation
           m(i)=ceiling(1.0_sp*n/k_seg(i))
           if (n-(m(i)-1)*k_seg(i) .eq. 1) then
              print *,'Segmentation problem: Last segment has only one point.'
              print *,'Please adjust k_seg values in module parameters.'
              stop
           end if
        end do
end subroutine chsett3
!
!=======================================================================================================================
!
subroutine deallocate1
        use datatx
        implicit none
!       Deallocates t, x.
        deallocate(t)
        deallocate(x)
end subroutine deallocate1
!
!=======================================================================================================================
!
subroutine deallocate2
        use residuals
        use trends
        implicit none
!       Deallocates x_trend_rm, e_rm.
        deallocate(x_trend_rm)
        deallocate(e_rm)
end subroutine deallocate2
!
!=======================================================================================================================
!
subroutine deallocate3
        use setting, only: k_seg, m
        use statistics
        implicit none
!       Deallocates k_seg, m, r_rm, false-alarm levels.
        deallocate(k_seg)
        deallocate(m)
        deallocate(r_rm)
        deallocate(r_rm_900)
        deallocate(r_rm_910)
        deallocate(r_rm_920)
        deallocate(r_rm_930)
        deallocate(r_rm_940)
        deallocate(r_rm_950)
        deallocate(r_rm_960)
        deallocate(r_rm_970)
        deallocate(r_rm_980)
        deallocate(r_rm_990)
        deallocate(r_rm_991)
        deallocate(r_rm_992)
        deallocate(r_rm_993)
        deallocate(r_rm_994)
        deallocate(r_rm_995)
        deallocate(r_rm_996)
        deallocate(r_rm_997)
        deallocate(r_rm_998)
        deallocate(r_rm_999)
end subroutine deallocate3
!
!=======================================================================================================================
!
subroutine deallocate4
        use segmented
        implicit none
!       Deallocates x_trend_rm_seg, s_rm_seg.
        deallocate(x_trend_rm_seg)
        deallocate(s_rm_seg)
end subroutine deallocate4
!
!=======================================================================================================================
!
subroutine deallocate5
        use datatx
        implicit none
!       Deallocates t_rm_sim, x_rm_sim.
        deallocate(t_rm_sim)
        deallocate(x_rm_sim)
end subroutine deallocate5
!
!=======================================================================================================================
!
subroutine deallocate6
        use trends
        implicit none
!       Deallocates x_trend_rm_sim.
        deallocate(x_trend_rm_sim)
end subroutine deallocate6
!
!=======================================================================================================================
!
subroutine deallocate7
        use residuals
        implicit none
!       Deallocates e_rm_sim.
        deallocate(e_rm_sim)
end subroutine deallocate7
!
!=======================================================================================================================
!
subroutine deallocate8
        use statistics
        implicit none
!       Deallocates r_rm_sim.
        deallocate(r_rm_sim)
end subroutine deallocate8
!
!=======================================================================================================================
!
subroutine deallocate9
        use segmented
        implicit none
!       Deallocates x_trend_rm_sim_seg, s_rm_sim_seg.
        deallocate(x_trend_rm_sim_seg)
        deallocate(s_rm_sim_seg)
end subroutine deallocate9
!
!=======================================================================================================================
!
subroutine info1
        implicit none
!       Welcomes.
        integer :: i
        do i=1,100
           print *
        end do
        print '(a)', 'EDC-VAR:                                                          '
        print *
        print '(a)', 'Testing for a relation between trend and variability              '
        print *
        print '(a)', 'Manfred Mudelsee                                                  '
        print '(a)', 'www.mudelsee.com                                                  '
        print '(a)', 'April 2011                                                        '
        print *
        print *
end subroutine info1
!
!=======================================================================================================================
!
function lstau(a,t,x,n)
        use nrtype
        implicit none
        integer, intent(in) :: n
        real(sp), dimension(n), intent(in) :: t
        real(sp), dimension(n), intent(in) :: x
        real(sp), intent(in) :: a
!       Least-squares function for tau estimation.
        real(sp) :: lstau
        integer :: i
        lstau=0.0_sp
        do i=2,n
           lstau=lstau+(x(i)-x(i-1)*sign(1.0_sp,a)*abs(a)**(t(i)-t(i-1)))**2.0_sp
        end do
end function lstau
!
!=======================================================================================================================
!
subroutine minls(a_ar1,n,tol,tol2,t,x,amin,nmu_)
        use nrtype
        implicit none
        integer, intent(in) :: n
        real(sp), dimension(n), intent(in) :: t
        real(sp), dimension(n), intent(in) :: x
        real(sp), intent(in)  :: a_ar1
        real(sp), intent(in)  :: tol
        real(sp), intent(in)  :: tol2
        real(sp), intent(out) :: amin
        integer, intent(out)  :: nmu_
!       Minimizes least-squares function s using Brent's method.
        real(sp) :: dum1=-999.0_sp
        real(sp) :: dum2=-999.0_sp
        real(sp) :: dum3=-999.0_sp
        real(sp) :: dum4=-999.0_sp
        real(sp) :: a_ar11=-999.0_sp
        real(sp) :: a_ar12=-999.0_sp
        real(sp) :: a_ar13=-999.0_sp
        real(sp) :: brent
        interface
        function lstau(a,t,x,n)
               use nrtype
               implicit none
               integer, intent(in) :: n
               real(sp), dimension(n), intent(in) :: t
               real(sp), dimension(n), intent(in) :: x
               real(sp), intent(in) :: a
               real(sp) :: lstau
        end function lstau
        end interface
        nmu_=0
        dum1=brent(-2.0_sp,                a_ar1, +2.0_sp, lstau, tol, a_ar11, t,x,n)
        dum2=brent( a_ar1, 0.5_sp*(a_ar1+1.0_sp), +2.0_sp, lstau, tol, a_ar12, t,x,n)
        dum3=brent(-2.0_sp, 0.5_sp*(a_ar1-1.0_sp),  a_ar1, lstau, tol, a_ar13, t,x,n)
        if ((abs(a_ar12-a_ar11) > tol2 .and. abs(a_ar12-a_ar1) > tol2) .or. &
            (abs(a_ar13-a_ar11) > tol2 .and. abs(a_ar13-a_ar1) > tol2))     &
           nmu_=1
        dum4=min(dum1,dum2,dum3)
        if (dum4 == dum2) then
           amin=a_ar12
        else if (dum4 == dum3) then
           amin=a_ar13
        else
           amin=a_ar11
        end if
end subroutine minls
!
!=======================================================================================================================
!
subroutine output
        use parameters
        use setting
        use statistics
        implicit none
!       Outputs.
        integer :: i
        open (unit=1, file=outfile, status='unknown', form='formatted', action='write')
        write (unit=1, fmt='(a)')         'Output from /EDC-VAR/test/t1.exe                                        '
        write (unit=1, fmt='(a)')         'Author: Manfred Mudelsee (www.mudelsee.com)                             '
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a)')         'EDC ice core dD                                                         '
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a,a)')       'Isotope stage:                                         ',mis
        if (mis .eq. mis1) then
           write (unit=1, fmt='(a)')      'Timescale after Lemieux-Dudon et al. (2010)                             '
        end if
        if (mis .eq. mis5) then
           write (unit=1, fmt='(a)')      'EDC3 timescale shifted                                                  '
           write (unit=1, fmt='(a)')      '(at 112.4 ka by +2 ka and at 134 ka by +1 ka)                           '
        end if
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a,a)')       'Input file:                                            ',infile
        write (unit=1, fmt='(a,i6)')      'Number of data points:                                 ',n
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a)')         'Trend estimation #1: Running mean (RM)                                  '
        write (unit=1, fmt='(a,i6)')      'Number of window points:                               ',k_rm
        write (unit=1, fmt='(a,f0.6)')    'Persistence time of trend residuals (in time units):   ',tau_rm
        write (unit=1, fmt='(a,f0.6)')    'Equivalent autocorrelation coefficient:                ',a_rm
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a)')         'Segmentation for calculating                                            '
        write (unit=1, fmt='(a)')         'test statistic as correlation between                                   '
        write (unit=1, fmt='(a)')         'average trend (within segment) and                                      '
        write (unit=1, fmt='(a)')         'standard deviation (within segment)                                     '
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a,i6)')      'Number of segmentations:                               ',n_segmentation
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a)')         'Calculation of null distribution of the                                 '
        write (unit=1, fmt='(a)')         'test statistic and false-alarm levels:                                  '
        write (unit=1, fmt='(a)')         'MBB resampling                                                          '
        write (unit=1, fmt='(a)')         '(index denotes trend estimation)                                        '
        write (unit=1, fmt='(a,i6)')      'Block length l_RM:                                     ',lopt_rm
        write (unit=1, fmt='(a,i6)')      'Number of simulations:                                 ',n_sim
        write (unit=1, fmt='(a)')         '                                                                        '
        write (unit=1, fmt='(a)')         'Test results                                                            '
        write (unit=1, fmt='(a)')         '(index denotes trend estimation)                                        '
        write (unit=1, fmt='(a,a)')       'Test statistic:                                        ','r'
        write (unit=1, fmt='(a,a)')       'False-alarm levels at P = 90% etc.                     ','r_900 etc.'
        write (unit=1, fmt='(a)') '                                                                                '
        write (unit=1, fmt='(23(a))')                                                                              &
                                        '#_segmentation  ',                                                        &
                                        '#_points_per_segment  ',                                                  &
                                        '#_segments  ',                                                            &
                                        'r_RM  ',                                                                  &
                                        'r_RM_900  ',                                                              &
                                        'r_RM_910  ',                                                              &
                                        'r_RM_920  ',                                                              &
                                        'r_RM_930  ',                                                              &
                                        'r_RM_940  ',                                                              &
                                        'r_RM_950  ',                                                              &
                                        'r_RM_960  ',                                                              &
                                        'r_RM_970  ',                                                              &
                                        'r_RM_980  ',                                                              &
                                        'r_RM_990  ',                                                              &
                                        'r_RM_991  ',                                                              &
                                        'r_RM_992  ',                                                              &
                                        'r_RM_993  ',                                                              &
                                        'r_RM_994  ',                                                              &
                                        'r_RM_995  ',                                                              &
                                        'r_RM_996  ',                                                              &
                                        'r_RM_997  ',                                                              &
                                        'r_RM_998  ',                                                              &
                                        'r_RM_999  '
        do i=1,n_segmentation
           write (unit=1, fmt='(i4,12x,i4,18x,i4,7x,f6.3,1x,18(f8.3,2x),f8.3)')                                    &
                                        i,k_seg(i),m(i),r_rm(i),                                                   &
                                        r_rm_900(i),r_rm_910(i),r_rm_920(i),r_rm_930(i),r_rm_940(i),               &
                                        r_rm_950(i),r_rm_960(i),r_rm_970(i),r_rm_980(i),r_rm_990(i),               &
                                        r_rm_991(i),r_rm_992(i),r_rm_993(i),r_rm_994(i),r_rm_995(i),               &
                                        r_rm_996(i),r_rm_997(i),r_rm_998(i),r_rm_999(i)
        end do
        close (unit=1, status='keep')
end subroutine output
!
!=======================================================================================================================
!
subroutine pearsn(x,y,n,r)
        use nrtype
        use nrutil
        implicit none
!       Pearson's correlation coefficient (Numerical Recipes, modified).
        integer, intent(in) :: n
        real(sp), dimension(n), intent(in) :: x,y
        real(sp), intent(out) :: r
        real(sp), parameter :: tiny=1.0e-08_sp
        real(sp), dimension(n) :: xt,yt
        real(sp) :: ax,ay,sxx,sxy,syy
        ax=sum(x)/n
        ay=sum(y)/n
        xt(:)=x(:)-ax
        yt(:)=y(:)-ay
        sxx=dot_product(xt,xt)
        syy=dot_product(yt,yt)
        sxy=dot_product(xt,yt)
        r=sxy/(sqrt(sxx*syy)+tiny)
end subroutine pearsn
!
!===================================================================================================
!
subroutine ranseed
        use nrtype
        use nr, only: ran
        implicit none
!       Seeds random number generator.
        integer :: seed=0
        real(sp) :: dummy=-999.0_sp
!       print *,'seed random number generator (negative integer): '
!       read (5,*) seed
        seed=-610622867
        dummy=ran(seed)
end subroutine ranseed
!
!=======================================================================================================================
!
subroutine read1
        use nrtype
        use datatx
        use setting, only: infile,n
        implicit none
!       Reads t, x.
        character (len=1) c1
        real(sp) :: tedc3
        real(sp) :: z
        integer :: i
        open (unit=1, file=infile, status='old', form='formatted', action='read')
        read (unit=1, fmt='(a1)') c1 ! header line
        do i=1,n
           read (unit=1, fmt=*) z,tedc3,t(i),x(i) ! depth, EDC3 age, t = EDC3 age shifted, x = dD
        end do
        close (unit=1, status='keep')
end subroutine read1
!
!=======================================================================================================================
!
subroutine rhoest(n,x,rho)
        use nrtype
        implicit none
        integer, intent(in) :: n
        real(sp), dimension(n), intent(in) :: x
        real(sp), intent(out) :: rho
!       Estimates autocorrelation coefficient (equidistant data).
        integer i
        real(sp) :: sum1=0.0_sp
        real(sp) :: sum2=0.0_sp
        do i=2,n
           sum1=sum1+x(i)*x(i-1)
           sum2=sum2+x(i)**2.0_sp
        end do
        rho=sum1/sum2
end subroutine rhoest
!
!=======================================================================================================================
!
subroutine simulation_t(n,t,t_sim)
        use nrtype
        implicit none
        integer, intent(in) :: n ! data size
        real(sp), dimension(n), intent(in) :: t ! time
        real(sp), dimension(n), intent(out) :: t_sim ! simulated t
!
!       Simulates t.
!       Currently: Trend estimation #1 running mean does not require t simulation because
!       it works with the points and not with the times.
        t_sim=t
end subroutine simulation_t
!
!=======================================================================================================================
!
subroutine simulation_x(n,x_trend,e,lopt,x_sim)
        use nrtype
        implicit none
        integer, intent(in) :: n ! data size
        real(sp), dimension(n), intent(in) :: x_trend ! trend
        real(sp), dimension(n), intent(in) :: e ! residuals
        integer, intent(in) :: lopt ! block length
        real(sp), dimension(n), intent(out) :: x_sim ! simulated x (block bootstrap resampling)
!       Simulates x (MBB resampling).
        real(sp), dimension(n) :: e_mbb
        call bootstrap1(n,e,lopt,e_mbb)
        x_sim=x_trend+e_mbb
end subroutine simulation_x
!
!=======================================================================================================================
!
subroutine tauest(t_in,x_in,n,c,tau)
        use nrtype
        use nr, only: avevar
        implicit none
        character (len=2), intent(in) :: c
        integer, intent(in) :: n
        real(sp), dimension(n), intent(in) :: t_in
        real(sp), dimension(n), intent(in) :: x_in
        real(sp), intent(out) :: tau
!
!       Estimates AR(1) persistence time tau using
!       the least-squares algorithm of Mudelsee (2002).
!
!       Notes:
!       (1) assumes t = age (see Point 1)
!       (2) automatic bias correction (see Point 5)
!
        real(sp), dimension(n) :: x
        real(sp), dimension(n) :: t
        real(sp) :: avex=-999.0_sp
        real(sp) :: varx=-999.0_sp
        real(sp) :: delta=-999.0_sp
        real(sp) :: scalt=-999.0_sp
        real(sp) :: rho=-999.0_sp
        real(sp) :: rho_non=-999.0_sp
        real(sp), parameter :: e_inv=0.3678794_sp
        real(sp), parameter :: tol =3.0e-08_sp
        real(sp), parameter :: tol2=1.0e-06_sp
        real(sp) :: amin=-999.0_sp
        integer :: mult=-999
        integer :: i=-999
!
! 1.    Rename and change time direction
!       ================================
!
        do i=n,1,-1
           t(i)=-1.0_sp*t_in(n+1-i)
           x(i)=+1.0_sp*x_in(n+1-i)
        end do
!
! 2.    Scaling (x)
!       ===========
!
        call avevar(x,avex,varx)
        x=x/sqrt(varx)
!
! 3.    Scaling (t)
!       ===========
!
!       => start value of a = 1/e
!
        delta=abs(t(n)-t(1))/(n-1)
        call rhoest(n,x,rho)
        if (rho <= 0.0_sp) then
           rho=0.05_sp
        else if (rho >= 1.0_sp) then
           rho=0.95_sp
        end if
        scalt=-1.0_sp*log(rho)/delta
        t=t*scalt
!
! 4.    Estimation (tau)
!       ================
!
        call minls(e_inv,n,tol,tol2,t,x,amin,mult)
!
! 5.    Result
!       ======
!
        if (mult == 1 .or. amin <= 0.0_sp .or. amin >= 1.0_sp) then
           if (mult == 1) then
              print '(a,a)', 'data ',c
              print *, 'Estimation problem (tau): LS function has > 1 minima'
              read (*,'()')
           end if
           if (amin <= 0.0_sp) then
              print '(a,a)', 'data ',c
              print *, 'Estimation problem (tau): a_min =< 0'
              read (*,'()')
           else if (amin >= 1.0_sp) then
              print '(a,a)', 'data ',c
              print *, 'Estimation problem (tau): a_min >= 1'
              read (*,'()')
           end if
           print *, 'Input tau per hand:'
           read (5,*) tau
           tau=tau*1.0_sp
        else
!
!       Bias correction and rescaling
!
           rho_non=1.0_sp*amin**(delta*scalt)
           rho_non=rho_non+(1.0_sp+3.0_sp*rho_non)/(n-1)
           amin=1.0_sp*exp(log(rho_non)/(delta*scalt))
           tau=-1.0_sp/(scalt*log(amin))
        end if
end subroutine tauest
!
!=======================================================================================================================
!
subroutine trend_rm(x,n,k,y)
        use nrtype
        implicit none
        integer, intent(in) :: n ! number of data points
        real(sp), dimension(n), intent(in) :: x ! x-value
        integer, intent(in) :: k ! number of window points
        real(sp), dimension(n), intent(out) :: y ! running mean
!
!       Trend estimation via running-mean smoothing:
!       k window points and boundary kernel.
!
        integer :: i
!
! 1.    Left boundary
!       =============
!
        do i=1,(k-1)/2
! test
!          print *, 'counter left boundary: ',i
           y(i)=1.0_sp*sum(x(1:2*i-1))/(2*i-1)
        end do
!
! 2.    Centre
!       ======
!
        do i=(k+1)/2,n-(k-1)/2
! test
!          print *, 'counter centre: ',i
           y(i)=1.0_sp*sum(x(i-(k-1)/2:i+(k-1)/2))/k
        end do
!
! 3.    Right boundary
!       ==============
!
        do i=n-(k-1)/2+1,n
! test
!          print *, 'counter right boundary: ',i
           y(i)=1.0_sp*sum(x(2*i-n:n))/(2*n-2*i+1)
        end do
!
! test
!       do i=1,n
!          print '(a,i4,a,i4,a,f8.2,a,f8.2)', 'i = ',i,'    k = ',k,'    x(i) = ',x(i),'    y(i) = ',y(i)
!       end do
!       read (*,'()')
end subroutine trend_rm
!
!=======================================================================================================================
!
subroutine trend_seg(n,x_trend,k,m,x_trend_seg)
        use nrtype
        implicit none
        integer, intent(in) :: n ! number of data points
        real(sp), dimension(n), intent(in) :: x_trend ! trend
        integer, intent(in) :: k ! number of points per segment
        integer, intent(in) :: m ! number of segments
        real(sp), dimension(m), intent(out) :: x_trend_seg ! segmented trend
!       Segmentation, trend.
        integer :: j
        do j=1,m-1
           x_trend_seg(j)=sum(x_trend((j-1)*k+1:j*k))/k
        end do
        x_trend_seg(m)=sum(x_trend((m-1)*k+1:n))/(n-(m-1)*k)
!
! test
!       do j=1,m
!          print '(a,i4,a,i4,a,i4,a,f8.2)', 'j = ',j,'    m = ',m,'    k_seg = ',k,'    x_trend_seg(j) = ',x_trend_seg(j)
!       end do
!       read (*,'()')
end subroutine trend_seg
!
!=======================================================================================================================
!
subroutine variability_seg(n,e,k,m,s_seg)
        use nrtype
        use nr, only: avevar
        implicit none
        integer, intent(in) :: n ! number of data points
        real(sp), dimension(n), intent(in) :: e ! residual
        integer, intent(in) :: k ! number of points per segment
        integer, intent(in) :: m ! number of segments
        real(sp), dimension(m), intent(out) :: s_seg ! variability, calculated as standard deviation over
                                                     ! residuals within segment
!       Segmentation, variability.
        integer :: j
        real(sp) :: avex,varx
        do j=1,m-1
           call avevar(e((j-1)*k+1:j*k),avex,varx)
           s_seg(j)=sqrt(varx)
        end do
        call avevar(e((m-1)*k+1:n),avex,varx)
        s_seg(m)=sqrt(varx)
!
! test
!       do j=1,m
!          print '(a,i4,a,i4,a,i4,a,f8.2)', 'j = ',j,'    m = ',m,'    k_seg = ',k,'    s_seg(j) = ',s_seg(j)
!       end do
!       read (*,'()')
end subroutine variability_seg
!
!=======================================================================================================================
!
! Used Numerical Recipes files on the computer:
!
include 'd:\nr7790\avevar.f90'
include 'd:\nr7790\ran.f90'
include 'd:\nr7790\sort.f90'

