!> @file database.F90
!! @brief Program which loads the hydrologic database generated by "database.py" (e.g. "hydrologic_database.nc")
!! and completes it with the lake volume/area/elevation curves of every depression into "hydrologic_database_final.nc".
!! @author Alexandre Gauvain
!! @date 10/08/2023

module LeaveNodeModule
!-----------------------------------------------------------------------
! NAME
!     LeaveNodeModule
!
! DESCRIPTION
!     Linked list of the leaf depressions found below a node of the depression
!     tree, with their bounding box.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  implicit none
  type :: LeaveNode
    integer :: id, imin, imax, jmin, jmax
    type(LeaveNode), pointer :: next => null()
  end type LeaveNode
end module LeaveNodeModule

module TreeNodeModule
!-----------------------------------------------------------------------
! NAME
!     TreeNodeModule
!
! DESCRIPTION
!     Node of the depression tree used to build the hydrologic database.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  use, intrinsic :: iso_fortran_env
  implicit none
  integer, parameter :: qp = real128
  type :: TreeNode
    integer :: id, imin, imax, jmin, jmax, index
    logical :: is_leaf
    integer, allocatable :: leaves(:)
    double precision :: Zmax
    double precision :: Vmax, Vmin, Aw, Amax
    logical(kind=1), allocatable :: label_mask(:,:)
    type(TreeNode), pointer :: parent => null()
    type(TreeNode), pointer :: left => null()
    type(TreeNode), pointer :: right => null()
  end type TreeNode
end module TreeNodeModule

program database
!-----------------------------------------------------------------------
! NAME
!     database
!
! DESCRIPTION
!     Completes the hydrologic database written by "database.py" with the lake
!     volume, area and elevation profiles of every depression, and writes the
!     final database "hydrologic_database_final.nc".
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!     Usage: database input_file output_folder. Without arguments, both
!     are asked.
!     The profiles of a leaf depression are computed from the cells of its
!     label. The ones of a parent depression come from a union-find on the
!     cells sorted by elevation.
!     The planet radius is the global attribute 'R' of the input file (3389.5
!     km if missing). The depression table must be sorted by 'Zmax', parents
!     after their children.
!     The final database and "hydrologic_database.log" are written only when
!     all the depressions are completed.
!-----------------------------------------------------------------------

use :: netcdf
use :: omp_lib
use :: TreeNodeModule
use :: LeaveNodeModule
implicit none

! Variables

! Files path
character(len=1000) :: infile
character(len=1000) :: outfolder
character(10) :: format_string 

! Longitudes:
integer :: lonlen
double precision, allocatable :: longitude(:)

! Latitudes
integer :: latlen
double precision, allocatable :: latitude(:)

! Depressions
integer :: deplen, IDdep, size_leaf, leaf

! Leaves informations
integer :: leaflen
integer, allocatable :: leaves_len(:), leaves(:)

! Function between water volume and water surface elevation
double precision, allocatable :: lake_vol(:,:), lake_area(:,:), lake_elev(:,:)
integer :: lake_dislen
double precision, allocatable :: lake_perc(:), elev(:,:)
double precision :: h, hmax, z

! Graph table
integer, allocatable :: child_list(:), lchild_list(:), rchild_list(:), parent_list(:), check_list(:)
integer, allocatable :: pit_cell_list(:), out_cell_list(:), sibling_list(:), geolink_list(:) 
integer, allocatable :: imin_list(:), imax_list(:), jmin_list(:), jmax_list(:), i_out_list(:), j_out_list(:)
double precision, allocatable :: Vmax_list(:), Vmin_list(:), pit_elev_list(:), lake_lat_list(:), lake_long_list(:)
double precision, allocatable :: Zmin_list(:), Zmax_list(:), watershed_area_list(:)
integer :: child, parent,  pit_cell, out_cell, sibling, geolink, imin, imax, jmin, jmax
integer :: loc(2)
double precision :: Zmin, Zmax, lake_long, lake_lat, watershed_area
double precision :: Vmax, Vmin, Vmin_childs, volume_lake, area_lake

! Maps informations
double precision, allocatable :: elevation(:,:), area(:,:), depth(:), cells_area(:)
integer,allocatable :: label(:,:), cell_number(:,:), flowdirs(:,:)

! NetCDF stuff
integer :: status ! NetCDF return code
integer :: ncid ! NetCDF file ID
integer :: dimid ! NetCDF dimension ID
integer :: varid ! NetCDF variable ID
integer :: lon_dimid,lat_dimid ! NetCDF dimension IDs
integer :: lon_varid,lat_varid ! NetCDF variable IDs
integer :: dimid_lon, dimid_lat, dimid_dep, dimid_lakedis, dimid_maxleafid, dimid_maxdepid,dimid_nleaf, dimid_ndep
integer :: varid_lon, varid_lat, varid_child, varid_rchild, varid_lchild, varid_parent
integer :: varid_sibling, varid_geolink, varid_Vmax, varid_Vmin, varid_check
integer :: varid_out_cell, varid_pit_cell, varid_lake_lon, varid_lake_lat, varid_Zmin, varid_Zmax, varid_watershed_area
integer :: varid_elevation, varid_label, varid_area, varid_flowdirs, varid_cell_number
integer :: varid_imin, varid_imax, varid_jmin, varid_jmax, varid_i_out, varid_j_out
integer :: varid_lake_perc, varid_lake_elev, varid_lake_area, varid_lake_vol

! Tree stuff 
type(TreeNode), target :: root
type(TreeNode), target, dimension(:), allocatable :: nodes

! Complete Database stuff
integer, allocatable :: label_mask_box(:,:), elevation_mask_box(:,:), label_mask(:,:), mask_indices(:,:)
real, allocatable :: elevation_box(:,:), area_box(:,:), lake_depth(:,:)
integer, allocatable :: leaves_store(:)
integer :: nleaf, ndep, i, j, k
integer :: total_sum, count, count_compute
integer :: print_count(10)

! Time stuff
real(8) :: cpu_start_time
! Computation time
real(kind=8) :: temps, t_cpu_0, t_cpu_1, t_cpu
integer :: ir, t0, t1

! Optimised version (cell pass + union-find) variables
logical, allocatable :: is_leaf_id(:), computed(:)
real, allocatable :: zmin4(:), cells_area4(:), lake_area4(:,:), watershed_area4(:)
integer, allocatable :: loc_i(:), loc_j(:), dsu(:), cell_order(:)
integer(8), allocatable :: bstart(:), bfill(:)
double precision, allocatable :: Tsum(:), Lsum(:), Sub(:)
real :: e4, a4
double precision :: emin, emax, ev, a, Alow
double precision :: radius ! Planet radius (km), global attribute "R" of the input file
integer :: r, nviol, nb, b, c, x, p
integer(8) :: q

