! ----------------------------------------------------------------------------- ! 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 hydro_core !=============================================================================== ! MODULE: Hydrology !=============================================================================== ! Hydrology module for simulating water flow, storage and evaporation in a ! depression hierarchy. ! ! List of subroutines: ! - OverFlowInto: Handles overflow of water from one depression to another. ! - InterpolateFromZ: Interpolates area and volume from elevation. ! - InterpolateFromV: Interpolates area and elevation from volume. ! - FixeWater: Fixes water in a depression and its descendants. ! - EvapLake: Removes water from a depression tree due to evaporation. ! - AvailableWater: Calculates available water for evaporation in a depression tree. !=============================================================================== use hydro_tree use hydro_database use hydro_diaghydro !$ use, intrinsic :: omp_lib use numerics, only: di, dp, qp, k4 implicit none real(real128) :: Z_interp, A_interp real(real128) :: extra_water, Vin, A_lac, Vin_store real(real128) :: Evap_real, Evap_potential, Evap_real_store, Evap_potential_store, area real(real128) :: convergence_indicator_A, convergence_indicator_V real(real128) :: precip_save(1000) = 0 real(real128) :: planet_area_active, water_area_active real(real128) :: alpha = 1 integer(di) :: active_bypass = 1 real(real128) :: lon_val = 0, lat_val = 0 real(real128) :: fixed_area = 0. logical(k4) :: fixed = .false. integer :: status integer(di), parameter :: HOMOGENEOUS_GEL = 1_di ! GEL state is homogeneous integer(di), parameter :: LONLAT_GEL = 2_di ! GEL state is defined by a lon/lat point integer(di), parameter :: FIXED_OCEAN = 3_di ! GEL state is fixed to ocean contains !------------------------------------------------------------------------------- ! SUBROUTINE: OverFlowInto !------------------------------------------------------------------------------- ! Description: ! Recursively handles the overflow of excess water (Vexcess) from a given child ! depression to its bypass, sibling, geolink, or parent, updating the water ! volumes and bypass links accordingly. Uses OpenMP locks for thread safety.! !------------------------------------------------------------------------------- recursive subroutine OverFlowInto(nleaf, child, Vexcess, bypass, active_bypass, verbose) ! Arguments integer , intent(in) :: child ! Index of the current child node integer , intent(in) :: nleaf ! Number of leaf nodes integer , intent(in) :: active_bypass ! Flag for bypass activation (0 or 1)) real(real128), intent(inout) :: Vexcess ! Excess water volume to be distributed integer , intent(inout) :: bypass ! Index of the bypass node logical , intent(in) :: verbose ! Verbose flag ! Local variables real(real128) :: Vfree=0 ! Available volume in the depression integer :: child_tmp ! Temporary variable for child index !$ call omp_set_lock(nodes(child)%lock) ! Set lock for the current node if (verbose) then ! Debugging output write(*,*) nodes(child)%child, nodes(child)%bypass write(*,*) nodes(child)%V, '/', nodes(child)%Vmax, nodes(child)%fixed, Vexcess end if if (nodes(child)%fixed == 1) then ! Check if the node is fixed Vexcess = 0 ! No excess water to distribute bypass = nodes(child)%child !$ call omp_unset_lock(nodes(child)%lock) ! Release lock of the current node return ! Exit the subroutine end if ! Vexcess is upstream flow which comes from upstream depressions ! nodes(child)%Vin = nodes(child)%Vin + Vexcess if (nodes(child)%V >= nodes(child)%Vmax) then ! If depression is full(V>Vmax) nodes(child)%Vout = nodes(child)%Vout + Vexcess ! Add excess water to outflow else ! If depression is not full (V0) then ! If sibling depression is not full but contains water call OverFlowInto(nleaf, nodes(nodes(child)%sibling)%child, Vexcess, bypass, active_bypass, verbose=verbose) !$ call omp_set_lock(nodes(child)%lock) ! Set lock for the current node if (active_bypass == 1) then nodes(child)%bypass = bypass ! Update bypass else if (active_bypass == 0) then if (nodes(child)%active == 1) then bypass = nodes(child)%child ! Set bypass to the current child endif nodes(child)%bypass = bypass ! Update bypass end if !$ call omp_unset_lock(nodes(child)%lock) ! Release lock of the current node return ! Exit the subroutine end if !4th : Transfer excess water to the downstream depression if (nodes(nodes(child)%sibling)%V == 0.) then ! If sibling depression is not activate (V equal to 0) call OverFlowInto(nleaf, nodes(nodes(child)%geolink)%child, Vexcess, bypass, active_bypass, verbose=verbose) !$ call omp_set_lock(nodes(child)%lock) ! Set lock for the current node if (active_bypass == 1) then nodes(child)%bypass = bypass ! Update bypass else if (active_bypass == 0) then if (nodes(child)%active == 1) then bypass = nodes(child)%child ! Set bypass to the current child endif nodes(child)%bypass = bypass ! Update bypass end if !$ call omp_unset_lock(nodes(child)%lock) ! Release lock of the current node return ! Exit the subroutine end if ! 5th : Transfer excess water to the parent depression if sibling depression is full if (nodes(nodes(child)%sibling)%V >= nodes(nodes(child)%sibling)%Vmax) then !$ call omp_set_lock(nodes(nodes(child)%parent)%lock) ! Set lock for the parent node if (nodes(child)%parent_pointer%V < nodes(child)%parent_pointer%Vmax) then !$ call omp_set_lock(nodes(child)%lock) ! Set lock for the current node nodes(child)%active = 0 ! Set current child as inactive !$ call omp_unset_lock(nodes(child)%lock) ! Release lock of the current node !$ call omp_set_lock(nodes(nodes(child)%sibling)%lock) ! Set lock for the sibling node nodes(nodes(child)%sibling)%active = 0 ! Set sibling as inactive !$ call omp_unset_lock(nodes(nodes(child)%sibling)%lock)! Release lock of the sibling node nodes(nodes(child)%parent)%active = 1 end if !$ call omp_unset_lock(nodes(nodes(child)%parent)%lock) ! Release lock of the parent node child_tmp = nodes(child)%parent_pointer%child ! Get the child index of the parent node call OverFlowInto(nleaf, child_tmp, Vexcess, bypass, active_bypass, verbose=verbose) !$ call omp_set_lock(nodes(child)%lock) ! Set lock for the current node if (active_bypass == 1) then nodes(child)%bypass = bypass ! Update bypass else if (active_bypass == 0) then if (nodes(child)%active == 1) then bypass = nodes(child)%child ! Set bypass to the current child endif nodes(child)%bypass = bypass ! Update bypass end if !$ call omp_unset_lock(nodes(child)%lock) ! Release lock of the current node return ! Exit the subroutine end if write(*,*) 'Warning: No overflow path found for child:', child, 'Vexcess:', Vexcess end subroutine OverFlowInto subroutine ComputeDrainageArea() implicit none end subroutine ComputeDrainageArea subroutine OverflowAllDep() implicit none !Local variables integer :: i,j !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(extra_water, j) do i = 1,ndep if (nodes(child_list(i))%V > nodes(child_list(i))%Vmax) then !$ call omp_set_lock(nodes(child_list(i))%lock) extra_water = nodes(child_list(i))%V - nodes(child_list(i))%Vmax !extra_water = 0 nodes(child_list(i))%V = nodes(child_list(i))%Vmax j = child_list(i) !$ call omp_unset_lock(nodes(child_list(i))%lock) call OverFlowInto(nleaf, child_list(i), extra_water, j, active_bypass, verbose=.false.) end if end do !$OMP END DO end subroutine OverflowAllDep !------------------------------------------------------------------------------- ! SUBROUTINE: InterpolateFromZ !------------------------------------------------------------------------------- ! Description: ! Performs linear interpolation to compute area (A) and volume (V) of a ! depression given an elevation (Z), using provided lookup tables. !------------------------------------------------------------------------------- subroutine InterpolateFromZ(dis_len,Z_list, V_list, A_list, Z, A, V) ! Arguments integer, intent(in) :: dis_len ! Length of lookup tables double precision, intent(in) :: A_list(dis_len) ! Area values [m^2] double precision, intent(in) :: Z_list(dis_len) ! Elevation values [m] double precision, intent(in) :: V_list(dis_len) ! Volume values [m^3] double precision, intent(in) :: Z ! Elevation to interpolate [m] double precision, intent(out) :: A ! Interpolated area [m^2] double precision, intent(out) :: V ! Interpolated volume [m^3] ! Local variables integer :: i ! Index for the loop real(8) :: t ! Temporary variable for interpolation ! Find location for the linear interpolation do i = 1, dis_len - 1 if (Z >= Z_list(i) .and. Z <= Z_list(i + 1)) then exit end if end do ! linear interpolation t = (Z - Z_list(i)) / (Z_list(i + 1) - Z_list(i)) A = A_list(i) + t * (A_list(i + 1) - A_list(i)) V = V_list(i) + t * (V_list(i + 1) - V_list(i)) end subroutine InterpolateFromZ !------------------------------------------------------------------------------- ! SUBROUTINE: InterpolateFromV !------------------------------------------------------------------------------- ! Description: ! Performs linear interpolation to compute area (A) and elevation (Z) of a ! depression given a volume (V), using provided lookup tables. !------------------------------------------------------------------------------- subroutine InterpolateFromV(dis_len,Z_list, V_list, A_list, V, A, Z) ! Arguments integer, intent(in) :: dis_len ! Length of lookup tables real(real128), intent(in) :: A_list(dis_len) ! Area values [m^2] real(real128), intent(in) :: Z_list(dis_len) ! Elevation values [m] real(real128), intent(in) :: V_list(dis_len) ! Volume values [m^3] real(real128), intent(in) :: V ! Volume to interpolate [m^3] real(real128), intent(out) :: A ! Interpolated area [m^2] real(real128), intent(out) :: Z ! Interpolated elevation [m] ! Local variables integer :: i ! Index for the loop real(real128) :: t ! Temporary variable for interpolation ! Check if V is zero if (V == 0.) then A = A_list(1) Z = Z_list(1) return end if ! Check if V is within bounds if (V < V_list(1) .or. V > V_list(dis_len)) then print *, "Error in InterpolateFromV subroutine: V out of bounds" write(*,*) V, V_list(1), V_list(dis_len) stop end if ! find location for the linear interpolation do i = 1, dis_len - 1 if (V >= V_list(i) .and. V <= V_list(i + 1)) then exit end if end do ! linear interpolation t = (V - V_list(i)) / (V_list(i + 1) - V_list(i)) A = A_list(i) + (t * (A_list(i + 1) - A_list(i))) Z = Z_list(i) + (t * (Z_list(i + 1) - Z_list(i))) end subroutine InterpolateFromV !------------------------------------------------------------------------------- ! SUBROUTINE: FixeWater !------------------------------------------------------------------------------- ! Description: ! Recursively sets the 'fixed' flag for a given child node and all its descendants, ! indicating that water in these depressions is fixed and cannot be changed. !------------------------------------------------------------------------------- recursive subroutine FixeChild(child) ! Arguments integer, intent(in) :: child ! Index of the current child node nodes(child)%fixed = 1 if (nodes(child)%is_leaf .eqv. .true.) then return else call FixeChild(nodes(child)%right_pointer%child) call FixeChild(nodes(child)%left_pointer%child) return end if end subroutine FixeChild recursive subroutine FixeParent(parent) ! Arguments integer, intent(in) :: parent ! Index of the parent node nodes(parent)%fixed = 1 if (nodes(parent)%parent == 0) then return else call FixeParent(nodes(parent)%parent_pointer%child) return end if end subroutine FixeParent subroutine FixeWaterAll() integer :: i !!$OMP DO do i = 1,ndep if (nodes(child_list(i))%active == 1 .and. nodes(child_list(i))%V > 0) then call FixeChild(child_list(i)) end if end do !!$OMP END DO end subroutine FixeWaterAll !------------------------------------------------------------------------------- ! SUBROUTINE: EvapLake !------------------------------------------------------------------------------- ! Description: ! Removes a specified volume of water (Vevap) from a depression and its descendants ! due to evaporation, updating the actual evaporated volume (Vevap_real). ! Traverses the depression tree in breadth-first order. !------------------------------------------------------------------------------- subroutine EvapLake(child, Vevap, Vevap_real) ! Arguments integer, intent(in) :: child ! Index of the current child node real(real128), intent(inout) :: Vevap ! Volume to evaporate [m^3] real(real128), intent(inout) :: Vevap_real ! Real evaporated volume [m^3] ! Local variables integer :: current_child ! Current child index integer :: back ! Back index for queue integer :: front ! Front index for queue real(real128) :: Vevap_tmp ! Temporary evaporated volume [m^3] integer, allocatable :: queue(:) ! Queue for breadth-first traversal allocate(queue(nodes(child)%n_down+1)) ! Allocate queue for breadth-first traversal queue = 0 ! Initialize queue queue(1) = child ! Set the first element of the queue to the child index back = 1 ! Initialize back index front = 1 ! Initialize front index do while (front <= back .and. Vevap > 0) current_child = queue(front) front = front + 1 !$ call omp_set_lock(nodes(current_child)%lock) ! Set lock for the current node if (nodes(current_child)%V >= Vevap) then Vevap_real = Vevap_real + Vevap nodes(current_child)%V = nodes(current_child)%V - Vevap Vevap = 0 else Vevap_real = Vevap_real + nodes(current_child)%V Vevap_tmp = Vevap - nodes(current_child)%V nodes(current_child)%V = 0 if (nodes(current_child)%is_leaf .eqv. .false.) then nodes(current_child)%active = 0 ! Set current child as inactive nodes(current_child)%processed = 1 ! Set current child as processed nodes(current_child)%right_pointer%active = 1 ! Set right child as active nodes(current_child)%left_pointer%active = 1 ! Set left child as active end if Vevap = Vevap_tmp end if !$ call omp_unset_lock(nodes(current_child)%lock) ! Release lock of the current node if (nodes(current_child)%is_leaf .eqv. .false.) then back = back + 1 ! Add left child to the queue queue(back) = nodes(current_child)%left_pointer%child back = back + 1 ! Add right child to the queue queue(back) = nodes(current_child)%right_pointer%child end if end do deallocate(queue) end subroutine EvapLake !------------------------------------------------------------------------------- ! SUBROUTINE: ResetUpperActiveDep !------------------------------------------------------------------------------- ! Description: ! Resets the upper active depression for each child node in the hydrology model tree, ! setting it to 0 (indicating that the upper active depression is the child itself). !------------------------------------------------------------------------------- subroutine ResetUpperActiveDep() ! Local variables integer :: i !$OMP PARALLEL DO do i = 1,ndep nodes(child_list(i))%upper_active = 0 ! Reset upper active depression to the child itself end do !$OMP END PARALLEL DO end subroutine ResetUpperActiveDep !------------------------------------------------------------------------------- ! SUBROUTINE: RetrieveActiveUpperDep !------------------------------------------------------------------------------- ! Description: ! Retrieves the active upper depression for each lead node in the hydrology model tree. !------------------------------------------------------------------------------- subroutine ActiveUpperDep() ! Local variables integer :: i !$OMP PARALLEL DO do i = 1,ndep if (nodes(child_list(i))%active == 1) then call UpdateActiveUpperDep2LowerDep(child_list(i), child_list(i)) ! Update active upper depression for the current child and its descendants end if end do !$OMP END PARALLEL DO end subroutine ActiveUpperDep !------------------------------------------------------------------------------- ! SUBROUTINE: UpdateActiveUpperDep2LowerDep !------------------------------------------------------------------------------- ! Description: ! Recursively updates the active upper depression for a given child node and all ! its descendants in the hydrology model tree, ensuring that the active upper ! depression is correctly propagated down the tree. Uses OpenMP locks for thread safety. !------------------------------------------------------------------------------- recursive subroutine UpdateActiveUpperDep2LowerDep(child, upper_active) ! Arguments integer, intent(in) :: child ! Index of the current child node integer, intent(in) :: upper_active ! Index of the active upper depression !$ call omp_set_lock(nodes(child)%lock) ! Set lock for the current node nodes(child)%upper_active = upper_active ! Set the active upper depression to the current child !$ call omp_unset_lock(nodes(child)%lock) ! Release lock for the current node if (nodes(child)%is_leaf .eqv. .false.) then ! If the current child is not a leaf call UpdateActiveUpperDep2LowerDep(nodes(child)%left_pointer%child, upper_active) ! Update left child call UpdateActiveUpperDep2LowerDep(nodes(child)%right_pointer%child, upper_active) ! Update right child end if end subroutine UpdateActiveUpperDep2LowerDep !------------------------------------------------------------------------------- ! SUBROUTINE: AvailableWater !------------------------------------------------------------------------------- ! Description: ! Computes the available water that can be evaporated from a given node and ! its descendants in the hydrology model tree structure, without modifying ! the actual water content of the nodes. Traverses the depression tree in ! breadth-first order. !------------------------------------------------------------------------------- subroutine AvailableWater(child, Vevap, Vevap_real) ! Arguments integer, intent(in) :: child ! Index of the current child node real(real128), intent(inout) :: Vevap ! Volume requested for evaporation [m^3] real(real128), intent(inout) :: Vevap_real ! Real available volume for evaporation [m^3] ! Local variables integer :: current_child ! Current child index integer :: back ! Back index for queue integer :: front ! Front index for queue real(real128) :: Vevap_tmp ! Temporary evaporated volume [m^3] integer, allocatable :: queue(:) ! Queue for breadth-first traversal allocate(queue(nodes(child)%n_down+1)) queue = 0 queue(1) = child back = 1 front = 1 do while (front <= back .and. Vevap > 0) current_child = queue(front) front = front + 1 !$ call omp_set_lock(nodes(current_child)%lock) if (nodes(current_child)%V >= Vevap) then Vevap_real = Vevap_real + Vevap Vevap = 0 else Vevap_real = Vevap_real + nodes(current_child)%V Vevap_tmp = Vevap - nodes(current_child)%V Vevap = Vevap_tmp end if !$ call omp_unset_lock(nodes(current_child)%lock) if (nodes(current_child)%is_leaf .eqv. .false.) then back = back + 1 queue(back) = nodes(current_child)%left_pointer%child back = back + 1 queue(back) = nodes(current_child)%right_pointer%child end if end do deallocate(queue) end subroutine AvailableWater !------------------------------------------------------------------------------- ! SUBROUTINE: InitializeModel !------------------------------------------------------------------------------- ! Description: !------------------------------------------------------------------------------- subroutine initialize_model(init_file, init_gel_state, gel_init_lon, gel_init_lat, global_eq_lay, active_bypass_in) use numerics, only:qp, di use display, only: print_msg, LVL_NFO implicit none character(*), intent(in) :: init_file integer(di), intent(in) :: init_gel_state real(qp), intent(in) :: gel_init_lon real(qp), intent(in) :: gel_init_lat real(qp), intent(in) :: global_eq_lay integer(di), intent(in) :: active_bypass_in ! Local variables integer :: i, j, status logical :: here ! This subroutine initializes the hydrology model. call print_msg('> Initializing hydrology',LVL_NFO) active_bypass = active_bypass_in !1. Initialize from file if provided inquire(file = init_file,exist = here) if (here) then status = load_diaghydro_active(init_file, deplen, Z_list, V_list, A_list, Vout_list, Vin_list, active_list) !status = load_diaghydro_iteration(init_file, deplen, Z_list, V_list, A_list, Vout_list, Vin_list, active_list) do i = 1,ndep if (nodes(child_list(i))%Vmax /= V_trans_dep) then nodes(child_list(i))%A = A_list(i) end if nodes(child_list(i))%V = V_list(i) nodes(child_list(i))%Z = Z_list(i) nodes(child_list(i))%Vout = Vout_list(i) nodes(child_list(i))%Vin = Vin_list(i) nodes(child_list(i))%active = active_list(i) if (active_list(i) == 1) then call set_children_volume_to_max(nodes, child_list(i),child_list(i)) end if end do else write(*,*) "No initialization file provided. Initializing with method: ", init_gel_state select case (init_gel_state) !2. Initialize according to the chosen method case(HOMOGENEOUS_GEL) !2.1 Initialize on the whole planet with a Global Equivalent Layer (GEL) !$OMP PARALLEL !$OMP DO PRIVATE(Vin) SCHEDULE(STATIC) do i = 1,ndep if (nodes(child_list(i))%is_leaf .eqv. .true.) then nodes(child_list(i))%active = 1 Vin = nodes(child_list(i))%watershed_area * global_eq_lay nodes(child_list(i))%V = nodes(child_list(i))%V + Vin end if end do !$OMP END DO call OverflowAllDep() if (fixed .eqv. .true.) then call FixeWaterAll() end if !$OMP END PARALLEL case(LONLAT_GEL) !2.2. Initializing at a given lon-lat coordinate with a GEL lon_val = gel_init_lon lat_val = gel_init_lat do i = 1, lonlen-1 if (lon_val >= longitude(i) .and. lon_val <= longitude(i+1)) then exit end if end do do j = 1, latlen-1 if (lat_val <= latitude(j) .and. lat_val >= latitude(j+1)) then exit end if end do write(*,*) "water will be injected in depression No.", label(i,j), "(coordinates:", lon_val, lat_val,")" Vin = planet_area * global_eq_lay nodes(label(i,j))%active = 1 nodes(label(i,j))%V = nodes(label(i,j))%V + Vin !$OMP PARALLEL call OverflowAllDep() if (fixed .eqv. .true.) then call FixeWaterAll() end if !$OMP END PARALLEL case(FIXED_OCEAN) !2.3. Initializing with fixed ocean depressions !$OMP PARALLEL !$OMP DO SCHEDULE(STATIC) do i = 1,ndep if (nodes(child_list(i))%is_leaf .eqv. .false.) then if (nodes(child_list(i))%lake_area(1) < fixed_area .and. nodes(child_list(i))%lake_area(lake_dislen) > fixed_area) then nodes(child_list(i))%V = 1 nodes(child_list(i))%active = 1 call set_children_volume_to_max(nodes, child_list(i),child_list(i)) else if (nodes(child_list(i))%lake_area(1) > fixed_area .and. nodes(child_list(i))%left_pointer%lake_area(lake_dislen) < fixed_area) then if (nodes(child_list(i))%lake_area(1) > fixed_area .and. nodes(child_list(i))%right_pointer%lake_area(lake_dislen) < fixed_area) then nodes(child_list(i))%V = 1 nodes(child_list(i))%active = 1 call set_children_volume_to_max(nodes, child_list(i), child_list(i)) end if end if end if end do !$OMP END DO if (fixed .eqv. .true.) then !$OMP SINGLE call FixeWaterAll() write(*,*) "Fixed ocean depressions:", sum(nodes(:)%fixed) !$OMP END SINGLE end if call PrintVariableInGraph('Water volume in graph = ', 'volume') !$OMP DO PRIVATE(Vin) SCHEDULE(STATIC) do i = 1,ndep if (nodes(child_list(i))%active == 1 .and. nodes(child_list(i))%fixed == 0) then A_lac = ((nodes(child_list(i))%P * nodes(child_list(i))%watershed_area) + (nodes(child_list(i))%P * 100000000000.))/ nodes(child_list(i))%E if (A_lac > nodes(child_list(i))%lake_area(lake_dislen)) then nodes(child_list(i))%V = nodes(child_list(i))%Vmax nodes(child_list(i))%active = 1 call set_children_volume_to_max(nodes, child_list(i), child_list(i)) end if !Vin = nodes(child_list(i))%watershed_area * global_eq_lay !nodes(child_list(i))%V = nodes(child_list(i))%V + Vin end if end do !$OMP END DO call PrintVariableInGraph('Water volume in graph = ', 'volume') call OverflowAllDep() !$OMP END PARALLEL case default write(*,*) "Error: Invalid initialization method. Please choose 1, 2, or 3." stop end select !$OMP PARALLEL call PrintVariableInGraph('Water volume in graph = ', 'volume') call PrintVariableInGraph('GEL on the planet = ', 'gel') call PrintVariableInGraph('Area of the planet = ', 'planet_area') call PrintVariableInGraph('Area of active depression = ', 'area') call PrintVariableInGraph('Number of active depressions = ', 'active') call InterpolateAllDepFromV() !$OMP END PARALLEL status = save_output_active(init_file) end if call ResetUpperActiveDep() call ActiveUpperDep() end subroutine initialize_model subroutine run_simulation(init_file) implicit none character(*), intent(in) :: init_file integer :: i, it, count, count_dt, count_save, iter_idx, pos_record integer :: status real(real128) :: dt, dt_init, dt_min, dt_max, dt_factor, slope_save real(real128) :: oscillation_analysis(1000) call print_msg('> Evolving hydrology',LVL_NFO) oscillation_analysis = 0.0 V_list_prev = 2 A_list_prev = 2 V_list = 1 Vout_list = 1 Vout_list_prev = 2 convergence_indicator_V = 1 count = 0 count_dt = 0 count_save = 0 iter_idx = 1 pos_record = 0 dt_init = 1.0_real128 dt = dt_init dt_min = 1e-5 dt_max = 1.0_real128 dt_factor = 0.5_real128 slope_save = 0 do it = 1, 1 call ResetUpperActiveDep() call ActiveUpperDep() call update_rate_upper_active() count = count + 1 count_save = count_save + 1 !$OMP PARALLEL call PrintVariableInGraph('Water volume at the begining = ', 'volume') call PrintVariableInGraph('Area of active depression = ', 'area') call PrintVariableInGraph('Number of active depressions = ', 'active') !$OMP SINGLE !call PrintHeader("Iteration n°", count) V_list = 0.0_real128 Vout_list = 0.0_real128 Vin_list = 0.0_real128 active_list = 0.0_real128 Z_list = 0.0_real128 A_list = 0.0_real128 Vevap_list = 0.0_real128 Vinj_list = 0.0_real128 Qin = 0.0_real128 Qout = 0.0_real128 planet_area_active = 0.0_real128 do i = 1,ndep if (nodes(child_list(i))%active == 1) then planet_area_active = planet_area_active + nodes(child_list(i))%watershed_area !write(*,*) 'Depression No.', i, 'Area = ', nodes(child_list(i))%watershed_area end if end do print*, 'dt =', dt !$OMP END SINGLE !$OMP DO SCHEDULE(STATIC) PRIVATE(Evap_potential, Evap_real) do i = 1,ndep nodes(child_list(i))%Vevap = 0 if (nodes(child_list(i))%active == 1 .and. nodes(child_list(i))%fixed == 0) then Evap_potential = nodes(child_list(i))%A * nodes(child_list(i))%E * dt Evap_real = 0 call AvailableWater(child_list(i), Evap_potential, Evap_real) nodes(child_list(i))%Vevap = Evap_real end if end do !$OMP END DO !$OMP SINGLE alpha = sum(nodes(:)%Vevap * nodes(:)%active) / (sum(nodes(:)%P * nodes(:)%A * nodes(:)%active) + sum(nodes(:)%Runoff * (nodes(:)%watershed_area - nodes(:)%A) * nodes(:)%active)) write(*,*) 'Alpha (Evap/Precip) =', alpha !$OMP END SINGLE !$OMP DO SCHEDULE(STATIC) do i = 1,ndep nodes(child_list(i))%dV = 0 if (nodes(child_list(i))%active== 1 .and. nodes(child_list(i))%fixed == 0) then nodes(child_list(i))%dV = (nodes(child_list(i))%A * nodes(child_list(i))%P * alpha * dt) + (nodes(child_list(i))%Runoff * (nodes(child_list(i))%watershed_area - nodes(child_list(i))%A) * alpha * dt) - nodes(child_list(i))%Vevap end if end do !$OMP END DO !$OMP SINGLE area = 0.0_real128 Evap_real_store = 0.0_real128 !$OMP END SINGLE !$OMP DO SCHEDULE(STATIC) REDUCTION(+:Evap_real_store, area) PRIVATE(Evap_real, Evap_potential) do i = 1,ndep if (nodes(child_list(i))%active == 1 .and. nodes(child_list(i))%fixed == 0) then if (nodes(child_list(i))%dV < 0) then Evap_potential = abs(nodes(child_list(i))%dV) Evap_real = 0 call EvapLake(child_list(i), Evap_potential, Evap_real) Evap_real_store = Evap_real_store + Evap_real Vevap_list(i) = Evap_real area = area + nodes(child_list(i))%watershed_area end if end if end do !$OMP END DO !$OMP SINGLE V_list = 0.0_real128 write(*,*) 'Water volume in graph after evap =', sum(nodes(:)%V) write(*,*) 'Water volume in atm after evap =', Evap_real_store write(*,*) 'Water volume =', Evap_real_store + sum(nodes(:)%V) !$OMP END SINGLE !$OMP DO do i = 1,ndep V_list(i) = nodes(child_list(i))%V end do !$OMP END DO call PrintVariableInGraph('Area of active depression =', 'area') !$OMP SINGLE !write(*,*) 'Precip Rate =', sum(gcm_grid(:,:)%P*gcm_grid(:,:)%cell_area)/sum(gcm_grid(:,:)%cell_area) * dt, 'm' write(*,*) 'Evap_real_store =', Evap_real_store, 'm3' ! Inject Water Vin_store = 0.0_real128 planet_area_active = 0 !$OMP END SINGLE !$OMP DO SCHEDULE(GUIDED) PRIVATE(Vin) REDUCTION(+:Vin_store, planet_area_active) do i = 1,ndep if (nodes(child_list(i))%active == 1 .and. nodes(child_list(i))%fixed == 0) then if (nodes(child_list(i))%dV > 0) then Vin = nodes(child_list(i))%dV nodes(child_list(i))%V = nodes(child_list(i))%V + Vin Vin_store = Vin_store + Vin Vinj_list(i) = Vin planet_area_active = planet_area_active + nodes(child_list(i))%watershed_area end if end if end do !$OMP END DO ! Move overflow into graph depressions ! intialize bypass !$OMP DO do i = 1,ndep nodes(child_list(i))%bypass = nodes(child_list(i))%child end do !$OMP END DO call OverflowAllDep() ! Calculate elevation Z and area A for each lake call InterpolateAllDepFromV() !$OMP END PARALLEL call ResetUpperActiveDep() call ActiveUpperDep() call PrintVariableInGraph('GEL in the planet = ', 'gel') call PrintVariableInGraph('Water volume at the end = ', 'volume') end do status = save_output_active(init_file) end subroutine run_simulation subroutine update_rate_upper_active() implicit none integer :: i do i = 1,ndep if (nodes(child_list(i))%is_leaf) then nodes(nodes(child_list(i))%upper_active)%P = nodes(nodes(child_list(i))%upper_active)%P + (nodes(child_list(i))%watershed_area * nodes(child_list(i))%P / nodes(nodes(child_list(i))%upper_active)%watershed_area) nodes(nodes(child_list(i))%upper_active)%E = nodes(nodes(child_list(i))%upper_active)%E + (nodes(child_list(i))%watershed_area * nodes(child_list(i))%E / nodes(nodes(child_list(i))%upper_active)%watershed_area) nodes(nodes(child_list(i))%upper_active)%Runoff = nodes(nodes(child_list(i))%upper_active)%Runoff + (nodes(child_list(i))%watershed_area * nodes(child_list(i))%Runoff / nodes(nodes(child_list(i))%upper_active)%watershed_area) end if end do end subroutine update_rate_upper_active recursive subroutine set_children_volume_to_max(tree, child, active_child) ! Use when the model is initialized from a file. The subroutine fills volume below the active depression use :: hydro_tree type(TreeNode), target, intent(inout) :: tree(IDdep) integer, intent(in) :: child, active_child if (tree(child)%is_leaf) then return end if if (associated(tree(child)%left_pointer)) then tree(child)%left_pointer%V = tree(child)%left_pointer%Vmax tree(child)%left_pointer%active = 0 tree(child)%left_pointer%upper_active = active_child call set_children_volume_to_max(tree, tree(child)%left_pointer%child, active_child) !write(*,*) "set_left_children_volume_to_max: child" end if if (associated(tree(child)%right_pointer)) then tree(child)%right_pointer%V = tree(child)%right_pointer%Vmax tree(child)%right_pointer%active = 0 tree(child)%right_pointer%upper_active = active_child call set_children_volume_to_max(tree, tree(child)%right_pointer%child, active_child) !write(*,*) "set_right_children_volume_to_max: child" end if end subroutine set_children_volume_to_max subroutine InterpolateAllDepFromV() implicit none !Local variables integer :: i !$OMP SINGLE Z_list = 0.0_real128 V_list = 0.0_real128 A_list = 0.0_real128 Vout_list = 0.0_real128 Vin_list = 0.0_real128 active_list = 0 !$OMP END SINGLE !$OMP DO SCHEDULE(GUIDED) PRIVATE(Z_interp, A_interp) do i = 1,ndep if (nodes(child_list(i))%Vmax /= V_trans_dep .and. nodes(child_list(i))%active == 1.) then call InterpolateFromV(lake_dislen, nodes(child_list(i))%lake_elev, nodes(child_list(i))%lake_vol, nodes(child_list(i))%lake_area, nodes(child_list(i))%V, A_interp, Z_interp) nodes(child_list(i))%Z = Z_interp nodes(child_list(i))%A = A_interp end if Z_list(i) = nodes(child_list(i))%Z V_list(i) = nodes(child_list(i))%V A_list(i) = nodes(child_list(i))%A Vout_list(i) = nodes(child_list(i))%Vout nodes(child_list(i))%Vout = 0 Vin_list(i) = nodes(child_list(i))%Vin nodes(child_list(i))%Vin = 0 active_list(i) = nodes(child_list(i))%active end do !$OMP END DO end subroutine InterpolateAllDepFromV subroutine PrintVariableInGraph(message, varname) implicit none character(len=*), intent(in) :: message, varname ! Message to print !$OMP SINGLE select case (trim(varname)) case ("volume") write(*,*) trim(message), sum(nodes(:)%V), trim('m³') case ("area") write(*,*) trim(message), sum(nodes(:)%active * nodes(:)%watershed_area), trim('m²') case ("active") write(*,*) trim(message), sum(nodes(:)%active), trim('depressions') case ("planet_area") write(*,*) trim(message), planet_area, trim('m²') case ("gel") write(*,*) trim(message), sum(nodes(:)%V)/planet_area, trim('mGEL') end select !$OMP END SINGLE end subroutine PrintVariableInGraph function detect_oscillation(x, n) result(is_oscillating) implicit none integer, intent(in) :: n real(real128), intent(in) :: x(n) logical :: is_oscillating integer :: i, sign_changes real(real128) :: d_prev, d_curr sign_changes = 0 d_prev = x(2) - x(1) do i = 2, n-1 d_curr = x(i+1) - x(i) if (d_curr * d_prev < 0.0d0) then sign_changes = sign_changes + 1 end if d_prev = d_curr end do ! Critère : plus de 20% de changements de signe → oscillation is_oscillating = (sign_changes > 0.2d0 * real(n)) write(*,*) 'Oscillation detection: ', is_oscillating, ' (Sign changes: ', sign_changes, ')', 0.2d0 * real(n) end function detect_oscillation recursive subroutine get_volume_below_node(child, total_volume) ! Compute total volume below a given node integer, intent(in) :: child real(real128), intent(inout) :: total_volume total_volume = total_volume + nodes(child)%V if (nodes(child)%is_leaf) then return end if if (associated(nodes(child)%left_pointer)) then call get_volume_below_node(nodes(child)%left_pointer%child, total_volume) end if if (associated(nodes(child)%right_pointer)) then call get_volume_below_node(nodes(child)%right_pointer%child, total_volume) end if end subroutine get_volume_below_node recursive subroutine compute_volume_above_depression(child, volume, lake_area) ! Compute total volume above depression integer, intent(in) :: child real(real128), intent(in) :: volume real(real128), intent(in) :: lake_area ! local variables real(real128) :: area ! initialize variables area = 0.0_real128 if (nodes(child)%active == 1) then nodes(child)%Vtop = nodes(child)%V else if (nodes(child)%Vmax == V_trans_dep) then nodes(child)%Vtop = nodes(child)%Vmax + (volume * nodes(child)%A / lake_area) else nodes(child)%Vtop = nodes(child)%Vmax + (volume * nodes(child)%lake_area(lake_dislen) / lake_area) end if end if if (nodes(child)%is_leaf) then return !else !if (nodes(child)%Vmax /= V_trans_dep) then !write(*,*) '----------------------------------------------' !write(*,*) 'right :',nodes(child)%left_pointer%lake_area(lake_dislen) !write(*,*) 'left :',nodes(child)%right_pointer%lake_area(lake_dislen) !write(*,*) 'current:',nodes(child)%lake_area(1) !write(*,*) 'sum :',nodes(child)%left_pointer%lake_area(lake_dislen) + nodes(child)%right_pointer%lake_area(lake_dislen) ! write(*,*) 'right :',nodes(nodes(child)%rchild)%lake_elev(lake_dislen) !write(*,*) 'left :',nodes(nodes(child)%lchild)%lake_elev(lake_dislen) !write(*,*) 'current:',nodes(child)%lake_elev(1) ! end if end if if (associated(nodes(child)%left_pointer)) then if (nodes(child)%left_pointer%Vmax == V_trans_dep) then area = area + nodes(child)%left_pointer%A else area = area + nodes(child)%left_pointer%lake_area(lake_dislen) end if end if if (associated(nodes(child)%right_pointer)) then if (nodes(child)%right_pointer%Vmax == V_trans_dep) then area = area + nodes(child)%right_pointer%A else area = area + nodes(child)%right_pointer%lake_area(lake_dislen) end if end if if (associated(nodes(child)%left_pointer)) then call compute_volume_above_depression(nodes(child)%left_pointer%child, nodes(child)%Vtop, area) end if if (associated(nodes(child)%right_pointer)) then call compute_volume_above_depression(nodes(child)%right_pointer%child, nodes(child)%Vtop, area) end if end subroutine compute_volume_above_depression subroutine FillDepression(id, nAffected, dh) implicit none integer, intent(in) :: id logical, allocatable :: dep_mask(:,:) integer, allocatable :: dep_elev(:,:), dep_label(:,:) real(real128), allocatable :: dep_area(:) real :: min_val integer :: big ! A large number representing "infinity" real(real128) :: WaterVol ! real(8) :: wtd(nCells) !!------------------------------------------------------------------ !! Local variables !!------------------------------------------------------------------ ! logical :: visited(nCells) ! cells already inserted into the queue ! logical :: affected(nCells) ! cells already flooded by the lake logical, allocatable :: visited(:,:), affected(:,:) integer, allocatable :: PQ(:,:) ! priority queue (min-heap of cell IDs) integer :: PQsize ! number of elements in the queue !integer :: c, n, k ! loop indices and current cell integer :: nx, ny, i, j, i_cell, j_cell, ii ,jj real(real128), intent(out) :: nAffected ! number of flooded cells real(real128) :: Te ! current lake water level (elevation) real(real128) :: Vreq ! volume required to raise water level real(real128), intent(out) :: dh ! final partial water level increment integer :: nleafs integer, allocatable :: leaf_ids(:) logical, allocatable :: is_leaf(:) integer :: leaf_id_max !------------------------------------------------------------------ ! Initialization !------------------------------------------------------------------ ! Initial lake level is the pit elevation dh = 0 Te = 0 !dep_elev(i_cell, j_cell) * dep_area(j_cell) nAffected = 0 !dep_area(j_cell) ! visited = .false. ! affected = .false. ! wtd = 0.0d0 big = huge(dep_elev) ! Compute water volume in the depression if (nodes(id)%is_leaf) then WaterVol = nodes(id)%V else WaterVol = 0 call get_volume_below_node(nodes(id)%child, WaterVol) end if if (WaterVol <= 0.0d0) return ! No water → nothing to do call load_ID_map(nodes(id)%imin, nodes(id)%imax, nodes(id)%jmin, nodes(id)%jmax, dep_label, dep_elev, nx, ny) ! hydro_inout.F90 - Database Module allocate(dep_mask(nx,ny),stat=status) if (nodes(id)%is_leaf) then dep_mask = (dep_label == nodes(id)%child) else nleafs = 0 call get_leaves_number_below_node(nodes(id)%child, nleafs) allocate(leaf_ids(nleafs)) nleafs = 0 call get_leaves_below_node(nodes(id)%child,leaf_ids, nleafs) dep_mask = .false. leaf_id_max = maxval(dep_label) allocate(is_leaf(leaf_id_max)) is_leaf = .false. is_leaf(leaf_ids) = .true. do j = 1, size(dep_label,2) do i = 1, size(dep_label,1) dep_mask(i,j) = is_leaf(dep_label(i,j)) end do end do deallocate(leaf_ids) end if deallocate(dep_label) ! Clip areas cells allocate(dep_area(ny)) dep_area = cell_area(nodes(id)%jmin:nodes(id)%jmax) ! Find pit cell min_val = minval(merge(dep_elev, big, dep_mask)) outer: do i = 1, nx do j = 1, ny if (dep_mask(i,j) .and. dep_elev(i,j) == min_val) then i_cell = i j_cell = j exit outer end if end do end do outer !Allocatate vairables allocate(PQ(2, count(dep_mask))) ! 1st row: is i index, 2nd row: j index allocate(visited(nx,ny)) allocate(affected(nx,ny)) !Initialize priority queue with the pit (lowest cell) PQsize = 0 PQ = 0 visited = .false. affected = .false. call PQ_push(PQ, PQsize, i_cell, j_cell, dep_elev) visited(i_cell, j_cell) = .true. !do i = 1, size(dep_mask,2) ! write(*,'(*(L1))') affected(:,i) !end do !!------------------------------------------------------------------ !! Main filling loop !! The lake level is raised step by step following topography !!------------------------------------------------------------------ do while (PQsize > 0) ! Get the next lowest cell surrounding the flooded area call PQ_pop(PQ, PQsize, i_cell,j_cell, dep_elev) ! Volume needed to raise the lake level from Te up to elev(c) Vreq = nAffected * dep_elev(i_cell, j_cell) - Te !write(*,*) Vreq, Watervol !!-------------------------------------------------------------- !! Case 1: not enough water to reach the next elevation step !!-------------------------------------------------------------- if (WaterVol < Vreq) then !! Distribute remaining water uniformly over flooded cells dh = (WaterVol + Te) / nAffected !do n = 1, nCells ! if (affected(n)) wtd(n) = wtd(n) + dh !end do !write(*,*) 'Area =', nAffected, dh return ! final lake level reached (closed basin) end if !!-------------------------------------------------------------- !! Case 2: enough water to reach elev(c) !!-------------------------------------------------------------- !! If this is not the spill point, the cell becomes flooded if (.not. affected(i_cell, j_cell)) then affected(i_cell, j_cell) = .true. nAffected = nAffected + dep_area(j_cell) Te = Te + (dep_elev(i_cell, j_cell) * dep_area(j_cell)) end if !!-------------------------------------------------------------- !! Expand laterally: add neighbors belonging to the depression !!-------------------------------------------------------------- do ii = max(1, i_cell-1), min(nx, i_cell+1) do jj = max(1, j_cell-1), min(ny, j_cell+1) if (ii == i_cell .and. jj == j_cell) cycle if (dep_mask(ii, jj) .and. .not. visited(ii, jj)) then call PQ_push(PQ, PQsize, ii, jj, dep_elev) visited(ii, jj) = .true. end if end do end do !!-------------------------------------------------------------- !! Safety: ensure the spill point is eventually processed !!-------------------------------------------------------------- !if (PQsize == 0 .and. .not. visited(OutCell)) then ! call PQ_push(PQ, PQsize, OutCell, elev) ! visited(OutCell) = .true. !end if end do end subroutine FillDepression subroutine PQ_push(PQ, size, i, j, elev) implicit none ! Insert a cell into the priority queue ordered by elevation integer, intent(inout) :: PQ(:,:) integer, intent(inout) :: size integer, intent(in) :: i, j integer, intent(in) :: elev(:,:) integer :: k, parent integer :: tmp(2) size = size + 1 k = size PQ(1,k) = i PQ(2,k) = j ! Bubble up while heap property is violated do while (k > 1) parent = k / 2 if (elev(PQ(1,parent),PQ(2,parent)) <= elev(PQ(1,k),PQ(2,k))) exit tmp = PQ(:,parent) PQ(:,parent) = PQ(:,k) PQ(:,k) = tmp k = parent end do end subroutine PQ_push subroutine PQ_pop(PQ, size, i, j, elev) implicit none ! Remove and return the cell(i,j) with the lowest elevation integer, intent(inout) :: PQ(:,:) integer, intent(inout) :: size integer, intent(out) :: i, j integer, intent(in) :: elev(:,:) integer :: k, left, right, smallest integer :: tmp(2) i = PQ(1,1) j = PQ(2,1) PQ(:,1) = PQ(:,size) size = size - 1 k = 1 ! Restore heap property by pushing down do left = 2*k right = 2*k + 1 smallest = k if (left <= size) then if (elev(PQ(1,left),PQ(2,left)) < elev(PQ(1,smallest),PQ(2,smallest))) smallest = left end if if (right <= size) then if (elev(PQ(1,right),PQ(2,right)) < elev(PQ(1,smallest),PQ(2,smallest))) smallest = right end if if (smallest == k) exit tmp = PQ(:,k) PQ(:,k) = PQ(:,smallest) PQ(:,smallest) = tmp k = smallest end do end subroutine PQ_pop end module hydro_core