Changes relative to floodplains height

F90subroutines/routing_interface.f90

@@ -435,10 +435,10 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
      & basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, &
      & 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_mean, &
+     & inflow_number, inflow_grid, inflow_basin, floodcri, routing_area, routing_orog_mean, &
      & routing_orog_min, routing_orog_max, routing_floodp, &
      & routing_cg, topo_resid, route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
-     & global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
+     & global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch,routing_floodcri)
   !
   USE ioipsl
   USE grid
@@ -478,6 +478,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout)          :: inflow_number !!
   INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_basin !!
   INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid !!
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: floodcri !!
   !
   ! Output variables
   !
@@ -497,6 +498,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
   INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: global_basinid !! ID of basin (unitless)
   REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out)  :: route_outlet !! Coordinate of outlet (-)
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: route_type !! Coordinate of outlet (-)
+  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_floodcri
   !
   !
   !
@@ -504,7 +506,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
   CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, gridcenters, nwbas, inflowmax, num_largest, &
        & basin_count, basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, 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)
+       & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
 
   routing_area(:,:) = routing_area_glo(:,:)
   routing_orog_mean(:,:) = routing_orog_mean_glo(:,:)
@@ -516,7 +518,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
   route_nbbasin(:) = route_count_glo(:)
   route_togrid(:,:) = route_togrid_glo(:,:)
   route_tobasin(:,:) = route_tobasin_glo(:,:)
-
+  routing_floodcri(:,:) = routing_floodcri_glo(:,:)
   routing_fetch(:,:) = route_fetch_glo(:,:)
 
   route_nbintobas(:) = route_nbintobas_glo(:)
@@ -769,7 +771,6 @@ SUBROUTINE correct_outflows(nbpt, nwbas,nbhalo, outflow_grid, outflow_basin, &
 
 END SUBROUTINE correct_outflows
 
-
 SUBROUTINE correct_inflows(nbpt, nwbas, inflowmax, outflow_grid,&
             & outflow_basin, basin_count, inflow_number, inflow_grid, &
             & inflow_basin)
@@ -819,6 +820,57 @@ SUBROUTINE correct_inflows(nbpt, nwbas, inflowmax, outflow_grid,&
 
 END SUBROUTINE correct_inflows
 
+
+SUBROUTINE get_floodcri(nbpt, nwbas,inflowmax, inflow_number, inflow_basin, inflow_grid, &
+                 & basin_count, basin_floodp, basin_orog_min,floodcri)
+  !
+  USE ioipsl
+  USE grid
+  USE routing_tools
+  USE routing_reg
+  !
+  !! INPUT VARIABLES
+  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
+  INTEGER(i_std), INTENT (in) :: nwbas !!
+  INTEGER(i_std), INTENT (in) :: inflowmax
+
+  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in)           :: inflow_number
+  INTEGER(i_std), DIMENSION(nbpt), INTENT(in)           :: basin_count
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in)              :: basin_floodp, basin_orog_min
+  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(in) :: inflow_basin
+  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(in) :: inflow_grid
+
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out)              :: floodcri
+
+  INTEGER :: ig, ib, ug, ub, i_inf, nbas
+  REAL    :: diff, orog, d
+
+  floodcri(:,:) = 0
+
+  DO ig=1,nbpt
+    nbas = basin_count(ig)
+    DO ib=1,nbas
+      IF (basin_floodp(ig,ib) .GT. 0) THEN
+        IF (inflow_number(ig,ib) .EQ. 0) THEN
+          diff = 2
+        ELSE
+          diff = 999
+          orog = basin_orog_min(ig,ib)
+          DO i_inf=1,inflow_number(ig,ib)
+            ug = inflow_grid(ig,ib,i_inf)
+            ub = inflow_basin(ig,ib,i_inf)
+            d = basin_orog_min(ug,ub) - orog
+            IF (d<diff) diff = d
+          END DO
+          IF (diff < 2) diff = 2
+          floodcri(ig, ib) = diff
+        END IF      
+      END IF
+    END DO
+  END DO
+
+END SUBROUTINE get_floodcri
+
 SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, outflow_grid, outflow_basin, basin_count)
   !
   USE ioipsl
