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

Test the dz limiter. It must now be between the minimum dz in the original...

Test the dz limiter. It must now be between the minimum dz in the original data and the difference orogÃ_max and orog_min of the HTU.
parent 60c91757
...@@ -44,7 +44,7 @@ SUBROUTINE closeinterface() ...@@ -44,7 +44,7 @@ SUBROUTINE closeinterface()
END SUBROUTINE closeinterface END SUBROUTINE closeinterface
SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_index, sub_area, max_basins, & SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_index, sub_area, max_basins, &
& min_topoind, meanrlen, meanrdz, lon_rel, lat_rel, trip, basins, topoindex, rlen, rdz, fac, hierarchy, & & min_topoind, meanrlen, meanrdz, minrdz, lon_rel, lat_rel, trip, basins, topoindex, rlen, rdz, fac, hierarchy, &
& orog, floodp, nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, rlen_bx, rdz_bx, rweight_bx, fac_bx, & & orog, floodp, nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, rlen_bx, rdz_bx, rweight_bx, fac_bx, &
& hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx) & hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
! !
...@@ -62,7 +62,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i ...@@ -62,7 +62,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
INTEGER, INTENT(in) :: sub_index(nbpt,nbvmax_in,2) !! Indices of the points we need on the fine grid (unitless) INTEGER, INTENT(in) :: sub_index(nbpt,nbvmax_in,2) !! Indices of the points we need on the fine grid (unitless)
REAL, INTENT(inout) :: max_basins !! The current maximum of basins REAL, INTENT(inout) :: max_basins !! The current maximum of basins
REAL, INTENT(in) :: min_topoind !! The current minimum of topographic index (km) REAL, INTENT(in) :: min_topoind !! The current minimum of topographic index (km)
REAL, INTENT(in) :: meanrlen, meanrdz !! Mean length and altitude change in the original grid (m) REAL, INTENT(in) :: meanrlen, meanrdz, minrdz !! Mean length and altitude change in the original grid (m)
REAL, INTENT(in) :: sub_area(nbpt,nbvmax_in) !! Area on the fine grid (m^2) REAL, INTENT(in) :: sub_area(nbpt,nbvmax_in) !! Area on the fine grid (m^2)
! !
REAL, INTENT(in) :: lon_rel(nbpt,nbvmax_in) !! REAL, INTENT(in) :: lon_rel(nbpt,nbvmax_in) !!
...@@ -114,7 +114,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i ...@@ -114,7 +114,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
DO ib=1,nbpt DO ib=1,nbpt
CALL routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, & CALL routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& meanrlen, meanrdz, lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, & & meanrlen, meanrdz, minrdz, lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, &
& rlen, rdz, fac, hierarchy, orog, floodp, nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), & & rlen, rdz, fac, hierarchy, orog, floodp, nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), &
& basin_bx(ib,:,:), topoind_bx(ib,:,:), rlen_bx(ib,:,:), rdz_bx(ib,:,:), fac_bx(ib,:,:), & & basin_bx(ib,:,:), topoind_bx(ib,:,:), rlen_bx(ib,:,:), rdz_bx(ib,:,:), fac_bx(ib,:,:), &
& hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), & & hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), &
......
...@@ -15,6 +15,7 @@ MODULE routing_reg ...@@ -15,6 +15,7 @@ MODULE routing_reg
REAL(r_std), SAVE :: topoindexmin = 0.0 REAL(r_std), SAVE :: topoindexmin = 0.0
REAL(r_std), SAVE :: hydro_meanlen = 0.0 REAL(r_std), SAVE :: hydro_meanlen = 0.0
REAL(r_std), SAVE :: hydro_meandz = 0.0 REAL(r_std), SAVE :: hydro_meandz = 0.0
REAL(r_std), SAVE :: hydro_mindz = 0.0
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_orog_mean_glo !! Mean topography (m) REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_orog_mean_glo !! Mean topography (m)
...@@ -46,6 +47,7 @@ MODULE routing_reg ...@@ -46,6 +47,7 @@ MODULE routing_reg
! values are then averaged over the HTU. ! values are then averaged over the HTU.
! topo_option = 4 : The length and elevation change of all rivers within the HTU are computed. ! 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. ! These values are then used to compute the topoindex of the HTU.
! topo_option = 5 : The average of the dz have to be larger than min(dz) and smaller than orog_max-orog_min
! !
INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex INTEGER(i_std), SAVE :: topo_option = 4 !! Option to calculate topoindex
! !
...@@ -54,10 +56,14 @@ MODULE routing_reg ...@@ -54,10 +56,14 @@ MODULE routing_reg
! Kill_option = 2 : When the difference of average elevation then we assume HTUs cascade into each other ! Kill_option = 2 : When the difference of average elevation then we assume HTUs cascade into each other
INTEGER(i_std), SAVE :: kill_option = 2 INTEGER(i_std), SAVE :: kill_option = 2
! !
! Option to limit the dz when computing topoindex
!
INTEGER(i_std), SAVE :: dzrang_limiter = 1
! Tested combinations : ! Tested combinations :
! Asum : topo_option = 3, kill_option = 1 ! Asum : topo_option = 3, kill_option = 1, dzrang_limiter = 0
! Bsum : topo_option = 4, kill_option = 1 ! Bsum : topo_option = 4, kill_option = 1, dzrang_limiter = 0
! Csum : topo_option = 4, kill_option = 2 ! Csum : topo_option = 4, kill_option = 2, dzrang_limiter = 0
! Dsum : topo_option = 4, kill_option = 2, dzrang_limiter = 1
! !
CONTAINS CONTAINS
!! ================================================================================================================================ !! ================================================================================================================================
...@@ -104,7 +110,7 @@ CONTAINS ...@@ -104,7 +110,7 @@ CONTAINS
!_ ================================================================================================================================ !_ ================================================================================================================================
SUBROUTINE routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, & SUBROUTINE routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& meanrlen, meanrdz, lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, rlen, rdz, & & meanrlen, meanrdz, mindz, lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, rlen, rdz, &
& fac, hierarchy, orog, floodp, nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, rlen_bx, rdz_bx, fac_bx, & & fac, hierarchy, orog, floodp, nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, rlen_bx, rdz_bx, fac_bx, &
& hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx) & hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
! !
...@@ -118,7 +124,7 @@ CONTAINS ...@@ -118,7 +124,7 @@ CONTAINS
INTEGER(i_std), INTENT(in) :: sub_index(nbpt,nbvmax,2) !! Indices of the points we need on the fine grid (unitless) INTEGER(i_std), INTENT(in) :: sub_index(nbpt,nbvmax,2) !! Indices of the points we need on the fine grid (unitless)
REAL(r_std), INTENT(inout) :: max_basins !! The current maximum of basins REAL(r_std), INTENT(inout) :: max_basins !! The current maximum of basins
REAL(r_std), INTENT(in) :: min_topoind !! The current minimum of topographic index (m) REAL(r_std), INTENT(in) :: min_topoind !! The current minimum of topographic index (m)
REAL(r_std), INTENT(in) :: meanrlen, meanrdz !! Mean length and altitude change in the original grid (m) REAL(r_std), INTENT(in) :: meanrlen, meanrdz, mindz !! Mean length and altitude change in the original grid (m)
REAL(r_std), INTENT(in) :: sub_area(nbpt,nbvmax) !! Area on the fine grid (m^2) REAL(r_std), INTENT(in) :: sub_area(nbpt,nbvmax) !! Area on the fine grid (m^2)
! !
REAL(r_std), INTENT(in) :: lon_rel(nbpt,nbvmax) !! REAL(r_std), INTENT(in) :: lon_rel(nbpt,nbvmax) !!
...@@ -187,6 +193,7 @@ CONTAINS ...@@ -187,6 +193,7 @@ CONTAINS
topoindexmin = min_topoind topoindexmin = min_topoind
hydro_meanlen = meanrlen hydro_meanlen = meanrlen
hydro_meandz = meanrdz hydro_meandz = meanrdz
hydro_mindz = mindz
! !
! Set everything to undef to locate easily empty points ! Set everything to undef to locate easily empty points
! !
...@@ -1680,7 +1687,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar ...@@ -1680,7 +1687,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!_ ================================================================================================================================ !_ ================================================================================================================================
! !
! !
IF ( topo_option .LT. 1 .OR. topo_option .GT. 4 ) THEN IF ( topo_option .LT. 1 .OR. topo_option .GT. 5 ) THEN
WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option WRITE(numout,*) ' You chose wrong topo_option for calculating topoindex: ', topo_option
WRITE(numout,*) ' It should be: 1, 2, 3 or 4' WRITE(numout,*) ' It should be: 1, 2, 3 or 4'
STOP 'routing_reg_globalize' STOP 'routing_reg_globalize'
...@@ -1867,7 +1874,12 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar ...@@ -1867,7 +1874,12 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
ENDDO ENDDO
IF ( cnt > 0 ) THEN IF ( cnt > 0 ) THEN
basin_rlen(ib,ij) = basin_rlen(ib,ij)/cnt basin_rlen(ib,ij) = basin_rlen(ib,ij)/cnt
basin_rdz(ib,ij) = basin_rdz(ib,ij)/cnt IF ( dzrang_limiter == 0 ) THEN
basin_rdz(ib,ij) = basin_rdz(ib,ij)/cnt
ELSE
! dZ has to be between hydro_mindz and the maximum height difference within the HTU
basin_rdz(ib,ij) = MAX(MIN(basin_rdz(ib,ij)/cnt, basin_orog_max(ib,ij)-basin_orog_min(ib,ij)), hydro_mindz)
ENDIF
ELSE ELSE
! Mean values of oringal hydro date are used to fill-in missing points. ! Mean values of oringal hydro date are used to fill-in missing points.
basin_rlen(ib,ij) = hydro_meanlen basin_rlen(ib,ij) = hydro_meanlen
......
...@@ -147,9 +147,10 @@ def calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) : ...@@ -147,9 +147,10 @@ def calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) :
# #
meanrlen = np.nanmean(rlen) meanrlen = np.nanmean(rlen)
meanrdz = np.nanmean(rdz) meanrdz = np.nanmean(rdz)
mindz=np.nanmin(rdz)
# #
return gather(topoind, index, default = 10), gather(rlen, index, default = 10), gather(rdz, index, default = 10), \ return gather(topoind, index, default = 10), gather(rlen, index, default = 10), gather(rdz, index, default = 10), \
"Topographic index which serves for the residence time", mintopo, maxtopo, meanrlen, meanrdz "Topographic index which serves for the residence time", mintopo, maxtopo, meanrlen, meanrdz, mindz
# #
class HydroGrid : class HydroGrid :
def __init__(self, lolacorners, wfile) : def __init__(self, lolacorners, wfile) :
...@@ -264,12 +265,13 @@ class HydroData : ...@@ -264,12 +265,13 @@ class HydroData :
self.topoinddesc = "desc" self.topoinddesc = "desc"
mintopo = 1 mintopo = 1
else: else:
self.topoind, self.rlen, self.rdz, self.topoinddesc, mintopo, maxtopo, mrlen, mrdz = \ self.topoind, self.rlen, self.rdz, self.topoinddesc, mintopo, maxtopo, mrlen, mrdz, minrdz = \
calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid)
self.topoindmin = part.domainmin(mintopo) self.topoindmin = part.domainmin(mintopo)
self.topoindmax = part.domainmax(maxtopo) self.topoindmax = part.domainmax(maxtopo)
self.meanlength = part.domainmean(mrlen) self.meanlength = part.domainmean(mrlen)
self.meandz = part.domainmean(mrdz) self.meandz = part.domainmean(mrdz)
self.mindz = part.domainmin(minrdz)
elif "topoind" in nf.variables.keys() and "orog" in nf.variables.keys() : elif "topoind" in nf.variables.keys() and "orog" in nf.variables.keys() :
self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index) self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, default = 10) self.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, default = 10)
......
...@@ -138,7 +138,7 @@ class HydroOverlap : ...@@ -138,7 +138,7 @@ class HydroOverlap :
self.rweight_bx, self.fac_bx, self.hierarchy_bx, self.orog_bx, self.floodp_bx, self.lon_bx, self.lat_bx, \ self.rweight_bx, self.fac_bx, self.hierarchy_bx, self.orog_bx, self.floodp_bx, self.lon_bx, self.lat_bx, \
self.lshead_bx = routing_interface.gethydrogrid(ijdimmax, maxpercent, sub_pts, sub_index, sub_area, \ self.lshead_bx = routing_interface.gethydrogrid(ijdimmax, maxpercent, sub_pts, sub_index, sub_area, \
hydrodata.basinsmax, hydrodata.topoindmin, hydrodata.meanlength, \ hydrodata.basinsmax, hydrodata.topoindmin, hydrodata.meanlength, \
hydrodata.meandz, sub_lon, sub_lat, trip_tmp, \ hydrodata.meandz, hydrodata.mindz, sub_lon, sub_lat, trip_tmp, \
basins_tmp, topoind_tmp, rlen_tmp, rdz_tmp, fac_tmp,\ basins_tmp, topoind_tmp, rlen_tmp, rdz_tmp, fac_tmp,\
hierarchy_tmp, orog_tmp, floodp_tmp) hierarchy_tmp, orog_tmp, floodp_tmp)
# #
......
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