From 80eb806263ec871f7bbaf24db75108db8fe0e4f0 Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Fri, 29 May 2020 21:25:14 +0200
Subject: [PATCH] Solving land points in ModelGrid with no Hydrogrid land
points
---
F90subroutines/routing_interface.f90 | 5 +++--
F90subroutines/routing_reg.f90 | 24 +++++++++++++++++++++++-
Interface.py | 5 +++--
3 files changed, 29 insertions(+), 5 deletions(-)
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index e741f22..e59e47a 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -425,7 +425,7 @@ SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_
END SUBROUTINE rivclassification
-SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, &
+SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridarea, cfrac, gridcenters, basin_count, &
& basin_notrun, basin_area, basin_orog, basin_floodp, basin_cg, basin_topoind, fetch_basin, &
& basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin, routing_area, routing_orog, routing_floodp, &
@@ -446,6 +446,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
INTEGER(i_std), INTENT (in) :: num_largest
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: cfrac
+ REAL(r_std), DIMENSION(nbpt,2), INTENT(in) :: gridcenters
!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_notrun !!
@@ -488,7 +489,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
!
!
- CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, inflowmax, num_largest, &
+ CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, gridcenters, nwbas, inflowmax, num_largest, &
& basin_count, basin_notrun, basin_area, basin_orog, basin_floodp, basin_cg, &
& basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 832fbfe..13bb650 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2603,7 +2603,7 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, inflowmax, num_largest, &
+SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcenters, nwbas, inflowmax, num_largest, &
& basin_count, basin_notrun, basin_area, basin_orog, basin_floodp, &
& basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
@@ -2620,6 +2620,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, in
!
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea !! The size of each grid box in X and Y (m)
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box (unitless;0-1)
+ REAL(r_std), DIMENSION(nbpt, 2), INTENT(in) :: gridcenters
!
INTEGER(i_std), INTENT(in) :: nwbas, inflowmax !!
INTEGER(i_std), INTENT (in) :: num_largest
@@ -2689,6 +2690,27 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, in
route_count_glo(ib) = basin_count(ib)
route_nbintobas_glo(ib) = basin_count(ib)
origin_nbintobas_glo(ib) = basin_notrun(ib)
+ ! If the land point doesn't overlap with any hydrogrid point
+ ! (cf 0.5° trip on coast or some particular grid points )
+ IF ( basin_count(ib) .EQ. 0) THEN
+ basin_area(ib,1) = contfrac(ib)*gridarea(ib)
+ basin_orog(ib,1) = 0
+ basin_floodp(ib,1) = 0
+ basin_cg(ib,1,1) = gridcenters(ib,1)
+ basin_cg(ib,1,2) = gridcenters(ib,2)
+ ! To define in function of the minimum topoindex on the map
+ ! (to revise)
+ basin_topoind(ib,1) = 1
+ basin_id(ib,1) = 0
+ basin_coor(ib,1,1) = basin_cg(ib,1,1)
+ basin_coor(ib,1,2) = basin_cg(ib,1,2)
+ fetch_basin(ib,1) = basin_area(ib,1)
+ basin_type(ib,1) = -4
+ outflow_grid(ib,1) = -2
+ basin_count(ib) = 1
+ route_count_glo(ib) = 1
+ route_nbintobas_glo(ib) = 1
+ END IF
!
ENDDO
!
diff --git a/Interface.py b/Interface.py
index 3e1daaa..d954856 100644
--- a/Interface.py
+++ b/Interface.py
@@ -328,7 +328,7 @@ class HydroSuper :
#
self.nwbas = part.domainmax(np.max(nb_basin))
# Set the number of inflows per basin. For the moment twice the maximum number of basins.
- self.inflowmax = self.nwbas*2
+ self.inflowmax = max(10, self.nwbas*2)
print("Maximum number of basin created : {0}".format(self.nwbas))
ijdim=[]
for i in range(self.nbpt) :
@@ -590,7 +590,8 @@ class HydroGraph :
self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
routing_interface.finish_truncate(nbpt = self.nbpt, inflowmax = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, \
- cfrac = modelgrid.contfrac, basin_count = hydrosuper.basin_count, \
+ cfrac = modelgrid.contfrac, gridcenters = np.array(modelgrid.centers), \
+ basin_count = hydrosuper.basin_count, \
basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
basin_orog = hydrosuper.basin_orog, basin_floodp = hydrosuper.basin_floodp, \
basin_cg = hydrosuper.basin_cg, \
--
GitLab