Solving land points in ModelGrid with no Hydrogrid land points

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)

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
     !

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, \