MODULE SALINITY

use NUMERICS, only : PROGONKA, CHECK_PROGONKA, IND_STAB_FACT_DB
use NUMERIC_PARAMS, only : vector_length
use LAKE_DATATYPES, only : ireals, iintegers
use DZETA_MOD, only : VARMEAN

contains
SUBROUTINE S_DIFF(dt,dhilow)

use ARRAYS, only : &
& deepice, water, &
& water_salinity_indic, soil_salinity_indic, &
& nstep
use ARRAYS_WATERSTATE, only : lamsal, Sal1, Sal2
use ARRAYS_BATHYM, only : h1, l1, dhw0, dhw, bathymwater
use ARRAYS_GRID, only : ddz, nsoilcols
use ARRAYS_SOIL, only : wsoil, Sals1, Sals2,rosdry,por,dzs

use DRIVING_PARAMS, only : sedim, soilswitch, M, ns, salsoil

use ATMOS, only: &
& Sflux0

use PHYS_FUNC, only: &
& w_sedim

use PHYS_CONSTANTS, only : &
& roi_d_row0, &
& salice0, &
& row0

use T_SOLVER_MOD, only : &
& DIFF_COEF


implicit none
      
! S_diff solves salinity (mineralization) diffusion equation  

! Input variables
real(kind=ireals), intent(in) :: dt 
real(kind=ireals), intent(in) :: dhilow

! Local variables
real(kind=ireals) :: Sflux1,Sflux_soil_bot,x,y,z
real(kind=ireals), dimension(1:vector_length):: a,b,c,d,Sal
      
a(:) = 0.; b(:) = 0.; c(:) = 0.; d(:) = 0

! Sflux1 --- salinity flux at the bottom boundary       
! Sflux0 --- salinity flux at the top boundary       

Sflux1 = 0.d0
Sflux_soil_bot = 0.d0
if (l1 > 0) then
! Salinity flux to water = flux from atmosphere + flux from melting/freezing of saline ice
!  if (dhilow < 0.) then
!    x = max(dhilow,-l1)
!  else
!    x = dhilow
!  endif
  Sflux0 = Sflux0 - dhilow/dt*roi_d_row0*salice0
endif

! It is allowed for atmospheric aerosol to immediately go to 
! water instead of first contaminating the snow cover.
! if (ice == 1) Sflux0 = 0.d0

! It is assumed, that the water freezing to ice, melting from ice and snow,
! the rain are all freshwater sinks/sources. It is not true in fact, as far as
! the acid rain may occur and contaminated snow thaw. For these cases
! the present scheme must be updated.

! Defines, if gravitational sedimentation of tracer is taken into account
! sedim = 1 it is taken into account
! sedim = 0 it is neglected
! sedim = 0

if (water == 1) then
  call DIFF_COEF(a,b,c,d,2,M,2,M,water_salinity_indic,dt)
  x = 0.5*( - bathymwater(1)%area_half*lamsal(1)/(h1*ddz(1)* &
  & bathymwater(1)%area_int) + dhw0/(2.d0*dt) )
  c(1)   = x - ddz(1)*h1/(2.d0*dt) ! - dhw0/(2*dt) is wrong sign!
  b(1)   = x
  d(1)   = - ddz(1)*h1*Sal1(1)/(2.d0*dt) - Sflux0 + x*(Sal1(2) - Sal1(1))
  if (deepice == 1) then
!   Case water,deepice and soil; upper ice and snow are allowed       
!   Salinity diffusion in water       
    x = 0.5*( - bathymwater(M)%area_half*lamsal(M)/(ddz(M)*h1* &
    & bathymwater(M+1)%area_int) + (dhw-dhw0)/(2.d0*dt) )
    c(M+1) = x - ddz(M)*h1/(2.d0*dt) 
    a(M+1) = x
    d(M+1) = - Sal1(M+1)*ddz(M)*h1/(2.d0*dt) + Sflux1 - x*(Sal1(M+1) - Sal1(M))
    call PROGONKA (a,b,c,d,Sal,1,M+1)
    Sal (1:M+1) = max(Sal(1:M+1),0._ireals)
    Sal2(1:M+1) = Sal(1:M+1)
    if (soilswitch%par == 1 .and. salsoil%par == 1) then
