Initialization of beta variable calculation (floodplains) + parameter maxpercent

F90subroutines/routing_interface.f90

@@ -43,9 +43,9 @@ SUBROUTINE closeinterface()
   !
 END SUBROUTINE closeinterface
 
-SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, sub_pts, sub_index, sub_area,  max_basins, min_topoind, &
-     & lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy,orog,floodp, nbi, nbj, area_bx, trip_bx, basin_bx, &
-     & topoind_bx, fac_bx, hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
+SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_index, sub_area,  max_basins, &
+     & min_topoind, lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy,orog,floodp, nbi, nbj, area_bx, &
+     & trip_bx, basin_bx, topoind_bx, fac_bx, hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
   !
   USE ioipsl
   USE grid
@@ -54,7 +54,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, sub_pts, sub_index, sub_area,
   !
   IMPLICIT NONE
   !
-  INTEGER, INTENT(in) :: nbvmax_in, ijdimmax
+  INTEGER, INTENT(in) :: nbvmax_in, ijdimmax, maxpercent_in
   !
   INTEGER, INTENT(in) :: nbpt !! Domain size (unitless)
   INTEGER, INTENT(in) :: sub_pts(nbpt) !! Number of high resolution points on this grid (unitless)
@@ -96,6 +96,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, sub_pts, sub_index, sub_area,
   ! It is transfered here but should be argument to the various subroutines.
   !
   nbvmax = nbvmax_in
+  maxpercent = maxpercent_in
   !
   IF ( debug ) WRITE(numout,*) "Memory Mgt getgrid : nbvmax, ijdimmax = ", nbvmax, ijdimmax
   !
@@ -179,7 +180,7 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
           & basin_bbout, lshead, coast_pts, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, &
           & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, &
           & basin_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
-          & outflow_grid, outflow_basin, nbcoastal, coastal_basin)
+          & basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
   !
   USE ioipsl
   USE grid
@@ -228,6 +229,7 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
   REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_fac  !!
   REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_topoind    !! Topographic index of the residence time for a basin (m)
   REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_lshead     !! Large scale heading out of the grid box.
+  REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas)    :: basin_beta_fp  !! Beta parameter for floodplains
   INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: outflow_grid  !! Type of outflow on the grid box (unitless)
   INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: outflow_basin !!
   INTEGER(i_std), INTENT (out), DIMENSION(nbpt)       :: nbcoastal           !!
@@ -244,7 +246,7 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
           & basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy,&
           & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac,&
           & basin_topoind, basin_id, basin_coor, &
-          & basin_type, basin_flowdir, basin_lshead, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
+          & basin_type, basin_flowdir, basin_lshead, basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
   ENDDO
   !
 END SUBROUTINE globalize
@@ -433,10 +435,10 @@ END SUBROUTINE rivclassification
 
 SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridarea, cfrac, gridcenters, basin_count, &
      & basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, &
-     & basin_floodp, basin_cg, basin_topoind, fetch_basin, &
+     & basin_floodp, basin_cg, basin_topoind, basin_beta_fp, fetch_basin, &
      & basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
      & inflow_number, inflow_grid, inflow_basin, floodcri, routing_area, routing_orog_mean, &
-     & routing_orog_min, routing_orog_max, routing_floodp, &
+     & routing_orog_min, routing_orog_max, routing_floodp, routing_beta,&
      & routing_cg, topo_resid, route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
      & global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch,routing_floodcri)
   !
@@ -469,6 +471,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_floodp !!
   REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout)  :: basin_cg
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_topoind !! Topographic index of the residence time for a basin (m)
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_beta_fp
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: fetch_basin !!
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
@@ -487,6 +490,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_orog_min !! Mean Orography (m)
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_orog_max !! Mean Orography (m)
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_floodp !! Surface of Floodplains (m^2)
+  REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_beta
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: routing_fetch !! Upstream are flowing into HTU (m^2)
   REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out)  :: routing_cg !! Centre of gravity of HTU (Lat, Lon)
   REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out)    :: topo_resid !! Topographic index of the retention time (m)
@@ -505,7 +509,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, &
+       & basin_cg, basin_topoind, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
        & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
 
   routing_area(:,:) = routing_area_glo(:,:)
@@ -520,6 +524,7 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
   route_tobasin(:,:) = route_tobasin_glo(:,:)
   routing_floodcri(:,:) = routing_floodcri_glo(:,:)
   routing_fetch(:,:) = route_fetch_glo(:,:)