! 0. getargument 
  ! Check if argument is the path of database.nc is found
  if (command_argument_count() == 2) then
    call get_command_argument(1, infile)
    call get_command_argument(2, outfolder)
    write(*,*) "Input file: ", infile
    write(*,*) "Save folder: ", outfolder
  else
    ! If not, ask the odatabase.nc 
    write(*,*) "Input file: "
    read(*,*) infile
    write(*,*) "Save folder: "
    read(*,*) outfolder
  end if

! 1. Load input datasets
  write(*,*) "###################################"
  write(*,*) "#### Load Hydrologic Database #####"
  write(*,*) "###################################"
  ! 1.1. Open input datafile
  status=nf90_open(infile,NF90_NOWRITE,ncid)
  write(*,*) "Reading input file: ",trim(infile)

  ! Planet radius (km) given in database_v9.py
  status=nf90_get_att(ncid,NF90_GLOBAL,"R",radius)
  if (status /= NF90_NOERR) then
    radius = 3389.5
    write(*,*) "Warning: no attribute R in the input file, default value used"
  end if
  write(*,*) "R (km) =", radius

  !  1.2 Identify/load dimensions
  ! Longitude
  status=nf90_inq_dimid(ncid,"lon",dimid)
  status=nf90_inquire_dimension(ncid,dimid,len=lonlen)
  write(*,*) "lonlen=",lonlen
  ! Latitude
  status=nf90_inq_dimid(ncid,"lat",dimid)
  status=nf90_inquire_dimension(ncid,dimid,len=latlen)
  write(*,*) "latlen=",latlen
  ! Depressions
  status=nf90_inq_dimid(ncid,"depression",dimid) ! Find depression dimension
  status=nf90_inquire_dimension(ncid,dimid,len=deplen) ! Find depression length
  write(*,*) "deplen=",deplen

  ! lake_dis
  status=nf90_inq_dimid(ncid,"lake",dimid) ! Find vol2cell dimension
  status=nf90_inquire_dimension(ncid,dimid,len=lake_dislen) ! Find vol2cell length
  write(*,*) "lake_dislen=",lake_dislen
  
  ! 1.3 Allocate dimensions and read variables
  ! 1.3.1 Latitudes and Longitudes
  ! Longitude
  allocate(longitude(lonlen),stat=status) ! Allocate longitude()
  status=nf90_inq_varid(ncid,"longitude",varid) ! Find longitude varid
  status=nf90_get_var(ncid,varid,longitude)! read longitude
  write(*,*) 'lonlen=',lonlen
  write(*,*) 'longitude(1)=',longitude(1)
  write(*,*) 'longitude(lonlen)=',longitude(lonlen)
  ! Latitude
  allocate(latitude(latlen),stat=status) ! Allocate latitude()
  status=nf90_inq_varid(ncid,"latitude",varid) ! Find latitude varid
  status=nf90_get_var(ncid,varid,latitude) ! Read latitude
  write(*,*) 'latlen=',latlen
  write(*,*) 'latitude(1)=',latitude(1)
  write(*,*) 'latitude(latlen)=',latitude(latlen)
  
  ! 1.3.2 Graph table  
  ! Child
  allocate(child_list(deplen),stat=status) ! Allocate child()
  status=nf90_inq_varid(ncid,"child",varid) ! Find child varid
  status=nf90_get_var(ncid,varid,child_list) ! Load and read child
  write(*,*) 'child_list(1)=',child_list(1)
  write(*,*) 'child_list(deplen)=',child_list(deplen)
  
  ! Parent
  allocate(parent_list(deplen),stat=status) ! Allocate parent()
  status=nf90_inq_varid(ncid,"parent",varid) ! Find parent varid
  status=nf90_get_var(ncid,varid,parent_list) ! Load and read parent
  write(*,*) 'parent_list(1)=',parent_list(1)
  write(*,*) 'parent_list(deplen)=',parent_list(deplen)
  
  ! lchild
  allocate(lchild_list(deplen),stat=status) ! Allocate lchild()
  status=nf90_inq_varid(ncid,"lchild",varid) ! Find lchild varid
  status=nf90_get_var(ncid,varid,lchild_list) ! Load and read lchild
  write(*,*) 'lchild_list(1)=',lchild_list(1)
  write(*,*) 'lchild_list(deplen)=',lchild_list(deplen)
  
  ! rchild
  allocate(rchild_list(deplen),stat=status) ! Allocate rchild()
  status=nf90_inq_varid(ncid,"rchild",varid) ! Find rchild varid
  status=nf90_get_var(ncid,varid,rchild_list) ! Load and read rchild
  write(*,*) 'rchild_list(1)=',rchild_list(1)
  write(*,*) 'rchild_list(deplen)=',rchild_list(deplen)
  
  
  ! Vmax
  allocate(Vmax_list(deplen),stat=status) ! Allocate Vmax()
  status=nf90_inq_varid(ncid,"Vmax",varid) ! Find Vmax varid
  status=nf90_get_var(ncid,varid,Vmax_list) ! Load and read Vmax
  write(*,*) 'Vmax_list(1)=',Vmax_list(1)
  write(*,*) 'Vmax_list(deplen)=',Vmax_list(deplen)
  
  ! Vmin
  allocate(Vmin_list(deplen),stat=status) ! Allocate Vmax()
  status=nf90_inq_varid(ncid,"Vmin",varid) ! Find Vmax varid
  status=nf90_get_var(ncid,varid,Vmin_list) ! Load and read Vmax
  write(*,*) 'Vmin_list(1)=',Vmin_list(1)
  write(*,*) 'Vmin_list(deplen)=',Vmin_list(deplen)
  
  ! pit_cell
  allocate(pit_cell_list(deplen),stat=status) ! Allocate pit_cell()
  status=nf90_inq_varid(ncid,"pit_cell",varid) ! Find pit_cell varid
  status=nf90_get_var(ncid,varid,pit_cell_list) ! Load and read pit_cell
  write(*,*) 'pit_cell_list(1)=',pit_cell_list(1)
  write(*,*) 'pit_cell_list(deplen)=',pit_cell_list(deplen)
  
  ! out_cell
  allocate(out_cell_list(deplen),stat=status) ! Allocate out_cell()
  status=nf90_inq_varid(ncid,"out_cell",varid) ! Find out_cell varid
  status=nf90_get_var(ncid,varid,out_cell_list) ! Load and pread out_cell
  write(*,*) 'out_cell_list(1)=',out_cell_list(1)
  write(*,*) 'out_cell_list(deplen)=',out_cell_list(deplen)
  
  ! i_out
  allocate(i_out_list(deplen),stat=status) ! Allocate i_out()
  status=nf90_inq_varid(ncid,"i_out",varid) ! Find i_out varid
  status=nf90_get_var(ncid,varid,i_out_list) ! Load and pread i_out
  write(*,*) 'i_out_list(1)=',i_out_list(1)
  write(*,*) 'i_out_list(deplen)=',i_out_list(deplen)
  
  ! j_out
  allocate(j_out_list(deplen),stat=status) ! Allocate j_out()
  status=nf90_inq_varid(ncid,"j_out",varid) ! Find j_out varid
  status=nf90_get_var(ncid,varid,j_out_list) ! Load and pread j_out
  write(*,*) 'j_out_list(1)=',j_out_list(1)
  write(*,*) 'j_out_list(deplen)=',j_out_list(deplen)
  
  ! lake_long
  allocate(lake_long_list(deplen),stat=status) ! Allocate lake_long()
  status=nf90_inq_varid(ncid,"lake_long",varid) ! Find lake_long varid
  status=nf90_get_var(ncid,varid,lake_long_list) ! Load and read lake_long
  write(*,*) 'lake_long_list(1)=',lake_long_list(1)
  write(*,*) 'lake_long_list(deplen)=',lake_long_list(deplen)
  
  ! lake_lat
  allocate(lake_lat_list(deplen),stat=status) ! Allocate lake_lat()
  status=nf90_inq_varid(ncid,"lake_lat",varid) ! Find lake_lat varid
  status=nf90_get_var(ncid,varid,lake_lat_list) ! Load and read lake_lat
  write(*,*) 'lake_lat_list(1)=',lake_lat_list(1)
  write(*,*) 'lake_lat_list(deplen)=',lake_lat_list(deplen)
  
  ! Zmin
  allocate(Zmin_list(deplen),stat=status) ! Allocate Zmin()
  status=nf90_inq_varid(ncid,"Zmin",varid) ! Find Zmin varid
  status=nf90_get_var(ncid,varid,Zmin_list) ! Load and read Zmin
  write(*,*) 'Zmin_list(1)=',Zmin_list(1)
  write(*,*) 'Zmin_list(deplen)=',Zmin_list(deplen)
  
  ! Zmax
  allocate(Zmax_list(deplen),stat=status) ! Allocate Zmax()
  status=nf90_inq_varid(ncid,"Zmax",varid) ! Find Zmax varid
  status=nf90_get_var(ncid,varid,Zmax_list) ! Load and read Zmax
  write(*,*) 'Zmax_list(1)=',Zmax_list(1)
  write(*,*) 'Zmax_list(deplen)=',Zmax_list(deplen)
  
  ! sibling
  allocate(sibling_list(deplen),stat=status) ! Allocate sibling()
  status=nf90_inq_varid(ncid,"sibling",varid) ! Find sibling varid
  status=nf90_get_var(ncid,varid,sibling_list) ! Load and read sibling
  write(*,*) 'sibling_list(1)=',sibling_list(1)
  write(*,*) 'sibling_list(deplen)=',sibling_list(deplen)
  
  ! geolink
  allocate(geolink_list(deplen),stat=status) ! Allocate geolink()
  status=nf90_inq_varid(ncid,"geolink",varid) ! Find geolink varid
  status=nf90_get_var(ncid,varid,geolink_list) ! Load and read geolink
  write(*,*) 'geolink_list(1)=',geolink_list(1)
  write(*,*) 'geolink_list(deplen)=',geolink_list(deplen)
  
  ! watershed_area
  allocate(watershed_area_list(deplen),stat=status) ! Allocate watershed_area()
  status=nf90_inq_varid(ncid,"watershed_area",varid) ! Find watershed_area varid
  status=nf90_get_var(ncid,varid,watershed_area_list) ! Load and read watershed_area
  write(*,*) 'watershed_area_list(1)=',watershed_area_list(1)
  write(*,*) 'watershed_area_list(deplen)=',watershed_area_list(deplen)
  watershed_area_list = 0

  ! check
  allocate(check_list(deplen),stat=status) ! Allocate check_list()
  status=nf90_inq_varid(ncid,"check",varid) ! Find check varid
  status=nf90_get_var(ncid,varid,check_list) ! Load and read check
  write(*,*) 'check_list(1)=',check_list(1)
  write(*,*) 'check_list(deplen)=',check_list(deplen)

  ! 1.3.4 Index to zoom on leaf into maps informations
  ! Imin
  allocate(imin_list(deplen),stat=status) ! Allocate imin
  status=nf90_inq_varid(ncid,"imin",varid) ! Find imin varid
  status=nf90_get_var(ncid,varid,imin_list) ! Load and read imin
  write(*,*) 'imin_list(1)=',imin_list(1)
  write(*,*) 'imin_list(deplen)=',imin_list(deplen)
  
  ! Imax
  allocate(imax_list(deplen),stat=status) ! Allocate imax
  status=nf90_inq_varid(ncid,"imax",varid) ! Find imax varid
  status=nf90_get_var(ncid,varid,imax_list) ! Load and read imax
  write(*,*) 'imax_list(1)=',imax_list(1)
  write(*,*) 'imax_list(deplen)=',imax_list(deplen)
  
  ! Jmin
  allocate(jmin_list(deplen),stat=status) ! Allocate jmin
  status=nf90_inq_varid(ncid,"jmin",varid) ! Find jmin varid
  status=nf90_get_var(ncid,varid,jmin_list) ! Load and read jmin
  write(*,*) 'jmin_list(1)=',jmin_list(1)
  write(*,*) 'jmin_list(deplen)=',jmin_list(deplen)
  
  ! Jmax
  allocate(jmax_list(deplen),stat=status) ! Allocate jmax
  status=nf90_inq_varid(ncid,"jmax",varid) ! Find jmax varid
  status=nf90_get_var(ncid,varid,jmax_list) ! Load and read jmax
  write(*,*) 'jmax_list(1)=',jmax_list(1)
  write(*,*) 'jmax_list(deplen)=',jmax_list(deplen)

  !Load ID Label of depression map
  allocate(label(lonlen,latlen),stat=status) ! Allocate label()
  status=nf90_inq_varid(ncid,"label",varid) ! Find label varid
  status=nf90_get_var(ncid,varid,label) ! Load and read label
  write(*,*) 'shape(label)=',shape(label)

  ! Elevation
  allocate(elevation(lonlen,latlen),stat=status) ! Allocate elevation()
  status=nf90_inq_varid(ncid,"elevation",varid) ! Find elevation varid
  status=nf90_get_var(ncid,varid,elevation) ! Load and read elevation
  write(*,*) 'shape(elevation)=',shape(elevation)
  write(*,*) 'elevation(1,1)=',elevation(1,1)
  write(*,*) 'elevation(lonlen,latlen)=',elevation(lonlen,latlen)

  write(*,*) "Load elevation map"
  ! 1.3.5 Read and load function between water volume and water surface elevation
  ! Lake volume discretisation
  allocate(lake_perc(lake_dislen),stat=status) ! Allocate lake_perc()
  status=nf90_inq_varid(ncid,"lake_perc",varid) ! Find lake_perc varid
  status=nf90_get_var(ncid,varid,lake_perc) ! Load and read lake_perc
  write(*,*) 'shape(lake_perc)=',shape(lake_perc)

  ! Nleaf
  nleaf = maxval(label)
  write(*,*) 'nleaf=',nleaf
  
  ! Size_leaf
  size_leaf = sum(merge(1,0,child_list<=nleaf))
  write(*,*) 'size_leaf=',size_leaf
  
  ! Ndep
  ndep = size(child_list)
  write(*,*) 'ndep=',ndep
  
  ! IDdep
  IDdep = maxval(child_list)+1
  write(*,*) 'IDdep=',IDdep

  ! Nombre de Threads utilisé
  !$OMP PARALLEL
    write(*,*) 'Num THREADS =', OMP_GET_MAX_THREADS()
  !$OMP END PARALLEL
  
