Commit 56e82d69 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Propose various metods to compute topoindex on the HTUs. The river length and...

Propose various metods to compute topoindex on the HTUs. The river length and elevation change within the HTU are not archived in the routing-graph file.
parent c8b9a77d
......@@ -454,11 +454,11 @@ END SUBROUTINE rivclassification
SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridarea, cfrac, gridcenters, basin_count, &
& basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, &
& basin_floodp, basin_cg, basin_topoind, basin_beta_fp, fetch_basin, &
& basin_floodp, basin_cg, basin_topoind, basin_rlen, basin_rdz, basin_beta_fp, fetch_basin, &
& basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin, floodcri, routing_area, routing_orog_mean, &
& routing_orog_min, routing_orog_max, routing_floodp, routing_beta,&
& routing_cg, topo_resid, route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
& routing_cg, topo_resid, topo_rlen, topo_rdz, route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
& global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch,routing_floodcri, &
& rfillval, ifillval)
!
......@@ -490,7 +490,9 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_orog_max !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_cg
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (km)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_rlen
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_rdz
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_beta_fp
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
......@@ -516,7 +518,9 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_beta
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_fetch !! Upstream are flowing into HTU (m^2)
REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out) :: routing_cg !! Centre of gravity of HTU (Lat, Lon)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: topo_resid !! Topographic index of the retention time (m)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: topo_resid !! Topographic index of the retention time (km)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: topo_rlen !!
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: topo_rdz !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(out) :: route_nbbasin !!
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_togrid !! Grid into which the basin flows (unitless)
INTEGER(i_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_tobasin !! Basin in to which the water goes (unitless)
......@@ -532,8 +536,8 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
!
CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, gridcenters, nwbas, inflowmax, num_largest, &
& basin_count, basin_notrun, basin_area, basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp,&
& basin_cg, basin_topoind, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
& basin_cg, basin_topoind, basin_rlen, basin_rdz, basin_beta_fp, fetch_basin, basin_id, basin_coor, &
& basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
routing_area(:,:) = rfillval
routing_orog_mean(:,:) = rfillval
......@@ -542,6 +546,8 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
routing_floodp(:,:) = rfillval
routing_cg(:,:,:) = rfillval
topo_resid(:,:) = rfillval
topo_rlen(:,:) = rfillval
topo_rdz(:,:) = rfillval
route_nbbasin(:) = ifillval
route_togrid(:,:) = ifillval
route_tobasin(:,:) = ifillval
......@@ -562,7 +568,11 @@ SUBROUTINE finish_truncate(nbpt, inflowmax, nbasmax, nwbas, num_largest, gridare
routing_orog_max(ij,ib) = routing_orog_max_glo(ij,ib)
routing_floodp(ij,ib) = routing_floodp_glo(ij,ib)
routing_cg(ij,ib,:) = routing_cg_glo(ij,ib,:)
!
topo_resid(ij,ib) = topo_resid_glo(ij,ib)
topo_rlen(ij,ib) = topo_rlen_glo(ij,ib)
topo_rdz(ij,ib) = topo_rdz_glo(ij,ib)
!
route_nbbasin(:) = route_count_glo(:)
route_togrid(ij,ib) = route_togrid_glo(ij,ib)
route_tobasin(ij,ib) = route_tobasin_glo(ij,ib)
......@@ -630,9 +640,9 @@ END SUBROUTINE finish_inflows
SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, numops, basin_count, basin_area, &
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, &
& basin_cg, basin_topoind, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir,&
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, &
& basin_topoind, basin_rlen, basin_rdz, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, &
& basin_flowdir, outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
!
!
USE ioipsl
......@@ -660,7 +670,9 @@ SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, num
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_cg
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (km)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_rlen
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_rdz
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_beta_fp
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
......@@ -711,9 +723,9 @@ SUBROUTINE killbas(nbpt, inflowmax, nbasmax, nwbas, ops, tokill, totakeover, num
END IF
END DO
CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, inflowmax, basin_count, basin_area, &
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind, basin_beta_fp,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
& inflow_grid, inflow_basin)
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind, basin_rlen, &
& basin_rdz, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
END IF
END IF
END DO
......
......@@ -24,7 +24,10 @@ MODULE routing_reg
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_beta_glo !!
!
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: routing_cg_glo !! Centre of gravity of HTU (Lat, Lon)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (km)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_rlen_glo !! HTU river length (m)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_rdz_glo !! HTU river elevation change (m)
!
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: route_count_glo !! Number of basins per grid
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_togrid_glo !! Grid into which the basin flows (unitless)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: route_tobasin_glo !! Basin in to which the water goes (unitless)
......@@ -39,11 +42,17 @@ MODULE routing_reg
! Options to compute the topographic index based on the information available in the Hydrological files.
! topo_option = 1 : Constant topographic index
! topo_option = 2 : Simple average without any weighting
! topo_option = 3 : The sum over all rivers of the HTU are computed in _findbasin. These
! values are then averaged.
! topo_option = 4 : As above but the average topoindex is computed for each river.
! topo_option = 3 : The sum over all rivers within the HTU are computed in _findbasin. These
! values are then averaged over the HTU.
! topo_option = 4 : The length and elevation change of all rivers within the HTU are computed.
! These values are then used to compute the topoindex of the HTU.
!
INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex
!
INTEGER(i_std), SAVE :: topo_option = 3 !! Option to calculate topoindex
! Options to compute the properties when merging HTUs in routing_reg_killbas
! kill_option = 1 : The old way of simpling building the weighted average.
! Kill_option = 2 : When the difference of average elevation then we assume HTUs cascade into each other
INTEGER(i_std), SAVE :: kill_option = 2
!
CONTAINS
!! ================================================================================================================================
......@@ -116,9 +125,9 @@ CONTAINS
!
REAL(r_std), INTENT(inout) :: trip(nbpt,nbvmax) !! The trip field (unitless)
REAL(r_std), INTENT(inout) :: basins(nbpt,nbvmax) !! data on the fine grid
REAL(r_std), INTENT(inout) :: topoindex(nbpt,nbvmax) !! Topographic index of the residence time (m)
REAL(r_std), INTENT(inout) :: rlen(nbpt,nbvmax) !!
REAL(r_std), INTENT(inout) :: rdz(nbpt,nbvmax) !!
REAL(r_std), INTENT(inout) :: topoindex(nbpt,nbvmax) !! Topographic index of the residence time (km)
REAL(r_std), INTENT(inout) :: rlen(nbpt,nbvmax) !! River length within the hydrological grid (m)
REAL(r_std), INTENT(inout) :: rdz(nbpt,nbvmax) !! Elevation change within the hydrological grid (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
......@@ -132,9 +141,9 @@ CONTAINS
REAL(r_std), INTENT(out) :: lon_bx(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(out) :: lat_bx(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(out) :: lshead_bx(ijdimmax,ijdimmax) !! Large scale heading for outflow points.
REAL(r_std), INTENT(out) :: topoind_bx(ijdimmax,ijdimmax) !! Topographic index of the residence time for each of the smaller boxes (m)
REAL(r_std), INTENT(out) :: rlen_bx(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(out) :: rdz_bx(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(out) :: topoind_bx(ijdimmax,ijdimmax) !! Topographic index of the residence time for each of the smaller boxes (km)
REAL(r_std), INTENT(out) :: rlen_bx(ijdimmax,ijdimmax) !! River length within the smaller boxes (m)
REAL(r_std), INTENT(out) :: rdz_bx(ijdimmax,ijdimmax) !! Elevation change within the smaller boxes (m)
INTEGER(i_std), INTENT(out) :: trip_bx(ijdimmax,ijdimmax) !! The trip field for each of the smaller boxes (unitless)
INTEGER(i_std), INTENT(out) :: basin_bx(ijdimmax,ijdimmax) !!
REAL(i_std), INTENT(out) :: orog_bx(ijdimmax,ijdimmax) !!
......@@ -421,8 +430,8 @@ CONTAINS
INTEGER(i_std), INTENT(inout) :: trip(ijdimmax,ijdimmax) !! The trip field (unitless)
INTEGER(i_std), INTENT(inout) :: basin(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(inout) :: topoind(ijdimmax,ijdimmax) !! Topographic index of the residence time (km)
REAL(r_std), INTENT(inout) :: rlen(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(inout) :: rdz(ijdimmax,ijdimmax) !!
REAL(r_std), INTENT(inout) :: rlen(ijdimmax,ijdimmax) !! River length (m)
REAL(r_std), INTENT(inout) :: rdz(ijdimmax,ijdimmax) !! Elevation change (m)
REAL(r_std), INTENT(out) :: rweight(ijdimmax,ijdimmax) !!
!
! For debug and get coordinate of river outlet
......@@ -1010,7 +1019,7 @@ CONTAINS
!
! Work only for topo_option number 3 or 4
!
IF ( topo_option .EQ. 3 ) THEN
IF ( topo_option .EQ. 3 .OR. topo_option .EQ. 4) THEN
rivpas(:,:) = 0
! Loop for all sub-basins and each point belongs to that sub-basin
DO ij=1,nb_basin
......@@ -1603,10 +1612,10 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: orog_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: floodp_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: fac_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: topoind_bx !! Topographic index of the residence time for each of the smaller boxes (m)
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: rlen_bx !!
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: rdz_bx !!
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: rweight_bx !!
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: topoind_bx !! Topographic index for each of the smaller boxes (km)
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: rlen_bx !! River length within the HTU (m)
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: rdz_bx !! Elevation change of rivers within the HTU (m)
REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: rweight_bx !! Area covered by each river within the HTU
REAL(r_std) :: min_topoind !! The current minimum of topographic index (m)
INTEGER(i_std) :: nb_basin !! Number of sub-basins (unitless)
INTEGER(i_std), DIMENSION(nb_htu) :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
......@@ -1629,16 +1638,16 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
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_mean !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_min !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_max !!
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 (km)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_rlen
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_rdz
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_min !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_max !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_fac !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Mean topographic index for HTU (km)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_rlen !! Mean river length within HTU (m)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_rdz !! Mean elevation change within HTU (m)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_lshead !! Large scale heading out of the grid box.
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_std !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_beta_fp !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_std !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_beta_fp !!
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !!
INTEGER(i_std), DIMENSION(nbpt) :: nbcoastal !!
......@@ -1666,9 +1675,9 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!_ ================================================================================================================================
!
!
IF ( topo_option .LT. 1 .OR. topo_option .GT. 3 ) THEN
IF ( topo_option .LT. 1 .OR. topo_option .GT. 4 ) THEN
WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
WRITE(numout,*) ' It should be: 1, 2, or 3'
WRITE(numout,*) ' It should be: 1, 2, 3 or 4'
STOP 'routing_reg_globalize'
ENDIF
!
......@@ -1855,11 +1864,16 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
basin_rlen(ib,ij) = basin_rlen(ib,ij)/cnt
basin_rdz(ib,ij) = basin_rdz(ib,ij)/cnt
ELSE
! More reasonable variables are needed here.
! Mean values of oringal hydro date are used to fill-in missing points.
basin_rlen(ib,ij) = hydro_meanlen
basin_rdz(ib,ij) = hydro_meandz
ENDIF
!
!
IF ( topo_option .EQ. 4 ) THEN
basin_topoind(ib,ij) = SQRT(basin_rlen(ib,ij)**3./(basin_rdz(ib,ij)*1.0e6))
ENDIF
!
! To make sure that it has the lowest number if this is an outflow point we reset basin_hierarchy
!
IF (basin_bxout(ij) .LT. 0 .AND. basin_bxout(ij) .NE. -4) THEN
......@@ -2668,8 +2682,8 @@ SUBROUTINE isin_halo(ig, nbhalo, nbpt, halopts, isinhalo)
SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcenters, nwbas,&
& inflowmax, num_largest, basin_count, basin_notrun, basin_area, basin_orog_mean,&
& basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind,basin_beta_fp, fetch_basin,&
& basin_id, basin_coor, basin_type, basin_flowdir, &
& basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind, basin_rlen, basin_rdz, &
& basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin, floodcri)
!
IMPLICIT NONE
......@@ -2701,7 +2715,9 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_orog_max !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_cg !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Mean Topographic index of HTU (km)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_rlen
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_rdz
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_beta_fp !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
......@@ -2759,8 +2775,10 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
routing_floodp_glo(ib,:) = 0
routing_beta_glo(ib,:) = 0
routing_cg_glo(ib,:,:) = 0.0
!
! Perhaps here we should put a NaN so we can filter it out later !!!
topo_resid_glo(ib,:) = topoindexmin
topo_rlen_glo(ib,:) = hydro_meanlen
topo_rdz_glo(ib,:) = hydro_meandz
global_basinid_glo(ib,:) = zero
route_outlet_glo(ib,:,:) = zero
!
......@@ -2788,6 +2806,8 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
! To define in function of the minimum topoindex on the map
! (to revise)
basin_topoind(ib,1) = topoindexmin
basin_rlen(ib,1) = hydro_meanlen
basin_rdz(ib,1) = hydro_meandz
basin_id(ib,1) = 0
basin_coor(ib,1,1) = basin_cg(ib,1,1)
basin_coor(ib,1,2) = basin_cg(ib,1,2)
......@@ -2816,6 +2836,9 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
routing_cg_glo(ib,ij,:) = basin_cg(ib,ij,:)
!
topo_resid_glo(ib,ij) = basin_topoind(ib,ij)
topo_rlen_glo(ib,ij) = basin_rlen(ib,ij)
topo_rdz_glo(ib,ij) = basin_rdz(ib,ij)
!
global_basinid_glo(ib,ij) = basin_id(ib,ij)
route_outlet_glo(ib,ij,1) = basin_coor(ib,ij,1)
route_outlet_glo(ib,ij,2) = basin_coor(ib,ij,2)
......@@ -2923,9 +2946,9 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
!_ ================================================================================================================================
SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, basin_count, basin_area, &
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind, basin_beta_fp,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
& inflow_grid, inflow_basin)
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_cg, basin_topoind, basin_rlen, &
& basin_rdz, basin_beta_fp, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, &
& outflow_basin, inflow_number, inflow_grid, inflow_basin)
!
!
IMPLICIT NONE
......@@ -2947,7 +2970,9 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog_max !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_cg !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (km)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_rlen
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_rdz
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_beta_fp !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin !!
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
......@@ -2957,6 +2982,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number !!
INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_basin !!
INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_grid !!
REAL(r_std) :: meandz
!
!! LOCAL VARIABLES
INTEGER(i_std) :: inf, ibs, ing, inb, ibasf, igrif, it !! Indices (unitless)
......@@ -2968,20 +2994,37 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
! For the moment only area
!
!
!basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1)*basin_area(ib, totakeover) &
!& + basin_cg(ib, tokill, 1)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
!basin_cg(ib, totakeover, 2) = (basin_cg(ib, totakeover, 2)*basin_area(ib, totakeover)&
!& + basin_cg(ib, tokill, 2)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
!!$ IF (outflow_basin(ib,tokill) == totakeover .OR. outflow_basin(ib,totakeover) == tokill ) THEN
!!$ basin_topoind(ib, totakeover) = basin_topoind(ib, totakeover)+basin_topoind(ib, tokill)
!!$ ELSE
!!$ basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
!!$ & + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!!$ ENDIF
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
& + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
IF ( kill_option .EQ. 1 ) THEN
!
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
& + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
basin_rlen(ib, totakeover) = (basin_rlen(ib, totakeover)+basin_rlen(ib, tokill))*0.5
basin_rdz(ib, totakeover) = (basin_rdz(ib, totakeover)+basin_rdz(ib, tokill))*0.5
!
ELSE IF ( kill_option .EQ. 2 ) THEN
IF (outflow_basin(ib,tokill) == totakeover .OR. outflow_basin(ib,totakeover) == tokill ) THEN
meandz = 0.5*(basin_rdz(ib,tokill)+basin_rdz(ib,totakeover))
IF ( ABS(basin_orog_mean(ib,tokill)-basin_orog_mean(ib,totakeover)) > meandz ) THEN
! HTUs flow into each other
basin_topoind(ib, totakeover) = basin_topoind(ib, totakeover)+basin_topoind(ib, tokill)
basin_rlen(ib, totakeover) = basin_rlen(ib, totakeover)+basin_rlen(ib, tokill)
basin_rdz(ib, totakeover) = basin_rdz(ib, totakeover)+basin_rdz(ib, tokill)
ELSE
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
& + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
basin_rlen(ib, totakeover) = (basin_rlen(ib, totakeover)+basin_rlen(ib, tokill))*0.5
basin_rdz(ib, totakeover) = (basin_rdz(ib, totakeover)+basin_rdz(ib, tokill))*0.5
ENDIF
ELSE
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) &
& + basin_topoind(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
basin_rlen(ib, totakeover) = (basin_rlen(ib, totakeover)+basin_rlen(ib, tokill))*0.5
basin_rdz(ib, totakeover) = (basin_rdz(ib, totakeover)+basin_rdz(ib, tokill))*0.5
ENDIF
ELSE
CALL ipslerr_p(3,'routing_reg_killbas','The selected kill_option does not exist.','','')
ENDIF
!
basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib,tokill)
!
......@@ -3089,6 +3132,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
basin_cg(ib, basin_count(ib):nwbas,:) = zero
basin_topoind(ib, tokill:basin_count(ib)-1) = basin_topoind(ib, tokill+1:basin_count(ib))
basin_topoind(ib, basin_count(ib):nwbas) = zero
basin_rlen(ib, basin_count(ib):nwbas) = zero
basin_rdz(ib, basin_count(ib):nwbas) = zero
fetch_basin(ib, tokill:basin_count(ib)-1) = fetch_basin(ib, tokill+1:basin_count(ib))
fetch_basin(ib, basin_count(ib):nwbas) = zero
!
......@@ -3166,6 +3211,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
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)
ALLOCATE (topo_rlen_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (topo_rdz_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (route_count_glo(nbpt), stat=ier)
ALLOCATE (route_togrid_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (route_tobasin_glo(nbpt,nbasmax), stat=ier)
......
......@@ -68,9 +68,9 @@ class HydroGraph :
#
#
self.routing_area, self.routing_orog_mean, self.routing_orog_min,self.routing_orog_max, \
self.routing_floodp,self.routing_beta, self.routing_cg, self.topo_resid, self.route_nbbasin,\
self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch, self.floodcri = \
self.routing_floodp,self.routing_beta, self.routing_cg, self.topo_resid, self.topo_rlen, \
self.topo_rdz, self.route_nbbasin, self.route_togrid, self.route_tobasin, self.route_nbintobas, \
self.global_basinid, self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch, self.floodcri = \
routing_interface.finish_truncate(nbpt = self.nbpt, inflowmax = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, \
cfrac = modelgrid.contfrac, gridcenters = np.array(modelgrid.centers), \
......@@ -79,6 +79,7 @@ class HydroGraph :
basin_orog_mean = hydrosuper.basin_orog_mean, basin_orog_min = hydrosuper.basin_orog_min,\
basin_orog_max = hydrosuper.basin_orog_max, basin_floodp = hydrosuper.basin_floodp, \
basin_cg = hydrosuper.basin_cg, basin_topoind = hydrosuper.basin_topoind, \
basin_rlen = hydrosuper.basin_rlen, basin_rdz = hydrosuper.basin_rdz, \
basin_beta_fp = hydrosuper.basin_beta_fp , fetch_basin = hydrosuper.fetch_basin, \
basin_id = hydrosuper.basin_id, \
basin_coor = hydrosuper.basin_outcoor, basin_type = hydrosuper.basin_type, \
......@@ -312,7 +313,11 @@ class HydroGraph :
self.add_tstepdistrib(outnf, procgrid, NCFillValue, part, self.topo_resid[:,:], self.routing_area[:,:], \
vtyp, nbbins, tcst)
#
# Geometric properties of HTU
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
"riverlen", "Mean river length within HTU", "m", self.topo_rlen[:,:], vtyp)
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
"riverdz", "Mean river elevation change within HTU", "m", self.topo_rdz[:,:], vtyp)
# Inflow number
self.add_variable(outnf, procgrid, NCFillValue, part, ('z','y','x'), \
"route_innum", "Number of inflow", "-", self.route_innum[:,:], vtyp, orig_type="int")
......
......@@ -362,7 +362,8 @@ class HydroSuper :
totakeover = totakeover, numops = numops, basin_count = self.basin_count,\
basin_area = self.basin_area, basin_orog_mean = self.basin_orog_mean, \
basin_orog_min = self.basin_orog_min, basin_orog_max = self.basin_orog_max, basin_floodp = self.basin_floodp, \
basin_cg = self.basin_cg, basin_topoind = self.basin_topoind, basin_beta_fp = self.basin_beta_fp, fetch_basin = self.fetch_basin,\
basin_cg = self.basin_cg, basin_topoind = self.basin_topoind, basin_rlen = self.basin_rlen, \
basin_rdz = self.basin_rdz, basin_beta_fp = self.basin_beta_fp, fetch_basin = self.fetch_basin,\
basin_id = self.basin_id, basin_coor = self.basin_outcoor, basin_type = self.basin_type,\
basin_flowdir = self.basin_flowdir, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, \
inflow_number = self.inflow_number, inflow_grid = self.inflow_grid, inflow_basin = self.inflow_basin)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment