module base_cstret_mod

  implicit none
  integer                           , private :: nz, nt
  real                              , private :: t0, dt, eps
  integer, dimension(:), allocatable, private :: x
  logical, dimension(:), allocatable, private :: m
  real   , dimension(:), allocatable, private :: w, diag, offd, d, o
contains
  
  integer function base_cstret_init(nz_in, nt_in, t0_in, dt_in, eps_in) result(stat)
    integer :: nz_in, nt_in
    real    :: t0_in, dt_in, eps_in
    stat = 1
    nz = nz_in; nt = nt_in
    t0 = t0_in; dt = dt_in; eps = eps_in

write(0,*)"nz, nt, t0, dt, eps", nz, nt, t0, dt, eps

    allocate (x(nz), m(nz), w(nz), diag(nt), offd(nt), d(nt), o(nt))
    stat = 0
  end function

  integer function base_cstret_clos() result(stat)
    stat = 1
    deallocate(x, m, w, diag, offd, d, o)
    stat = 0
  end function
  
  integer function base_cstret_oper(coord, ord, mod) result(stat)
    real   , dimension(:) :: coord
    complex, dimension(:) :: ord
    complex, dimension(:) :: mod

    integer               :: iz, ix, i1, i2, k
    real                  :: rx, w1, w2
    complex               :: t

    stat = 1

    !!.. stretch define
    diag = 2.*eps; offd = -eps  !!.. regularization

    !!.. stretch define
    do iz = 1, nz
       rx = (coord(iz) - t0)/dt; ix = rx; rx = rx - ix
       if (ix < 0 .or. ix > nt - 2) then
          m(iz) = .true.; cycle
       end if
       x(iz) = ix; m(iz) = .false.; w(iz) = rx
       i1 = ix + 1; i2 = i1 + 1
       w2 = rx; w1 = 1. - w2
       diag(i1) = diag(i1) + w1 * w1
       diag(i2) = diag(i2) + w2 * w2
       offd(i1) = offd(i1) + w1 * w2
    end do

    !!.. stretch apply    
    mod = 0.
    do iz = 1, nz; if (m(iz)) cycle
       i1 = x(iz) + 1; i2 = i1 + 1
       w2 = w(iz); w1 = 1. - w2
       mod(i1) = mod(i1) + w1 * ord(iz)
       mod(i2) = mod(i2) + w2 * ord(iz)
    end do

    !!.. tridiagonal solver
    d(1) = diag(1)
    do k = 2, nt
       t = offd(k-1) ; o(k-1) = t/diag(k-1); d(k) = diag(k) - t*o(k-1)
    end do
    do k = 2, nt
       mod(k) = mod(k) - o(k-1)*mod(k-1)
    end do
    mod(nt) = mod(nt)/d(nt)
    do k = nt-1, 1, -1
       mod(k) = mod(k)/d(k)-o(k)*mod(k+1)
    end do

    stat = 0
  end function



end module