! 2. Build tree
  ! Build nodes
  allocate(nodes(IDdep))
  do i = 1, IDdep
    nodes(i)%index = 0
    nodes(i)%id = i
    nodes(i)%is_leaf = .false.
    nodes(i)%Zmax = 0
    nodes(i)%imin = lonlen + 1
    nodes(i)%imax = 0
    nodes(i)%jmin = latlen + 1
    nodes(i)%jmax = 0
    nodes(i)%Aw = 0
    nodes(i)%Amax = 0
    nodes(i)%Vmax = 0
    nodes(i)%Vmin = 0
  end do
  
  do i = 1, ndep
    nodes(child_list(i))%index = i
    nodes(child_list(i))%id = child_list(i)
    nodes(child_list(i))%Zmax = Zmax_list(i)
    ! Initialisation des indices
    nodes(child_list(i))%imin = lonlen + 1
    nodes(child_list(i))%imax = 0
    nodes(child_list(i))%jmin = latlen + 1
    nodes(child_list(i))%jmax = 0
    if (child_list(i) <= nleaf) then
      nodes(child_list(i))%is_leaf = .true.
      allocate(nodes(child_list(i))%leaves(1))
      nodes(child_list(i))%leaves(1) = child_list(i)
    else
      nodes(child_list(i))%is_leaf = .false.
    end if
  end do
  
  ! Build tree
  root%id = maxval(child_list) ! Racine
  root%is_leaf = .false.
  root%parent => null()
  
  do i = 1, ndep
    if (child_list(i) <= nleaf) then
      nodes(child_list(i))%parent => nodes(parent_list(i))
    end if
    if (child_list(i) > nleaf) then
      nodes(child_list(i))%parent => nodes(parent_list(i))
    end if
    if (associated(nodes(parent_list(i))%left)) then
      nodes(parent_list(i))%right => nodes(child_list(i))
    else
      nodes(parent_list(i))%left => nodes(child_list(i))
    end if
  end do

