MODULE BATHYMSUBR

contains
!> Subroutine BATHYMETRY calculates lake bathymetry characteristics
SUBROUTINE BATHYMETRY(M,Mice,ns,ix,iy,ndatamax,month,day,lakeform,hour,dt, &
                     & depth_area,area_lake,cellipt, &
                     & multsoil,trib_inflow,dhwtrib,vol,botar)

! Subroutine updates bathymetry variables

use LAKE_DATATYPES, only : &
& ireals, iintegers

use DRIVING_PARAMS, only : soilcolconjtype

use PHYS_CONSTANTS, only : g, row0

use ARRAYS, only : init

use ARRAYS_WATERSTATE, only: preswat

use ARRAYS_BATHYM, only : &
& area_int, area_half, &
& bathymwater, bathymice, &
& bathymdice, bathymsoil, &
& h1,l1,ls1,Lx,Ly, &
& form_ellipse, form_rectangle

use ARRAYS_GRID, only : &
& ddzi, nsoilcols, zsoilcols, &
& isoilcolc, ksoilcol, &
& z_full, z_half, &
& dzeta_int, dzeta_05int

use ARRAYS_SOIL, only : dzs

use NUMERIC_PARAMS, only : &
& pi, small_value

use ATMOS, only : pressure

use TIMEVAR, only : hour_sec

implicit none

! Input variables
integer(kind=iintegers), intent(in) :: M, Mice, ns
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: ndatamax
integer(kind=iintegers), intent(in) :: month, day
integer(kind=iintegers), intent(in) :: lakeform

real(kind=ireals),       intent(in) :: hour

real(kind=ireals),       intent(in) :: dt

real(kind=ireals),       intent(in) :: depth_area(1:ndatamax,1:2)
real(kind=ireals),       intent(in) :: area_lake, cellipt

logical, intent(in) :: multsoil

real(kind=ireals), intent(in) :: trib_inflow
real(kind=ireals), intent(inout) :: dhwtrib

real(kind=ireals), intent(out) :: vol, botar

! Local variables
real(kind=ireals), allocatable :: work1(:), work2(:), work3(:), work4(:), work5(:), work6(:)

integer(kind=iintegers) :: i, j

real(kind=ireals) :: aellipt, bellipt, z, dz, dadz, da

real(kind=ireals) :: watabstr = 0. ! Water abstraction, currently implemented for single lake 

logical, save :: water_abstraction = .false.
logical :: flag

real(kind=ireals), external :: ABSTR

do i = 1, M+1
  z_full(i) = dzeta_int(i)*h1
  preswat(i) = pressure + row0*g*z_full(i)
enddo
do i = 1, M    
  z_half(i) = dzeta_05int(i)*h1
enddo  

! Specification of profile output levels, if not given in setup file
! In this case these levls coincide with numerical dzeta-levels
! if (.not. allocated(zgrid_out)) then
!    ngrid_out%par = M+1
!    allocate (zgrid_out(1:ngrid_out%par))    
! endif
! zgrid_out(:) = z_full(:)

! Constant cross-section
if (depth_area(1,2) < 0.) then ! No data on morphometry
    area_int            (1:M+1)    = area_lake
    area_half           (1:M)      = area_lake
    bathymice (1:Mice+1)%area_int  = area_lake
    bathymice (1:Mice)  %area_half = area_lake
    bathymdice(1:Mice+1)%area_int  = area_lake
    bathymdice(1:Mice)  %area_half = area_lake
    !if (init(ix,iy) == 0 .and. multsoil) then
    !  bathymsoil(1:nsoilcols+1,ix,iy)%area_int  = area_lake
    !  bathymsoil(1:nsoilcols,  ix,iy)%area_half = area_lake     
    !endif