@@ -1002,7 +1054,6 @@ SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, nbhalo, in_core, halopts,  corepts,
   !
 END SUBROUTINE finalfetch
 
-
 SUBROUTINE finalrivclass(nbpt, nbasmax, nbcore, corepts, route_togrid, route_tobasin, routing_fetch, largest_rivarea, num_largest)
   !
   USE routing_reg

F90subroutines/routing_reg.f90

@@ -29,6 +29,7 @@ MODULE routing_reg
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:)  :: route_outlet_glo !! Coordinate of outlet (-)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: route_type_glo !! Coordinate of outlet (-)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: route_fetch_glo !! Upstream area at each HTU
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: routing_floodcri_glo !! Upstream area at each HTU
 
 CONTAINS
 !! ================================================================================================================================
@@ -2587,7 +2588,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
        & inflowmax, num_largest, basin_count, basin_notrun, basin_area, basin_orog_mean,&
        & basin_orog_min, basin_orog_max, 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)
+       & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
     !
     IMPLICIT NONE
     !
@@ -2627,6 +2628,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
     INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout)          :: inflow_number !!
     INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_basin !!
     INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid !!
+    REAL(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: floodcri !!
     !
 !! LOCAL VARIABLES
     INTEGER(i_std) :: ib, ij, ic, ibf, ijf, igrif, ibasf, cnt, pold, bold, ff(2) !! Indices (unitless)
@@ -2677,6 +2679,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
        routing_orog_mean_glo(ib,:) = 0
        routing_orog_min_glo(ib,:) = 0
        routing_orog_max_glo(ib,:) = 0
+       routing_floodcri_glo(ib,:) = 0
        ! 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
@@ -2726,6 +2729,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
           !
           route_togrid_glo(ib,ij) = outflow_grid(ib,ij)
           route_tobasin_glo(ib,ij) = outflow_basin(ib,ij)
+          routing_floodcri_glo(ib,ij) = floodcri(ib,ij)*2000
           !
           !
        ENDDO
@@ -3060,6 +3064,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
     ALLOCATE (route_outlet_glo(nbpt,nbasmax,2), stat=ier)
     ALLOCATE (route_type_glo(nbpt,nbasmax), stat=ier)
     ALLOCATE (origin_nbintobas_glo(nbpt), stat=ier)
+    ALLOCATE (routing_floodcri_glo(nbpt,nbasmax), stat=ier)
 
   END SUBROUTINE allocategraph
 END MODULE routing_reg

GraphHydro.py

@@ -64,7 +64,7 @@ class HydroGraph :
         self.routing_area, self.routing_orog_mean, self.routing_orog_min,self.routing_orog_max, \
             self.routing_floodp, self.routing_cg, self.topo_resid, self.route_nbbasin,\
             self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
-            self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
+            self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch, self.floodcri = \
                                     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, gridcenters = np.array(modelgrid.centers), \
@@ -79,7 +79,7 @@ class HydroGraph :
                                                                       basin_flowdir = hydrosuper.basin_flowdir, \
                                                                       outflow_grid = hydrosuper.outflow_grid, outflow_basin = hydrosuper.outflow_basin, \
                                                                       inflow_number = hydrosuper.inflow_number, inflow_grid = hydrosuper.inflow_grid, \
-                                                                      inflow_basin = hydrosuper.inflow_basin)
+                                                                      inflow_basin = hydrosuper.inflow_basin, floodcri = hydrosuper.floodcri)
 
         #
         #
@@ -197,6 +197,8 @@ class HydroGraph :
 
         self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_floodp", "Fraction of floodplains", "-", self.routing_floodp[:,:], vtyp)
 
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "floodcri", "Height at which all basin is flooded", "mm", self.floodcri[:,:], vtyp)
+
         # route number into basin
         self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, orig_type="int")
         #