+  routing_beta(:,:) = routing_beta_glo(:,:)
 
   route_nbintobas(:) = route_nbintobas_glo(:)
   global_basinid(:,:) = global_basinid_glo(:,:)
@@ -580,8 +585,8 @@ END SUBROUTINE finish_inflows
 
 SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, numops, basin_count, 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)
+     & basin_cg, basin_topoind, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir,&
+     & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
   !
   !
   USE ioipsl
@@ -610,6 +615,7 @@ SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, num
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_floodp !!
   REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout)  :: basin_cg
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_topoind !! Topographic index of the residence time for a basin (m)
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_beta_fp
   REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: fetch_basin !!
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
@@ -633,7 +639,14 @@ SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, num
          totak = totakeover(ib,op)
          IF (tok .GT. 0) THEN
             IF ( debug ) THEN
-               WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak
+               WRITE(numout,*) " "
+               WRITE(numout,*) "numop", op
+               WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak, "basin_count", basin_count(ib)
+               WRITE(numout,*) "outflow tokill", outflow_grid(ib,tok), outflow_basin(ib,tok)
+               WRITE(numout,*) "inflows totakeover", inflow_number(ib,totak)
+               DO it=1,inflow_number(ib,totak)
+                 WRITE(numout,*) "->", inflow_grid(ib,totak,it), inflow_basin(ib,totak,it)
+               END DO
             ENDIF
             ! Test if tokill is downstream of totakeover (avoid loop)
             igrif = outflow_grid(ib,totak)
@@ -652,7 +665,7 @@ SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, num
               END IF
             END DO
             CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, inflowmax, basin_count, basin_area, &
-               & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind,&
+               & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind, basin_beta_fp,&
                & fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
                & inflow_grid, inflow_basin)
          END IF

F90subroutines/routing_reg.f90

@@ -18,6 +18,8 @@ MODULE routing_reg
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: routing_orog_min_glo !! Mean topography (m)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: routing_orog_max_glo !! Mean topography (m)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: routing_floodp_glo !! Surface of floodplains (m^2)
+  REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: routing_beta_glo !!
+  !
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:)  :: routing_cg_glo !! Centre of gravity of HTU (Lat, Lon)
   REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:)    :: topo_resid_glo !! Topographic index of the retention time (m)
   INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:)   :: route_count_glo !! Number of basins per grid
@@ -1528,7 +1530,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        & basin_bbout, lshead, coast_pts, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, &
        & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, &
        & basin_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
-       & outflow_grid, outflow_basin, nbcoastal, coastal_basin)
+       & basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
     !
     IMPLICIT NONE
     !
@@ -1576,6 +1578,8 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
     REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_fac  !!
     REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind    !! Topographic index of the residence time for a basin (m)
     REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_lshead     !! Large scale heading out of the grid box.
+    REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_std  !!
+    REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_beta_fp  !! 
     INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid  !! Type of outflow on the grid box (unitless)
     INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !!
     INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal           !!
@@ -1588,6 +1592,15 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
     !
     INTEGER(i_std) :: option !! Option to calculate topoindex
     REAL(r_std) :: area_frac !! Area fraction for calculate topoindex
+    REAL(r_std) :: temp !! local variable
+    !
+    ! PARAMETERs beta calculation
+    ! limits for orog_std normalization
+    REAL(r_std) :: uplim_std = 20 !! 
+    REAL(r_std) :: lowlim_std = 0.1 !! 
+    ! limits for beta normalization
+    REAL(r_std) :: uplim_beta = 2 !! 
+    REAL(r_std) :: lowlim_beta = 0.2 !! 
     !
     LOGICAL, PARAMETER :: debug = .FALSE. !! Add new logical variable to debug
     !
@@ -1648,6 +1661,8 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        basin_orog_min(ib,ij) = 99999
        basin_orog_max(ib,ij) = zero
        basin_floodp(ib,ij) = zero
+       basin_orog_std(ib,ij) = 0
+       basin_beta_fp(ib,ij) = 0
        !
        SELECT CASE (hierar_method)
           !
@@ -1720,6 +1735,22 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        ! Finish the calculation of mean orography
        basin_orog_mean(ib,ij) = basin_orog_mean(ib,ij)/basin_area(ib,ij)
        !