! 3. Hydrologic Database
  write(*,*) "#######################################"
  write(*,*) "#### Complete Hydrologic Database #####"
  write(*,*) "#######################################"
  call system_clock(t0, ir)

  ! Cell area (m2) for each latitude: double precision (parent depressions, area(:,:) in v1)
  ! and simple precision (leaf depressions, area_box in v1)
  allocate(cells_area(latlen), cells_area4(latlen))
  cells_area = (4 * atan(1.0_8) * radius * 1000.0 / 180.0)**2 * &
              abs(longitude(1) - longitude(2)) * &
              abs(latitude(1) - latitude(2)) * &
              cos(latitude * 4 * atan(1.0_8) / 180.0)
  cells_area4 = real(cells_area)
  write(*,*) "Planet area = ", sum(cells_area) * lonlen

  ! Leaf depressions of the graph table
  allocate(is_leaf_id(0:IDdep))
  is_leaf_id = .false.
  do i = 1, ndep
    if (child_list(i) <= nleaf) is_leaf_id(child_list(i)) = .true.
  end do

  write(*,*) "1. Find bounding box, minimum elevation of each leaf depression"
  ! Column-major loop (same order as minval/maxloc on the bounding box in v1)
  allocate(zmin4(IDdep), loc_i(IDdep), loc_j(IDdep))
  zmin4 = huge(1.0)
  loc_i = 0
  loc_j = 0
  do j = 1, latlen
    do i = 1, lonlen
      leaf = label(i,j)
      if (leaf < 1 .or. leaf > IDdep) cycle
      nodes(leaf)%imin = min(nodes(leaf)%imin, i)
      nodes(leaf)%imax = max(nodes(leaf)%imax, i)
      nodes(leaf)%jmin = min(nodes(leaf)%jmin, j)
      nodes(leaf)%jmax = max(nodes(leaf)%jmax, j)
      if (.not. is_leaf_id(leaf)) cycle
      e4 = real(elevation(i,j))
      if (e4 < zmin4(leaf)) then
        zmin4(leaf) = e4
        loc_i(leaf) = i
        loc_j(leaf) = j
      end if
    end do
  end do

  ! Lake elevation
  allocate(lake_elev(lake_dislen, deplen),stat=status) ! Allocate lake_elev()
  lake_elev = 0
  ! Lake volume
  allocate(lake_vol(lake_dislen, deplen),stat=status) ! Allocate lake_vol()
  lake_vol = 0
  ! lake surface
  allocate(lake_area(lake_dislen, deplen),stat=status) ! Allocate lake_area()
  lake_area = 0
  status = nf90_close(ncid) ! Close NetCDF

  write(*,*) "2. Complete Database for leaf depressions"
  do i = 1, ndep
    child = child_list(i)
    if (child <= nleaf) then
      Zmin = zmin4(child)
      hmax = Zmax_list(i) - Zmin
      do j = 1, lake_dislen
        h = lake_perc(j) * hmax
        lake_elev(j,i) = h + Zmin
      end do
    end if
  end do

  ! Sums in the same order as v1 (column-major on the bounding box): same results
  allocate(lake_area4(lake_dislen, ndep), watershed_area4(ndep))
  lake_area4 = 0
  watershed_area4 = 0
  do j = 1, latlen
    a4 = cells_area4(j)
    do i = 1, lonlen
      leaf = label(i,j)
      if (leaf < 1 .or. leaf > IDdep) cycle
      if (.not. is_leaf_id(leaf)) cycle
      r = nodes(leaf)%index
      e4 = real(elevation(i,j))
      watershed_area4(r) = watershed_area4(r) + a4
      do k = 1, lake_dislen
        z = lake_elev(k,r)
        if (e4 < z) then
          lake_vol(k,r) = lake_vol(k,r) + (z - e4) * a4
          lake_area4(k,r) = lake_area4(k,r) + a4
        end if
      end do
    end do
  end do

  do i = 1, ndep
    child = child_list(i)
    if (child <= nleaf) then
      imin = nodes(child)%imin
      imax = nodes(child)%imax
      jmin = nodes(child)%jmin
      jmax = nodes(child)%jmax
      Zmin = zmin4(child)
      if (loc_i(child) > 0) then
        lake_long = longitude(loc_i(child))
        lake_lat = latitude(loc_j(child))
      else ! leaf without cell on the map (maxloc = 0 in v1)
        lake_long = longitude(max(1, min(lonlen, imin - 1)))
        lake_lat = latitude(max(1, min(latlen, jmin - 1)))
      end if
      lake_area(:,i) = lake_area4(:,i)
      nodes(child)%Vmin = 0
      nodes(child)%Vmax = lake_vol(lake_dislen,i)
      ! Map Zoom
      imin_list(i) = imin
      imax_list(i) = imax
      jmin_list(i) = jmin
      jmax_list(i) = jmax
      ! Hydrology
      Vmax_list(i) = lake_vol(lake_dislen,i)
      Vmin_list(i) = 0
      nodes(child)%Amax = lake_area(lake_dislen,i)
      watershed_area_list(i) = watershed_area4(i)
      pit_cell_list(i) = 0
      lake_long_list(i) = lake_long
      lake_lat_list(i) = lake_lat
      Zmin_list(i) = Zmin
      check_list(i) = 1
    end if
  end do
  deallocate(lake_area4, watershed_area4, zmin4, loc_i, loc_j)
  call system_clock(t1)
  write(*,*) "   time (s):", real(t1-t0)/real(ir)

