diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 26581492fb497debf2d979b91782f3fb0419f780..89f67adda14e1604b5878baeda2225bd027b4ab2 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -87,7 +87,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
!
INTEGER :: ii, ib
REAL :: resolution(nbpt,2)
-
+
!
! These two parameters are part of the routing_reg module saved variables. They are transfered here.
!
@@ -98,9 +98,9 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
resolution(ii,1) = SQRT(area(ii))
resolution(ii,2) = SQRT(area(ii))
ENDDO
-
+
nwbas = MAXVAL(sub_pts)
-
+
DO ib=1,nbpt
CALL routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, fac, hierarchy, &
@@ -137,7 +137,7 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
INTEGER, INTENT(out) :: nb_basin(nbpt) !! Number of sub-basins (unitless)
INTEGER, INTENT(out) :: basin_inbxid(nbpt,nbvmax_in) !!
REAL, INTENT(out) :: basin_outlet(nbpt,nbvmax_in,2) !!
- REAL, INTENT(out) :: basin_outtp(nbpt,nbvmax_in) !!
+ REAL, INTENT(out) :: basin_outtp(nbpt,nbvmax_in) !!
INTEGER, INTENT(out) :: basin_sz(nbpt,nbvmax_in) !!
INTEGER, INTENT(out) :: basin_bxout(nbpt,nbvmax_in) !!
INTEGER, INTENT(out) :: basin_bbout(nbpt,nbvmax_in) !!
@@ -160,7 +160,7 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
WRITE(*,*) "nbvmax or nbxmax have changed !!"
STOP
ENDIF
-
+
DO ib=1,nbpt
CALL routing_reg_findbasins(ib, nbi(ib), nbj(ib), trip_bx(ib,:,:), basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), &
& topoind_bx(ib,:,:), lshead_bx(ib,:,:), diaglalo, &
@@ -197,7 +197,7 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
INTEGER(i_std), INTENT (in), DIMENSION(nbpt) :: nb_basin !! Number of sub-basins (unitless)
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in,2) :: basin_outlet !! Outlet coordinate of each subgrid basin
- REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_outtp !!
+ REAL(r_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_outtp !!
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in,nbvmax_in,2) :: basin_pts !! Points in each basin
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_bxout !! outflow direction
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbvmax_in) :: basin_bbout !! outflow sub-basin
@@ -218,7 +218,7 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
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.
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,nwbas) :: outflow_basin !!
INTEGER(i_std), INTENT (out), DIMENSION(nbpt) :: nbcoastal !!
INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: coastal_basin !!
!!
@@ -279,7 +279,7 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
WRITE(*,*) "LINKUP : nbxmax has changed !!"
STOP
ENDIF
-
+
CALL routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, &
& basin_lshead, basin_hierarchy, basin_fac, diaglalo, outflow_grid, outflow_basin, inflow_number, inflow_grid, &
& inflow_basin, nbcoastal, coastal_basin, invented_basins)
@@ -327,7 +327,7 @@ SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, outflow_grid, basin_ar
ENDDO
!
END SUBROUTINE areanorm
-
+
SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac, &
& outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
!
@@ -369,7 +369,7 @@ SUBROUTINE fetch(nbpt, nwbas, nbcore, corepts, basin_count, basin_area, basin_id
!
CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, nbcore, fcorepts, basin_count, basin_area, basin_id, &
& basin_hierarchy, basin_fac, outflow_grid, outflow_basin, partial_sum, fetch_basin, outflow_uparea)
-
+
!
END SUBROUTINE fetch
@@ -419,7 +419,8 @@ SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_
!
END SUBROUTINE rivclassification
-SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, nbcore, corepts, basin_count, basin_notrun, basin_area, &
+
+SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, basin_notrun, basin_area, &
& 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_cg, topo_resid, route_nbbasin, route_togrid, &
& route_tobasin, route_nbintobas, global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
@@ -436,8 +437,8 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, nbcore, corept
INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax
INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), INTENT (in) :: num_largest
- INTEGER(i_std), INTENT (in) :: nbcore
- INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
+ REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea
+ REAL(r_std), DIMENSION(nbpt), INTENT(in) :: cfrac
!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_notrun !!
@@ -454,9 +455,9 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, nbcore, corept
!
! Changed nwbas to nbxmax*4 and *8 to save the memory
!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in*8), INTENT(inout) :: inflow_number !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in*8,nbxmax_in*8), INTENT(inout) :: inflow_basin !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in*8,nbxmax_in*8), INTENT(inout) :: inflow_grid !!
+ INTEGER(i_std), DIMENSION(nbpt,nbxmax_in), INTENT(inout) :: inflow_number !!
+ INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_basin !!
+ INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_grid !!
!
! Output variables
!
@@ -464,7 +465,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, nbcore, corept
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)
- INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: route_nbbasin !!
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: route_nbbasin !!
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_togrid !! Grid into which the basin flows (unitless)
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_tobasin !! Basin in to which the water goes (unitless)
INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: route_nbintobas !! Number of basin into current one (unitless)
@@ -473,19 +474,10 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, nbcore, corept
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 (-)
!
- INTEGER(i_std) :: ic
- INTEGER(i_std), DIMENSION(nbcore) :: fcorepts
!
- IF ( nbxmax_in .NE. nbxmax ) THEN
- WRITE(*,*) "TRUNCATE : nbxmax has changed !!"
- STOP
- ENDIF
!
- DO ic=1,nbcore
- fcorepts(ic) = corepts(ic)+1
- ENDDO
-
- CALL routing_reg_truncate(nbpt, nbasmax, area, contfrac, nwbas, num_largest, nbcore, fcorepts, basin_count, &
+
+ CALL routing_reg_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, basin_count, &
& basin_notrun, basin_area, basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
@@ -497,7 +489,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, nbcore, corept
route_tobasin(:,:) = route_tobasin_glo(:,:)
routing_fetch(:,:) = route_fetch_glo(:,:)
-
+
route_nbintobas(:) = route_nbintobas_glo(:)
global_basinid(:,:) = global_basinid_glo(:,:)
route_outlet(:,:,:) = route_outlet_glo(:,:,:)
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 215e5346e2ada071862cd71f45508fb910bba9c1..9ffbac79d6fe6757c76ab1e5cce91913e050f1c7 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2537,7 +2537,7 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, nbcore, corepts, basin_count, &
+SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, basin_count, &
& basin_notrun, basin_area, basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
!
@@ -2556,8 +2556,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
!
INTEGER(i_std), INTENT(in) :: nwbas !!
INTEGER(i_std), INTENT (in) :: num_largest
- INTEGER(i_std), INTENT (in) :: nbcore
- INTEGER(i_std), DIMENSION(nbcore), INTENT(in) :: corepts
+
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_notrun !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_id !!
@@ -2624,8 +2623,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
! core area of the domain.
!
nbtruncate = 0
- DO ic = 1, nbcore
- ib = corepts(ic)
+ DO ib = 1, nbpt
IF ( basin_count(ib) .GT. nbasmax ) THEN
nbtruncate = nbtruncate + 1
indextrunc(nbtruncate) = ib
@@ -3122,13 +3120,16 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
! Update the information needed in the basin "totakeover"
! For the moment only area
!
+ basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1)*basin_area(ib, totakeover) &
+ & + basin_cg(ib, tokill, 1)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
- basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib, tokill)
+ basin_cg(ib, totakeover, 2) = (basin_cg(ib, totakeover, 2)*basin_area(ib, totakeover)&
+ & + basin_cg(ib, tokill, 2)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
- basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1) + basin_cg(ib, tokill, 1))/2.0
- basin_cg(ib, totakeover, 2) = (basin_cg(ib, totakeover, 2) + basin_cg(ib, tokill, 2))/2.0
+ basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
+ & + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
- basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover) + basin_topoind(ib, tokill))/2.0
+ basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib,tokill)
!
! Add the fetch of the basin will kill to the one which gets the water
!
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index 9513a1d2009aa0f54e19054df9305d2a5b71f5c8..668a11257d2a60559fd71e9f90765f061cecd79d 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -31,6 +31,7 @@ import RPPtools as RPP
import HydroGrid as HG
import diagplots as DP
import Interface as IF
+from Truncate import truncation as TR
from spherical_geometry import vector
#
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
@@ -194,8 +195,14 @@ print("=================== Compute fetch ====================")
hsuper.fetch(part)
hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part)
-hgraph = IF.HydroGraph(nbasmax, hsuper, part)
-
-hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part)
+#hgraph = IF.HydroGraph(nbasmax, hsuper, part)
+#hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part)
+# For now we use a one proc truncation
+if rank ==0:
+ hs = TR(OutGraphFile.replace(".nc","_HydroSuper.nc"), nbasmax)
+ hs.get_fetch()
+ hs.truncate()
+ hs.finalrivclass(hsuper.largest_rivarea)
+ hs.dumpnetcdf(OutGraphFile)
IF.closeinterface(comm)