+       ! Standard deviation of orography
+       IF (basin_floodp(ib,ij) .GT. 0) THEN
+         temp = 0
+         DO iz=1,basin_sz(ij)
+           temp =  temp + (orog_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) -basin_orog_mean(ib,ij))**2
+         END DO
+         IF (basin_sz(ij) .GT. 0) basin_orog_std(ib,ij) = SQRT(temp)/(basin_sz(ij))
+         IF (basin_orog_std(ib,ij) .LT. lowlim_std) basin_orog_std(ib,ij) = lowlim_std
+         IF (basin_orog_std(ib,ij) .GT. uplim_std) basin_orog_std(ib,ij) = uplim_std
+       END IF
+       ! Calculation of beta floodplains
+       IF (basin_floodp(ib,ij) .GT. 0) THEN
+         temp = (basin_orog_std(ib,ij)-lowlim_std)/(uplim_std-lowlim_std)
+         basin_beta_fp(ib,ij) = temp*(uplim_beta-lowlim_beta)+0.2
+       END IF
+       !
        SELECT CASE (hierar_method)
           !
           CASE("MEAN")
@@ -2586,7 +2617,7 @@ SUBROUTINE isin_halo(ig, nbhalo, nbpt, halopts, isinhalo)
 
 SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, 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_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind,basin_beta_fp, fetch_basin,&
        & basin_id, basin_coor, basin_type, basin_flowdir, &
        & outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
     !
@@ -2620,6 +2651,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
     REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_floodp !!
     REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout)  :: basin_cg !!
     REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_topoind !! Topographic index of the residence time for a basin (m)
+    REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: basin_beta_fp !!
     REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout)    :: fetch_basin !!
     INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
     INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
@@ -2680,6 +2712,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
        routing_orog_min_glo(ib,:) = 0
        routing_orog_max_glo(ib,:) = 0
        routing_floodcri_glo(ib,:) = 0
+       routing_beta_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
@@ -2688,6 +2721,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
          basin_orog_max(ib,1) = 0
          basin_orog_mean(ib,1) = 0
          basin_floodp(ib,1) = 0
+         basin_beta_fp(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
@@ -2717,6 +2751,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
           routing_orog_max_glo(ib,ij) = basin_orog_max(ib,ij)
           ! Now the final floodplains is converted into fraction
           routing_floodp_glo(ib,ij) = basin_floodp(ib,ij)/basin_area(ib,ij)
+          routing_beta_glo(ib,ij) = basin_beta_fp(ib,ij)
           routing_cg_glo(ib,ij,:) = basin_cg(ib,ij,:)
           topo_resid_glo(ib,ij) = basin_topoind(ib,ij)
           global_basinid_glo(ib,ij) = basin_id(ib,ij)
@@ -2828,7 +2863,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
 !_ ================================================================================================================================
 
 SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, basin_count, basin_area, &
-       & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind,&
+       & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind, basin_beta_fp,&
        & fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
        & inflow_grid, inflow_basin)
     !
@@ -2853,6 +2888,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
     REAL(r_std), DIMENSION(nbpt,nwbas)    :: basin_floodp !!
     REAL(r_std), DIMENSION(nbpt,nwbas,2)  :: basin_cg !!
     REAL(r_std), DIMENSION(nbpt,nwbas)    :: basin_topoind !! Topographic index of the residence time for a basin (m)
+    REAL(r_std), DIMENSION(nbpt,nwbas)    :: basin_beta_fp !!
     REAL(r_std), DIMENSION(nbpt,nwbas)    :: fetch_basin !!
     INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
     INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !!
@@ -3065,6 +3101,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
     ALLOCATE (route_type_glo(nbpt,nbasmax), stat=ier)
     ALLOCATE (origin_nbintobas_glo(nbpt), stat=ier)
     ALLOCATE (routing_floodcri_glo(nbpt,nbasmax), stat=ier)
+    ALLOCATE (routing_beta_glo(nbpt,nbasmax), stat=ier)
 
   END SUBROUTINE allocategraph
 END MODULE routing_reg

GraphHydro.py