! 4. Complete Database for parent depressions
  ! For a parent depression P (Zmin = spillover elevation of its children, Zmax = its own),
  ! v1 fills the cells of P below Zmin and computes, for each level z:
  !   area(z) = sum_{c in P, max(e_c,Zmin) <= z} a_c ,  vol(z) = sum (z - max(e_c,Zmin)) a_c
  ! = Alow*(z-Zmin) + sum_{c in P, Zmin < e_c <= z} (z - e_c) a_c   (Alow = area of P below Zmin)
  ! Since Zmax increases from a leaf to the root, each cell has only one depression with
  ! Zmin < e_c <= Zmax on its branch ("entry" depression). It is found for all the cells
  ! with a union-find on the cells sorted by elevation: O(ncell) instead of a scan of the
  ! bounding box of each parent depression.
  write(*,*) "3. Complete Database for parent depressions"
  call system_clock(t0, ir)

  allocate(computed(ndep))
  computed = .false.
  count = 0
  nviol = 0
  do i = 1, ndep
    child = child_list(i)
    ! The depression hierarchy must be sorted by Zmax, parents after their children
    if (nodes(parent_list(i))%index > 0) then
      if (nodes(parent_list(i))%index < i) nviol = nviol + 1
    end if
    if (child > nleaf) then
      Zmin = max(nodes(child)%left%Zmax, nodes(child)%right%Zmax)
      if (nodes(child)%left%Zmax /= nodes(child)%right%Zmax) then
        count = count + 1
        if (count <= 10) write(*,*) 'Problème avec :', child
      end if
      Zmin_list(i) = Zmin
      if (Zmin /= Zmax_list(i) .and. check_list(i) == 0) then
        computed(i) = .true.
        hmax = Zmax_list(i) - Zmin
        do j = 2, lake_dislen
          lake_elev(j,i) = lake_perc(j) * hmax + Zmin
        end do
      end if
    end if
  end do
  write(*,*) 'Number of parent depressions with left%Zmax /= right%Zmax:', count
  if (nviol > 0) then
    write(*,*) 'ERROR: depression table not sorted by Zmax (parent before child):', nviol
    stop 1
  end if

  ! Sort the cells by elevation (counting sort, elevation is an integer value)
  emin = minval(elevation)
  emax = maxval(elevation)
  if (any(elevation /= anint(elevation))) then
    write(*,*) 'ERROR: elevation must be integer values'
    stop 1
  end if
  if (int(lonlen,8) * int(latlen,8) > int(huge(1),8)) then
    write(*,*) 'ERROR: too many cells for 32 bits indices'
    stop 1
  end if
  nb = int(emax - emin) + 1
  allocate(bstart(nb+1))
  bstart = 0
  do j = 1, latlen
    do i = 1, lonlen
      leaf = label(i,j)
      if (leaf < 1 .or. leaf > IDdep) cycle
      if (.not. is_leaf_id(leaf)) cycle
      b = int(elevation(i,j) - emin) + 2
      bstart(b) = bstart(b) + 1
    end do
  end do
  bstart(1) = 1
  do b = 2, nb+1
    bstart(b) = bstart(b) + bstart(b-1)
  end do
  allocate(cell_order(bstart(nb+1)-1))
  allocate(bfill(nb))
  bfill = bstart(1:nb)
  do j = 1, latlen
    do i = 1, lonlen
      leaf = label(i,j)
      if (leaf < 1 .or. leaf > IDdep) cycle
      if (.not. is_leaf_id(leaf)) cycle
      b = int(elevation(i,j) - emin) + 1
      cell_order(bfill(b)) = i + (j-1)*lonlen
      bfill(b) = bfill(b) + 1
    end do
  end do
  deallocate(bfill)

  ! Union-find: a depression is merged into its parent when the elevation exceeds its Zmax
  allocate(dsu(IDdep), Tsum(IDdep), Lsum(IDdep), Sub(IDdep))
  do k = 1, IDdep
    dsu(k) = k
  end do
  Tsum = 0  ! area of the cells whose entry depression is k
  Lsum = 0  ! area of the cells whose entry depression is k and below its Zmin
  p = 1
  do b = 1, nb
    ev = emin + (b-1)
    do while (p <= ndep)
      if (Zmax_list(p) >= ev) exit
      dsu(child_list(p)) = parent_list(p)
      p = p + 1
    end do
    do q = bstart(b), bstart(b+1)-1
      c = cell_order(q)
      j = (c-1)/lonlen + 1
      i = c - (j-1)*lonlen
      x = label(i,j)
      do while (dsu(x) /= x)
        dsu(x) = dsu(dsu(x))
        x = dsu(x)
      end do
      a = cells_area(j)
      Tsum(x) = Tsum(x) + a
      if (x > nleaf) then
        r = nodes(x)%index
        if (r > 0) then
          if (ev <= Zmin_list(r)) then
            Lsum(x) = Lsum(x) + a
          else if (computed(r) .and. ev <= Zmax_list(r)) then
            do k = 2, lake_dislen
              z = lake_elev(k,r)
              if (ev <= z) then
                lake_vol(k,r) = lake_vol(k,r) + (z - ev) * a
                lake_area(k,r) = lake_area(k,r) + a
              end if
            end do
          end if
        end if
      end if
    end do
  end do
  deallocate(cell_order, bstart, dsu)

  ! Parent depressions (children before parents)
  Sub = 0   ! area of the cells whose entry depression is in the sub-tree of k (below Zmin of its parent)
  do i = 1, ndep
    child = child_list(i)
    if (child <= nleaf) then
      Sub(child) = Tsum(child)
    else
      Zmin = Zmin_list(i)
      imin = min(nodes(child)%left%imin, nodes(child)%right%imin)
      imax = max(nodes(child)%left%imax, nodes(child)%right%imax)
      jmin = min(nodes(child)%left%jmin, nodes(child)%right%jmin)
      jmax = max(nodes(child)%left%jmax, nodes(child)%right%jmax)
      nodes(child)%imin = imin
      nodes(child)%imax = imax
      nodes(child)%jmin = jmin
      nodes(child)%jmax = jmax

      imin_list(i) = imin
      imax_list(i) = imax
      jmin_list(i) = jmin
      jmax_list(i) = jmax
      pit_cell_list(i) = 0
      lake_long_list(i) = 0
      lake_lat_list(i) = 0

      nodes(child)%Aw = nodes(child)%left%Aw + nodes(child)%right%Aw
      watershed_area_list(i) = nodes(child)%Aw

      if (computed(i)) then
        Alow = Sub(nodes(child)%left%id) + Sub(nodes(child)%right%id) + Lsum(child)
        lake_vol(1,i) = 0
        lake_area(1,i) = nodes(child)%left%Amax + nodes(child)%right%Amax
        lake_elev(1,i) = Zmin
        do j = 2, lake_dislen
          lake_vol(j,i) = (lake_elev(j,i) - Zmin) * Alow + lake_vol(j,i)
          lake_area(j,i) = Alow + lake_area(j,i)
        end do
        Vmax_list(i) = lake_vol(lake_dislen,i)
        nodes(child)%Vmax = lake_vol(lake_dislen,i)
        nodes(child)%Amax = lake_area(lake_dislen,i)
      else
        Vmax_list(i) = 0
        watershed_area_list(i) = nodes(child)%left%Aw + nodes(child)%right%Aw
        lake_vol(:,i) = 0
        lake_area(:,i) = nodes(child)%left%Amax + nodes(child)%right%Amax
        nodes(child)%Amax = lake_area(1,i)
      end if
      Sub(child) = Sub(nodes(child)%left%id) + Sub(nodes(child)%right%id) + Tsum(child)
      check_list(i) = 1
    end if
  end do
  deallocate(Tsum, Lsum, Sub, computed)
  call system_clock(t1)
  write(*,*) "   time (s):", real(t1-t0)/real(ir)


  print_count = 0
  write(*,*) "5. Complete Database from tree informations"
  do i = 1, ndep
    child = child_list(i)
    if (child > nleaf) then
      rchild_list(i) = nodes(child)%right%id
      lchild_list(i) = nodes(child)%left%id
      Vmin_list(i) = 0
    end if 
  end do

  if (sum(check_list) == size(check_list)) then
    ! label and elevation are not modified: no need to reload them

    call SaveNETCDF(trim(outfolder) // trim("/hydrologic_database_final.nc"), lonlen, latlen, deplen, lake_dislen, &
    latitude, longitude, child_list, lchild_list, rchild_list, parent_list, check_list, &
    pit_cell_list, out_cell_list, sibling_list, geolink_list, &
    imin_list, imax_list, jmin_list, jmax_list, i_out_list, j_out_list, &
    Vmax_list, Vmin_list, lake_lat_list, lake_long_list, &
    Zmin_list, Zmax_list, watershed_area_list, &
    lake_vol, lake_area, lake_elev, label, elevation)
    ! Open file
    open(unit=10, file=trim(outfolder) // trim("/hydrologic_database.log"), status='replace', action='write', &
        form='formatted')
    ! Write in file
    write(10,'(A)')  "Construction of the hydrological database is completed!"
    ! Close file
    close(10)
  end if 

  contains

