! ----------------------------------------------------------------------------- ! Copyright © 2025 Alexandre Gauvain ! Contributor(s): Alexandre Gauvain ! Email: alexandre.gauvain@lmd.ipsl.fr or alexandre.gauvain.ag@gmail.com ! Affiliation: LMD/IPSL, CNRS, Sorbonne Université, Paris, France ! Created: 2023-09-03 ! Modified: 2025-06-10 ! ! This software is a computer program whose purpose is to model hydrological ! processes at planetary scale, including water flow, storage, climate ! interactions, and groundwater. ! ! This software is governed by the CeCILL license under ! French law and abiding by the rules of distribution of free software. ! You can use, modify and/or redistribute the software under the terms of ! the CeCILL license as circulated by CEA, CNRS and INRIA ! at the following URL "http://www.cecill.info". ! ! As a counterpart to the access to the source code and rights to copy, ! modify and redistribute granted by the license, users are provided only ! with a limited warranty and the software's author, the holder of the ! economic rights, and the successive licensors have only limited liability. ! ! In this respect, the user's attention is drawn to the risks associated ! with loading, using, modifying and/or developing or reproducing the ! software by the user in light of its specific status of free software, ! that may mean that it is complicated to manipulate, and that also therefore ! means that it is reserved for developers and experienced professionals ! having in-depth computer knowledge. Users are therefore encouraged to load ! and test the software's suitability as regards their requirements in ! conditions enabling the security of their systems and/or data to be ensured ! and, more generally, to use and operate it in the same conditions as ! regards security. ! ! The fact that you are presently reading this means that you have had ! knowledge of the CeCILL license and that you accept its terms. ! ! ----------------------------------------------------------------------------- !=============================================================================== ! MODULE: TreeNodeModule !=============================================================================== ! DESCRIPTION: ! Defines the TreeNode type and the BuildTree subroutine for representing and ! building a tree of hydrological depressions and their relationships. !=============================================================================== module hydro_tree use, intrinsic :: iso_fortran_env, only : real64, real128 use, intrinsic :: omp_lib use hydro_database, only : IDdep, ndep, lake_dislen, nleaf, child_list, parent_list, & lchild_list, rchild_list, sibling_list, geolink_list, & imin_list, imax_list, jmin_list, jmax_list, & Vmax_list, lake_lat_list, lake_long_list, Zmin_list, Zmax_list, & watershed_area_list, lake_elev, lake_vol, lake_area, & latitude, longitude implicit none real(real128) :: V_trans_dep=1e-3_real128 ! Volume of transient depressions (m^3) type :: TreeNode integer :: child ! Index of this node in the tree (dimensionless) integer :: parent ! Index of parent node (dimensionless) integer :: lchild ! Index of left child node (dimensionless) integer :: rchild ! Index of right child node (dimensionless) integer :: geolink ! Index of geolinked node (dimensionless) integer :: sibling ! Index of sibling node (dimensionless) integer :: imin ! Minimum i index (grid index) integer :: imax ! Maximum i index (grid index) integer :: jmin ! Minimum j index (grid index) integer :: jmax ! Maximum j index (grid index) integer :: bypass ! Index for bypass routing (dimensionless) integer :: top_active ! Index of upper active depression (dimensionless) logical :: is_leaf ! True if node is a leaf (dimensionless) integer :: gcm_idx ! Index in GCM grid (dimensionless) !$ integer(omp_lock_kind) :: lock ! OpenMP lock (dimensionless) integer :: active = 0 ! Active flag (1 or 0) integer :: upper_active = 0 ! ID of upper active depression (dimensionless) integer :: fixed = 0 ! Fixed flag (1 or 0) integer :: processed = 0 ! Processed flag (1 or 0) integer :: n_down = 0 ! Number of downstream nodes (dimensionless) double precision :: lake_lon ! Longitude of lake center (degrees_east) double precision :: lake_lat ! Latitude of lake center (degrees_north) double precision :: Zmin ! Minimum elevation (m) double precision :: Zmax ! Maximum elevation (m) double precision :: lat_min ! Minimum latitude (degrees) double precision :: lat_max ! Maximum latitude (degrees) double precision :: lon_min ! Minimum longitude (degrees) double precision :: lon_max ! Maximum longitude (degrees) real(real128) :: Z = 0 ! Current elevation (m) real(real128) :: V = 0 ! Current volume (m^3) real(real128) :: A = 0 ! Current area (m^2) real(real128) :: Vout = 0 ! Outflow volume (m^3) real(real128) :: Vtop = 0 ! Total volume (m^3) real(real128) :: Vin = 0 ! Inflow volume (m^3) real(real128) :: dV = 0 ! Change in volume (m^3) real(real128) :: Tsurf = 0 ! Surface temperature (K) real(real128) :: E = 0 ! Evaporation (m/y, specify in context) real(real128) :: P = 0 ! Precipitation (m/y, specify in context) real(real128) :: Runoff = 0 ! Runoff (m/y, specify in context) real(real128) :: Vmax = 0 ! Maximum volume (m^3) real(real128) :: Vevap = 0 ! Maximum volume (m^3) real(real128) :: watershed_area = 0 ! Watershed area (m^2) real(real128), allocatable :: lake_elev(:) ! Lake elevation profile (m) real(real128), allocatable :: lake_vol(:) ! Lake volume profile (m^3) real(real128), allocatable :: lake_area(:) ! Lake area profile (m^2) ! Pointers to other nodes type(TreeNode), pointer :: self_pointer => null() ! Pointer to self (TreeNode) type(TreeNode), pointer :: parent_pointer => null() ! Pointer to parent (TreeNode) type(TreeNode), pointer :: sibling_pointer => null() ! Pointer to sibling (TreeNode) type(TreeNode), pointer :: geolink_pointer => null() ! Pointer to geolinked node (TreeNode) type(TreeNode), pointer :: left_pointer => null() ! Pointer to left child (TreeNode) type(TreeNode), pointer :: right_pointer => null() ! Pointer to right child (TreeNode) end type TreeNode ! Array of TreeNode structures type(TreeNode), target, allocatable :: nodes(:) contains !------------------------------------------------------------------------------- ! SUBROUTINE: BuildTree !------------------------------------------------------------------------------- ! DESCRIPTION: ! Initializes and links the nodes of the tree based on input arrays describing ! the topology and properties of each depression. ! Sets up pointers and properties for each TreeNode in the model. ! Handles allocation of arrays and OpenMP locks for thread safety. ! Computes derived properties such as watershed area and node relationships. ! Issues warnings for missing watershed areas. !------------------------------------------------------------------------------- ! ARGUMENTS: ! - IDdep [in] : Number of nodes in the tree (integer) ! - ndep [in] : Number of depressions (integer) ! - lake_dislen [in] : Length of lake discretization arrays (integer) ! - nleaf [in] : Number of leaf nodes (integer) ! - lake_vol [in] : Lake volumes (lake_dislen, ndep) (double precision) ! - lake_area [in] : Lake areas (lake_dislen, ndep) (double precision) ! - lake_elev [in] : Lake elevations (lake_dislen, ndep) (double precision) ! - child_list [in] : Child node indices (ndep) (integer) ! - lchild_list [in] : Left child node indices (ndep) (integer) ! - rchild_list [in] : Right child node indices (ndep) (integer) ! - parent_list [in] : Parent node indices (ndep) (integer) ! - sibling_list [in] : Sibling node indices (ndep) (integer) ! - geolink_list [in] : Geolink node indices (ndep) (integer) ! - imin_list [in] : Minimum i indices (ndep) (integer) ! - imax_list [in] : Maximum i indices (ndep) (integer) ! - jmin_list [in] : Minimum j indices (ndep) (integer) ! - jmax_list [in] : Maximum j indices (ndep) (integer) ! - Vmax_list [in] : Maximum volumes (ndep) (double precision) ! - lake_lat_list [in] : Lake latitudes (ndep) (double precision) ! - lake_long_list [in] : Lake longitudes (ndep) (double precision) ! - Zmin_list [in] : Minimum elevations (ndep) (double precision) ! - Zmax_list [in] : Maximum elevations (ndep) (double precision) ! - watershed_area_list[in] : Watershed areas (ndep) (double precision) ! - nodes [out] : Array of TreeNode structures (IDdep) (type(TreeNode)) !------------------------------------------------------------------------------- subroutine allocate_tree(this) implicit none type(TreeNode), target, allocatable, intent(inout) :: this(:) allocate(this(IDdep)) end subroutine allocate_tree subroutine build_tree(this) implicit none type(TreeNode), target, intent(inout) :: this(:) ! Out tree integer :: i, count !$OMP PARALLEL !$OMP DO do i = 1, ndep this(child_list(i))%child = child_list(i) this(child_list(i))%bypass = child_list(i) this(child_list(i))%parent = parent_list(i) this(child_list(i))%lchild = lchild_list(i) this(child_list(i))%rchild = rchild_list(i) this(child_list(i))%geolink = geolink_list(i) this(child_list(i))%sibling = sibling_list(i) this(child_list(i))%imin = imin_list(i) this(child_list(i))%imax = imax_list(i) this(child_list(i))%jmin = jmin_list(i) this(child_list(i))%jmax = jmax_list(i) this(child_list(i))%lake_lon = lake_long_list(i) this(child_list(i))%lake_lat = lake_lat_list(i) this(child_list(i))%Zmax = Zmax_list(i) this(child_list(i))%Zmin = Zmin_list(i) this(child_list(i))%Vmax = Vmax_list(i) this(child_list(i))%watershed_area = watershed_area_list(i) this(child_list(i))%lat_min = latitude(jmax_list(i)) this(child_list(i))%lat_max = latitude(jmin_list(i)) this(child_list(i))%lon_min = longitude(imin_list(i)) this(child_list(i))%lon_max = longitude(imax_list(i)) + 1 !write(*,*) 'lon_min, lon_max, lat_min, lat_max :', this(child_list(i))%lon_min, this(child_list(i))%lon_max, this(child_list(i))%lat_min, this(child_list(i))%lat_max this(child_list(i))%V = 0.0_real128 allocate(this(child_list(i))%lake_elev(lake_dislen)) allocate(this(child_list(i))%lake_vol(lake_dislen)) allocate(this(child_list(i))%lake_area(lake_dislen)) this(child_list(i))%lake_elev = lake_elev(:,i) this(child_list(i))%lake_vol = lake_vol(:,i) !if (lake_area(1,i) == lake_area(lake_dislen,i)) then !this(child_list(i))%lake_area = [(lake_area(1,i) + (lake_area(lake_dislen,i)+1 - lake_area(1,i)) * (j-1) / (lake_dislen-1), j = 1, lake_dislen)] !else this(child_list(i))%lake_area = lake_area(:,i) !end if ! Pointers this(child_list(i))%self_pointer => this(child_list(i)) this(child_list(i))%parent_pointer => this(parent_list(i)) this(child_list(i))%sibling_pointer => this(sibling_list(i)) this(child_list(i))%geolink_pointer => this(geolink_list(i)) if (child_list(i) > nleaf) then this(child_list(i))%right_pointer => this(rchild_list(i)) this(child_list(i))%left_pointer => this(lchild_list(i)) end if if (child_list(i) <= nleaf) then this(child_list(i))%active = 1 this(child_list(i))%is_leaf = .true. else this(child_list(i))%is_leaf = .false. end if this(child_list(i))%n_down = 1 !$ call omp_init_lock(this(child_list(i))%lock) end do !$OMP END DO !$OMP DO do i = 1, ndep if (this(child_list(i))%Vmax == 0.) then !.and. child_list(i) > nleaf) then this(child_list(i))%Vmax = V_trans_dep !AG end if end do !$OMP END DO !$OMP END PARALLEL count = 0 do i = 1, ndep if (this(child_list(i))%Vmax == V_trans_dep .and. child_list(i) > nleaf) then !AG !write(*,*) child_list(i),'/', ndep, nleaf, this(child_list(i))%Vmax, V_trans_dep,this(child_list(i))%is_leaf this(child_list(i))%A = 0 this(child_list(i))%Z = this(child_list(i))%right_pointer%Zmax if (this(child_list(i))%right_pointer%Vmax == V_trans_dep) then !AG this(child_list(i))%A = this(child_list(i))%A + this(child_list(i))%right_pointer%A else this(child_list(i))%A = this(child_list(i))%A + this(child_list(i))%right_pointer%lake_area(lake_dislen) end if if (this(child_list(i))%left_pointer%Vmax == V_trans_dep) then !AG this(child_list(i))%A = this(child_list(i))%A + this(child_list(i))%left_pointer%A else this(child_list(i))%A = this(child_list(i))%A + this(child_list(i))%left_pointer%lake_area(lake_dislen) end if end if if (child_list(i) > nleaf) then if (this(child_list(i))%left_pointer%watershed_area == 0.) then write(*,*) 'Warning: lchild Aw equal to 0. Depression :', child_list(i), this(child_list(i))%left_pointer%child, this(child_list(i))%left_pointer%watershed_area end if if (this(child_list(i))%right_pointer%watershed_area == 0.) then write(*,*) 'Warning: rchild Aw equal to 0. Depression :', child_list(i), this(child_list(i))%right_pointer%child end if this(child_list(i))%watershed_area = this(child_list(i))%left_pointer%watershed_area + this(child_list(i))%right_pointer%watershed_area this(child_list(i))%n_down = this(child_list(i))%n_down + this(child_list(i))%left_pointer%n_down + this(child_list(i))%right_pointer%n_down end if end do end subroutine build_tree subroutine deallocate_nodes(this) implicit none type(TreeNode), target, allocatable, intent(inout) :: this(:) if (allocated(this)) then deallocate(this) end if end subroutine deallocate_nodes recursive subroutine get_leaves_below_node(child, leaf_ids, nleaf) implicit none integer, intent(in) :: child integer, allocatable, intent(inout) :: leaf_ids(:) integer, intent(inout) :: nleaf if (nodes(child)%is_leaf) then nleaf = nleaf + 1 leaf_ids(nleaf) = nodes(child)%child else call get_leaves_below_node(nodes(child)%left_pointer%child, leaf_ids, nleaf) call get_leaves_below_node(nodes(child)%right_pointer%child, leaf_ids, nleaf) end if end subroutine get_leaves_below_node recursive subroutine get_leaves_number_below_node(child, nleafs) implicit none integer, intent(in) :: child integer, intent(inout) :: nleafs if (nodes(child)%is_leaf) then nleafs = nleafs + 1 else call get_leaves_number_below_node(nodes(child)%left_pointer%child, nleafs) call get_leaves_number_below_node(nodes(child)%right_pointer%child, nleafs) end if end subroutine get_leaves_number_below_node end module hydro_tree