else

  ! Multiple soil columns bathymetry (constant in time)
  if_multsoil : if (multsoil) then
    allocate(work3(1:nsoilcols) )
    do i = 1, nsoilcols
      work3(i) = 0.5*( zsoilcols(i,ix,iy) + zsoilcols(i+1,ix,iy) )
    enddo
    if (init(ix,iy) == 0) then
      allocate(work5(1:nsoilcols+1), work6(1:nsoilcols))
      call LININTERPOL (depth_area(1,1),depth_area(1,2),ndatamax, &
                       & zsoilcols(1,ix,iy),work5,nsoilcols+1,flag)
      call LININTERPOL (depth_area(1,1),depth_area(1,2),ndatamax, &
                       & work3,work6,nsoilcols,flag) 
      bathymsoil(1:nsoilcols+1,ix,iy)%area_int  = work5(1:nsoilcols+1)
      bathymsoil(1:nsoilcols,  ix,iy)%area_half = work6(1:nsoilcols)
      deallocate (work5, work6)
      !! Horizontal cross-section area of soil columns at soil numerical grid levels
      !do i = 1, nsoilcols-1
      !  bathymsoil(i,ix,iy)%sarea_int(1) = small_value
      !  z = 0.; j = 1
      !  dz = zsoilcols(i+1,ix,iy) - zsoilcols(i,ix,iy)
      !  da = bathymsoil(i,ix,iy)%area_int - bathymsoil(i+1,ix,iy)%area_int
      !  dadz = da / dz 
      !  do while (z < dz)
      !    j = j + 1
      !    z = z + dzs(j-1)
      !    bathymsoil(i,ix,iy)%sarea_int(j) = bathymsoil(i,ix,iy)%sarea_int(j-1) + dadz*dzs(j-1)
      !    bathymsoil(i,ix,iy)%sarea_half(j-1) = &
      !    & 0.5*(bathymsoil(i,ix,iy)%sarea_int(j) + bathymsoil(i,ix,iy)%sarea_int(j-1))
      !  enddo
      !  bathymsoil(i,ix,iy)%sarea_int(j)    = min(bathymsoil(i,ix,iy)%sarea_int(j)   ,da)
      !  bathymsoil(i,ix,iy)%sarea_half(j-1) = min(bathymsoil(i,ix,iy)%sarea_half(j-1),da)
      !  bathymsoil(i,ix,iy)%sarea_int(j+1:ns)  = bathymsoil(i,ix,iy)%sarea_int(j)
      !  bathymsoil(i,ix,iy)%sarea_half(j:ns-1) = bathymsoil(i,ix,iy)%sarea_half(j-1)
      !enddo
      !!Lowest soil column
      !i = nsoilcols
      !bathymsoil(i,ix,iy)%sarea_int(1:ns)    = bathymsoil(nsoilcols,ix,iy)%area_int
      !bathymsoil(i,ix,iy)%sarea_half(1:ns-1) = bathymsoil(nsoilcols,ix,iy)%area_int
    endif
    dz = - maxval(depth_area(:,1)) + h1 + ls1 ! Relating coordinate systems before interpolation
    do i = 2, nsoilcols
      do j = M, 1, -1
        if (zsoilcols(i,ix,iy) + dz >  z_full(j) .and. &
            zsoilcols(i,ix,iy) + dz <= z_full(j+1)) then
          bathymsoil(i,ix,iy)%itop = j
          bathymsoil(i,ix,iy)%dzSL = work3(i) + dz - z_full(j)
        endif
      enddo
    enddo
    bathymsoil(1,ix,iy)%itop = 1
    bathymsoil(1,ix,iy)%dzSL = work3(1) + dz - z_full(1)
    do i = 2, nsoilcols
      do j = 1, M
        if (zsoilcols(i,ix,iy) + dz >  z_full(j) .and. &
            zsoilcols(i,ix,iy) + dz <= z_full(j+1)) then
          bathymsoil(i-1,ix,iy)%ibot = j
        endif
      enddo
    enddo
    bathymsoil(nsoilcols,ix,iy)%ibot = M+1
    do i = 1, nsoilcols-1
      do j = 1, M-1
        if (work3(i) + dz >  z_half(j) .and. &
            work3(i) + dz <= z_half(j+1)) then
          bathymsoil(i,ix,iy)%icent = j+1
          bathymsoil(i,ix,iy)%dzSLc = work3(i) + dz - z_half(j)
        endif
      enddo
    enddo
    deallocate(work3)
  else
    !i = nsoilcols
    !bathymsoil(i,ix,iy)%sarea_int(1:ns)    = area_lake
    !bathymsoil(i,ix,iy)%sarea_half(1:ns-1) = area_lake
  endif if_multsoil


  ! Interpolation of the predefined area-depth dependence to the current grid...
  ! ...in water layer
  if (h1 > small_value) then
    if     (soilcolconjtype == 1) then
      allocate(work1(1:nsoilcols+1),work2(1:nsoilcols+1))
      work2(:) = bathymsoil(:,ix,iy)%area_int
      work1(:) = zsoilcols(:,ix,iy) - maxval(zsoilcols(:,ix,iy)) + &
      & h1 + ls1 ! Relating coordinate systems before interpolation
      call PIECEWISECONSTINT (work1,work2,nsoilcols+1,z_full,area_int,M+1,1_iintegers)
      call PIECEWISECONSTINT (work1,work2,nsoilcols+1,z_half,area_half,M,1_iintegers)
      area_int(M+1) = bathymsoil(nsoilcols+1,ix,iy)%area_int
    elseif (soilcolconjtype == 2) then
      allocate(work1(1:ndatamax),work2(1:ndatamax))
      work2(:) = depth_area(:,2)
      work1(:) = depth_area(:,1) - maxval(depth_area(:,1)) + &
      & h1 + ls1 ! Relating coordinate systems before interpolation
      call LININTERPOL (work1,work2,ndatamax,z_full,area_int,M+1,flag)
      call LININTERPOL (work1,work2,ndatamax,z_half,area_half,M,flag)
    endif
  endif
  !... in ice layer
  if (l1 > small_value) then
    allocate (work3(1:Mice+1), work4(1:Mice))
    work3(1) = 0.
    do i = 2, Mice+1
      work3(i) = work3(i-1) + ddzi(i-1)*l1
    enddo
    work4(1) = 0.5*ddzi(1)*l1
    do i = 2, Mice
      work4(i) = work4(i-1) + 0.5*(ddzi(i-1) + ddzi(i))*l1
    enddo
    allocate(work5(1:Mice+1), work6(1:Mice))
    if     (soilcolconjtype == 1) then
      work1(:) = zsoilcols(:,ix,iy) - maxval(zsoilcols(:,ix,iy)) + h1 + l1 + ls1 
      call PIECEWISECONSTINT (work1,work2,nsoilcols,work3,work5,Mice+1,1_iintegers)
      call PIECEWISECONSTINT (work1,work2,nsoilcols,work4,work6,Mice,1_iintegers)
    elseif (soilcolconjtype == 2) then
      work1(:) = depth_area(:,1) - maxval(depth_area(:,1)) + &
      & h1 + l1 + ls1 
      call LININTERPOL (work1,work2,ndatamax,work3,work5,Mice+1,flag)
      call LININTERPOL (work1,work2,ndatamax,work4,work6,Mice,flag) 
    endif
    bathymice(1:Mice+1)%area_int  = work5(1:Mice+1)
    bathymice(1:Mice)  %area_half = work6(1:Mice)
    deallocate(work3, work4, work5, work6)
  endif
  ! ... in deep ice layer
  if (ls1 > small_value) then
    allocate (work3(1:Mice+1), work4(1:Mice))
    work3(1) = 0.
    do i = 2, Mice+1
      work3(i) = work3(i-1) + ddzi(i-1)*ls1
    enddo
    work4(1) = 0.5*ddzi(1)*ls1
    do i = 2, Mice
      work4(i) = work4(i-1) + 0.5*(ddzi(i-1) + ddzi(i))*ls1
    enddo
    allocate(work5(1:Mice+1), work6(1:Mice))
    if     (soilcolconjtype == 1) then
      work1(:) = zsoilcols(:,ix,iy) - maxval(zsoilcols(:,ix,iy)) + ls1 
      call PIECEWISECONSTINT (work1,work2,nsoilcols,work3,work5,Mice+1,1_iintegers)
      call PIECEWISECONSTINT (work1,work2,nsoilcols,work4,work6,Mice,1_iintegers)
      work5(Mice+1) = bathymsoil(nsoilcols+1,ix,iy)%area_int
    elseif (soilcolconjtype == 2) then
      work1(:) = depth_area(:,1) - maxval(depth_area(:,1)) + ls1 
      call LININTERPOL (work1,work2,ndatamax,work3,work5,Mice+1,flag)
      call LININTERPOL (work1,work2,ndatamax,work4,work6,Mice,flag) 
    endif
    bathymdice(1:Mice+1)%area_int  = work5(1:Mice+1)
    bathymdice(1:Mice)  %area_half = work6(1:Mice)
    deallocate(work3, work4, work5, work6)
  endif

  deallocate (work1,work2)