!     Salinity diffusion in soil
      call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
      x = 0.5*dt*wsoil(1)/dzs(1)
      c(1) = - 1.d0 - x
      b(1) = x
      d(1) = - Sals1(1,nsoilcols) + x*(Sals1(2,nsoilcols) + Sals1(1,nsoilcols))
      x = 0.5*wsoil(ns-1)*dt/dzs(ns-1)
      c(ns) = - 1.d0 + x
      a(ns) = - x
      d(ns) = - Sals1(ns,nsoilcols) + 2.d0*dt*Sflux_soil_bot/dzs(ns-1) - &
      & x*(Sals1(ns,nsoilcols) + Sals1(ns-1,nsoilcols))
      !print*, CHECK_PROGONKA(ns,a,b,c,d,Sal)
      call PROGONKA (a,b,c,d,Sal,1,ns)
      Sal(1:ns) = max(Sal(1:ns),0._ireals)
      Sals2(1:ns,nsoilcols) = Sal(1:ns)
    else
      Sals2(:,nsoilcols) = Sals1(:,nsoilcols)
    endif
  else
!   Case water, soil; upper ice and snow are allowed       
    if (soilswitch%par == 1 .and. salsoil%par == 1) then
!------WATER-SOIL INTERFACE-------------------------
      x = 0.5*( 0.5*(dhw-dhw0)/dt - bathymwater(M)%area_half*lamsal(M)/ &
      & (ddz(M)*h1*bathymwater(M+1)%area_int) )
      y = rosdry(1)*(1 - por(1))/row0
      z = 0.25*wsoil(1)
      a(M+1) = x
      b(M+1) = z*y
      c(M+1) = - ( 0.5*(dzs(1) + ddz(M)*h1)/dt + z - x ) 
      d(M+1) = - Sal1(M+1)*0.5*(dzs(1) + ddz(M)*h1)/dt - &
      & x*(Sal1(M+1) - Sal1(M)) + z*Sal1(1) + z*y*Sals1(2,nsoilcols)
      call DIFF_COEF(a,b,c,d,2,ns-1,M+2,M+ns-1,soil_salinity_indic,dt)
      x = 0.5*wsoil(ns-1)*dt/(dzs(ns-1))
      c(M+ns) = - 1.d0 + x
      a(M+ns) = - x
      d(M+ns) = - Sals1(ns,nsoilcols) + 2*dt*Sflux_soil_bot/dzs(ns-1) - &
      & x*(Sals1(ns,nsoilcols) + Sals1(ns-1,nsoilcols))
      !print*, CHECK_PROGONKA(M+ns,a,b,c,d,Sal)
      !print*, 'a', a(1:M+ns)
      !read*
      !print*, 'b', b(1:M+ns)
      !read*
      !print*, 'c', c(1:M+ns)
      !read*
      !print*, 'd', d(1:M+ns)
      !read*
      call PROGONKA (a,b,c,d,Sal,1,M+ns)
      Sal(1:M+ns) = max(Sal(1:M+ns),0._ireals)
      Sals2(2:ns,nsoilcols) = Sal(M+2:M+ns)
      Sals2(1,nsoilcols) = Sal(M+1)/y !mind y!
      Sal2(1:M+1) = Sal(1:M+1)
    else ! Zero salinity flux at the bottom
      x = 0.5*( - bathymwater(M)%area_half*lamsal(M)/(ddz(M)*h1* &
      & bathymwater(M+1)%area_int) + (dhw-dhw0)/(2.d0*dt) )
      c(M+1) = x - ddz(M)*h1/(2.d0*dt) 
      a(M+1) = x
      d(M+1) = - Sal1(M+1)*ddz(M)*h1/(2.d0*dt) + Sflux1 - x*(Sal1(M+1) - Sal1(M))
      call PROGONKA (a,b,c,d,Sal,1,M+1)
      Sal (1:M+1) = max(Sal(1:M+1),0._ireals)
      Sal2(1:M+1) = Sal(1:M+1)
    endif
  endif
  !if (sedim%par == 1) call SAL_SEDIM(ddz,h1,dt,Sal2)