!=======================================================================
subroutine SaveNETCDF(file, lonlen, latlen, deplen, lake_dislen, latitude, longitude,&
  child_list, lchild_list, rchild_list, parent_list, check_list, &
  pit_cell_list, out_cell_list, sibling_list, geolink_list, &
  imin_list, imax_list, jmin_list, jmax_list, i_out_list, j_out_list, &
  Vmax_list, Vmin_list, lake_lat_list, lake_long_list, &
  Zmin_list, Zmax_list, watershed_area_list, &
  lake_vol, lake_area, lake_elev, label, elevation)
!-----------------------------------------------------------------------
! NAME
!     SaveNETCDF
!
! DESCRIPTION
!     Write the completed hydrologic database in a NetCDF file.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  implicit none
  ! Arguments
  character(len=*), intent(in) :: file ! Output NetCDF file
  ! Dimensions
  integer, intent(in) :: lonlen      ! Number of longitudes
  integer, intent(in) :: latlen      ! Number of latitudes
  integer, intent(in) :: deplen      ! Number of depressions
  integer, intent(in) :: lake_dislen ! Number of filling levels of the lake profiles
  ! Variables
  double precision, intent(in) :: longitude(lonlen)              ! Longitudes [deg]
  double precision, intent(in) :: latitude(latlen)               ! Latitudes [deg]
  integer,          intent(in) :: child_list(deplen)             ! Identifiers of the depressions
  integer,          intent(in) :: lchild_list(deplen)            ! Left children of the depressions
  integer,          intent(in) :: rchild_list(deplen)            ! Right children of the depressions
  integer,          intent(in) :: parent_list(deplen)            ! Parents of the depressions
  integer,          intent(in) :: check_list(deplen)             ! Completion flag of the depressions (1 or 0)
  integer,          intent(in) :: pit_cell_list(deplen)          ! Pit cells of the depressions
  integer,          intent(in) :: out_cell_list(deplen)          ! Outlet cells of the depressions
  integer,          intent(in) :: sibling_list(deplen)           ! Siblings of the depressions
  integer,          intent(in) :: geolink_list(deplen)           ! Downstream depressions (geolink)
  integer,          intent(in) :: imin_list(deplen)              ! First longitude index of the bounding boxes
  integer,          intent(in) :: imax_list(deplen)              ! Last longitude index of the bounding boxes
  integer,          intent(in) :: jmin_list(deplen)              ! First latitude index of the bounding boxes
  integer,          intent(in) :: jmax_list(deplen)              ! Last latitude index of the bounding boxes
  integer,          intent(in) :: i_out_list(deplen)             ! Longitude index of the outlet cells
  integer,          intent(in) :: j_out_list(deplen)             ! Latitude index of the outlet cells
  double precision, intent(in) :: Vmax_list(deplen)              ! Maximum volume of the depressions [m^3]
  double precision, intent(in) :: Vmin_list(deplen)              ! Minimum volume of the depressions [m^3]
  double precision, intent(in) :: lake_lat_list(deplen)          ! Latitude of the lakes [deg]
  double precision, intent(in) :: lake_long_list(deplen)         ! Longitude of the lakes [deg]
  double precision, intent(in) :: Zmin_list(deplen)              ! Minimum elevation of the depressions [m]
  double precision, intent(in) :: Zmax_list(deplen)              ! Maximum elevation of the depressions [m]
  double precision, intent(in) :: watershed_area_list(deplen)    ! Watershed area of the depressions [m^2]
  double precision, intent(in) :: lake_vol(lake_dislen, deplen)  ! Lake volume profiles [m^3]
  double precision, intent(in) :: lake_area(lake_dislen, deplen) ! Lake area profiles [m^2]
  double precision, intent(in) :: lake_elev(lake_dislen, deplen) ! Lake elevation profiles [m]
  ! Maps
  double precision, intent(in) :: elevation(lonlen,latlen) ! Elevation map [m]
  integer,          intent(in) :: label(lonlen,latlen)     ! Label map of the leaf depressions
  ! Netcdf stuff
  integer :: status ! NetCDF return code
  integer :: ncid ! NetCDF file ID
  integer :: dimid ! NetCDF dimension ID
  integer :: varid ! NetCDF variable ID
  integer :: lon_dimid,lat_dimid ! NetCDF dimension IDs
  integer :: lon_varid,lat_varid ! NetCDF variable IDs
  integer :: dimid_lon, dimid_lat, dimid_dep, dimid_lakedis, dimid_maxleafid, dimid_maxdepid,dimid_nleaf, dimid_ndep
  integer :: varid_lon, varid_lat, varid_child, varid_rchild, varid_lchild, varid_parent
  integer :: varid_sibling, varid_geolink, varid_Vmax, varid_Vmin, varid_check
  integer :: varid_out_cell, varid_pit_cell, varid_lake_lon, varid_lake_lat, varid_Zmin, varid_Zmax, varid_watershed_area
  integer :: varid_elevation, varid_label, varid_area, varid_flowdirs, varid_cell_number
  integer :: varid_imin, varid_imax, varid_jmin, varid_jmax, varid_i_out, varid_j_out
  integer :: varid_lake_perc, varid_lake_elev, varid_lake_area, varid_lake_vol

  write(*,*) "#######################################"
  write(*,*) "###### Save Hydrologic Database #######"
  write(*,*) "#######################################"

  write(*,*) file

  !Create file.nc
  status = nf90_create(file, NF90_NETCDF4, ncid)
  
  ! Define dimensions
  status = nf90_put_att(ncid, NF90_GLOBAL, "R", radius)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing attribute 'R': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "lon", lonlen, dimid_lon)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lon': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "lat", latlen, dimid_lat)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lat': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "depression", deplen, dimid_dep)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'depression': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "lake", lake_dislen, dimid_lakedis)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lake': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "maxleafid", maxval(label), dimid_maxleafid)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'maxleafid': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "maxdepid", maxval(child_list)+1, dimid_maxdepid)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'maxdepid': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "nleaf", sum(merge(1,0,child_list<=nleaf)), dimid_nleaf)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'nleaf': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_dim(ncid, "ndep", size(child_list), dimid_ndep)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'ndep': ", trim(nf90_strerror(status))
  end if

  !Define variables
  status = nf90_def_var(ncid, "longitude", NF90_FLOAT, dimid_lon, varid_lon)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining variable 'longitude': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "latitude", NF90_FLOAT, dimid_lat, varid_lat)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'latitude': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "child", NF90_INT, dimid_dep, varid_child)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'child': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "rchild", NF90_INT, dimid_dep, varid_rchild)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'rchild': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "lchild", NF90_INT, dimid_dep, varid_lchild)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lchild': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "parent", NF90_INT, dimid_dep, varid_parent)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'parent': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "sibling", NF90_INT, dimid_dep, varid_sibling)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'sibling': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "geolink", NF90_INT, dimid_dep, varid_geolink)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'geolink': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "Vmax", NF90_FLOAT, dimid_dep, varid_Vmax)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'Vmax': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "Vmin", NF90_FLOAT, dimid_dep, varid_Vmin)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'Vmin': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "pit_cell", NF90_INT, dimid_dep, varid_pit_cell)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'pit_cell': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "out_cell", NF90_INT, dimid_dep, varid_out_cell)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'out_cell': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "lake_lon", NF90_FLOAT, dimid_dep, varid_lake_lon)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lake_lon': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "lake_lat", NF90_FLOAT, dimid_dep, varid_lake_lat)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lake_lat': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "Zmin", NF90_FLOAT, dimid_dep, varid_Zmin)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'Zmin': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "Zmax", NF90_FLOAT, dimid_dep, varid_Zmax)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'Zmax': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "watershed_area", NF90_FLOAT, dimid_dep, varid_watershed_area)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'watershed_area': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "elevation", NF90_FLOAT, [dimid_lon, dimid_lat], varid_elevation)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'elevation': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "label", NF90_INT, [dimid_lon, dimid_lat], varid_label)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'label': ", trim(nf90_strerror(status))
  end if
  !status = nf90_def_var(ncid, "area", NF90_FLOAT, [dimid_lon, dimid_lat], varid_area)
  !status = nf90_def_var(ncid, "flowdirs", NF90_INT, [dimid_lon, dimid_lat], varid_flowdirs)
  !status = nf90_def_var(ncid, "cell_number", NF90_INT, [dimid_lon, dimid_lat], varid_cell_number)
  status = nf90_def_var(ncid, "imin", NF90_INT, dimid_dep, varid_imin)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'imin': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "imax", NF90_INT, dimid_dep, varid_imax)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'imax': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "jmin", NF90_INT, dimid_dep, varid_jmin)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'jmin': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "jmax", NF90_INT, dimid_dep, varid_jmax)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'jmax': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "i_out", NF90_INT, dimid_dep, varid_i_out)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'i_out': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "j_out", NF90_INT, dimid_dep, varid_j_out)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'j_out': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "lake_perc", NF90_FLOAT, dimid_lakedis, varid_lake_perc)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lake_perc': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "lake_elev", NF90_FLOAT, [dimid_lakedis, dimid_dep], varid_lake_elev)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lake_elev': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "lake_vol", NF90_FLOAT, [dimid_lakedis, dimid_dep], varid_lake_vol)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lake_vol': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "lake_area", NF90_FLOAT, [dimid_lakedis, dimid_dep], varid_lake_area)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'lake_area': ", trim(nf90_strerror(status))
  end if
  status = nf90_def_var(ncid, "check", NF90_INT, dimid_dep, varid_check)
  if (status /= NF90_NOERR) then
    write(*, *) "Error defining dimension 'check': ", trim(nf90_strerror(status))
  end if
  status = nf90_enddef(ncid)
  if (status /= NF90_NOERR) then
    write(*, *) "Error ending define mode: ", trim(nf90_strerror(status))
  end if

  
  status = nf90_put_var(ncid, varid_lon, longitude)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'longitude' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lat, latitude)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'latitude' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_child, child_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'child' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_rchild, rchild_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'rchild' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lchild, lchild_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'lchild' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_parent, parent_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'parent' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_sibling, sibling_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'sibling' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_geolink, geolink_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'geolink' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_Vmax, Vmax_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'Vmax' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_Vmin, Vmin_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'Vmin' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_pit_cell, pit_cell_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'pit_cell' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_out_cell, out_cell_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'out_cell' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lake_lon, lake_long_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'lake_lon' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lake_lat, lake_lat_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'lake_lat' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_Zmin, Zmin_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'Zmin' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_Zmax, Zmax_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'Zmax' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_watershed_area, watershed_area_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'watershed_area' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_elevation, elevation)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'elevation' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_label, label)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'label' data: ", trim(nf90_strerror(status))
  end if
  !status = nf90_put_var(ncid, varid_area, area)
  !status = nf90_put_var(ncid, varid_flowdirs, flowdirs)
  !status = nf90_put_var(ncid, varid_cell_number, cell_number)
  status = nf90_put_var(ncid, varid_imin, imin_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'imin' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_imax, imax_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'imax' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_jmin, jmin_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'jmin' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_jmax, jmax_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'jmax' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_i_out, i_out_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'i_out' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_j_out, j_out_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'j_out' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lake_perc, lake_perc)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'lake_perc' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lake_elev, real(lake_elev,8))
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'lake_elev' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lake_vol, real(lake_vol,8))
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'lake_vol' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_lake_area, real(lake_area,8))
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'lake_area' data: ", trim(nf90_strerror(status))
  end if
  status = nf90_put_var(ncid, varid_check, check_list)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'check' data: ", trim(nf90_strerror(status))
  end if
  
  status = nf90_close(ncid)
  if (status /= NF90_NOERR) then
    write(*, *) "Error writing 'Vmax' data: ", trim(nf90_strerror(status))
  end if
  
  write(*,*) "NETCDF saved !"

