module math_mod
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: The module contains all the mathematical SUBROUTINE used in the PEM
!!!           
!!! Author: Adapted from Schorgofer MSIM (N.S, Icarus 2010), impletented here by LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

implicit none

!=======================================================================
  contains
!=======================================================================

SUBROUTINE deriv1(z,nz,y,y0,ybot,dzY)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the first derivative of a function y(z) on an irregular grid 
!!!           
!!! Author: From N.S (N.S, Icarus 2010), impletented here by LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! first derivative of a function y(z) on irregular grid
! upper boundary conditions: y(0) = y0
! lower boundary condition.: yp = ybottom

implicit none

integer, intent(in) :: nz      ! number of layer
real,    intent(in) :: z(nz)   ! depth layer
real,    intent(in) :: y(nz)   ! function which needs to be derived
real,    intent(in) :: y0,ybot ! boundary conditions
real, intent(out) :: dzY(nz) ! derivative of y w.r.t depth
! local
integer :: j
real    :: hm, hp, c1, c2, c3

hp = z(2) - z(1)
c1 = z(1)/(hp*z(2))
c2 = 1/z(1) - 1/(z(2) - z(1))
c3 = -hp/(z(1)*z(2))
dzY(1) = c1*y(2) + c2*y(1) + c3*y0
do j = 2,nz - 1
    hp = z(j + 1) - z(j)
    hm = z(j) - z(j - 1)
    c1 = +hm/(hp*(z(j + 1) - z(j - 1)))
    c2 = 1/hm - 1/hp
    c3 = -hp/(hm*(z(j + 1) - z(j - 1)))
    dzY(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
enddo
dzY(nz) = (ybot - y(nz - 1))/(2.*(z(nz) - z(nz - 1)))

END SUBROUTINE deriv1

!=======================================================================

SUBROUTINE deriv2_simple(z,nz,y,y0,yNp1,yp2)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the second derivative of a function y(z) on an irregular grid 
!!!           
!!! Author: N.S (raw copy/paste from MSIM)
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! second derivative y_zz on irregular grid
! boundary conditions: y(0) = y0, y(nz + 1) = yNp1

implicit none

integer, intent(in) :: nz
real,    intent(in) :: z(nz), y(nz), y0, yNp1
real, intent(out) :: yp2(nz)
integer :: j
real    :: hm, hp, c1, c2, c3

c1 = +2./((z(2) - z(1))*z(2))
c2 = -2./((z(2) - z(1))*z(1))
c3 = +2./(z(1)*z(2))
yp2(1) = c1*y(2) + c2*y(1) + c3*y0
do j = 2,nz - 1
    hp = z(j + 1) - z(j)
    hm = z(j) - z(j - 1)
    c1 = +2./(hp*(z(j + 1) - z(j - 1)))
    c2 = -2./(hp*hm)
    c3 = +2./(hm*(z(j + 1) - z(j-1)))
    yp2(j) = c1*y(j + 1) + c2*y(j) + c3*y(j - 1)
enddo
yp2(nz) = (yNp1 - 2*y(nz) + y(nz - 1))/(z(nz) - z(nz - 1))**2

END SUBROUTINE deriv2_simple

!=======================================================================

SUBROUTINE  deriv1_onesided(j,z,nz,y,dy_zj)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: First derivative of function y(z) at z(j)  one-sided derivative on irregular grid
!!!           
!!! Author: N.S (raw copy/paste from MSIM)
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  implicit none

integer, intent(in) :: nz, j
real,    intent(in) :: z(nz), y(nz)
real, intent(out) :: dy_zj
real :: h1, h2, c1, c2, c3

if (j < 1 .or. j > nz - 2) then
    dy_zj = -1.
else
    h1 = z(j + 1) - z(j)
    h2 = z(j + 2)- z(j + 1)
    c1 = -(2*h1 + h2)/(h1*(h1 + h2))
    c2 = (h1 + h2)/(h1*h2)
    c3 = -h1/(h2*(h1 + h2))
    dy_zj = c1*y(j) + c2*y(j + 1) + c3*y(j + 2)
endif

END SUBROUTINE deriv1_onesided

!=======================================================================

PURE SUBROUTINE colint(y,z,nz,i1,i2,integral)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose:  Column integrates y on irregular grid
!!!           
!!! Author: N.S (raw copy/paste from MSIM)
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
implicit none
  
integer, intent(in) :: nz, i1, i2
real,    intent(in) :: y(nz), z(nz)
real, intent(out) :: integral
integer i
real dz(nz)

dz(1) = (z(2) - 0.)/2
do i = 2,nz - 1
    dz(i) = (z(i + 1) - z(i - 1))/2.
enddo
dz(nz) = z(nz) - z(nz - 1)
integral = sum(y(i1:i2)*dz(i1:i2))

END SUBROUTINE colint

!=======================================================================

SUBROUTINE findroot(y1,y2,z1,z2,zr)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!! Purpose: Compute the root zr, between two values y1 and y2 at depth z1,z2
!!! 
!!! Author: LL
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

implicit none

real, intent(in) :: y1, y2 ! difference between surface water density and at depth  [kg/m^3]
real, intent(in) :: z1, z2 ! depth [m]
real, intent(out) :: zr    ! depth at which we have zero

zr = (y1*z2 - y2*z1)/(y1 - y2)

END SUBROUTINE findroot

end module math_mod