else
! Case of the soil and the ice above     
  if (soilswitch%par == 1 .and. salsoil%par == 1) then
    call DIFF_COEF(a,b,c,d,2,ns-1,2,ns-1,soil_salinity_indic,dt)
    x = 0.5*dt*wsoil(1)/dzs(1)
    c(1) = - 1.d0 - x
    b(1) = x
    d(1) = - Sals1(1,nsoilcols) + x*(Sals1(1,nsoilcols) + Sals1(2,nsoilcols))
    x = 0.5*wsoil(ns-1)*dt/(dzs(ns-1))
    c(ns) = - 1.d0 + x
    a(ns) = - x
    d(ns) = - Sals1(ns,nsoilcols) + 2.d0*dt*Sflux_soil_bot/dzs(ns-1) - &
    & x*(Sals1(ns,nsoilcols) + Sals1(ns-1,nsoilcols))
    !print*, CHECK_PROGONKA(ns,a,b,c,d,Sal)
    call PROGONKA (a,b,c,d,Sal,1,ns)
    Sal(1:ns) = max(Sal(1:ns),0._ireals)
    Sals2(1:ns,nsoilcols) = Sal(1:ns)
  else
    Sals2(:,nsoilcols) = Sals1(:,nsoilcols)
  endif
endif 

!print*, VARMEAN(Sal2,bathymwater,1_iintegers)

!print*, 'Sals2', Sals2(:,nsoilcols)
!read*
!print*, 'Sals1', Sals1
!read*
!print*, 'Sal2', Sal2
!read*
!print*, 'Sal1', Sal1
!read*

END SUBROUTINE S_DIFF


!SUBROUTINE SAL_SEDIM(ddz,h1,dt,Sal)
!
!! The subroutine SAL_SEDIM updates the salinity profile
!! due to gravitational sedimentation
!
!use DRIVING_PARAMS, only : M
!use PHYS_FUNC!, only: &
!!& W_SEDIM
!
!implicit none
!
!! Input variables
!real(kind=ireals), intent(in) :: ddz  (M) ! Spacing of dzeta-coordinate grid
!      
!real(kind=ireals), intent(in) :: h1 ! Lake depth, m
!real(kind=ireals), intent(in) :: dt ! Timestep,   sec
!
!! Input/output variables
!real(kind=ireals), intent(inout) :: Sal (M+1) ! Salinity at main levels, kg/kg
!
!! Local variables
!real(kind=ireals) :: a(M+1) 
!real(kind=ireals) :: b(M+1) 
!real(kind=ireals) :: c(M+1) 
!real(kind=ireals) :: d(M+1) 
!
!real(kind=ireals) :: w_sediment (M) !Speed of gravitational sedimentation, m/s
!
!real(kind=ireals) x ! Help variable
!
!integer i ! loop index
!
!logical indstab
!logical ind_bound
!
!! The speed of gravitational sedimentation, positive downwards
!do i = 1, M
!  w_sediment(i) = W_SEDIM()
!enddo
!
!! Boundary conditions at the top boundary (dzeta = 0)
!x = 0.5d0*ddz(1)*h1/dt
!c(1) = - 0.5d0*w_sediment(1) - x
!b(1) =   0.5d0*w_sediment(1)
!d(1) = - Sal(1)*x
!
!! Boundary conditions at the bottom boundary (dzeta = 1)
!x = 0.5d0*ddz(M)*h1/dt
!a(M+1) = - 0.5d0*w_sediment(M)
!c(M+1) =   0.5d0*w_sediment(M) - x
!d(M+1) = - Sal(M+1)*x
!
!! The coefficients of tridiagonal matrix
!do i = 2, M
!  x = 0.5d0*(ddz(i-1) + ddz(i))*h1/dt
!  a(i) = - 0.5d0*w_sediment(i-1)
!  c(i) =   0.5d0*w_sediment(i-1) - 0.5d0*w_sediment(i) - x
!  b(i) =   0.5d0*w_sediment(i)
!  d(i) = - Sal(i)*x
!enddo
!  
!ind_bound = .true.
!call IND_STAB_FACT_DB(a,b,c,1,M+1,indstab,ind_bound)
!call PROGONKA(a,b,c,d,Sal,1,M+1)
!if (.not.indstab) then
!  print*,'Info: Unstable factorization method in SAL_SEDIM'
!  print*,'The accuracy flag is', CHECK_PROGONKA(M+1,a,b,c,d,Sal)
!  if (.not.CHECK_PROGONKA(M+1,a,b,c,d,Sal) ) STOP
!endif
!
!RETURN
!END SUBROUTINE SAL_SEDIM

END MODULE SALINITY