endif


! Cross-section parameters for water
if (h1 > small_value) then
  do i = 1, M+1
    ! Assuming horizontal cross-section to be an ellipse at every level
    call LXLY(area_int(i),Lx(i),Ly(i))
    call LXLY(area_half(i),bathymwater(i)%Lx_half,bathymwater(i)%Ly_half)
  enddo
  bathymwater(1:M+1)%area_int  = area_int (1:M+1)
  bathymwater(1:M+1)%area_half = area_half(1:M+1)
  bathymwater(1:M+1)%Lx        = Lx       (1:M+1)
  bathymwater(1:M+1)%Ly        = Ly       (1:M+1)
  bathymwater(1:M+1)%rad_int  =  sqrt(area_int (1:M+1)/pi)
endif
! Cross-section parameters for ice
if (l1 > small_value) then
  do i = 1, Mice+1
    ! Assuming horizontal cross-section to be an ellipse at every level
    call LXLY(bathymice(i)%area_int,bathymice(i)%Lx,bathymice(i)%Ly)
    call LXLY(bathymice(i)%area_half,bathymice(i)%Lx_half,bathymice(i)%Ly_half)
  enddo
endif
! Cross-section parameters for deep ice
if (ls1 > small_value) then
  do i = 1, Mice+1
    call LXLY(bathymdice(i)%area_int,bathymdice(i)%Lx,bathymdice(i)%Ly)
    call LXLY(bathymdice(i)%area_half,bathymdice(i)%Lx_half,bathymdice(i)%Ly_half)
  enddo
