diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 1f5b8492bc0bff10f85a116a7c8753f86b6980c2..d9cef5eff8b566d4b1bf63f2b487963d22968074 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -417,7 +417,7 @@ SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_
END SUBROUTINE rivclassification
-SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, basin_notrun, basin_area, &
+SUBROUTINE finish_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)
@@ -474,7 +474,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfra
!
!
- CALL routing_reg_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, basin_count, &
+ CALL routing_reg_end_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)
@@ -493,7 +493,188 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfra
route_type(:,:) = route_type_glo(:,:)
origin_nbintobas(:) = origin_nbintobas_glo(:)
-END SUBROUTINE truncate
+END SUBROUTINE finish_truncate
+
+
+SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numops, basin_count, 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)
+ !
+ !
+ 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) :: nbxmax_in, nbasmax
+ INTEGER(i_std), INTENT (in) :: nwbas !!
+ INTEGER(i_std), INTENT (in) :: ops !!
+ INTEGER(i_std), DIMENSION(nbpt, ops), INTENT (in) :: tokill, totakeover
+ INTEGER(i_std), DIMENSION(nbpt), INTENT (in) :: numops
+ !
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_id !!
+ REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_coor !!
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_type !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_flowdir !! Water flow directions in the basin (unitless)
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
+ 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) :: 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 !!
+ !
+ ! Changed nwbas to nbxmax*4 and *8 to save the memory
+ !
+ 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 !!
+
+ !! LOCAL
+ INTEGER(i_std) :: ib, op, tok, totak
+
+
+ DO ib=1,nbpt
+ DO op=1,numops(ib)
+ IF (basin_count(ib) > nbasmax) THEN
+ tok = tokill(ib,op)
+ totak = totakeover(ib,op)
+ IF (tok .GT. 0) THEN
+ WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak
+ CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, 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)
+ END IF
+ END IF
+ END DO
+ END DO
+
+ ! REPRENDRE LA FIN DE TRUNCATE
+
+END SUBROUTINE killbas
+
+
+SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
+ !
+ 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), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_grid !! Type of outflow on the grid box (unitless)
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
+
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
+ ! Local
+ INTEGER(i_std), DIMENSION(nbpt,nwbas) :: flag !!
+ INTEGER(i_std) :: ig, ib, it, igrif, ibasf, basnum, test
+
+
+ flag(:,:) = 0
+
+ WRITE(numout,*) "Checking routing"
+ test = 0
+
+ DO ig=1,nbpt
+ basnum = basin_count(ig)
+ DO ib=1,basnum
+ IF (flag(ig,ib) .EQ. 0) THEN
+ flag(ig,ib) = -1
+ igrif = outflow_grid(ig,ib)
+ ibasf = outflow_basin(ig,ib)
+ !
+ DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -1) .AND. (flag(igrif,ibasf) .NE. -99))
+ flag(igrif, ibasf) = -1
+ it = outflow_grid(igrif,ibasf)
+ ibasf = outflow_basin(igrif,ibasf)
+ igrif = it
+ END DO
+
+
+ IF ((flag(igrif,ibasf) .EQ. -99) .OR. (igrif .LE. 0)) THEN
+ flag(ig,ib) = -99
+ igrif = outflow_grid(ig,ib)
+ ibasf = outflow_basin(ig,ib)
+
+ DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -99))
+ flag(igrif, ibasf) = -99
+ it = outflow_grid(igrif,ibasf)
+ ibasf = outflow_basin(igrif,ibasf)
+ igrif = it
+ END DO
+
+ ELSE IF (flag(igrif,ibasf) .EQ. -1) THEN
+ test = test + 1
+
+ flag(ig,ib) = -99
+ igrif = outflow_grid(ig,ib)
+ ibasf = outflow_basin(ig,ib)
+ DO WHILE (flag(igrif,ibasf) .EQ. -1)
+ flag(igrif, ibasf) = -99
+ it = outflow_grid(igrif,ibasf)
+ ibasf = outflow_basin(igrif,ibasf)
+ igrif = it
+ END DO
+ END IF
+ END IF
+ END DO
+ END DO
+
+ WRITE(numout,*) "**** Fetch has",test, "loop errors"
+
+
+END SUBROUTINE checkrouting
+
+
+SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, outflow_grid, outflow_basin, basin_count)
+ !
+ 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 !!
+
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: fetch_basin !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_grid !! Type of outflow on the grid box (unitless)
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
+ ! Local
+ INTEGER(i_std) :: ig, ib, it, bt, igrif, ibasf, basnum, test
+
+ WRITE(numout,*) "Checking Fetch coherence"
+ test = 0
+
+ DO ig=1,nbpt
+ basnum = basin_count(ig)
+
+ DO ib=1,basnum
+ igrif = outflow_grid(ig,ib)
+ ibasf = outflow_basin(ig,ib)
+
+ IF (igrif .GT. 0) THEN
+ IF (fetch_basin(ig,ib) .GT. fetch_basin(igrif,ibasf)) THEN
+ WRITE(numout,*) "****ERROR Fetch, nbpt:",ig
+ WRITE(numout,*) igrif, ibasf
+ WRITE(numout,*) it, bt
+ test = 1
+ END IF
+ END IF
+ END DO
+ END DO
+ IF (test .EQ. 0) WRITE(numout,*) "**** Fetch has no error"
+
+END SUBROUTINE checkfetch
+
+
SUBROUTINE finalfetch(nbpt, nbasmax, nbcore, corepts, routing_area, basin_count, route_togrid, route_tobasin, partial_sum, &
& fetch_basin, outflow_uparea)
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 0fb14293419f4207a94abc08b25b991c0fe9a37e..b9d9d53241c53b496692680863a5e63466fe2f51 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2554,7 +2554,7 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, basin_count, &
+SUBROUTINE routing_reg_end_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)
!
@@ -2623,352 +2623,6 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
!
! Read option for truncation method
!
- trunc_equa = .FALSE.
-!!$ CALL getin('ROUTING_TRUNCEQUA', trunc_equa)
- !
- IF ( .NOT. ALLOCATED(indextrunc)) THEN
- ALLOCATE(indextrunc(nbpt))
- ENDIF
- !
- ! We have to go through the grid as least as often as we have to reduce the number of basins
- ! For good measure we add 3 more passages.
- !
- !
- DO iter = 1, MAXVAL(basin_count) - nbasmax +3
- !
- ! Get the points over which we wish to truncate. We only work with points of the
- ! core area of the domain.
- !
- nbtruncate = 0
- DO ib = 1, nbpt
- IF ( basin_count(ib) .GT. nbasmax ) THEN
- nbtruncate = nbtruncate + 1
- indextrunc(nbtruncate) = ib
- ENDIF
- ENDDO
- !
- ! Go through the basins which need to be truncated.
- !
- DO ibt=1,nbtruncate
- !
- ib = indextrunc(ibt)
- !
- ! Check if we have basin which flows into a basin in the same grid
- ! kbas = basin we will have to kill
- ! sbas = basin which takes over kbas
- !
- kbas = 0
- sbas = 0
- !
- ! 1) Merge two basins which flow into the ocean as coastal or return flow
- ! (outflow_grid = -2 or -3). Well obviously only if we have more than 1 and
- ! have not found anything yet!
- !
- IF ( (COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 .OR.&
- & COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -3) .GT. 1) .AND.&
- & kbas*sbas .EQ. 0) THEN
- !
- multbas = 0
- multbas_sz(:) = 0
- !
- IF ( COUNT(outflow_grid(ib,1:basin_count(ib)) .EQ. -2) .GT. 1 ) THEN
- obj = -2
- ELSE
- obj = -3
- ENDIF
- !
- ! First we get the list of all basins which go out as coastal or return flow (obj)
- !
- DO ij=1,basin_count(ib)
- IF ( outflow_grid(ib,ij) .EQ. obj ) THEN
- multbas = multbas + 1
- multbas_sz(multbas) = ij
- tmp_area(multbas) = fetch_basin(ib,ij)
- ENDIF
- ENDDO
- !
- ! Test the new truncation priority
- !
- IF ( trunc_equa ) THEN
- !
- ! Take the smallest to be transfered to the second smallest
- !
- iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
- kbas = multbas_sz(iml(1))
- tmp_area(iml(1)) = -1
- iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
- sbas = multbas_sz(iml(1))
- !
- ELSE
- !
- ! Take the smallest to be transfered to the largest
- !
- iml = MAXLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
- sbas = multbas_sz(iml(1))
- iml = MINLOC(tmp_area(1:multbas), MASK = tmp_area(1:multbas) .GT. zero)
- kbas = multbas_sz(iml(1))
- !
- ENDIF
- ENDIF
- !
- !
- ! 2) If we have basins flowing into the same grid but different basins then we put them
- ! together. Obviously we first work with the grid which has most streams running into it
- ! and putting the smallest in the largests catchments.
- !
- IF ( kbas*sbas .EQ. 0) THEN
- !
- tmp_ids(1:basin_count(ib)) = outflow_grid(ib,1:basin_count(ib))
- multbas = 0
- multbas_sz(:) = 0
- !
- ! First obtain the list of basins which flow into the same basin
- !
- DO ij=1,basin_count(ib)
- IF ( outflow_grid(ib,ij) .GT. 0 .AND.&
- & COUNT(tmp_ids(1:basin_count(ib)) .EQ. outflow_grid(ib,ij)) .GT. 1) THEN
- multbas = multbas + 1
- DO ii=1,basin_count(ib)
- IF ( tmp_ids(ii) .EQ. outflow_grid(ib,ij)) THEN
- multbas_sz(multbas) = multbas_sz(multbas) + 1
- multbas_list(multbas,multbas_sz(multbas)) = ii
- tmp_ids(ii) = -99
- ENDIF
- ENDDO
- ELSE
- tmp_ids(ij) = -99
- ENDIF
- ENDDO
- !
- ! Did we come up with any basins to deal with this way ?
- !
- IF ( multbas .GT. 0 ) THEN
- !
- iml = MAXLOC(multbas_sz(1:multbas))
- ik = iml(1)
- !
- ! Take the smallest and largest of these basins !
- !
- DO ii=1,multbas_sz(ik)
- tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
- ENDDO
- !
- IF ( trunc_equa ) THEN
- iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- kbas = multbas_list(ik,iml(1))
- tmp_area(iml(1)) = -1
- iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- sbas = multbas_list(ik,iml(1))
- ELSE
- iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- sbas = multbas_list(ik,iml(1))
- iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- kbas = multbas_list(ik,iml(1))
- ENDIF
- !
- ENDIF
- !
- ENDIF
- !
- ! 3) If we have twice the same basin we put them together even if they flow into different
- ! directions. If one of them goes to the ocean it takes the advantage.
- !
- IF ( kbas*sbas .EQ. 0) THEN
- !
- tmp_ids(1:basin_count(ib)) = basin_id(ib,1:basin_count(ib))
- multbas = 0
- multbas_sz(:) = 0
- !
- ! First obtain the list of basins which have sub-basins in this grid box.
- ! (these are identified by their IDs)
- !
- DO ij=1,basin_count(ib)
- IF ( COUNT(tmp_ids(1:basin_count(ib)) .EQ. basin_id(ib,ij)) .GT. 1) THEN
- multbas = multbas + 1
- DO ii=1,basin_count(ib)
- IF ( tmp_ids(ii) .EQ. basin_id(ib,ij)) THEN
- multbas_sz(multbas) = multbas_sz(multbas) + 1
- multbas_list(multbas,multbas_sz(multbas)) = ii
- tmp_ids(ii) = -99
- ENDIF
- ENDDO
- ELSE
- tmp_ids(ij) = -99
- ENDIF
- ENDDO
- !
- ! We are going to work on the basin with the largest number of sub-basins.
- ! (IF we have a basin which has subbasins !)
- !
- IF ( multbas .GT. 0 ) THEN
- !
- iml = MAXLOC(multbas_sz(1:multbas))
- ik = iml(1)
- !
- ! If one of the basins goes to the ocean then it is going to have the priority
- !
- tmp_area(:) = zero
- IF ( COUNT(outflow_grid(ib,multbas_list(ik,1:multbas_sz(ik))) .LT. 0) .GT. 0) THEN
- DO ii=1,multbas_sz(ik)
- IF ( outflow_grid(ib,multbas_list(ik,ii)) .LT. 0 .AND. sbas .EQ. 0 ) THEN
- sbas = multbas_list(ik,ii)
- ELSE
- tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
- ENDIF
- ENDDO
- ! take the smallest of the subbasins
- iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- kbas = multbas_list(ik,iml(1))
- ELSE
- !
- ! Else we take simply the largest and smallest
- !
- DO ii=1,multbas_sz(ik)
- tmp_area(ii) = fetch_basin(ib, multbas_list(ik,ii))
- ENDDO
- !
- IF ( trunc_equa ) THEN
- iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- kbas = multbas_list(ik,iml(1))
- tmp_area(iml(1)) = -1
- iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- sbas = multbas_list(ik,iml(1))
- ELSE
- iml = MAXLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- sbas = multbas_list(ik,iml(1))
- iml = MINLOC(tmp_area(1:multbas_sz(ik)), MASK = tmp_area(1:multbas_sz(ik)) .GT. zero)
- kbas = multbas_list(ik,iml(1))
- ENDIF
- !
- ENDIF
- !
- !
- ENDIF
- ENDIF
- !
- ! Then we call routing_reg_killbas to clean up the basins in this grid
- !
- IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
- CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, 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)
- ENDIF
- !
- ENDDO
- !
- ENDDO
- !
- ! We only merge a basin which flows into a basin of the same grid if it's
- ! necessary.
- !
- IF ( COUNT(basin_count .GT. nbasmax) .GT. 0 ) THEN
- !
- ! Get the points over which we wish to truncate
- !
- nbtruncate = 0
- DO ib = 1, nbpt
- IF ( basin_count(ib) .GT. nbasmax ) THEN
- nbtruncate = nbtruncate + 1
- indextrunc(nbtruncate) = ib
- ENDIF
- ENDDO
- !
- ! Go through the basins which need to be truncated.
- !
- DO ibt=1,nbtruncate
- !
- ib = indextrunc(ibt)
- !
- ! Check if we have basin which flows into a basin in the same grid
- ! kbas = basin we will have to kill
- ! sbas = basin which takes over kbas
- !
- kbas = 0
- sbas = 0
- !
- ! 4) Can we find a basin which flows into a basin of the same grid ?
- !
- DO ij=1,basin_count(ib)
- DO ii=1,basin_count(ib)
- IF ( outflow_grid(ib,ii) .EQ. ib .AND. outflow_basin(ib, ii) .EQ. ij .AND. kbas*sbas .NE. 0) THEN
- kbas = ii
- sbas = ij
- ENDIF
- ENDDO
- ENDDO
- !
- ! Then we call routing_reg_killbas to clean up the basins in this grid
- !
- IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
- CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, 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)
- ENDIF
- !
- ENDDO
- !
- ENDIF
- !
- ! If there are any grids left with too many basins we need to take out the big hammer !
- ! We will only do it if this represents less than 5% of all points.
- !
- IF ( COUNT(basin_count .GT. nbasmax) .GT. 0 ) THEN
- !
- !
- IF ( COUNT(basin_count .GT. nbasmax)/nbpt*100 .GT. 5 ) THEN
- WRITE(numout,*) 'We have ', COUNT(basin_count .GT. nbasmax)/nbpt*100, '% of all points which do not yet'
- WRITE(numout,*) 'have the right trunctaction. That is too much to apply a brutal method'
- DO ib = 1, nbpt
- IF ( basin_count(ib) .GT. nbasmax ) THEN
- !
- WRITE(numout,*) 'We did not find a basin which could be supressed. We will'
- WRITE(numout,*) 'not be able to reduce the truncation in grid ', ib
- DO ij=1,basin_count(ib)
- WRITE(numout,*) 'grid, basin nb and id :', ib, ij, basin_id(ib,ij)
- WRITE(numout,*) 'Outflow grid and basin ->', outflow_grid(ib,ij), outflow_basin(ib, ij)
- ENDDO
- ENDIF
- ENDDO
- CALL ipslerr_p(3,'routing_reg_truncate','No basin found which could be supressed.','','')
- ELSE
- !
- !
- DO ib = 1,nbpt
- DO WHILE ( basin_count(ib) .GT. nbasmax )
- !
- ! Here we simply put the smallest basins into the largest ones. It is really a brute force
- ! method but it will only be applied if everything has failed.
- !
- DO ii = 1,basin_count(ib)
- tmp_area(ii) = fetch_basin(ib, ii)
- ENDDO
- !
- IF ( trunc_equa ) THEN
- iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
- kbas = iml(1)
- tmp_area(iml(1)) = -1
- iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
- sbas =iml(1)
- ELSE
- iml = MAXLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
- sbas =iml(1)
- iml = MINLOC(tmp_area(1:basin_count(ib)), MASK = tmp_area(1:basin_count(ib)) .GT. 0.)
- kbas = iml(1)
- ENDIF
- !
- IF ( kbas .GT. 0 .AND. sbas .GT. 0 ) THEN
- CALL routing_reg_killbas(nbpt, ib, kbas, sbas, nwbas, basin_count, 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)
- ENDIF
- ENDDO
- ENDDO
- !
- ENDIF
- !
- !
- ENDIF
- !
! Now that we have reached the right truncation (resolution) we will start
! to produce the variables we will use to route the water.
!
@@ -3077,7 +2731,7 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
WRITE(numout,*) 'Surface of all continents :', SUM(gridarea_tmp)
!
- END SUBROUTINE routing_reg_truncate
+ END SUBROUTINE routing_reg_end_truncate
!
!! ================================================================================================================================
!! SUBROUTINE : routing_reg_killbas
@@ -3175,6 +2829,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
igrif = it
ENDDO
!
+ !
+ !
! Redirect the flows which went into the basin to be killed before we change everything
!
DO inf = 1, inflow_number(ib, tokill)
diff --git a/Interface.py b/Interface.py
index 471ca5dbc17f8a00b23c4f096c16dfe37bce7f63..827855a7d40bd6e96a149f225c8d0afe272208cc 100644
--- a/Interface.py
+++ b/Interface.py
@@ -51,7 +51,10 @@ if gendoc.lower() == "true" :
docwrapper.write("====================================================================\n")
docwrapper.write(routing_interface.fetch.__doc__)
docwrapper.write("====================================================================\n")
- docwrapper.write(routing_interface.truncate.__doc__)
+ docwrapper.write(routing_interface.finish_truncate.__doc__)
+ docwrapper.write("====================================================================\n")
+ docwrapper.write(routing_interface.killbas.__doc__)
+
docwrapper.close
#
# Functions to access the interfaces
@@ -204,6 +207,7 @@ def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fet
maxdiff_sorted = np.max(np.abs(sorted_outareas[0:largest_pos]-old_sorted))
old_sorted[:] = sorted_outareas[0:largest_pos]
iter_count += 1
+
#
fetch_error = np.sum(np.abs(fetch_out[part.landcorelist,:]-fetch_in[part.landcorelist,:]), axis=1)\
/np.sum(routing_area[part.landcorelist,:], axis=1)
@@ -233,7 +237,7 @@ class HydroOverlap :
for ip in range(sub_pts[ib]) :
sub_index[ib,ip,:] = [sub_index_in[ib][0][ip],sub_index_in[ib][1][ip]]
#
- part.landsendtohalo(sub_area, order='F')
+ part.landsendtohalo(np.array(sub_area), order='F')
#
trip_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
basins_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
@@ -275,11 +279,13 @@ class HydroOverlap :
#
#
class HydroSuper :
- def __init__(self, nbvmax, hydrodata, hydrooverlap) :
+ def __init__(self, nbvmax, hydrodata, hydrooverlap, nbasmax) :
#
# Keep largest possible number of HTUs
#
+ self.nbasmax = nbasmax
self.nbhtuext = nbvmax
+ self.nwbas = hydrooverlap.nwbas
#
# Call findbasins
#
@@ -302,7 +308,9 @@ class HydroSuper :
hydrooverlap.hierarchy_bx, hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \
nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, self.basin_pts, basin_bxout, \
basin_bbout, basin_lshead, coast_pts, hydrooverlap.nwbas)
+
self.nbpt = self.basin_count.shape[0]
+
return
#
def linkup(self, hydrodata) :
@@ -314,8 +322,10 @@ class HydroSuper :
self.basin_flowdir, self.basin_lshead, self.basin_hierarchy, \
self.basin_fac, self.outflow_grid, self.outflow_basin, \
self.nbcoastal, self.coastal_basin, float(hydrodata.basinsmax))
+ self.nbxmax_in = self.inflow_number.shape[1]
return
#
+
def fetch(self, part) :
#
fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
@@ -357,11 +367,31 @@ class HydroSuper :
#
#
#
- self.num_largest = routing_interface.rivclassification(part.landcorelist, self.basin_count, self.outflow_grid, self.outflow_basin, \
- self.fetch_basin, self.largest_rivarea)
- print("Rank :"+str(part.rank)+" Area of smallest large rivers : ", self.largest_rivarea, " Nb of Large rivers on proc : ",self.num_largest)
+ self.num_largest = routing_interface.rivclassification( part.landcorelist, self.basin_count, self.outflow_grid, self.outflow_basin, \
+ self.fetch_basin, self.largest_rivarea)
+ #print("Rank :"+str(part.rank)+" Area of smallest large rivers : ", self.largest_rivarea, " Nb of Large rivers on proc : ",self.num_largest)
return
+
+ def check_fetch(self):
+
+ routing_interface.checkfetch(nbpt = self.nbpt, nwbas = self.nwbas, fetch_basin = self.fetch_basin, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, basin_count = self.basin_count)
+
+ return
+
+ def check_routing(self):
+
+ routing_interface.checkrouting(nbpt = self.nbpt, nwbas = self.nwbas, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, basin_count = self.basin_count)
+
+ return
#
+ #
+ def killbas(self, tokill, totakeover, numops):
+ ops = tokill.shape[1]
+
+ routing_interface.killbas(nbpt = self.nbpt, nbxmax_in = self.nbxmax_in, nbasmax = self.nbasmax, nwbas = self.nwbas, ops = ops, tokill = tokill, totakeover = totakeover, numops = numops, basin_count = self.basin_count, basin_area = self.basin_area, \
+ basin_cg = self.basin_cg, basin_topoind = self.basin_topoind, 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)
+
#
def add_variable(self,outnf, procgrid, NCFillValue, part, coord, name, title, units, data, vtyp):
var = procgrid.landscatter(data.astype(vtyp), order='F')
@@ -477,24 +507,27 @@ class HydroSuper :
#
#
class HydroGraph :
- def __init__(self, nbasmax, hydrosuper, part) :
+ def __init__(self, nbasmax, hydrosuper, part, modelgrid) :
#
self.nbasmax = nbasmax
self.nbpt = hydrosuper.basin_count.shape[0]
+ nwbas = hydrosuper.basin_topoind.shape[1]
+ nbxmax_in = hydrosuper.inflow_grid.shape[1]
#
self.routing_area, 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 = \
- routing_interface.truncate(nbasmax, hydrosuper.num_largest, part.landcorelist, hydrosuper.basin_count, \
- hydrosuper.basin_notrun, hydrosuper.basin_area, hydrosuper.basin_cg, \
- hydrosuper.basin_topoind, hydrosuper.fetch_basin, hydrosuper.basin_id, \
- hydrosuper.basin_outcoor, hydrosuper.basin_type, hydrosuper.basin_flowdir, \
- hydrosuper.outflow_grid, hydrosuper.outflow_basin, \
- hydrosuper.inflow_number,hydrosuper.inflow_grid,hydrosuper.inflow_basin)
+ routing_interface.finish_truncate(nbpt = self.nbpt, nbxmax_in = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, cfrac = modelgrid.contfrac, basin_count = hydrosuper.basin_count, \
+ basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, basin_cg = hydrosuper.basin_cg, \
+ basin_topoind = hydrosuper.basin_topoind, 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, \
+ 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)
+
#
self.routing_fetch = finalfetch(part, self.routing_area, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.routing_fetch)
#
self.num_largest = routing_interface.finalrivclass(part.landcorelist, self.route_togrid, self.route_tobasin, self.routing_fetch, \
- hydrosuper.largest_rivarea)
+ hydrosuper.largest_rivarea)
#
return
#
diff --git a/RoutingPreProc.py b/RoutingPreProc.py
index c1489d7fea1dd160cdc88a923f42b77f29b3b61e..e358847eabf431ec330eb79c7fa16dc36434e647 100644
--- a/RoutingPreProc.py
+++ b/RoutingPreProc.py
@@ -11,6 +11,7 @@ import matplotlib as mpl
mpl.use('Agg')
from matplotlib.backends.backend_pdf import PdfPages
import matplotlib.pyplot as plt
+import time
#
# Gert the information from the configuration file.
#
@@ -31,7 +32,7 @@ import RPPtools as RPP
import HydroGrid as HG
import diagplots as DP
import Interface as IF
-from Truncate import truncation as TR
+from Truncate import trunc as TR
from spherical_geometry import vector
#
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
@@ -153,10 +154,9 @@ print("nbasmax : ", nbasmax)
INFO("hydrodata")
hydrodata = HG.HydroData(hydrogrid.ncfile, hydrogrid.box, sub_index)
-INFO("nitiatmgrid")
+INFO("initiatmgrid")
IF.initatmgrid(rank, nbcore, nbpt, modelgrid)
-
INFO("hoverlap")
hoverlap = IF.HydroOverlap(nbpt, nbvmax, sub_pts, sub_index, sub_area, sub_lon, sub_lat, part, modelgrid, hydrodata)
@@ -174,7 +174,7 @@ if rank ==0 :
if lonint[0] != lonint[1] and latint[0] != latint[1] :
DP.MAPPLOT("MapHydroGrid", lonint, latint, hoverlap, hoverlap.hierarchy_bx, modelgrid.polyll, title="Distance to ocean")
-hsuper = IF.HydroSuper(nbvmax, hydrodata, hoverlap)
+hsuper = IF.HydroSuper(nbvmax, hydrodata, hoverlap, nbasmax)
#
# Plot the hydrological supermesh
#
@@ -187,26 +187,29 @@ hsuper.linkup(hydrodata)
if rank ==0 :
if lonint[0] != lonint[1] and latint[0] != latint[1] :
DP.SUPERMESHPLOT("MapSuperGrid_Afterlinkup", lonint, latint, hoverlap, hsuper, modelgrid.polyll)
-
#
# Do some memory management and synchronize procs.
#
+
comm.Barrier()
del hoverlap
gc.collect()
print("=================== Compute fetch ====================")
+t = time.time()
hsuper.fetch(part)
+comm.Barrier()
+t1 = time.time()
+print("Time for fetch: {:0.2f} s.".format(t1-t))
+comm.Barrier()
+
hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), 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)
+print("Rank : {0} - Basin_count Before Truncate : ".format(part.rank), hsuper.basin_count)
+hs = TR(hsuper, part, comm, modelgrid, numop = 200)
+print("Rank : {0} - Basin_count After Truncate : ".format(part.rank), hsuper.basin_count)
+
+hgraph = IF.HydroGraph(nbasmax, hsuper, part, modelgrid)
+hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part)
IF.closeinterface(comm)
diff --git a/Truncate.py b/Truncate.py
index 6115dd52754f57edb076ff41c592a77dd5249aa6..bdcdb61632372a0fc883795bc64d98349214ace5 100755
--- a/Truncate.py
+++ b/Truncate.py
@@ -13,276 +13,467 @@ import time
import sys
import os
from inspect import currentframe, getframeinfo
-
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
-
import routing_interface
######################
-class truncation:
+class trunc:
"""
- Indices in Fortran index : need to adapt when using numpy
+ __init__ : realize all the truncate using the methods from the class
+
+ This class is organized in 3 type of methods
+ 1. METHOD FOR EACH STEP
+ Each of the 5 steps has its method that organize the (tokill, totakeover)
+ 2. METHOD to GET THE LIST OF CHANGES
+ When needed (for step 1 to 4), these are the method to find the list of changes
+ 3. BASIC METHODS
+ Basic methods used
"""
- def __init__(self, dfile, nbasmax):
- self.nc = NetCDFFile(dfile, "r")
-
- self.nbasmax = nbasmax
- self.num_largest = 10
- self.FillValue=np.nan
- self.IntFillValue=999999
- #
- self.proc_num = self.nc.variables["proc_num"][:]
- self.nbpt_proc = self.nc.variables["nbpt_proc"][:]
- #
- self.nbpt = self.nc.variables["nbpt_glo"][:]
- self.nbpt_num = int(np.max(self.nbpt))
- self.nj, self.ni = self.nbpt.shape
- #
- self.il_ji = [[k[0][0], k[1][0]] for k in [np.where(self.nbpt == i) for i in range(1,self.nbpt_num+1)]]
- #
- # Get i,j variables
- self.contfrac_e = self.nc.variables["contfrac"][:]
- self.area_e = self.nc.variables["area"][:]
- self.basin_count_e = self.nc.variables["basin_count"][:]
- self.basin_notrun_e = self.nc.variables["basin_notrun"][:]
- self.basin_area_e = self.nc.variables["basin_area"][:]
- self.basin_cg_lon_e = self.nc.variables["CG_lon"][:]
- self.basin_cg_lat_e = self.nc.variables["CG_lat"][:]
- self.basin_topoind_e = self.nc.variables["basin_topoind"][:]
- self.fetch_basin_e = self.nc.variables["fetch_basin"][:]
- self.basin_id_e = self.nc.variables["basin_id"][:]
- self.basin_outcoor_lon_e = self.nc.variables["outcoor_lon"][:]
- self.basin_outcoor_lat_e = self.nc.variables["outcoor_lat"][:]
- self.basin_type_e = self.nc.variables["basin_type"][:]
- self.basin_flowdir_e = self.nc.variables["basin_flowdir"][:]
- self.outflow_grid_e = self.nc.variables["HTUoutgrid"][:] # name
- self.outflow_basin_e = self.nc.variables["HTUoutbasin"][:] # name
- self.inflow_grid_e = self.nc.variables["HTUingrid"][:] # name
- self.inflow_basin_e = self.nc.variables["HTUinbas"][:] # name
- self.inflow_number_e = self.nc.variables["HTUinnum"][:]
- #
- # Conversion
+
+ def __init__(self, hydrosuper, part, comm,modelgrid, numop = 100):
+ """
+ __init__ : is the main program, it realize all the truncate
+ """
#
- self.gridarea = self.convert_toland(self.area_e).astype(np.int32)
- self.contfrac = self.convert_toland(self.contfrac_e).astype(np.int32)
- self.basin_count = self.convert_toland(self.basin_count_e).astype(np.int32)
- self.basin_notrun = self.convert_toland(self.basin_notrun_e).astype(np.int32)
- self.basin_area = self.convert_toland(self.basin_area_e).astype(np.float32)
- self.basin_cg_lon = self.convert_toland(self.basin_cg_lon_e).astype(np.float32)
- self.basin_cg_lat = self.convert_toland(self.basin_cg_lat_e).astype(np.float32)
- self.basin_cg = np.zeros((self.basin_cg_lon.shape[0], self.basin_cg_lon.shape[1], 2), order = 'F').astype(np.float32)
- self.basin_cg[:,:,0] = self.basin_cg_lat
- self.basin_cg[:,:,1] = self.basin_cg_lon
-
- self.basin_topoind = self.convert_toland(self.basin_topoind_e).astype(np.float32)
- self.fetch_basin = self.convert_toland(self.fetch_basin_e).astype(np.float32)
- self.basin_id = self.convert_toland(self.basin_id_e).astype(np.int32)
- self.basin_outcoor_lon = self.convert_toland(self.basin_outcoor_lon_e)
- self.basin_outcoor_lat = self.convert_toland(self.basin_outcoor_lat_e)
- self.basin_outcoor = np.zeros((self.basin_outcoor_lon.shape[0], self.basin_outcoor_lon.shape[1], 2), order = 'F').astype(np.float32)
- self.basin_outcoor[:,:,0] = self.basin_outcoor_lat
- self.basin_outcoor[:,:,1] = self.basin_outcoor_lon
-
- self.basin_type = self.convert_toland(self.basin_type_e).astype(np.float32)
- self.basin_flowdir = self.convert_toland(self.basin_flowdir_e).astype(np.int32)
- self.outflow_grid = self.convert_toland(self.outflow_grid_e).astype(np.int32)
- self.outflow_basin = self.convert_toland(self.outflow_basin_e).astype(np.int32)
+ # Initialize some variables
+ #
+ self.nbpt = int(hydrosuper.nbpt)
+ self.nbasmax = hydrosuper.nbasmax
+ self.gridarea = modelgrid.area
+ self.landcorelist = part.landcorelist
+ self.landcorelist_f = [l+1 for l in part.landcorelist]
+ self.landhalolist_f=[l+1 for l in range(hydrosuper.nbpt) if l not in self.landcorelist]
+ self.numop = numop
+
+
+ # Step 1 : coastal flow
+ self.max_ops_num = numop
+ i = 1
+ while self.max_ops_num == numop:
+ self.step1(hydrosuper, part)
+ hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
+
+ max_bas = part.domainmax(np.max(self.basin_count))
+ if part.rank == 0:
+ print("Step 1 iter. {0}, maxops: {1}->max_bas : {2}".format(i, self.max_ops_num, max_bas))
+ i+=1
+ hydrosuper.fetch(part)
+ comm.Barrier()
+
+
+ # Step 2 : HTUs going to the same neighbour
+ i = 1
+ self.max_ops_num_halo = self.numop; self.max_ops_num_core = self.numop
+ while (self.max_ops_num_halo == self.numop) or (self.max_ops_num_core == self.numop):
+
+ # Points with outflow in the core
+ if self.max_ops_num_core == self.numop:
+ self.step2(hydrosuper, part, "core")
+ hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
+
+ max_bas = part.domainmax(np.max(self.basin_count))
+ if part.rank == 0:
+ print("Step 2 iter. {0} (core) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_core, max_bas))
+ hydrosuper.fetch(part)
+ comm.Barrier()
+ # Points with outflow in the halo
+ if self.max_ops_num_halo == self.numop:
+ self.step2(hydrosuper, part, "halo")
+ hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
+
+ max_bas = part.domainmax(np.max(self.basin_count))
+ if part.rank == 0:
+ print("Step 2 iter. {0} (halo) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_halo, max_bas))
+ hydrosuper.fetch(part)
+ comm.Barrier()
+ i += 1
+
- self.inflow_grid = self.convert_toland(self.inflow_grid_e).astype(np.int32)
- self.inflow_basin = self.convert_toland(self.inflow_basin_e).astype(np.int32)
- self.inflow_number = self.convert_toland(self.inflow_number_e).astype(np.int32)
+ # Step 3 : HTUs with the same ID
+ i = 1
+ self.max_ops_num_halo = self.numop; self.max_ops_num_core = self.numop
+ while (self.max_ops_num_halo == self.numop) or (self.max_ops_num_core == self.numop):
+ # points with outflow in the core
+ if self.max_ops_num_core == self.numop:
+ self.step3(hydrosuper, part, "core")
+ hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
+
+ max_bas = part.domainmax(np.max(self.basin_count))
+ if part.rank == 0:
+ print("Step 3 iter. {0} (core) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_core, max_bas))
+ hydrosuper.fetch(part)
+ comm.Barrier()
+
+ # points with outflow in the halo
+ if self.max_ops_num_halo == numop:
+ self.step3(hydrosuper, part, "halo")
+ hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
+ max_bas = part.domainmax(np.max(self.basin_count))
+ if part.rank == 0:
+ print("Step 3 iter. {0} (halo) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_halo, max_bas))
+ hydrosuper.fetch(part)
+ comm.Barrier()
+ i += 1
+
+ # Step 4 : HTUs going to the same point
+ i = 1
+ self.max_ops_num = self.numop
+ while (self.max_ops_num == numop):
+ self.step4(hydrosuper, part)
+ hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
+ max_bas = part.domainmax(np.max(self.basin_count))
+ if part.rank == 0:
+ print("Step 4 iter. {0} (core) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num, max_bas))
+ hydrosuper.fetch(part)
+ comm.Barrier()
+ i += 1
+
- def convert_toland(self, array):
- if array.ndim == 2:
- dim = (self.nbpt_num)
- if array.ndim == 3:
- dim = (self.nbpt_num, array.shape[0])
- if array.ndim == 4:
- dim = (self.nbpt_num, array.shape[1], array.shape[1]) # (inflow, htu, lat, lon) -> nbpt, htu, inflow
-
- # Ne pas mettre à 0, mettre à IntFillValue
- out = np.zeros(dim, order = 'F')
-
- for il in range(self.nbpt_num):
- j,i = self.il_ji[il]
- if array.ndim == 2:
- out[il] = array[j,i]
- if array.ndim == 3:
- out[il,:] = array[:,j,i]
- if array.ndim == 4:
- out[il,:,:array.shape[0]] = np.swapaxes(array[:,:,j,i],0,1)
- return out
-
- ###############
-
- def convert_toij(self, array):
- if array.ndim == 1:
- dim = (self.nj,self.ni)
- if array.ndim == 2:
- dim = (array.shape[1], self.nj, self.ni)
- if array.ndim == 3:
- dim = (array.shape[2], array.shape[1], self.nj, self.ni) # (grid, htu, inflow) -> inflow, htu, lat, lon
-
- # Ne pas mettre à 0, mettre à IntFillValue
- dtype = array.dtype
- out = np.zeros(dim, order = 'F', dtype = dtype)
- if dtype == np.int32:
- out[:] = self.IntFillValue
- elif dtype == np.float32:
- out[:] = self.FillValue
-
- for il in range(self.nbpt_num):
- j,i = self.il_ji[il]
- if array.ndim == 1:
- out[j,i] = array[il]
- if array.ndim == 2:
- out[:,j,i] = array[il,:]
- if array.ndim == 3:
- out[:,:,j,i] = np.swapaxes(array[il,:,:], 0, 1) # (inflow, htu, lat, lon)
- return out
-
- ###############
-
- def truncate(self):
- nwbas = self.basin_topoind.shape[1]
- nbxmax_in = self.inflow_grid.shape[1]
- inf = np.copy(self.inflow_number)
- self.routing_area, 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 = \
- routing_interface.truncate(nbpt = self.nbpt_num, nbxmax_in = nbxmax_in, nbasmax = self.nbasmax, nwbas = nwbas, num_largest = self.num_largest, gridarea = self.gridarea, cfrac = self.contfrac, basin_count = self.basin_count, \
- basin_notrun = self.basin_notrun, basin_area = self.basin_area, basin_cg = self.basin_cg, \
- basin_topoind = self.basin_topoind, fetch_basin = self.fetch_bis, 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)
-
- # Adjust inflows -> add this part in fortran
- for il in range(self.nbpt_num):
- for ib in range(self.nbasmax):
- self.inflow_grid[il,ib,self.inflow_number[il,ib]:] = self.IntFillValue
- self.inflow_basin[il,ib,self.inflow_number[il,ib]:] = self.IntFillValue
-
- self.inflow_grid[il,self.nbasmax:,:] = self.IntFillValue
- self.inflow_basin[il,self.nbasmax:,:] = self.IntFillValue
- self.inflow_number[il,self.nbasmax:] = self.IntFillValue
-
-
- def finalrivclass(self, largest_rivarea):
-
- A = np.where((self.route_tobasin > self.nbasmax) & (self.routing_fetch > largest_rivarea))
- A = [[A[0][l], A[1][l]] for l in range(len(A[0]))]
- num_largest = len(A)
- for il, ib in A:
- self.route_tobasin[il,ib] = self.nbasmax +3
+ # Step 5 : Brutal method
+ self.step5(hydrosuper, part)
+ hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
+ max_bas = part.domainmax(np.max(self.basin_count))
+ if part.rank == 0:
+ print("Step 5, maxops: {0}->max_bas : {1}".format(self.numopmax, max_bas))
+
+ hydrosuper.fetch(part)
+ comm.Barrier()
- print('NUMBER OF RIVERS :', len(A))
+ ########################
+ #
+ # METHODs FOR EACH STEP
+ #
- def get_downstream(self, L, ig, ib):
+ def step1(self, hydrosuper, part):
+ """
+ Step 1 : coastal outflow or river outflow
"""
- Useful to add / remove fetch
+ self.load_data(hydrosuper)
+ self.tokill = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.totakeover = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.numops = np.zeros((self.nbpt), dtype = np.int32)
+
+ for k in range(len(self.landcorelist)):
+ changes = self.get_changes_step1(k)
+ if len(changes)>0:
+ for ch in changes:
+ orig_ops = self.numops[self.landcorelist[k]]
+ if orig_ops < self.numop:
+ self.add_changes(ch, k)
+ self.reindex_changes(k)
+
+ self.max_ops_num = part.domainmax(max(self.numops))
+ part.landsendtohalo(self.numops, order='F')
+ part.landsendtohalo(self.tokill, order='F')
+ part.landsendtohalo(self.totakeover, order='F')
+
+ return
+ #
+ def step2(self, hydrosuper, part, neigh_domain):
+ """
+ Step 2 : neighbour
"""
- L.append([ig,ib])
- if self.outflow_grid[ig-1,ib-1]>0:
- self.get_downstream(L, int(self.outflow_grid[ig-1, ib-1]), int(self.outflow_basin[ig-1, ib-1]))
- else:
- L.append(int(self.outflow_grid[ig-1, ib-1]))#,int(self.outbasin[ib-1,j,i])])
-
- def get_fetch(self):
-
- L_area = []
- for ig, J in enumerate(self.il_ji):
- j,i = J
- L_area.append([self.basin_area[ig,k] for k in range(self.basin_count[ig])])
+ self.load_data(hydrosuper)
+ self.tokill = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.totakeover = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.numops = np.zeros((self.nbpt), dtype = np.int32)
- L_fetch_calc = []
- for ig, J in enumerate(self.il_ji):
- j,i = J
- L_fetch_calc.append([0 for k in range(self.basin_count[ig])])
-
- for ig in range(1, len(self.basin_count)+1): # plus propre nbpt_num
- for ib in range(1, self.basin_count[ig-1]+1):
- fe = L_area[ig-1][ib-1]
- L = []
- self.get_downstream(L, ig, ib)
- for igf,ibf in L[:-1]:
- L_fetch_calc[igf-1][ibf-1] += fe
-
- self.fetch_bis = np.zeros(self.fetch_basin.shape, order = 'F').astype(np.float32)
- for ig in range(1, len(self.basin_count)+1):# plus propre nbpt_num
- self.fetch_bis[ig-1,:self.basin_count[ig-1]] = np.array(L_fetch_calc[ig-1])
-
- return
+ for k in range(len(self.landcorelist)):
+ changes = self.get_changes_step2(k, neigh_domain)
+ if len(changes)>0:
+ for ch in changes:
+ orig_ops = self.numops[self.landcorelist[k]]
+ if orig_ops < self.numop:
+ self.add_changes(ch, k)
+
+ self.reindex_changes(k)
+
+ if neigh_domain == "core":
+ self.max_ops_num_core = part.domainmax(max(self.numops))
+ elif neigh_domain == "halo":
+ self.max_ops_num_halo = part.domainmax(max(self.numops))
- ####################################################
+ part.landsendtohalo(self.numops, order='F')
+ part.landsendtohalo(self.tokill, order='F')
+ part.landsendtohalo(self.totakeover, order='F')
+
+ return
#
- # Generation of the NetCDF output
+ def step3(self, hydrosuper, part, neigh_domain):
+ """
+ Step 3 : Same basin_id
+ """
+ self.load_data(hydrosuper)
+ self.tokill = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.totakeover = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.numops = np.zeros((self.nbpt), dtype = np.int32)
+
+ for k in range(len(self.landcorelist)):
+ changes = self.get_changes_step3(k, neigh_domain)
+ if len(changes)>0:
+ for ch in changes:
+ orig_ops = self.numops[self.landcorelist[k]]
+ if orig_ops < self.numop:
+ self.add_changes(ch, k)
+ self.reindex_changes(k)
+
+ if neigh_domain == "core":
+ self.max_ops_num_core = part.domainmax(max(self.numops))
+ elif neigh_domain == "halo":
+ self.max_ops_num_halo = part.domainmax(max(self.numops))
+
+ part.landsendtohalo(self.numops, order='F')
+ part.landsendtohalo(self.tokill, order='F')
+ part.landsendtohalo(self.totakeover, order='F')
+
+ return
#
+ def step4(self, hydrosuper, part):
+ """
+ Step 4 : htus inside point
+ """
+ self.load_data(hydrosuper)
+ self.tokill = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.totakeover = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.numops = np.zeros((self.nbpt), dtype = np.int32)
+
+ for k in range(len(self.landcorelist)):
+ tokk, totakk = self.get_changes_step4(k)
+
+ if len(tokk)>0:
+ ops = min(self.numop,len(tokk))
+
+ self.tokill[self.landcorelist[k],:ops] = tokk[:ops]
+ self.totakeover [self.landcorelist[k],:ops] = totakk[:ops]
+ self.numops[self.landcorelist[k]] = ops
+
+ self.reindex_changes(k)
+
+ self.max_ops_num = part.domainmax(max(self.numops))
+ part.landsendtohalo(self.numops, order='F')
+ part.landsendtohalo(self.tokill, order='F')
+ part.landsendtohalo(self.totakeover, order='F')
- def addcoordinates(self, outnf) :
- #
- for varn in ["longitude", "latitude", "lon_bnd", "land", "area", "proc_num", "nbpt_proc", "nbpt_glo"]:
- ovar = self.nc.variables[varn]
- nvar = outnf.createVariable(varn, ovar.dtype, ovar.dimensions)
- nvar[:] = ovar[:]
- for attrn in ovar.ncattrs():
- attrv = ovar.getncattr(attrn)
- nvar.setncattr(attrn,attrv)
- outnf.sync()
return
-
- def add_variable(self, outnf, data, NCFillValue, name, vtyp, dim, title, unit, type_data):
- dtype = data.dtype
- data_ij = self.convert_toij(data)
- data_ij = data_ij.astype(vtyp)
-
- if dtype == np.int32:
- data_ij[data_ij >= self.IntFillValue] = NCFillValue
- elif dtype == np.float32:
- data_ij[np.isnan(data_ij)] = NCFillValue
-
- ncvar = outnf.createVariable(name, vtyp, dim, fill_value=NCFillValue)
- ncvar.title = title
- ncvar.units = unit
- ncvar[:] = data_ij[:]
-
- def dumpnetcdf(self, filename) :
- #
- NCFillValue=1.0e20
- insize = 100
- vtyp=np.float64
- cornerind=[0,2,4,6]
- nbcorners = len(cornerind)
- #
- outnf=NetCDFFile(filename, 'w', format='NETCDF4_CLASSIC')
- # Dimensions
- outnf.createDimension('x', self.ni)
- outnf.createDimension('y', self.nj)
- outnf.createDimension('in', insize)
- outnf.createDimension('land', self.nbpt_num)
- outnf.createDimension('htu', self.nbasmax)
- outnf.createDimension('bnd', nbcorners)
-
- self.add_variable(outnf, self.route_togrid, NCFillValue, "routetogrid", vtyp, ('htu','y','x'), "Grid into which the basin flows", "-", "int")
- self.add_variable(outnf, self.route_tobasin, NCFillValue, "routetobasin", vtyp, ('htu','y','x'), "Basin in to which the water goes", "-", "int")
- self.add_variable(outnf, self.global_basinid, NCFillValue, "basinid", vtyp, ('htu','y','x'), "ID of basin", "-", "int")
- self.add_variable(outnf, self.route_nbintobas, NCFillValue, "routenbintobas", vtyp, ('y','x'), "Number of basin into current one", "-", "int")
- self.add_variable(outnf, self.origin_nbintobas, NCFillValue, "originnbintobas", vtyp, ('y','x'), "Number of sub-grid basin into current one before truncation", "-", "int")
- self.add_variable(outnf, self.route_outlet[:,:,0], NCFillValue, "outletlat", vtyp, ('htu','y','x'), "Latitude of Outlet","degrees north", "float")
- self.add_variable(outnf, self.route_outlet[:,:,1], NCFillValue, "outletlon", vtyp, ('htu','y','x'), "Longitude of outlet","degrees east", "float")
- self.add_variable(outnf, self.route_type[:,:], NCFillValue, "outlettype", vtyp, ('htu','y','x'), "Type of outlet","code", "int")
- self.add_variable(outnf, self.topo_resid[:,:], NCFillValue, "topoindex", vtyp, ('htu','y','x'), "Topographic index of the retention time","m", "float")
- self.add_variable(outnf, self.routing_cg[:,:,1], NCFillValue, "CG_lon", vtyp, ('htu','y','x'), "Longitude of centre of gravity of HTU","degrees east", "float")
- self.add_variable(outnf, self.routing_cg[:,:,0], NCFillValue, "CG_lat", vtyp, ('htu','y','x'), "Latitude of centre of gravity of HTU","degrees east", "float")
- self.add_variable(outnf, self.routing_fetch, NCFillValue, "fetch", vtyp, ('htu','y','x'), "Fetch contributing to each HTU","m^2", "float")
-
- self.add_variable(outnf, self.inflow_number[:,:self.nbasmax], NCFillValue, "inflow_number", vtyp, ('htu','y','x'), "Fetch contributing to each HTU","m^2", "int")
- self.add_variable(outnf, self.inflow_grid[:,:self.nbasmax,:100], NCFillValue, "inflow_grid", vtyp, ('in', 'htu','y','x'), "Grid of HTUs flowing into the basin","-", "int")
- self.add_variable(outnf, self.inflow_basin[:,:self.nbasmax,:100], NCFillValue, "inflow_basin", vtyp, ('in', 'htu','y','x'), "Basin of HTUs flowing into the basin","-", "int")
- # Revoir inflows
-
- outnf.close()
+ #
+ def step5(self, hydrosuper, part):
+ """
+ # Step 5 : htus inside point
+ """
+ # Limit of the area of basin to delete by step 5 (in %)
+ limit_step5 = 5
+ self.load_data(hydrosuper)
+
+ # Number of basin to delete in each point of the core
+ num2del = []
+ for k, k_land in enumerate(self.landcorelist):
+ num2del.append(self.basin_count[k_land]-self.nbasmax)
+ # Get the maximum over the whole domain
+ self.numopmax = part.domainmax(np.max(num2del))
+
+ if self.numopmax>self.numop:
+ print("Error - cannot finish to resolve truncate by step 5")
+ print("Try with an higher nbasmax")
+ sys.exit()
+
+
+ self.tokill = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.totakeover = np.zeros((self.nbpt,self.numop), dtype = np.int32)
+ self.numops = np.zeros((self.nbpt), dtype = np.int32)
+
+ for k, k_land in enumerate(self.landcorelist):
+ if num2del[k] > 0:
+ B = np.argsort(self.fetch[k])[:num2del[k]]
+ tokk = [k+1 for k in B]
+ totakk = [np.argsort(self.fetch[k])[-1] + 1] * num2del[k]
+
+ area_basin_to_truncate = np.sum([self.basin_area[k][l] for l in B])
+
+ if (area_basin_to_truncate/np.sum(self.basin_area[k]))*100 < limit_step5:
+ ops = len(tokk)
+ self.tokill[k_land,:ops] = tokk
+ self.totakeover[k_land,:ops] = totakk
+ self.numops[k_land] = ops
+
+ self.reindex_changes(k)
+
+ else:
+ print("Error - cannot finish to resolve truncate by step 5")
+ print("Try with an higher nbasmax")
+ sys.exit()
+
+ part.landsendtohalo(self.numops, order='F')
+ part.landsendtohalo(self.tokill, order='F')
+ part.landsendtohalo(self.totakeover, order='F')
+
return
+
+ ########################
+ #
+ # METHOD to GET THE LIST OF CHANGES
+ #
+
+ def get_changes_step1(self, npt):
+ changes = []
+
+ # outflow_grid = -2
+ A = np.where(self.outflow_grid[npt] == -2)[0]
+ if len(A)>1:
+ B = self.fetch[npt][A]
+ changes.append([A[k]+1 for k in np.argsort(B)])
+ # outflow_grid = -3
+ A = np.where(self.outflow_grid[npt] == -3)[0]
+ if len(A)>1:
+ B = self.fetch[npt][A]
+ changes.append([A[k]+1 for k in np.argsort(B)])
+
+ return changes
+
+ def get_changes_step2(self, npt, neigh_domain):
+ """
+ neigh_domain = "halo" or "core"
+ """
+ # List of neighbour which are outflow
+ outgrid = self.outflow_grid[npt]
+ outgrid = ma.masked_where(outgrid == self.landcorelist[npt], outgrid)
+ outgrid = ma.masked_where(outgrid <= 0, outgrid)
+ outgrid = ma.unique(outgrid)
+
+ # List of neighbour in core/halo
+ core = []; halo = []
+ for g in outgrid:
+ if g in self.landhalolist_f:
+ halo.append(g)
+ else:
+ core.append(g)
+
+ # Establishing the list of changes
+ changes = []
+
+ if neigh_domain == "core":
+ for k in core:
+ A = np.where(self.outflow_grid[npt] == k)[0]
+ if len(A)>1:
+ B = self.fetch[npt][A]
+ changes.append([A[k]+1 for k in np.argsort(B)])
+
+ elif neigh_domain == "halo":
+ for k in halo:
+ A = np.where(self.outflow_grid[npt] == k)[0]
+ if len(A)>1:
+ B = self.fetch[npt][A]
+ changes.append([A[k]+1 for k in np.argsort(B)])
+
+ return changes
+
+
+ def get_changes_step3(self, npt, neigh_domain):
+ """
+ neigh_domain = "halo" or "core"
+ """
+ # List of basin ID included in the grid point
+ basid = self.basin_id[npt]
+ basid = ma.unique(basid)
+
+ core = []; halo = []
+
+ # For each basin id, list of point with outflow in the core / halo
+ for bid in basid:
+ htus = np.where(basid == bid)[0]
+ outgrid = [self.outflow_grid[npt][htu] for htu in htus]
+ c = []; h = []
+
+ for g,htu in zip(outgrid,htus):
+ if g in self.landhalolist_f:
+ h.append(htu)
+ elif g in self.landcorelist_f:
+ c.append(htu)
+ if len(h)>1:
+ halo.append(h)
+ if len(c)>1:
+ core.append(c)
+
+ # Establishing the list of changes
+ changes = []
+
+ if neigh_domain == "core":
+ for c in core:
+ B = self.fetch[npt][c]
+ changes.append([c[k]+1 for k in np.argsort(B)])
+
+ elif neigh_domain == "halo":
+ for h in halo:
+ B = self.fetch[npt][h]
+ changes.append([h[k]+1 for k in np.argsort(B)])
+
+ return changes
+
+ def get_changes_step4(self, npt):
+ tok = []
+ totak = []
+
+ # List of HTUs flowing in the same grid
+ outgrid = self.outflow_grid[npt]
+ htus = np.where(outgrid == self.landcorelist[npt])[0]
+
+ # Establishing the lists of changes
+ if len(htus)>0:
+ B = self.fetch[npt][htus]
+ tok = [htus[k]+1 for k in np.argsort(B)]
+ totak = [self.outflow_basin[npt][htus[k]] for k in np.argsort(B)]
+
+ return tok, totak
+
+
+ ########################
+ #
+ # BASIC METHODS
+ #
+ def load_data(self, hydrosuper):
+ """
+ Load the core data from hydrosuper
+ """
+ self.basin_count = hydrosuper.basin_count
+ #
+ self.outflow_grid = [hydrosuper.outflow_grid[k,:self.basin_count[k]] for k in self.landcorelist]
+ self.outflow_basin = [hydrosuper.outflow_basin[k,:self.basin_count[k]] for k in self.landcorelist]
+ self.basin_id = [hydrosuper.basin_id[k,:self.basin_count[k]] for k in self.landcorelist]
+ self.basin_area = [hydrosuper.basin_area[k,:self.basin_count[k]] for k in self.landcorelist]
+ self.fetch = [hydrosuper.fetch_basin[k,:self.basin_count[k]] for k in self.landcorelist]
+
+ def add_changes(self, HTUs_to_truncate, k):
+ """
+ Add the changes to (tokill, totakeover, numops) from HTUs_to_truncate
+ HTUs_to_truncate : List of HTUs to truncate, ordered by fetch
+ """
+ orig_ops = self.numops[self.landcorelist[k]]
+ poss_ops = len(HTUs_to_truncate)-1
+ ops = min(poss_ops,self.numop-orig_ops)
+
+ self.tokill[self.landcorelist[k],orig_ops:orig_ops+ops] = HTUs_to_truncate[:ops]
+ self.totakeover[self.landcorelist[k],orig_ops:orig_ops+ops] = HTUs_to_truncate[1:ops+1]
+ self.numops[self.landcorelist[k]] += ops
+
+ def reindex_changes(self, k):
+ """
+ Reindex the (tokill, totakeover) to anticipate the changes made before
+ """
+ ops = self.numops[self.landcorelist[k]]
+ for l in range(ops-1):
+ tok = self.tokill[self.landcorelist[k],l]
+ tokill_up = self.tokill[self.landcorelist[k],l+1:]
+ totakeover_up = self.totakeover[self.landcorelist[k],l+1:]
+
+ tokill_up[tokill_up > tok] -= 1
+ totakeover_up[totakeover_up > tok] -= 1
+
+ self.tokill[self.landcorelist[k],l+1:] = tokill_up
+ self.totakeover[self.landcorelist[k],l+1:] = totakeover_up
+
+
+
+
diff --git a/tests/Iberia_n48/run.def b/tests/Iberia_n48/run.def
index 3e11dd5898bc1d3bc3a4322565ed469a7861dc97..e252d30a3064c468aa45316579c5e5125ee09113 100644
--- a/tests/Iberia_n48/run.def
+++ b/tests/Iberia_n48/run.def
@@ -15,7 +15,7 @@ Documentation = true
#
# Configuration for the graph to be generated
#
-nbasmax = 35
+nbasmax = 100
#
# Output
#