!====================================================================
! STL Decomposition (simplified version, additive model)
! Input: monthly time series Y(t), t = 1..n
! Output: T(t) - trend, S(t) - seasonal, R(t) - residual
! Note: This implementation assumes fixed seasonal period (e.g., 12)
! and uses simple moving average instead of full LOESS.
!====================================================================

program stl_decomposition
  Use NetCDF
  implicit none
  !integer, parameter :: n = 240            ! length of series (e.g., 20 years monthly)
  integer, parameter :: p = 12             ! seasonal period
  integer :: nt, n, i, j, m, iter, num_iter, nlon,nlat,ntime
  !real :: T(n), S(n), R(n)
  !real :: Y_deseason(n), Y_detrend(n)
  real :: seasonal_cycle(p)
  Character (Len=255) :: Fname
  Real, Allocatable :: X(:,:,:), Y(:), T(:), S(:), R(:), Y_deseason(:), Y_detrend(:)
  Real, Allocatable, dimension(:) :: lon,lat,time
  integer :: it1, it2
  Logical :: status

  !-------------------------------
  ! Load or initialize your time series
  !-------------------------------
  Fname="/home/tanpv/Data/SEACLID/ClimIDX_2022/Clim_IDX_Longtime_2022.nc"
  Call rd_ncdf_3d (Fname,lon,lat,time,nlon,nlat,ntime,X,"dmi")
  nt = (2022-1850+1)*12
  n  = (2021-1870)*12
  it1=241; it2=it1+n-1
  print*,"ntime,nt,n: ",ntime,nt,n
  Allocate (Y(n), T(n), S(n), R(n), Y_deseason(n), Y_detrend(n))
  Y(1:n) = X(1,1,it1:it2)
  !do i = 1, n
  !   Y(i) = sin(2.0 * 3.1416 * i / 12.0) + 0.01*i + 0.2 * cos(0.5*i)  ! example synthetic series
  !end do

  !-------------------------------
  ! STL Parameters
  !-------------------------------
  num_iter = 1000

  ! Initialize trend and seasonal
  T(:) = 0.0
  S(:) = 0.0

  !-------------------------------
  ! Begin outer iterations
  !-------------------------------
  do iter = 1, num_iter

     ! Step 1: Detrend
     do i = 1, n
        Y_detrend(i) = Y(i) - T(i)
     end do

     ! Step 2: Seasonal smoothing (monthly averages)
     seasonal_cycle(:) = 0.0
     do m = 1, p
        j = 0
        do i = m, n, p
           seasonal_cycle(m) = seasonal_cycle(m) + Y_detrend(i)
           j = j + 1
        end do
        if (j > 0) seasonal_cycle(m) = seasonal_cycle(m) / j
     end do

     do i = 1, n
        S(i) = seasonal_cycle(mod(i-1, p)+1)
     end do

     ! Step 3: Deseasonalize
     do i = 1, n
        Y_deseason(i) = Y(i) - S(i)
     end do

     ! Step 4: Trend smoothing (moving average)
     do i = 1, n
        T(i) = 0.0
        j = 0
        do m = max(1,i-6), min(n,i+6)  ! 13-point moving average
           T(i) = T(i) + Y_deseason(m)
           j = j + 1
        end do
        if (j > 0) T(i) = T(i) / j
     end do

     ! Step 5: Update residual
     do i = 1, n
        R(i) = Y(i) - T(i) - S(i)
     end do

  end do

  !-------------------------------
  ! Output results (can be written to file)
  !-------------------------------
  print *, "t", "Y", "T", "S", "R"
  open (3,file="Decomp.txt",status="unknown")
  do i = 1, n
     write (3,*) i, Y(i), T(i), S(i), R(i)
  end do
  close (3)
  CONTAINS
  Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90"
  Include "/home/tanpv/Mylibs/sub_wr_ncdf_Tan.F90"
end program stl_decomposition