Skip to content
GitLab
Projects
Groups
Snippets
Help
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in / Register
Toggle navigation
Open sidebar
IPSL
LMD
InTro
RoutingPP
Commits
ce3feb48
Commit
ce3feb48
authored
Apr 20, 2021
by
Anthony
Browse files
Improvement of rdz / rlen / topoindex for mainstream river
parent
4fc02c8b
Changes
3
Hide whitespace changes
Inline
Side-by-side
Showing
3 changed files
with
140 additions
and
97 deletions
+140
-97
F90subroutines/routing_interface.f90
F90subroutines/routing_interface.f90
+3
-2
F90subroutines/routing_reg.f90
F90subroutines/routing_reg.f90
+135
-94
Interface.py
Interface.py
+2
-1
No files found.
F90subroutines/routing_interface.f90
View file @
ce3feb48
...
...
@@ -129,7 +129,7 @@ END SUBROUTINE gethydrogrid
SUBROUTINE
findbasins
(
nbpt
,
nb_htu
,
nbv
,
ijdimmax
,
nbi
,
nbj
,
trip_bx
,
basin_bx
,
fac_bx
,
hierarchy_bx
,
topoind_bx
,
&
&
rlen_bx
,
rdz_bx
,
rweight_bx
,
lshead_bx
,
nb_basin
,
basin_inbxid
,
basin_outlet
,
basin_outtp
,
&
&
basin_sz
,
basin_bxout
,
basin_bbout
,
basin_pts
,
basin_lshead
,
coast_pts
,
lontmp
,
lattmp
,
orog_bx
,
&
&
basin_topoindex_stream
,
basin_rlen_stream
,
basin_dz_stream
)
&
area_bx
,
basin_topoindex_stream
,
basin_rlen_stream
,
basin_dz_stream
)
!
USE
ioipsl
USE
grid
...
...
@@ -151,6 +151,7 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
REAL
,
INTENT
(
inout
)
::
rweight_bx
(
nbpt
,
ijdimmax
,
ijdimmax
)
!! Weight of each river within the HTU
REAL
,
INTENT
(
in
)
::
lshead_bx
(
nbpt
,
ijdimmax
,
ijdimmax
)
!! Large scale heading for outflow points.
REAL
,
INTENT
(
in
)
::
orog_bx
(
nbpt
,
ijdimmax
,
ijdimmax
)
!!
REAL
,
INTENT
(
in
)
::
area_bx
(
nbpt
,
ijdimmax
,
ijdimmax
)
!!
!
!
! OUTPUT VARIABLES
...
...
@@ -188,7 +189,7 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
&
rlen_bx
(
ib
,:,:),
rdz_bx
(
ib
,:,:),
rweight_bx
(
ib
,:,:),
lshead_bx
(
ib
,:,:),
diaglalo
,
nb_basin
(
ib
),
&
&
basin_inbxid
(
ib
,:),
basin_outlet
(
ib
,:,:),
basin_outtp
(
ib
,:),
basin_sz
(
ib
,:),
basin_bxout
(
ib
,:),
&
&
basin_bbout
(
ib
,:),
basin_pts
(
ib
,:,:,:),
basin_lshead
(
ib
,:),
coast_pts
(
ib
,:),
lontmp
(
ib
,:,:),
lattmp
(
ib
,:,:),&
&
orog_bx
(
ib
,:,:),
basin_topoindex_stream
(
ib
,:),
basin_rlen_stream
(
ib
,:),
basin_dz_stream
(
ib
,:))
&
orog_bx
(
ib
,:,:),
area_bx
(
ib
,:,:),
basin_topoindex_stream
(
ib
,:),
basin_rlen_stream
(
ib
,:),
basin_dz_stream
(
ib
,:))
ENDDO
END
SUBROUTINE
findbasins
...
...
F90subroutines/routing_reg.f90
View file @
ce3feb48
...
...
@@ -53,7 +53,7 @@ MODULE routing_reg
! topo_option = 6 : Keep and average of all sums of topoindex which are build for all pixels
! within 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.
...
...
@@ -428,8 +428,8 @@ CONTAINS
SUBROUTINE
routing_reg_findbasins
(
nb_htu
,
nbv
,
ib
,
ijdimmax
,
nbi
,
nbj
,
trip
,
basin
,
fac
,
hierarchy
,
topoind
,
&
&
rlen
,
rdz
,
rweight
,
lshead
,
diaglalo
,
nb_basin
,
basin_inbxid
,
basin_outlet
,
basin_outtp
,
basin_sz
,
&
&
basin_bxout
,
basin_bbout
,
basin_pts
,
basin_lshead
,
coast_pts
,
lontmp
,
lattmp
,
orog
,
basin_topoindex_st
rea
m
,
&
&
basin_rlen_stream
,
basin_dz_stream
)
&
basin_bxout
,
basin_bbout
,
basin_pts
,
basin_lshead
,
coast_pts
,
lontmp
,
lattmp
,
orog
,
a
rea
,
&
&
basin_topoindex_stream
,
basin_rlen_stream
,
basin_dz_stream
)
!
IMPLICIT
NONE
!
...
...
@@ -443,6 +443,7 @@ CONTAINS
REAL
(
r_std
),
INTENT
(
in
)
::
fac
(
ijdimmax
,
ijdimmax
)
!!
REAL
(
r_std
),
INTENT
(
in
)
::
lshead
(
ijdimmax
,
ijdimmax
)
REAL
(
r_std
),
INTENT
(
in
)
::
orog
(
ijdimmax
,
ijdimmax
)
REAL
(
r_std
),
INTENT
(
in
)
::
area
(
ijdimmax
,
ijdimmax
)
REAL
(
r_std
),
DIMENSION
(:,:),
INTENT
(
in
)
::
diaglalo
!! Point (in Lat/Lon) where diagnostics will be printed.
!
! Modified
...
...
@@ -757,7 +758,7 @@ CONTAINS
DO
ibas
=
1
,
nbb
fac_out
(
ibas
)
=
fac
(
toutloc
(
ibas
,
1
),
toutloc
(
ibas
,
2
))
END
DO
DO
ibas
=
1
,
5
DO
ibas
=
1
,
4
ff
=
MAXLOC
(
fac_out
)
sortedtrifac
(
ibas
)
=
ff
(
1
)
fac_lim
=
fac_out
(
sortedtrifac
(
ibas
))
...
...
@@ -779,14 +780,15 @@ CONTAINS
!
! In this case nbdiv = 3 only if the upstream or downstream are not too
! small (>1% of the grid)
IF
((
fac_loc_upstr
.GT.
5
)
.AND.
(
fac_loc_main
-
fac_loc_trib
-
fac_loc_upstr
.GT.
5
))
THEN
IF
((
fac_loc_upstr
/
REAL
(
totsz
)
.GT.
1.
/
REAL
(
100
))
.AND.
&
&
((
fac_loc_main
-
fac_loc_trib
-
fac_loc_upstr
)/
REAL
(
totsz
)
.GT.
1.
/
REAL
(
100
)))
THEN
nbdiv
=
3
ELSE
nbdiv
=
2
END
IF
! Perform the operation if the tributary is large enough (global fac)
! and if the HTU is not too small -> more than 4 pixels
IF
((
fac_glo_trib
.GT.
fac_lim
)
.AND.
(
fac_loc_trib
.GT.
5
))
THEN
IF
((
fac_glo_trib
.GT.
fac_lim
)
.AND.
(
fac_loc_trib
/
REAL
(
totsz
)
.GT.
1.
/
REAL
(
100
)
))
THEN
changes
=
changes
+
1
CALL
routing_reg_divbas_cut
(
nb_htu
,
nbv
,
nbi
,
nbj
,
ibas
,
toutloc
(
ibas
,
1
),
toutloc
(
ibas
,
2
),&
&
tsz
(
ibas
),
toutbas
(
ibas
),
toutdir
(
ibas
),
&
...
...
@@ -819,7 +821,7 @@ CONTAINS
nbdiv
=
2
!
! Check with the local fac if the tributary is not too small (>2%)
IF
(
fac_loc_trib
/
REAL
(
totsz
)
.GT.
2
.
/
REAL
(
100
)
)
THEN
IF
(
fac_loc_trib
/
REAL
(
totsz
)
.GT.
1
.
/
REAL
(
100
)
)
THEN
changes
=
changes
+
1
CALL
routing_reg_divbas_cut
(
nb_htu
,
nbv
,
nbi
,
nbj
,
ibas
,
toutloc
(
ibas
,
1
),
toutloc
(
ibas
,
2
),&
&
tsz
(
ibas
),
toutbas
(
ibas
),
toutdir
(
ibas
),
&
...
...
@@ -1019,13 +1021,12 @@ CONTAINS
basin_rlen_stream
(:)
=
0
basin_dz_stream
(:)
=
0
DO
ip
=
1
,
nb_basin
! Another more efficient options would be
! to calculate topoindex_stream from disto
! WRITE(numout,*) "POINT", ib, ip
CALL
mainstream_topoindex
(
nb_htu
,
nbv
,
nbi
,
nbj
,
ip
,
basin_outloc
(
ip
,
1
),
basin_outloc
(
ip
,
2
),&
&
basin_sz
(
ip
),
basin_bbout
(
ip
),
basin_bxout
(
ip
),
basin_lshead
(
ip
),
&
&
basin_pts
,
trip
,
basin
,
fac
,
lontmp
,
lattmp
,
orog
,
rlen
,
rdz
,
&
&
basin_pts
,
trip
,
basin
,
fac
,
lontmp
,
lattmp
,
orog
,
rlen
,
rdz
,
area
,
&
&
basin_topoindex_stream
(
ip
),
basin_rlen_stream
(
ip
),
basin_dz_stream
(
ip
))
IF
(
basin_topoindex_stream
(
ip
)
.LE.
0
)
THEN
IF
(
basin_topoindex_stream
(
ip
)
.LE.
0
)
THEN
! We may have very small value (due to large slope)
! just send alert when its <= 0
WRITE
(
numout
,
*
)
basin_topoindex_stream
(
ip
)
...
...
@@ -1869,19 +1870,6 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
ELSE
toutbas
(
nbb
+
ip
-1
)
=
new_outbas
(
ip
)
ENDIF
! IF ( new_outbas(ip) .NE. undef_int )
!IF (new_nb .EQ. 3) THEN
! If we divide in 3 the second HTU is the upstream main
! it flows in the downstream main at nbb+2
! IF () THEN
! toutbas(nbb+ip-1) = nbb+2
! ELSE
! toutbas(nbb+ip-1) = new_outbas(ip)
! END IF
!ELSE IF (new_nb .EQ. 2) THEN
! toutbas(nbb+ip-1) = new_outbas(ip)
!END IF
!
toutloc
(
nbb
+
ip
-1
,:)
=
new_outloc
(
ip
,:)
tlat
(
nbb
+
ip
-1
)
=
new_lat
(
ip
)
...
...
@@ -1922,7 +1910,7 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
!_ ================================================================================================================================
!
SUBROUTINE
mainstream_topoindex
(
nb_htu
,
nbv
,
nbi
,
nbj
,
ibas
,
iloc
,
jloc
,
tsz
,
tout
,
toutd
,
headd
,
&
&
tpts
,
trip
,
basin
,
fac
,
lon
,
lat
,
orog
,
rlen
,
rdz
,
topoindex_stream
,
distance
,
dorog
)
&
tpts
,
trip
,
basin
,
fac
,
lon
,
lat
,
orog
,
rlen
,
rdz
,
area
,
topoindex_stream
,
distance
,
dorog
)
!
IMPLICIT
NONE
!
...
...
@@ -1941,6 +1929,7 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
INTEGER
(
i_std
),
INTENT
(
in
)
::
basin
(:,:)
!!
REAL
(
r_std
),
INTENT
(
in
)
::
fac
(:,:)
!!
REAL
(
r_std
),
INTENT
(
in
)
::
orog
(:,:),
rlen
(:,:),
rdz
(:,:)
!!
REAL
(
r_std
),
INTENT
(
in
)
::
area
(:,:)
REAL
(
r_std
),
INTENT
(
in
)
::
lon
(:,:),
lat
(:,:)
!!
INTEGER
(
i_std
),
INTENT
(
in
)
::
trip
(:,:)
!!
...
...
@@ -1951,23 +1940,25 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
!
!! LOCAL VARIABLES
REAL
(
r_std
),
PARAMETER
::
flag
=
-9999.
!!
LOGICAL
,
PARAMETER
::
debug
=
.FALSE.
LOGICAL
,
PARAMETER
::
checkgrid
=
.FALSE.
CHARACTER
(
LEN
=
7
)
::
fmt
!!
CHARACTER
(
LEN
=
9
)
::
fmtr
!!
CHARACTER
(
LEN
=
11
)
::
afmt
!!
CHARACTER
(
LEN
=
13
)
::
afmtr
!!
!
REAL
(
r_std
),
DIMENSION
(
nbi
,
nbj
)
::
factmp_glo
!!
REAL
(
r_std
),
DIMENSION
(
nbi
,
nbj
)
::
factmp_glo
REAL
(
r_std
),
DIMENSION
(
nbi
,
nbj
)
::
factmp
!!
REAL
(
r_std
)
::
rel_diff_fac
INTEGER
(
i_std
),
DIMENSION
(
nbi
,
nbj
)
::
triptmp
,
triptemp
!!
!
INTEGER
(
i_std
)
::
il
,
jl
,
ill
,
jll
,
jp
!!
INTEGER
(
i_std
)
::
ie
,
je
!!
INTEGER
(
i_std
)
::
ie
,
je
,
x
,
y
!!
INTEGER
(
i_std
)
::
p
,
cnt
,
k
,
l
,
ic
!!
INTEGER
(
i_std
)
::
ip
,
isz
!!
INTEGER
(
i_std
)
::
ff
(
1
)
!!
!
LOGICAL
::
okpoint
=
.TRUE.
LOGICAL
::
debug
=
.FALSE.
!
! Number of neighbours on the HydroShed grid (regular Lat/Lon)
!
...
...
@@ -1975,11 +1966,15 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
INTEGER
(
i_std
),
DIMENSION
(
nbne
,
2
)
::
inc
!!
!
INTEGER
(
i_std
)
::
main_loc
(
nbv
,
2
),
tri_loc
(
nbne
,
nbv
,
2
)
INTEGER
(
r_std
)
::
main_len
(
nbv
),
main_orog
(
nbv
),
main_dz
(
nbv
)
REAL
(
r_std
)
::
main_fac
(
nbv
),
tri_fac
(
nbne
,
nbv
)
REAL
(
r_std
)
::
main_len
(
nbv
),
main_orog
(
nbv
),
main_dz
(
nbv
)
REAL
(
r_std
)
::
main_fac
(
nbv
),
tri_fac
(
nbne
,
nbv
),
main_area
(
nbv
)
REAL
(
r_std
)
::
main_glofac
(
nbv
),
main_locfac
(
nbv
)
REAL
(
r_std
)
::
tmptri_fac
(
nbne
)
! Flow accumulation of all neighbour points
INTEGER
(
i_std
)
::
sortedtrifac
(
nbne
)
! And sort of tmptri_fac(nbne)
REAL
(
r_std
)
::
alltri_fac
(
nbv
)
! Sort flow accumulation of all tributaries
REAL
(
r_std
)
::
maxpixarea
,
rel_fac
REAL
(
r_std
),
DIMENSION
(
tsz
)
::
fac_glo_pts
,
fac_loc_pts
INTEGER
(
i_std
)
::
cn
,
cbeg
,
cend
!_ ================================================================================================================================
!
...
...
@@ -2018,12 +2013,39 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
triptmp
(
il
,
jl
)
=
trip
(
il
,
jl
)
triptemp
(
il
,
jl
)
=
trip
(
il
,
jl
)
ENDDO
! 1.0 If the first version of this subroutine doesn't work well,
! we need to write the validation step of input arrays here (for
! example, the input sub-basin should only has 1 outlet).
! (October 2016: it still works well :))
!
!
! Calculation of local flow accumulation
factmp
(:,:)
=
0
DO
isz
=
1
,
tsz
ip
=
tpts
(
ibas
,
isz
,
1
)
jp
=
tpts
(
ibas
,
isz
,
2
)
IF
(
trip
(
ip
,
jp
)
.GT.
0
.AND.
trip
(
ip
,
jp
)
.LT.
109
)
THEN
p
=
trip
(
ip
,
jp
)
il
=
ip
;
jl
=
jp
DO
WHILE
(
p
.GT.
0
.AND.
p
.LT.
9
.AND.
cnt
.LT.
nbi
*
nbj
)
ill
=
il
+
inc
(
p
,
1
)
jll
=
jl
+
inc
(
p
,
2
)
il
=
ill
;
jl
=
jll
! Here we compare factmp and factmp glo - factmp glo is a surface !
factmp
(
il
,
jl
)
=
factmp
(
il
,
jl
)
+
1
!area(ip,jp)
p
=
trip
(
il
,
jl
)
ENDDO
ENDIF
END
DO
!
! The inflow is from the max of diff_fac_pts (Jan suggestion)
!maxpixarea = 0
!DO x=1,nbi
! DO y=1,nbj
! IF (area(x,y) .LT. undef_sechiba) THEN
! IF (area(x,y) .GT. maxpixarea) THEN
! maxpixarea = area(x,y)
! END IF
! END IF
! END DO
!END DO
!
IF
(
triptmp
(
iloc
,
jloc
)
.LT.
-4
.OR.
triptmp
(
iloc
,
jloc
)
.GE.
109
)
THEN
CALL
ipslerr_p
(
3
,
'routing_reg_divbas'
,
'The TRIP of outlet is wrong'
,
''
,
''
)
ENDIF
...
...
@@ -2034,44 +2056,51 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
ENDIF
!
! 2.0
!
main_loc
(:,:)
=
0
main_fac
(:)
=
flag
main_len
(:)
=
0
main_orog
(:)
=
0
main_dz
(:)
=
0
tri_loc
(:,:,:)
=
0
main_fac
(:)
=
-1
main_glofac
(:)
=
-1
main_locfac
(:)
=
-1
main_len
(:)
=
-1
main_orog
(:)
=
-1
!main_area(:) = -1
main_dz
(:)
=
-1
tri_loc
(:,:,:)
=
-1
tri_fac
(:,:)
=
flag
!
! Start!
!
! Start!
cnt
=
0
okpoint
=
.TRUE.
il
=
iloc
;
jl
=
jloc
!
DO
WHILE
(
il
.GT.
0
.AND.
il
.LE.
nbi
.AND.
jl
.GT.
0
.AND.
jl
.LE.
nbj
&
&
.AND.
cnt
.LT.
nbi
*
nbj
.AND.
okpoint
)
!DO WHILE ( il .GT. 0 .AND. il .LE. nbi .AND. jl .GT. 0 .AND. jl .LE. nbj &
! & .AND. cnt .LT. nbi*nbj .AND. okpoint )
DO
WHILE
(
cnt
.LT.
nbi
*
nbj
.AND.
okpoint
)
!
! Count the length of the mainstream
!
cnt
=
cnt
+
1
main_loc
(
cnt
,
1
)
=
il
main_loc
(
cnt
,
2
)
=
jl
main_len
(
cnt
)
=
rlen
(
il
,
jl
)
! Vaut toujours 0 ?
main_len
(
cnt
)
=
rlen
(
il
,
jl
)
main_orog
(
cnt
)
=
orog
(
il
,
jl
)
main_dz
(
cnt
)
=
rdz
(
il
,
jl
)
main_fac
(
cnt
)
=
factmp_glo
(
il
,
jl
)
!
! Look for [nbne] neighbour points to find the tributaries
! October 2016: so far, nbne = 8 (for 8 directions of TRIP)
!
main_fac
(
cnt
)
=
factmp_glo
(
il
,
jl
)
-
factmp
(
il
,
jl
)
main_locfac
(
cnt
)
=
factmp
(
il
,
jl
)
main_glofac
(
cnt
)
=
factmp_glo
(
il
,
jl
)
!WRITE(numout,*) "cnt=", cnt, factmp_glo(il,jl)
!main_area(cnt) = area(il,jl)/maxpixarea
!
! Look for [nbne] neighbour points to find the tributaries
! October 2016: so far, nbne = 8 (for 8 directions of TRIP)
!
l
=
0
DO
k
=
1
,
nbne
ill
=
il
+
inc
(
k
,
1
)
jll
=
jl
+
inc
(
k
,
2
)
! If this point is still in the grid box
IF
(
ill
.GT.
0
.AND.
ill
.LE.
nbi
.AND.
jll
.GT.
0
.AND.
jll
.LE.
nbj
)
THEN
! If this point still belongs to the current sub-basin
! If this point still belongs to the current sub-basin
IF
(
triptmp
(
ill
,
jll
)
.GT.
0
.AND.
triptmp
(
ill
,
jll
)
.LE.
nbne
)
THEN
ie
=
ill
+
inc
(
INT
(
triptmp
(
ill
,
jll
)),
1
)
je
=
jll
+
inc
(
INT
(
triptmp
(
ill
,
jll
)),
2
)
...
...
@@ -2089,30 +2118,15 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
ENDIF
ENDDO
!
! Remettre la partie suivante
!
IF
(
l
>
0
)
THEN
!
! Sort the FAC of all neighbour points
tmptri_fac
(:)
=
flag
sortedtrifac
(:)
=
0
!
tmptri_fac
(:)
=
tri_fac
(:,
cnt
)
!tmptrifac -> local
DO
k
=
1
,
l
ff
=
MAXLOC
(
tmptri_fac
)
sortedtrifac
(
k
)
=
ff
(
1
)
tmptri_fac
(
ff
(
1
))
=
flag
ENDDO
! Continue trace upstream: the neighbour point with highest flow
! accumulation will belongs to main stream. We move to this point
! (change value of il, jl).
IF
(
tri_fac
(
sortedtrifac
(
1
),
cnt
)
.NE.
flag
)
THEN
il
=
tri_loc
(
sortedtrifac
(
1
),
cnt
,
1
)
jl
=
tri_loc
(
sortedtrifac
(
1
),
cnt
,
2
)
!
tri_fac
(
sortedtrifac
(
1
),
cnt
)
=
flag
!
ff
=
MAXLOC
(
tmptri_fac
)
k
=
ff
(
1
)
IF
(
tri_fac
(
k
,
cnt
)
.NE.
flag
)
THEN
il
=
tri_loc
(
k
,
cnt
,
1
)
jl
=
tri_loc
(
k
,
cnt
,
2
)
ELSE
okpoint
=
.FALSE.
END
IF
...
...
@@ -2121,21 +2135,47 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
END
IF
ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Main river output
! We stop where the fac_glo-fac_loc is minimum
!ff = MAXLOC(main_fac(:cnt))
cn
=
cnt
IF
(
cnt
.GT.
1
)
THEN
okpoint
=
.TRUE.
k
=
2
! Cas main_facglo almost equal to main_facloc
DO
WHILE
(
okpoint
)
rel_diff_fac
=
(
main_glofac
(
k
-1
)
-
main_glofac
(
k
))/
main_glofac
(
k
)
IF
(
rel_diff_fac
.GT.
0.15
)
THEN
cnt
=
k
-1
okpoint
=
.FALSE.
END
IF
IF
(
k
.EQ.
cnt
)
THEN
okpoint
=
.FALSE.
ELSE
k
=
k
+1
END
IF
END
DO
END
IF
IF
(
debug
)
THEN
WRITE
(
numout
,
*
)
"***"
WRITE
(
numout
,
*
)
"cnt_orig / cnt_new"
,
cn
,
cnt
WRITE
(
numout
,
*
)
"main_glofac"
,
main_glofac
(:
cn
)
WRITE
(
numout
,
*
)
"main_locfac"
,
main_locfac
(:
cn
)
WRITE
(
numout
,
*
)
"main_fac"
,
main_fac
(:
cn
)
WRITE
(
numout
,
*
)
"cnt_orig / cnt_new"
,
cn
,
cnt
WRITE
(
numout
,
*
)
"max/min(orog -> cn_init)"
,
MAXVAL
(
main_orog
(:
cn
)),
MINVAL
(
main_orog
(:
cn
))
WRITE
(
numout
,
*
)
"max/min(orog -> cn_red)"
,
MAXVAL
(
main_orog
(:
cnt
)),
MINVAL
(
main_orog
(:
cnt
))
WRITE
(
numout
,
*
)
"main_orog"
,
main_orog
(:
cn
)
WRITE
(
numout
,
*
)
"main_dz"
,
main_dz
(:
cn
)
WRITE
(
numout
,
*
)
" "
END
IF
il
=
main_loc
(
1
,
1
);
jl
=
main_loc
(
1
,
2
)
! Main river - most upstream pixel
ill
=
main_loc
(
cnt
,
1
);
jll
=
main_loc
(
cnt
,
2
)
IF
(
cnt
.GT.
1
)
THEN
dorog
=
MAXVAL
(
main_orog
(:
cnt
))
-
MINVAL
(
main_orog
(:
cnt
))
dorog
=
MAX
(
0.1
,
dorog
)
distance
=
SUM
(
main_len
(:
cnt
))
ELSE
! if the main river has only one pixel we use directly rdz
dorog
=
MAX
(
0.1
,
rdz
(
il
,
jl
))
distance
=
main_len
(
1
)
END
IF
dorog
=
main_orog
(
cnt
)
-
main_orog
(
1
)
IF
(
cnt
.EQ.
1
)
dorog
=
main_dz
(
1
)
dorog
=
MAX
(
0.1
,
dorog
)
distance
=
SUM
(
main_len
(:
cnt
))
! topoindex is in km while distance and dorog are in meters
topoindex_stream
=
SQRT
(
distance
**
3
/
dorog
)
/
1000
!
...
...
@@ -2473,6 +2513,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!
!
IF
(
topo_option
.EQ.
4
.OR.
topo_option
.EQ.
5
)
THEN
basin_topoind
(
ib
,
ij
)
=
SQRT
(
basin_rlen
(
ib
,
ij
)
**
3.
/(
basin_rdz
(
ib
,
ij
)
*
1.0e6
))
ENDIF
!
...
...
@@ -3643,18 +3684,18 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
!
! The two following conditions are if one basins flows into the other
! we sum the rlen and rdz and recalculate topoindex_stream
IF
(
outflow_basin
(
ib
,
tokill
)
==
totakeover
.AND.
outflow_grid
(
ib
,
tokill
)
==
ib
)
THEN
basin_rlen
(
ib
,
totakeover
)
=
basin_rlen
(
ib
,
totakeover
)
+
basin_rlen
(
ib
,
tokill
)
basin_rdz
(
ib
,
totakeover
)
=
basin_rdz
(
ib
,
totakeover
)
+
basin_rdz
(
ib
,
tokill
)
basin_topoindex_stream
(
ib
,
totakeover
)
=
SQRT
(
basin_rlen
(
ib
,
totakeover
)
**
3
/
basin_rdz
(
ib
,
totakeover
))
/
1000
ELSE
IF
(
outflow_basin
(
ib
,
totakeover
)
==
tokill
.AND.
outflow_grid
(
ib
,
totakeover
)
==
ib
)
THEN
basin_rlen
(
ib
,
totakeover
)
=
basin_rlen
(
ib
,
totakeover
)
+
basin_rlen
(
ib
,
tokill
)
basin_rdz
(
ib
,
totakeover
)
=
basin_rdz
(
ib
,
totakeover
)
+
basin_rdz
(
ib
,
tokill
)
basin_topoindex_stream
(
ib
,
totakeover
)
=
SQRT
(
basin_rlen
(
ib
,
totakeover
)
**
3
/
basin_rdz
(
ib
,
totakeover
))
/
1000
!
IF (outflow_basin(ib,tokill) == totakeover .AND. outflow_grid(ib,tokill) == ib ) THEN
!
basin_rlen(ib, totakeover) = basin_rlen(ib, totakeover) + basin_rlen(ib, tokill)
!
basin_rdz(ib, totakeover) = basin_rdz(ib, totakeover) + basin_rdz(ib,tokill)
!
basin_topoindex_stream(ib, totakeover) = SQRT(basin_rlen(ib, totakeover) **3 / basin_rdz(ib, totakeover)) / 1000
!
ELSE IF (outflow_basin(ib,totakeover) == tokill .AND. outflow_grid(ib,totakeover) == ib ) THEN
!
basin_rlen(ib, totakeover) = basin_rlen(ib, totakeover) + basin_rlen(ib, tokill)
!
basin_rdz(ib, totakeover) = basin_rdz(ib, totakeover) + basin_rdz(ib,tokill)
!
basin_topoindex_stream(ib, totakeover) = SQRT(basin_rlen(ib, totakeover) **3 / basin_rdz(ib, totakeover)) / 1000
! If they don't flow one into the other then take the value of the HTU with
! the higher fetch (naturally it keeps the totakeover value, so just treats
! this case)
ELSE
IF
(
fetch_basin
(
ib
,
tokill
)
.GT.
fetch_basin
(
ib
,
totakeover
))
THEN
IF
(
fetch_basin
(
ib
,
tokill
)
.GT.
fetch_basin
(
ib
,
totakeover
))
THEN
basin_rlen
(
ib
,
totakeover
)
=
basin_rlen
(
ib
,
tokill
)
basin_rdz
(
ib
,
totakeover
)
=
basin_rdz
(
ib
,
tokill
)
basin_topoindex_stream
(
ib
,
totakeover
)
=
basin_topoindex_stream
(
ib
,
tokill
)
...
...
Interface.py
View file @
ce3feb48
...
...
@@ -179,7 +179,8 @@ class HydroSuper :
hierarchy_bx
=
hydrooverlap
.
hierarchy_bx
,
topoind_bx
=
hydrooverlap
.
topoind_bx
,
\
rlen_bx
=
hydrooverlap
.
rlen_bx
,
rdz_bx
=
hydrooverlap
.
rdz_bx
,
\
rweight_bx
=
hydrooverlap
.
rweight_bx
,
lshead_bx
=
hydrooverlap
.
lshead_bx
,
\
lontmp
=
hydrooverlap
.
lon_bx
,
lattmp
=
hydrooverlap
.
lat_bx
,
orog_bx
=
hydrooverlap
.
orog_bx
)
lontmp
=
hydrooverlap
.
lon_bx
,
lattmp
=
hydrooverlap
.
lat_bx
,
orog_bx
=
hydrooverlap
.
orog_bx
,
\
area_bx
=
hydrooverlap
.
area_bx
)
#
# Adjust nwbas to the maximum found over the domain
#
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment