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
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        implicit none

! first derivative of a function y(z) on irregular grid
! upper boundary conditions: y(0)=y0
! lower boundary condition.: yp = ybottom
  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)))
 return
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
  return
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
  return
end subroutine deriv1_onesided


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)
        RETURN
        end subroutine findroot

end module math_mod