end subroutine SaveNETCDF

!=======================================================================
subroutine progress(i, N, printed_values)
!-----------------------------------------------------------------------
! NAME
!     progress
!
! DESCRIPTION
!     Print a progress bar every 10% of a loop.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  implicit none
  ! Arguments
  integer, intent(in) :: i                     ! Current iteration
  integer, intent(in) :: N                     ! Number of iterations
  integer, intent(inout) :: printed_values(10) ! Flags of the steps of 10% already printed (1 or 0)
  integer :: k
  real :: div
  character(len=17)::bar="???% |          |"
  integer, dimension(10) :: test=[10,20,30,40,50,60,70,80,90,100]
  logical, dimension(10) :: test_log
  integer, dimension(1) :: j
  div = real(i)/real(N)
  
  if (any(test == int(div * 100))) then
      j = findloc(test, int(div * 100))
      if (printed_values(j(1)) == 0) then
        write(unit=bar(1:3),fmt="(I3)") int(div * 100)
        test_log = test <= div*100
        do k=1, 10
          if (test_log(k)) then
            bar(6+k:6+k)="*"
          else
            bar(6+k:6+k)=" "
          end if
        enddo
        ! print the progress bar.
        write(unit=6,fmt="(a1,a1,a17)") '+',char(13), bar
        printed_values(j(1)) = 1
      end if 
  end if
  return
