From 8b2082f5e7a3182769fe9e4860f5f526ecf10899 Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Tue, 21 Apr 2020 18:07:03 +0200
Subject: [PATCH] Add the orography and the floodplains
---
F90subroutines/routing_interface.f90 | 77 ++++++++++++-------
F90subroutines/routing_reg.f90 | 59 ++++++++++++---
HydroGrid.py | 16 +++-
Interface.py | 109 ++++++++++++++++-----------
tests/Iberia_n48/run.def | 2 +-
tests/Iberia_n48_regular/run.def | 2 +-
6 files changed, 180 insertions(+), 85 deletions(-)
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 5f9593b..d2ac028 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -44,8 +44,8 @@ SUBROUTINE closeinterface()
END SUBROUTINE closeinterface
SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
- & lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy, nbi, nbj, area_bx, trip_bx, basin_bx, &
- & topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx)
+ & 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
@@ -71,6 +71,8 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
REAL, INTENT(inout) :: topoindex(nbpt,nbvmax_in) !! Topographic index of the residence time (m)
REAL, INTENT(inout) :: hierarchy(nbpt,nbvmax_in) !! data on the fine grid
REAL, INTENT(inout) :: fac(nbpt,nbvmax_in) !! data on the fine grid
+ REAL, INTENT(inout) :: orog(nbpt,nbvmax_in) !! data on the fine grid
+ REAL, INTENT(inout) :: floodp(nbpt,nbvmax_in) !! data on the fine grid
!
!! OUTPUT VARIABLES
INTEGER, INTENT(out) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless)
@@ -83,6 +85,8 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
REAL, INTENT(out) :: topoind_bx(nbpt,nbxmax_in,nbxmax_in) !! Topographic index of the residence time for each of the smaller boxes (m)
INTEGER, INTENT(out) :: trip_bx(nbpt,nbxmax_in,nbxmax_in) !! The trip field for each of the smaller boxes (unitless)
INTEGER, INTENT(out) :: basin_bx(nbpt,nbxmax_in,nbxmax_in) !!
+ REAL, INTENT(out) :: orog_bx(nbpt,nbxmax_in,nbxmax_in) !! Orography (m)
+ REAL, INTENT(out) :: floodp_bx(nbpt,nbxmax_in,nbxmax_in) !! floodplains (m^2)
!
INTEGER :: ii, ib
REAL :: resolution(nbpt,2)
@@ -101,8 +105,9 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
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, &
- & nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), basin_bx(ib,:,:), topoind_bx(ib,:,:), &
- & fac_bx(ib,:,:), hierarchy_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), lshead_bx(ib,:,:))
+ & orog, floodp, nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), basin_bx(ib,:,:), &
+ & topoind_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:),orog_bx(ib,:,:),floodp_bx(ib,:,:),&
+ & lon_bx(ib,:,:), lat_bx(ib,:,:), lshead_bx(ib,:,:))
ENDDO
!
END SUBROUTINE gethydrogrid
@@ -167,9 +172,11 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
END SUBROUTINE findbasins
-SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, fac_bx, topoind_bx, &
+SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, &
+ & fac_bx, topoind_bx, &
& min_topoind, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_pts, basin_bxout, &
& basin_bbout, lshead, coast_pts, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, &
+ & basin_orog, basin_floodp, &
& basin_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
& outflow_grid, outflow_basin, nbcoastal, coastal_basin)
!
@@ -188,6 +195,8 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: lat_bx !! Latitude of each small box in the grid box
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: trip_bx !! The trip field for each of the smaller boxes (unitless)
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: hierarchy_bx !! Level in the basin of the point
+ REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: orog_bx !! Orography (m)
+ REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: floodp_bx !! Floodplains area (m^2)
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: fac_bx !! Level in the basin of the point
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: topoind_bx !! Topographic index of the residence time for each of the smaller boxes (m)
REAL(r_std), INTENT (in) :: min_topoind !! The current minimum of topographic index (m)
@@ -211,6 +220,8 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_area !!
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas,2) :: basin_cg !! Centre of gravity of the HTU
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_hierarchy !!
+ REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_orog !!
+ REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_floodp !!
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.
@@ -228,10 +239,11 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
!!
DO ib=1,nbpt
CALL routing_reg_globalize(nbpt, ib, neighbours, area_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), trip_bx(ib,:,:), &
- & hierarchy_bx(ib,:,:), fac_bx(ib,:,:), &
+ & hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), fac_bx(ib,:,:), &
& topoind_bx(ib,:,:), min_topoind, nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), &
& basin_sz(ib,:), basin_pts(ib,:,:,:), basin_bxout(ib,:), basin_bbout(ib,:), lshead(ib,:), coast_pts(ib,:), nwbas, &
- & basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, basin_fac, basin_topoind, basin_id, basin_coor, &
+ & basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy,&
+ & basin_orog, basin_floodp, basin_fac, basin_topoind, basin_id, basin_coor, &
& basin_type, basin_flowdir, basin_lshead, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
ENDDO
!
@@ -417,10 +429,12 @@ SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_
END SUBROUTINE rivclassification
-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)
+SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, &
+ & basin_notrun, basin_area, basin_orog, basin_floodp, basin_cg, basin_topoind, fetch_basin, &
+ & basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
+ & inflow_number, inflow_grid, inflow_basin, routing_area, routing_orog, routing_floodp, &
+ & routing_cg, topo_resid, route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
+ & global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
!
USE ioipsl
USE grid
@@ -444,6 +458,8 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
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), INTENT(inout) :: basin_orog !!
+ 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) :: fetch_basin !!
@@ -459,6 +475,8 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
! Output variables
!
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_area !! Surface of basin (m^2)
+ REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_orog !! 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_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)
@@ -474,11 +492,14 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
!
!
- 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, &
+ CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, &
+ & basin_count, basin_notrun, basin_area, basin_orog, basin_floodp, basin_cg, &
+ & basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
routing_area(:,:) = routing_area_glo(:,:)
+ routing_orog(:,:) = routing_orog_glo(:,:)
+ routing_floodp(:,:) = routing_floodp_glo(:,:)
routing_cg(:,:,:) = routing_cg_glo(:,:,:)
topo_resid(:,:) = topo_resid_glo(:,:)
route_nbbasin(:) = route_count_glo(:)
@@ -545,6 +566,7 @@ END SUBROUTINE finish_inflows
SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numops, basin_count, basin_area, &
+ & basin_orog, basin_floodp, &
& basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
!
@@ -568,6 +590,8 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
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), INTENT(inout) :: basin_orog !!
+ 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) :: fetch_basin !!
@@ -589,7 +613,7 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
IF (basin_count(ib) > nbasmax) THEN
tok = tokill(ib,op)
totak = totakeover(ib,op)
- IF (tok .GT. 0) THEN
+ IF (tok .GT. 0) THEN
WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak
! Test if tokill is downstream of totakeover (avoid loop)
igrif = outflow_grid(ib,totak)
@@ -603,11 +627,12 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
tok = it
ELSE
it = outflow_grid(igrif,ibasf)
- ibasf = outflow_basin(igrif,ibasf)
+ ibasf = outflow_basin(igrif,ibasf)
igrif = it
END IF
END DO
- CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, basin_area, basin_cg, basin_topoind,&
+ CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, basin_area, &
+ & basin_orog, basin_floodp, basin_cg, basin_topoind,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
& inflow_grid, inflow_basin)
END IF
@@ -615,8 +640,6 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
END DO
END DO
- ! REPRENDRE LA FIN DE TRUNCATE
-
END SUBROUTINE killbas
@@ -633,7 +656,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
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 !!
@@ -641,13 +664,13 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
flag(:,:) = 0
-
+
WRITE(numout,*) "Checking routing"
test = 0
DO ig=1,nbpt
basnum = basin_count(ig)
- DO ib=1,basnum
+ DO ib=1,basnum
IF (flag(ig,ib) .EQ. 0) THEN
flag(ig,ib) = -1
igrif = outflow_grid(ig,ib)
@@ -655,7 +678,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
!
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)
+ it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
@@ -668,7 +691,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -99))
flag(igrif, ibasf) = -99
- it = outflow_grid(igrif,ibasf)
+ it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
@@ -681,7 +704,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
ibasf = outflow_basin(ig,ib)
DO WHILE (flag(igrif,ibasf) .EQ. -1)
flag(igrif, ibasf) = -99
- it = outflow_grid(igrif,ibasf)
+ it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
@@ -720,11 +743,11 @@ SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, outflow_grid, outflow_basin, bas
DO ig=1,nbpt
basnum = basin_count(ig)
- DO ib=1,basnum
+ DO ib=1,basnum
igrif = outflow_grid(ig,ib)
ibasf = outflow_basin(ig,ib)
-
- IF (igrif .GT. 0) THEN
+
+ 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
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 8fe2298..f38a21f 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -15,6 +15,8 @@ MODULE routing_reg
INTEGER(i_std), SAVE :: nbasmax_save
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
+ REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_orog_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_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
@@ -73,7 +75,8 @@ CONTAINS
SUBROUTINE 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, &
- & nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx)
+ & 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)
!
IMPLICIT NONE
!
@@ -98,6 +101,8 @@ CONTAINS
REAL(r_std), INTENT(inout) :: topoindex(nbpt,nbvmax) !! Topographic index of the residence time (m)
REAL(r_std), INTENT(inout) :: hierarchy(nbpt,nbvmax) !! data on the fine grid
REAL(r_std), INTENT(inout) :: fac(nbpt,nbvmax) !! data on the fine grid
+ REAL(r_std), INTENT(inout) :: orog(nbpt,nbvmax) !! data on the fine grid
+ REAL(r_std), INTENT(inout) :: floodp(nbpt,nbvmax) !! data on the fine grid
!
!! OUTPUT VARIABLES
INTEGER(i_std), INTENT(out) :: nbi, nbj !! Number of point in x and y within the grid (unitless)
@@ -110,6 +115,8 @@ CONTAINS
REAL(r_std), INTENT(out) :: topoind_bx(nbxmax,nbxmax) !! Topographic index of the residence time for each of the smaller boxes (m)
INTEGER(i_std), INTENT(out) :: trip_bx(nbxmax,nbxmax) !! The trip field for each of the smaller boxes (unitless)
INTEGER(i_std), INTENT(out) :: basin_bx(nbxmax,nbxmax) !!
+ REAL(i_std), INTENT(out) :: orog_bx(nbxmax,nbxmax) !!
+ REAL(i_std), INTENT(out) :: floodp_bx(nbxmax,nbxmax) !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc !! Indices (unitless)
@@ -152,6 +159,8 @@ CONTAINS
lon_bx(:,:) = undef_sechiba
lat_bx(:,:) = undef_sechiba
lshead_bx(:,:) = undef_sechiba
+ orog_bx(:,:) = undef_sechiba
+ floodp_bx(:,:) = undef_sechiba
cenlon = zero
cenlat = zero
!
@@ -193,6 +202,8 @@ CONTAINS
topoind_bx(iloc, jloc) = topoindex(ib, ip)
hierarchy_bx(iloc, jloc) = hierarchy(ib, ip)
fac_bx(iloc, jloc) = fac(ib, ip)
+ orog_bx(iloc, jloc) = orog(ib, ip)
+ floodp_bx(iloc, jloc) = floodp(ib, ip) * sub_area(ib, ip)
lon_bx(iloc, jloc) = lon_rel(ib, ip)
lat_bx(iloc, jloc) = lat_rel(ib, ip)
cenlon = cenlon + lon_rel(ib, ip)/sub_pts(ib)
@@ -1635,9 +1646,11 @@ CONTAINS
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, fac_bx, topoind_bx, &
+SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, &
+ & orog_bx, floodp_bx, fac_bx, topoind_bx, &
& min_topoind, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_pts, basin_bxout, &
& basin_bbout, lshead, coast_pts, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy, &
+ & basin_orog, basin_floodp, &
& basin_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
& outflow_grid, outflow_basin, nbcoastal, coastal_basin)
!
@@ -1654,6 +1667,8 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: lat_bx !! Latitude of each small box in the grid box
INTEGER(i_std), DIMENSION(nbxmax,nbxmax) :: trip_bx !! The trip field for each of the smaller boxes (unitless)
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: hierarchy_bx !! Level in the basin of the point
+ REAL(r_std), DIMENSION(nbxmax,nbxmax) :: orog_bx !! Level in the basin of the point
+ REAL(r_std), DIMENSION(nbxmax,nbxmax) :: floodp_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: fac_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: topoind_bx !! Topographic index of the residence time for each of the smaller boxes (m)
REAL(r_std) :: min_topoind !! The current minimum of topographic index (m)
@@ -1677,6 +1692,8 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area !!
REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_cg !! Centre of gravity of the HTU in latitude, longitude
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_hierarchy !!
+ REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog !!
+ REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_floodp !!
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.
@@ -1755,6 +1772,8 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
basin_cg(ib,ij,:) = zero
basin_hierarchy(ib,ij) = zero
basin_fac(ib,ij) = zero
+ basin_orog(ib,ij) = zero
+ basin_floodp(ib,ij) = zero
!
SELECT CASE (hierar_method)
!
@@ -1767,6 +1786,10 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
DO iz=1,basin_sz(ij)
!
basin_area(ib,ij) = basin_area(ib,ij) + area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
+ ! Numerator for mean orography, then need to be divided by the HTUs area
+ basin_orog(ib,ij) = basin_orog(ib,ij) + orog_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) &
+ & * area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
+ basin_floodp(ib,ij) = basin_floodp(ib,ij) + floodp_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
!
! Compute centre of gravity
!
@@ -1813,6 +1836,9 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
!
ENDDO
!
+ ! Finish the calculation of mean orography
+ basin_orog(ib,ij) = basin_orog(ib,ij)/basin_area(ib,ij)
+ !
SELECT CASE (hierar_method)
!
CASE("MEAN")
@@ -2054,7 +2080,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,1) = solved(sp,1) + 1
ENDIF
- WRITE(numout,*) sp,sb,"flow in the same grid, basin:",bop
+ WRITE(numout,*) sp,sb,"flow in the same grid, basin:",bop
!
ELSE IF ( outflow_grid(sp,sb) .EQ. -3 ) THEN
! Return flow
@@ -2296,20 +2322,20 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
& inflow_number, inflow_grid, inflow_basin)
! It is possible that the lowest hierarchy is in two grid cells
! In this case, if we have an error for this routing_updateflow
- ! We go to next step and make it a coastal flow
+ ! We go to next step and make it a coastal flow
IF (outflow_basin(sp, sb) < undef_int) THEN
found = 1
WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points"
ELSE
WRITE(numout,*) "Lowest hierarchy may be in two different grid points"
END IF
- WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb)
+ WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb)
WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop)
END IF
END IF
END DO
END IF
- ENDIF
+ ENDIF
!
! Last option : coastal outflow
IF (found .EQ. 0) THEN
@@ -2554,8 +2580,9 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!! \n
!_ ================================================================================================================================
-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, &
+SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, &
+ & basin_count, basin_notrun, basin_area, basin_orog, basin_floodp, &
+ & basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
!
IMPLICIT NONE
@@ -2581,6 +2608,8 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
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), INTENT(inout) :: basin_orog !!
+ 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) :: fetch_basin !!
@@ -2647,6 +2676,8 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
DO ij=1,basin_count(ib)
!
routing_area_glo(ib,ij) = basin_area(ib,ij)
+ routing_orog_glo(ib,ij) = basin_orog(ib,ij)
+ routing_floodp_glo(ib,ij) = basin_floodp(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)
@@ -2751,7 +2782,8 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, basin_cg, basin_topoind,&
+SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, &
+ & basin_orog, basin_floodp, basin_cg, basin_topoind,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
& inflow_grid, inflow_basin)
!
@@ -2770,6 +2802,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_type !!
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area !!
+ REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog !!
+ 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) :: fetch_basin !!
@@ -2802,6 +2836,11 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
!
basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib,tokill)
!
+ basin_floodp(ib, totakeover) = basin_floodp(ib, totakeover) + basin_floodp(ib,tokill)
+ !
+ basin_orog(ib, totakeover) = (basin_orog(ib, totakeover)*basin_area(ib, totakeover) &
+ & + basin_orog(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
+ !
! Add the fetch of the basin will kill to the one which gets the water
!
fetch_basin(ib, totakeover) = fetch_basin(ib, totakeover) + fetch_basin(ib, tokill)
@@ -2937,6 +2976,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
nbasmax_save = nbasmax
ALLOCATE (routing_area_glo(nbpt,nbasmax), stat=ier)
+ ALLOCATE (routing_orog_glo(nbpt,nbasmax), stat=ier)
+ ALLOCATE (routing_floodp_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (route_fetch_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (routing_cg_glo(nbpt,nbasmax,2), stat=ier)
ALLOCATE (topo_resid_glo(nbpt,nbasmax), stat=ier)
diff --git a/HydroGrid.py b/HydroGrid.py
index 70e75f9..66f8488 100644
--- a/HydroGrid.py
+++ b/HydroGrid.py
@@ -49,10 +49,10 @@ def corners(lon, lat) :
#
centersll.append([lon[j,i], lat[j,i]])
radiusll.append(RPP.maxradius([lon[j,i], lat[j,i]], lons, lats))
- #
+ #
index.append([j,i])
cornersll.append(boxll)
-
+
return cornersll, cornerspoly, centersll, radiusll, index
#
def gather(x, index) :
@@ -130,3 +130,15 @@ class HydroData :
#
self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index)
self.facdesc=nf.variables["fac"].long_name
+ #
+ if "orog" in nf.variables.keys():
+ self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index)
+ self.orogdesc=nf.variables["orog"].long_name
+ else:
+ self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
+ #
+ if "floodplains" in nf.variables.keys():
+ self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index)
+ self.floodplainsdesc=nf.variables["floodplains"].long_name
+ else:
+ self.floodplains = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
diff --git a/Interface.py b/Interface.py
index 086a88a..c5dfdeb 100644
--- a/Interface.py
+++ b/Interface.py
@@ -38,7 +38,7 @@ prec = 100.0
#
# Print the documentation for the FORTRAN interface
#
-if gendoc.lower() == "true" :
+if gendoc.lower() == "true" :
docwrapper = open('DocumentationInterface', 'w')
docwrapper.write(routing_interface.initatmgrid.__doc__)
docwrapper.write("====================================================================\n")
@@ -212,9 +212,9 @@ def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fet
#
fetch_error = np.sum(np.abs(fetch_out[part.landcorelist,:]-fetch_in[part.landcorelist,:]), axis=1)\
/np.sum(routing_area[part.landcorelist,:], axis=1)
- if np.max(fetch_error) > prec :
+ if np.max(fetch_error) > prec :
print("Rank :"+str(part.rank)+" Too large fetch error (fraction of greid area) : ", fetch_error)
-
+
print("Total fetch error in fraction of grid box : ", np.sum(fetch_error))
#
return fetch_out
@@ -245,6 +245,8 @@ class HydroOverlap :
topoind_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
fac_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
hierarchy_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
+ orog_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
+ floodp_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
#
trip_tmp[:,:] = np.nan
basins_tmp[:,:] = np.nan
@@ -254,6 +256,8 @@ class HydroOverlap :
topoind_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.topoind[ib][:])
fac_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.fac[ib][:])
hierarchy_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.disto[ib][:])
+ orog_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.orog[ib][:])
+ floodp_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.floodplains[ib][:])
#
trip_tmp[np.isnan(trip_tmp)] = undef_int
basins_tmp[np.isnan(trip_tmp)] = undef_int
@@ -264,9 +268,11 @@ class HydroOverlap :
print("GETHYDROGRID : nbvmax = ", nbvmax)
print("GETHYDROGRID : nbxmax = ", nbxmax)
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(nbxmax, sub_pts, sub_index, sub_area, \
- hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp, hierarchy_tmp)
+ hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp,\
+ hierarchy_tmp, orog_tmp, floodp_tmp)
#
# Plot some diagnostics for the hydrology grid within the atmospheric meshes.
#
@@ -302,16 +308,18 @@ class HydroSuper :
lon_bx_tmp[np.isnan(lon_bx_tmp)] = undef_int
lat_bx_tmp = hydrooverlap.lat_bx
lat_bx_tmp[np.isnan(lat_bx_tmp)] = undef_int
- self.basin_count, self.basin_notrun, self.basin_area, self.basin_cg, self.basin_hierarchy, self.basin_fac, self.basin_topoind, \
+ self.basin_count, self.basin_notrun, self.basin_area, self.basin_cg, self.basin_hierarchy, \
+ self.basin_orog, 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 = \
routing_interface.globalize(hydrooverlap.area_bx, lon_bx_tmp, lat_bx_tmp, hydrooverlap.trip_bx, \
- hydrooverlap.hierarchy_bx, hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \
+ hydrooverlap.hierarchy_bx, hydrooverlap.orog_bx, hydrooverlap.floodp_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) :
@@ -363,7 +371,7 @@ class HydroSuper :
self.fetch_basin = np.copy(fetch_basin)
#
- # Upstream area of the smalest river we call largest rivers.
+ # Upstream area of the smalest river we call largest rivers.
#
self.largest_rivarea = sorted_outareas[l-1]
#
@@ -378,21 +386,26 @@ class HydroSuper :
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
+ 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
+ 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)
+ 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_orog = self.basin_orog, basin_floodp = self.basin_floodp, \
+ 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):
@@ -405,7 +418,7 @@ class HydroSuper :
else :
ncvar = np.zeros((1,1,1))
ncvar[:] = part.gather(var)
-
+
#
def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
#
@@ -430,7 +443,7 @@ class HydroSuper :
#
addcoordinates(outnf, globalgrid, procgrid, part, vtyp, NCFillValue, nbcorners, cornerind)
addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
- #
+ #
# Variables
# Once gathered on root_proc we transform them into float64. Careful, Integer variables do not have NaN !
#
@@ -456,24 +469,30 @@ class HydroSuper :
# self.basin_area
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_area", "Basin area", "-", self.basin_area, vtyp)
#
+ # self.basin_orog
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_orog", "Basin orography", "-", self.basin_orog, vtyp)
+ #
+ # self.basin_floodp
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_floodp", "Basin floodplains", "-", self.basin_floodp, vtyp)
+ #
# self.basin_cg
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "CG_lon", "CG lon", "-", self.basin_cg[:,:,1], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "CG_lat", "CG lat", "-", self.basin_cg[:,:,0], vtyp)
#
# self.topoind
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_topoind", "Topoindex", "-", self.basin_topoind, vtyp)
- #
+ #
# outcoor
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "outcoor_lon", "outcoor lon", "-", self.basin_outcoor[:,:,1], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "outcoor_lat", "outcoor lat", "-", self.basin_outcoor[:,:,0], vtyp)
- #
+ #
# type
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_type", "type", "-", self.basin_type, vtyp)
#
# flowdir
self.add_variable(outnf, procgrid, NCFillValue, part, ('htuext','y','x'), "basin_flowdir", "flowdir", "-", self.basin_flowdir, vtyp)
#
- #
+ #
#self.outflow_grid
grgrid = part.l2glandindex(self.outflow_grid)
grgrid[self.outflow_grid == 0 ] = -2 # in case it flows out of the domain, the 0 should not remain
@@ -516,10 +535,12 @@ class HydroGraph :
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 = \
+ self.routing_area, self.routing_orog, self.routing_floodp, self.routing_cg, self.topo_resid, self.route_nbbasin,\
+ self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
+ self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
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_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
+ basin_orog = hydrosuper.basin_orog, basin_floodp = hydrosuper.basin_floodp, 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, \
@@ -527,7 +548,7 @@ class HydroGraph :
#
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)
#
@@ -536,10 +557,7 @@ class HydroGraph :
gingrid = part.l2glandindex( hydrosuper.inflow_grid[:,:,:self.max_inflow])
self.route_innum, self.route_ingrid, self.route_inbasin = routing_interface.finish_inflows(nbpt = self.nbpt, nbxmax_in = nbxmax_in, nbasmax = nbasmax, inf_max = self.max_inflow, basin_count = hydrosuper.basin_count, inflow_number = hydrosuper.inflow_number, inflow_grid = gingrid, inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
-
-
return
-
#
#
#
@@ -558,7 +576,7 @@ class HydroGraph :
else :
ncvar = np.zeros((1,1,1))
ncvar[:] = part.gather(var)
- #
+ #
#
#
def dumpnetcdf(self, filename, globalgrid, procgrid, part) :
@@ -584,14 +602,14 @@ class HydroGraph :
addenvironment(outnf, procgrid, part, vtyp, NCFillValue, self.nbpt)
#
# land grid index -> to facilitate the analyses of the routing
- #
+ #
nbpt_loc = np.zeros((self.nbpt,1)).astype(np.int32)
nbpt_loc[:,0] = np.arange(1, self.nbpt+1)
nbpt_glo = part.l2glandindex(nbpt_loc)
self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "nbpt_glo", "Grid point Global", "-", nbpt_glo[:,0], vtyp)
#
################
- #
+ #
# TEST: l2glandindex
itarget=-1
for il in range(procgrid.nbland) :
@@ -600,16 +618,16 @@ class HydroGraph :
d=np.sqrt((lo-3.13)**2+(la-39.70)**2)
if d < 0.05 :
itarget = il
-
+
if itarget >+ 0 :
print(part.rank, itarget, " Before route_togrid = ", self.route_togrid[itarget,:])
- # Conversion
+ # Conversion
grgrid = part.l2glandindex(self.route_togrid[:,:])
if itarget >+ 0 :
print(part.rank, itarget, " After route_togrid = ", self.route_togrid[itarget,:])
################
#
- # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
+ # The field route_togrid is with indices on the local grid. That needs to be converted to the global grid.
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetogrid", "Grid into which the basin flows", "-", grgrid, vtyp, "int")
# route to basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "routetobasin", "Basin in to which the water goes", "-", self.route_tobasin[:,:], vtyp, "int")
@@ -618,23 +636,28 @@ class HydroGraph :
#
# basin area
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_area", "area of basin", "m^2", self.routing_area[:,:], vtyp, "float")
+
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_orog", "Mean orography", "m", self.routing_orog[:,:], vtyp, "float")
+
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('htu', 'y','x'), "basin_floodp", "area of floodplains", "m^2", self.routing_floodp[:,:], vtyp, "float")
+
# route number into basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('y','x'), "routenbintobas", "Number of basin into current one", "-", self.route_nbintobas[:], vtyp, "int")
- #
+ #
# original number into basin
self.add_variable(outnf, procgrid, NCFillValue, part, ( 'y','x'), "originnbintobas", "Number of sub-grid basin into current one before truncation", "-", self.origin_nbintobas[:], vtyp, "int")
- #
+ #
# latitude of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlat", "Latitude of Outlet", "degrees north", self.route_outlet[:,:,0], vtyp, "float")
# longitude of outlet
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outletlon", "Longitude of Outlet", "degrees east", self.route_outlet[:,:,1], vtyp, "float")
# type of outlet
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "outlettype", "Type of outlet", "code", self.route_type[:,:], vtyp, "float")
#
# topographic index
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "topoindex", "Topographic index of the retention time", "m", self.topo_resid[:,:], vtyp, "float")
#
-
+
# Inflow number
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, "int")
#
@@ -644,14 +667,14 @@ class HydroGraph :
#
# Inflow basin
self.add_variable(outnf, procgrid, NCFillValue, part, ('inflow', 'htu','y','x'), "route_inbasin", "Basin from which the water flows", "-", self.route_inbasin[:,:,:], vtyp, "int")
-
+
#
# Save centre of gravity of HTU
#
# Check if it works
self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "CG_lon", "Longitude of centre of gravity of HTU", "degrees east", self.routing_cg[:,:,1], vtyp, "float")
-
- self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
+
+ self.add_variable(outnf, procgrid, NCFillValue, part, ('htu','y','x'), "CG_lat", "Latitude of centre of gravity of HTU", "degrees north", self.routing_cg[:,:,0], vtyp, "float")
#
# Save the fetch of each basin
#
@@ -662,7 +685,3 @@ class HydroGraph :
outnf.close()
#
return
-
-
-
-
diff --git a/tests/Iberia_n48/run.def b/tests/Iberia_n48/run.def
index 9ef68a4..6fe2932 100644
--- a/tests/Iberia_n48/run.def
+++ b/tests/Iberia_n48/run.def
@@ -19,7 +19,7 @@ nbasmax = 35
#
# Number of operation of simplification performed together
#
-numop = 200
+numop = 100
#
# Output
#
diff --git a/tests/Iberia_n48_regular/run.def b/tests/Iberia_n48_regular/run.def
index b0afac7..1cb8a04 100644
--- a/tests/Iberia_n48_regular/run.def
+++ b/tests/Iberia_n48_regular/run.def
@@ -7,7 +7,7 @@ EarthRadius = 6370000.
ModelGridFile = /bdd/MEDI/workspaces/polcher/NewRouting/EM_WFDEI_CRU_2000.nc
WEST_EAST = -9.75, 5.25
SOUTH_NORTH = 35.5, 43.5
-HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
+HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
#
# FORTRAN interface parameters
#
--
GitLab