endif
! Cross-section parameters for multiple soil columns
if (init(ix,iy) == 0 .and. multsoil) then
  ! Assuming horizontal cross-section to be an ellipse at every level
  do i = 1, nsoilcols+1
    call LXLY(bathymsoil(i,ix,iy)%area_int,bathymsoil(i,ix,iy)%Lx,bathymsoil(i,ix,iy)%Ly)
  enddo
endif

if (h1 > small_value .and. multsoil) then
  allocate(work1(1:nsoilcols+1))
  work1(:) = zsoilcols(:,ix,iy) - maxval(zsoilcols(:,ix,iy)) + &
  & h1 + ls1 ! Relating coordinate systems before interpolation
  ksoilcol(:) = nsoilcols
  do i = 1, nsoilcols-1
    z = 0.5*(work1(i) + work1(i+1))
    do j = 1, M
      if (z_full(j) <= z .and. z_full(j+1) > z) then
        isoilcolc(i) = j
        exit
      endif
    enddo
    do j = 1, M
      if (z_full(j) >= work1(i)   .and. &
        & z_full(j) <  work1(i+1)) then
        ksoilcol(j) = i
      endif
    enddo
  enddo
  deallocate(work1)
endif


! Lake volume and bottom area
vol = 0.
botar = area_int(M+1)
do i = M, 1, -1
  vol = vol + 0.5*(area_int(i+1) + area_int(i))* &
  & (z_full(i+1) - z_full(i))
  botar = botar + area_int(i) - area_int(i+1)
enddo


!TRIBUTARY INFLOW to a lake
if (water_abstraction .and. month == 12 .and. day == 31 &
& .and. (24. - hour)*hour_sec <= dt) then
! Water abstraction for the next year, m**3/s
  watabstr = ABSTR(vol)
  if (trib_inflow /= -9999.) then
    dhwtrib = dhwtrib - watabstr*dt/area_int(1)
  endif
endif

contains

!> Subroutine LXLY calculates the length and width of the lake
!! given its cross-section area and form
SUBROUTINE LXLY(area,Lx_,Ly_)

implicit none

real(kind=ireals), intent(in) :: area
real(kind=ireals), intent(out) :: Lx_, Ly_

if (lakeform == form_ellipse) then
  bellipt = sqrt(area/(pi*cellipt))
  aellipt = bellipt*cellipt
  Lx_ = 2.*aellipt
  Ly_ = 2.*bellipt
elseif (lakeform == form_rectangle) then
  Ly_ = sqrt(area/cellipt)
  Lx_ = cellipt*Ly_
else
  print*, 'Incorrect LAKEFORM identifier:', lakeform, ', STOP'
  STOP
endif

END SUBROUTINE LXLY

END SUBROUTINE BATHYMETRY

END MODULE BATHYMSUBR