end subroutine progress

!=======================================================================
subroutine FindLeaves(node, lonlen, latlen, num_leaves, id_leaves, imin, imax, jmin, jmax)
!-----------------------------------------------------------------------
! NAME
!     FindLeaves
!
! DESCRIPTION
!     Find the leaf depressions below a node of the tree and the bounding box
!     enclosing them.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  implicit none
  ! Arguments
  type(TreeNode), target, intent(in) :: node        ! Node of the tree
  integer, intent(in) :: lonlen                     ! Number of longitudes
  integer, intent(in) :: latlen                     ! Number of latitudes
  integer, intent(out) :: num_leaves                ! Number of leaves found
  integer, intent(out) :: imin, imax                ! Longitude indices of the bounding box of the leaves
  integer, intent(out) :: jmin, jmax                ! Latitude indices of the bounding box of the leaves
  integer, intent(out), allocatable :: id_leaves(:) ! Identifiers of the leaves found
  type(LeaveNode), pointer :: list_leaves, current
  type(LeaveNode), pointer :: temp
  integer :: i
  
  ! Initialize the count
  num_leaves = 0
  imin = lonlen + 1
  imax = 0
  jmin = latlen + 1
  jmax = 0

  ! Call the recursive function
  call RecursiveFindLeaves(node, num_leaves, list_leaves)
  allocate(temp)
  allocate(id_leaves(num_leaves))
  current => list_leaves
  do i = 1, num_leaves
    id_leaves(i) = current%id
    imin = min(imin, current%imin)
    imax = max(imax, current%imax)
    jmin = min(jmin, current%jmin)
    jmax = max(jmax, current%jmax)
    current => current%next
  end do
  
  ! Deallocate the memory for the linked list of leaves
  current => list_leaves
  do i = 1, num_leaves
    temp => current%next
    deallocate(current)
    current => temp
  end do
end subroutine FindLeaves

!=======================================================================
recursive subroutine RecursiveFindLeaves(node, num_leaves, list_leaves)
!-----------------------------------------------------------------------
! NAME
!     RecursiveFindLeaves
!
! DESCRIPTION
!     Recursively add the leaves below a node of the tree to a linked list.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  ! Arguments
  type(TreeNode), target, intent(in) :: node             ! Node of the tree
  integer, intent(inout) :: num_leaves                   ! Number of leaves, incremented by the ones found
  type(LeaveNode), pointer, intent(inout) :: list_leaves ! Linked list of the leaves, the ones found are added on top
  type(LeaveNode), pointer :: current
  
  ! Initialise current à NULL()
  current => NULL()
  
  ! Check if the node is a leaf
  if (associated(node%left) .or. associated(node%right)) then
    ! Recursively count leaves in left and right subtrees
    if (associated(node%left)) then
      call RecursiveFindLeaves(node%left, num_leaves, list_leaves)
    end if
    if (associated(node%right)) then
      call RecursiveFindLeaves(node%right, num_leaves, list_leaves)
    end if
  else
  ! If both left and right are not associated, it's a leaf
    allocate(current)
    current%id = node%id
    current%imin = node%imin
    current%imax = node%imax
    current%jmin = node%jmin
    current%jmax = node%jmax
    current%next => list_leaves
    list_leaves => current
    num_leaves = num_leaves + 1
  end if
end subroutine RecursiveFindLeaves
    

!=======================================================================
subroutine init_execution_timer(cpu_start_time)
!-----------------------------------------------------------------------
! NAME
!     init_execution_timer
!
! DESCRIPTION
!     Record the CPU start time.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  ! Arguments
  real(8), intent(out) :: cpu_start_time ! CPU start time [s]

  ! Record the CPU start time
  call CPU_TIME(cpu_start_time)
end subroutine init_execution_timer

!=======================================================================
subroutine print_execution_time(cpu_start_time)
!-----------------------------------------------------------------------
! NAME
!     print_execution_time
!
! DESCRIPTION
!     Print the CPU time elapsed since a start time.
!
! AUTHORS & DATE
!     A. Gauvain, 08/2023
!
! NOTES
!
!-----------------------------------------------------------------------
  ! Arguments
  real(8), intent(in) :: cpu_start_time ! CPU start time [s]
  real(8) :: cpu_end_time, cpu_elapsed_time

  ! Record the CPU end time
  call CPU_TIME(cpu_end_time)

  ! Calculate the CPU elapsed time in seconds
  cpu_elapsed_time = cpu_end_time - cpu_start_time

  ! Print the CPU execution time
  print '("Time = ",f9.6," seconds.")',cpu_elapsed_time
end subroutine print_execution_time


end program database