@@ -62,7 +62,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.routing_floodp,self.routing_beta, 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.floodcri = \
                                     routing_interface.finish_truncate(nbpt = self.nbpt, inflowmax = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
@@ -72,8 +72,8 @@ class HydroGraph :
                                                                       basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
                                                                       basin_orog_mean = hydrosuper.basin_orog_mean, basin_orog_min = hydrosuper.basin_orog_min,\
                                                                       basin_orog_max = hydrosuper.basin_orog_max, basin_floodp = hydrosuper.basin_floodp, \
-                                                                      basin_cg = hydrosuper.basin_cg, \
-                                                                      basin_topoind = hydrosuper.basin_topoind, fetch_basin = hydrosuper.fetch_basin, \
+                                                                      basin_cg = hydrosuper.basin_cg, basin_topoind = hydrosuper.basin_topoind, \
+                                                                      basin_beta_fp = hydrosuper.basin_beta_fp , fetch_basin = hydrosuper.fetch_basin, \
                                                                       basin_id = hydrosuper.basin_id, \
                                                                       basin_coor = hydrosuper.basin_outcoor, basin_type = hydrosuper.basin_type, \
                                                                       basin_flowdir = hydrosuper.basin_flowdir, \
@@ -217,6 +217,8 @@ class HydroGraph :
 
         self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "floodcri", "Height at which all basin is flooded", "mm", self.floodcri[:,:], vtyp)
 
+        self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), "basin_beta_fp", "Shape of the floodplins HTUs (concave / convex)", "-", self.routing_beta[:,:], 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")
         #

Interface.py

@@ -30,6 +30,7 @@ config = getargs.SetupConfig()
 gendoc = config.get("OverAll", "Documentation", fallback='false')
 nbxmax = config.getint("OverAll", "nbxmax", fallback=63)
 largest_pos = config.getint("OverAll", "ROUTING_RIVERS", fallback=200)
+maxpercent = config.getint("OverAll", "maxpercent", fallback=2)
 #
 undef_int = 999999999.9
 # Order of magnitude for the area precision in m^2.
@@ -140,7 +141,7 @@ class HydroOverlap :
         self.nbi, self.nbj, self.area_bx, self.trip_bx, self.basin_bx, self.topoind_bx, self.fac_bx, self.hierarchy_bx, \
             self.orog_bx, self.floodp_bx, \
             self.lon_bx, self.lat_bx, self.lshead_bx = \
-                    routing_interface.gethydrogrid(ijdimmax, sub_pts, sub_index, sub_area, \
+                    routing_interface.gethydrogrid(ijdimmax, maxpercent, sub_pts, sub_index, sub_area, \
                     hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp,\
                         hierarchy_tmp, orog_tmp, floodp_tmp)
         #
@@ -202,7 +203,8 @@ class HydroSuper :
         self.basin_count, self.basin_notrun, self.basin_area, self.basin_cg, self.basin_hierarchy, \
             self.basin_orog_mean, self.basin_orog_min, self.basin_orog_max, self.basin_floodp, \
             self.basin_fac, self.basin_topoind, self.basin_id, self.basin_outcoor, self.basin_type,\
-            self.basin_flowdir, self.basin_lshead, self.outflow_grid, self.outflow_basin, self.nbcoastal, self.coastal_basin = \
+            self.basin_flowdir, self.basin_lshead, self.basin_beta_fp, self.outflow_grid, self.outflow_basin,\
+            self.nbcoastal, self.coastal_basin = \
                     routing_interface.globalize(nbpt = self.nbpt, nb_htu = self.nbhtuext, nbv = nbv, ijdimmax = self.ijdimmax, \
                                                 area_bx = hydrooverlap.area_bx, lon_bx = lon_bx_tmp, lat_bx = lat_bx_tmp, trip_bx = hydrooverlap.trip_bx, \
                                                 hierarchy_bx = hydrooverlap.hierarchy_bx, orog_bx = hydrooverlap.orog_bx, floodp_bx =  hydrooverlap.floodp_bx,\
@@ -359,7 +361,7 @@ class HydroSuper :
                 totakeover = totakeover, numops = numops, basin_count = self.basin_count,\
                 basin_area = self.basin_area, basin_orog_mean = self.basin_orog_mean, \
                 basin_orog_min = self.basin_orog_min, basin_orog_max = self.basin_orog_max, basin_floodp = self.basin_floodp, \
-                basin_cg = self.basin_cg, basin_topoind = self.basin_topoind, fetch_basin = self.fetch_basin,\
+                basin_cg = self.basin_cg, basin_topoind = self.basin_topoind, basin_beta_fp = self.basin_beta_fp, fetch_basin = self.fetch_basin,\
                 basin_id = self.basin_id, basin_coor = self.basin_outcoor, basin_type = self.basin_type,\
                 basin_flowdir = self.basin_flowdir, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, \
                 inflow_number = self.inflow_number, inflow_grid = self.inflow_grid, inflow_basin = self.inflow_basin)