Skip to content

GitLab

Commit 166870ae authored by Anthony's avatar Anthony
Browse files

Perform the triple HTU division over large confluence + rlen / rdz in the... 

Perform the triple HTU division over large confluence + rlen / rdz in the output are now the value for mainstream.
parent 4e87c8f8
Hide whitespace changes
Inline Side-by-side
Showing with 135 additions and 55 deletions
F90subroutines/routing_interface.f90
View file @ 166870ae
...@@ -129,7 +129,7 @@ END SUBROUTINE gethydrogrid ...@@ -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, & 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, & & 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_sz, basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog_bx, &
& basin_topoindex_stream) & basin_topoindex_stream, basin_rlen_stream, basin_dz_stream)
! !
USE ioipsl USE ioipsl
USE grid USE grid
...@@ -165,6 +165,8 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, ...@@ -165,6 +165,8 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
REAL, INTENT(out) :: basin_lshead(nbpt,nb_htu) !! REAL, INTENT(out) :: basin_lshead(nbpt,nb_htu) !!
INTEGER, INTENT(out) :: coast_pts(nbpt,nb_htu) !! The coastal flow points (unitless) INTEGER, INTENT(out) :: coast_pts(nbpt,nb_htu) !! The coastal flow points (unitless)
REAL, INTENT(out) :: basin_topoindex_stream(nbpt,nb_htu) !! REAL, INTENT(out) :: basin_topoindex_stream(nbpt,nb_htu) !!
REAL, INTENT(out) :: basin_rlen_stream(nbpt,nb_htu) !!
REAL, INTENT(out) :: basin_dz_stream(nbpt,nb_htu) !!
! !
! For debug and get coordinate of river outlet ! For debug and get coordinate of river outlet
! !
...@@ -186,7 +188,7 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, ...@@ -186,7 +188,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), & & 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_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,:,:),& & basin_bbout(ib,:), basin_pts(ib,:,:,:), basin_lshead(ib,:), coast_pts(ib,:), lontmp(ib,:,:), lattmp(ib,:,:),&
& orog_bx(ib,:,:), basin_topoindex_stream(ib,:)) & orog_bx(ib,:,:), basin_topoindex_stream(ib,:), basin_rlen_stream(ib,:), basin_dz_stream(ib,:))
ENDDO ENDDO
END SUBROUTINE findbasins END SUBROUTINE findbasins
...@@ -194,7 +196,8 @@ END SUBROUTINE findbasins ...@@ -194,7 +196,8 @@ END SUBROUTINE findbasins
SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, & SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, &
& fac_bx, topoind_bx, rlen_bx, rdz_bx, rweight_bx, & & fac_bx, topoind_bx, rlen_bx, rdz_bx, rweight_bx, &
& min_topoind, nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_pts, basin_bxout, & & 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_bbout, lshead, coast_pts, nwbas, basin_rlen_stream, basin_dz_stream, basin_count, &
& basin_notrun, basin_area, basin_cg, basin_hierarchy, &
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac, basin_topoind, & & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac, basin_topoind, &
& basin_rlen, basin_rdz, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, & & basin_rlen, basin_rdz, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
& basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin) & basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
...@@ -231,6 +234,7 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_ ...@@ -231,6 +234,7 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: basin_bbout !! outflow sub-basin INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: basin_bbout !! outflow sub-basin
REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: lshead !! Large scale heading of outflow. REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: lshead !! Large scale heading of outflow.
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: coast_pts !! The coastal flow points (unitless) INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: coast_pts !! The coastal flow points (unitless)
REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: basin_rlen_stream, basin_dz_stream
!! Output VARIABLES !! Output VARIABLES
INTEGER(i_std), INTENT (in) :: nwbas !! INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), INTENT (out), DIMENSION(nbpt) :: basin_count !! INTEGER(i_std), INTENT (out), DIMENSION(nbpt) :: basin_count !!
...@@ -270,6 +274,10 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_ ...@@ -270,6 +274,10 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac,& & basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac,&
& basin_topoind, basin_rlen, basin_rdz, basin_id, basin_coor, basin_type, & & basin_topoind, basin_rlen, basin_rdz, basin_id, basin_coor, basin_type, &
& basin_flowdir, basin_lshead, basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin) & basin_flowdir, basin_lshead, basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
! Once all possible operation on topoindex have been performed, we keep
! rlen_stream and rdz_stream as rlen and rdz
basin_rlen(ib,:) = basin_rlen_stream(ib,:nwbas)
basin_rdz(ib,:) = basin_dz_stream(ib,:nwbas)
ENDDO ENDDO
! !
END SUBROUTINE globalize END SUBROUTINE globalize
......
F90subroutines/routing_reg.f90
View file @ 166870ae
...@@ -428,7 +428,8 @@ CONTAINS ...@@ -428,7 +428,8 @@ CONTAINS
SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, trip, basin, fac, hierarchy, topoind, & 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, & & 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_stream ) & basin_bxout, basin_bbout, basin_pts, basin_lshead, coast_pts, lontmp, lattmp, orog, basin_topoindex_stream, &
& basin_rlen_stream, basin_dz_stream )
! !
IMPLICIT NONE IMPLICIT NONE
! !
...@@ -471,6 +472,8 @@ CONTAINS ...@@ -471,6 +472,8 @@ CONTAINS
INTEGER(i_std), INTENT(out) :: basin_pts(nb_htu, nbv, 2) !! INTEGER(i_std), INTENT(out) :: basin_pts(nb_htu, nbv, 2) !!
INTEGER(i_std), INTENT(out) :: coast_pts(nb_htu) !! The coastal flow points (unitless) INTEGER(i_std), INTENT(out) :: coast_pts(nb_htu) !! The coastal flow points (unitless)
REAL(r_std), INTENT(out) :: basin_topoindex_stream(nb_htu) !! REAL(r_std), INTENT(out) :: basin_topoindex_stream(nb_htu) !!
REAL(r_std), INTENT(out) :: basin_rlen_stream(nb_htu) !!
REAL(r_std), INTENT(out) :: basin_dz_stream(nb_htu) !!
! !
!! LOCAL VARIABLES !! LOCAL VARIABLES
LOGICAL, PARAMETER :: debug=.FALSE. LOGICAL, PARAMETER :: debug=.FALSE.
...@@ -526,8 +529,9 @@ CONTAINS ...@@ -526,8 +529,9 @@ CONTAINS
INTEGER(i_std) :: sortedtrifac(nbne) ! INTEGER(i_std) :: sortedtrifac(nbne) !
REAL(r_std) :: fac_lim !! REAL(r_std) :: fac_lim !!
REAL(r_std), PARAMETER :: flag=-9999. !! REAL(r_std), PARAMETER :: flag=-9999. !!
INTEGER(i_std) :: option, il_trib, jl_trib, changes !! INTEGER(i_std) :: option, il_trib, jl_trib, changes !!
REAL(r_std) :: fac_glo_trib, fac_loc_trib, fac_loc_main INTEGER(i_std) :: nbdiv, il_upstr, jl_upstr
REAL(r_std) :: fac_glo_trib, fac_loc_trib, fac_loc_main, fac_loc_upstr
! !
!_ ================================================================================================================================ !_ ================================================================================================================================
! !
...@@ -770,15 +774,24 @@ CONTAINS ...@@ -770,15 +774,24 @@ CONTAINS
CALL routing_reg_divbas_search(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),& CALL routing_reg_divbas_search(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),&
& tsz(ibas), toutbas(ibas), toutdir(ibas), & & tsz(ibas), toutbas(ibas), toutdir(ibas), &
& toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, & & toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, &
& option, il_trib, jl_trib, fac_glo_trib, fac_loc_trib, fac_loc_main) & option, il_trib, jl_trib, fac_glo_trib, fac_loc_trib, fac_loc_main, &
& il_upstr, jl_upstr, fac_loc_upstr)
! !
! 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
nbdiv = 3
ELSE
nbdiv = 2
END IF
! Perform the operation if the tributary is large enough (global fac) ! Perform the operation if the tributary is large enough (global fac)
! and if the HTU is not too small -> more than 4 pixels ! and if the HTU is not too small -> more than 4 pixels
IF (( fac_glo_trib .GT. fac_lim ) .AND. ( fac_loc_trib .GT. 4 )) THEN IF (( fac_glo_trib .GT. fac_lim ) .AND. ( fac_loc_trib .GT. 5 )) THEN
changes = changes + 1 changes = changes + 1
CALL routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),& CALL routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),&
& tsz(ibas), toutbas(ibas), toutdir(ibas), & & tsz(ibas), toutbas(ibas), toutdir(ibas), &
& toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, il_trib, jl_trib, & & toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, &
& il_trib, jl_trib, nbdiv, il_upstr, jl_upstr, &
& new_nb, mp, new_bname, new_outdir, new_heading, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts) & new_nb, mp, new_bname, new_outdir, new_heading, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts)
! !
CALL routing_reg_divbas_divide(nb_htu, nbv, nbi, nbj, ijdimmax, tbname, tsz, tpts, toutdir, toutbas, & CALL routing_reg_divbas_divide(nb_htu, nbv, nbi, nbj, ijdimmax, tbname, tsz, tpts, toutdir, toutbas, &
...@@ -798,14 +811,20 @@ CONTAINS ...@@ -798,14 +811,20 @@ CONTAINS
CALL routing_reg_divbas_search(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),& CALL routing_reg_divbas_search(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),&
& tsz(ibas), toutbas(ibas), toutdir(ibas), & & tsz(ibas), toutbas(ibas), toutdir(ibas), &
& toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, & & toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, &
& option, il_trib, jl_trib, fac_glo_trib, fac_loc_trib,fac_loc_main) & option, il_trib, jl_trib, fac_glo_trib, fac_loc_trib,fac_loc_main, &
& il_upstr, jl_upstr, fac_loc_upstr)
!
! With this option nbdiv is always equal to 2 (no separtion between
! upstream / downstream of the HTU in rel. to the confluence).
nbdiv = 2
! !
! Check with the local fac if the tributary is not too small (>1%) ! Check with the local fac if the tributary is not too small (>2%)
IF ( fac_loc_trib / REAL(totsz) .GT. 1./REAL(100) ) THEN IF ( fac_loc_trib / REAL(totsz) .GT. 2./REAL(100) ) THEN
changes = changes + 1 changes = changes + 1
CALL routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),& CALL routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas,toutloc(ibas,1), toutloc(ibas,2),&
& tsz(ibas), toutbas(ibas), toutdir(ibas), & & tsz(ibas), toutbas(ibas), toutdir(ibas), &
& toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, il_trib, jl_trib, & & toutlshead(ibas), tpts, trip, basin, fac, lontmp, lattmp, &
& il_trib, jl_trib, nbdiv, il_upstr, jl_upstr, &
& new_nb, mp, new_bname, new_outdir, new_heading, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts) & new_nb, mp, new_bname, new_outdir, new_heading, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts)
CALL routing_reg_divbas_divide(nb_htu, nbv, nbi, nbj, ijdimmax, tbname, tsz, tpts, toutdir, toutbas, & CALL routing_reg_divbas_divide(nb_htu, nbv, nbi, nbj, ijdimmax, tbname, tsz, tpts, toutdir, toutbas, &
...@@ -997,14 +1016,22 @@ CONTAINS ...@@ -997,14 +1016,22 @@ CONTAINS
! because rlen/rdz may be changed ! because rlen/rdz may be changed
! !
basin_topoindex_stream(:) = 0 basin_topoindex_stream(:) = 0
basin_rlen_stream(:) = 0
basin_dz_stream(:) = 0
DO ip=1,nb_basin DO ip=1,nb_basin
! Another more efficient options would be ! Another more efficient options would be
! to calculate topoindex_stream from disto ! to calculate topoindex_stream from disto
CALL mainstream_topoindex(nb_htu, nbv, nbi, nbj, ip, basin_outloc(ip,1), basin_outloc(ip,2),& 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_sz(ip), basin_bbout(ip), basin_bxout(ip), basin_lshead(ip), &
& basin_pts, trip, basin, fac, lontmp, lattmp, orog, rlen, rdz, basin_topoindex_stream(ip)) & basin_pts, trip, basin, fac, lontmp, lattmp, orog, rlen, rdz, &
IF ( basin_topoindex_stream(ip) .LT. 1 ) THEN & basin_topoindex_stream(ip), basin_rlen_stream(ip), basin_dz_stream(ip))
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) WRITE(numout,*) basin_topoindex_stream(ip)
WRITE(numout,*) basin_rlen_stream(ip)
WRITE(numout,*) basin_sz(ip)
WRITE(numout,*) basin_dz_stream(ip)
CALL ipslerr_p(3,'routing_reg_findbasins','Error in the mainstream topoindex','','') CALL ipslerr_p(3,'routing_reg_findbasins','Error in the mainstream topoindex','','')
END IF END IF
END DO END DO
...@@ -1115,7 +1142,8 @@ CONTAINS ...@@ -1115,7 +1142,8 @@ CONTAINS
!_ ================================================================================================================================ !_ ================================================================================================================================
SUBROUTINE routing_reg_divbas_search(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, & SUBROUTINE routing_reg_divbas_search(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, &
& tpts, trip, basin, fac, lon, lat, option, il_trib, jl_trib, fac_glo_trib, fac_loc_trib, fac_loc_main) & tpts, trip, basin, fac, lon, lat, option, il_trib, jl_trib, fac_glo_trib, fac_loc_trib, fac_loc_main,&
& il_upstr, jl_upstr, fac_loc_upstr)
! !
IMPLICIT NONE IMPLICIT NONE
! !
...@@ -1139,7 +1167,9 @@ CONTAINS ...@@ -1139,7 +1167,9 @@ CONTAINS
! !
!! OUTPUT VARIABLES !! OUTPUT VARIABLES
INTEGER(i_std), INTENT(out) :: il_trib, jl_trib !! INTEGER(i_std), INTENT(out) :: il_trib, jl_trib !!
INTEGER(i_std), INTENT(out) :: il_upstr, jl_upstr !!
REAL(r_std), INTENT(out) :: fac_glo_trib, fac_loc_trib, fac_loc_main !! REAL(r_std), INTENT(out) :: fac_glo_trib, fac_loc_trib, fac_loc_main !!
REAL(r_std), INTENT(out) :: fac_loc_upstr
! !
!! LOCAL VARIABLES !! LOCAL VARIABLES
...@@ -1394,6 +1424,9 @@ CONTAINS ...@@ -1394,6 +1424,9 @@ CONTAINS
jl_trib = tri_loc(sortedtrifac(2),cnt,2) jl_trib = tri_loc(sortedtrifac(2),cnt,2)
fac_glo_trib = factmp_glo(il_trib, jl_trib) !tri_fac(sortedtrifac(2),cnt) fac_glo_trib = factmp_glo(il_trib, jl_trib) !tri_fac(sortedtrifac(2),cnt)
fac_loc_trib = factmp_loc(il_trib, jl_trib) fac_loc_trib = factmp_loc(il_trib, jl_trib)
il_upstr = il
jl_upstr = jl
fac_loc_upstr = factmp_loc(il, jl)
END IF END IF
END IF END IF
ELSE IF ( option .EQ. 2) THEN ELSE IF ( option .EQ. 2) THEN
...@@ -1410,6 +1443,11 @@ CONTAINS ...@@ -1410,6 +1443,11 @@ CONTAINS
jl_trib = tri_loc(sortedtrifac(2),cnt,2) jl_trib = tri_loc(sortedtrifac(2),cnt,2)
fac_glo_trib = factmp_glo(il_trib, jl_trib) !tri_fac(sortedtrifac(2),cnt) fac_glo_trib = factmp_glo(il_trib, jl_trib) !tri_fac(sortedtrifac(2),cnt)
fac_loc_trib = factmp_loc(il_trib, jl_trib) fac_loc_trib = factmp_loc(il_trib, jl_trib)
! Is not useful in this option, but this helps keep the
! structure of the subroutine
il_upstr = il
jl_upstr = jl
fac_loc_upstr = factmp_loc(il, jl)
END IF END IF
END IF END IF
END IF END IF
...@@ -1443,7 +1481,7 @@ END SUBROUTINE routing_reg_divbas_search ...@@ -1443,7 +1481,7 @@ END SUBROUTINE routing_reg_divbas_search
!_ ================================================================================================================================ !_ ================================================================================================================================
SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, & SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, tout, toutd, headd, &
& tpts, trip, basin, fac, lon, lat, il_trib, jl_trib, & & tpts, trip, basin, fac, lon, lat, il_trib, jl_trib, nbdiv, il_upstr, jl_upstr, &
& new_nb, oic, new_bname, new_outdir, new_head, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts) & new_nb, oic, new_bname, new_outdir, new_head, new_outbas, new_lat, new_lon, new_outloc, new_sz, new_pts)
! !
IMPLICIT NONE IMPLICIT NONE
...@@ -1464,6 +1502,8 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, ...@@ -1464,6 +1502,8 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz,
REAL(r_std), INTENT(in) :: fac(:,:) !! REAL(r_std), INTENT(in) :: fac(:,:) !!
REAL(r_std), INTENT(in) :: lon(:,:), lat(:,:) !! REAL(r_std), INTENT(in) :: lon(:,:), lat(:,:) !!
INTEGER(i_std), INTENT(in) :: il_trib, jl_trib !! INTEGER(i_std), INTENT(in) :: il_trib, jl_trib !!
INTEGER(i_std), INTENT(in) :: nbdiv
INTEGER(i_std), INTENT(in) :: il_upstr, jl_upstr !!
! !
!! MODIFIED VARIABLES !! MODIFIED VARIABLES
INTEGER(i_std), INTENT(inout) :: trip(:,:) !! INTEGER(i_std), INTENT(inout) :: trip(:,:) !!
...@@ -1591,19 +1631,35 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, ...@@ -1591,19 +1631,35 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz,
divloc(1,1) = il_trib divloc(1,1) = il_trib
divloc(1,2) = jl_trib divloc(1,2) = jl_trib
divloc(2,1) = iloc IF (nbdiv .EQ. 3) THEN
divloc(2,2) = jloc ! Downstream main tributary must come last
oic = 2 divloc(2,1) = il_upstr
! divloc(2,2) = jl_upstr
new_outbas(1) = 2 divloc(3,1) = iloc
new_outbas(2) = tout divloc(3,2) = jloc
ELSE IF (nbdiv .EQ. 2) THEN
divloc(2,1) = iloc
divloc(2,2) = jloc
END IF
!
! For what I understand oic is the index in divloc of the output HTU
IF (nbdiv .EQ. 2) THEN
oic = 2
new_outbas(1) = 2
new_outbas(2) = tout
ELSE IF (nbdiv .EQ. 3) THEN
oic = 3
new_outbas(1) = 3
new_outbas(2) = 3
new_outbas(3) = tout
END IF
! !
!IF ( debug ) WRITE (numout,*) 'Number of new sub-basin: ', l, oic !IF ( debug ) WRITE (numout,*) 'Number of new sub-basin: ', l, oic
! Here l always = 1 so we may simplify some variables ! Here l always = 1 so we may simplify some variables
! !
! Now, cut it ! Cut the sub-basin ! Release The Kraken ! ! Now, cut it ! Cut the sub-basin ! Release The Kraken !
! !
new_nb = 2 new_nb = nbdiv
! !
DO k = 1, new_nb DO k = 1, new_nb
! We need to start by the second (the tributary) elsewhere all the points goes to main ! We need to start by the second (the tributary) elsewhere all the points goes to main
...@@ -1669,7 +1725,7 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz, ...@@ -1669,7 +1725,7 @@ SUBROUTINE routing_reg_divbas_cut(nb_htu, nbv, nbi, nbj, ibas, iloc, jloc, tsz,
ENDDO ENDDO
! !
IF ( checksz .NE. tsz) THEN IF ( checksz .NE. tsz) THEN
WRITE(numout,*) ' Water got lost while I tried to divide basins' WRITE(numout,*) ' ater got lost while I tried to divide basins'
WRITE(numout,*) ' Number of new sub-basin : ', new_nb WRITE(numout,*) ' Number of new sub-basin : ', new_nb
ip = 2 !COUNT(divloc(:,1) .NE. 0) ip = 2 !COUNT(divloc(:,1) .NE. 0)
WRITE(numout,*) ' Number of mainstrem sub-basin : ', oic WRITE(numout,*) ' Number of mainstrem sub-basin : ', oic
...@@ -1774,17 +1830,13 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz, ...@@ -1774,17 +1830,13 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
! !
! toutbas define basin_bbout which is the predefined basin outflow (due to outflow in same grid) ! toutbas define basin_bbout which is the predefined basin outflow (due to outflow in same grid)
! Here we cut in two we just have to change the toutbas for the new basin ! Here we cut in two we just have to change the toutbas for the new basin
!IF ( new_outbas(1) .NE. undef_int ) THEN IF ( new_outbas(1) .NE. undef_int ) THEN
! IF ( mp .NE. 1 ) THEN IF ( mp .NE. 1 ) THEN
! toutbas(ibas) = new_outbas(1)+nbb-1 toutbas(ibas) = new_outbas(1)+nbb-1
! ELSE ELSE
! toutbas(ibas) = new_outbas(1) toutbas(ibas) = new_outbas(1)
! ENDIF ENDIF
!ENDIF ENDIF
! From changes : first one is always the tributary -> goes to the main HTU
! which is the new htu generated at nbb+1
toutbas(ibas) = nbb+1
! !
toutloc(ibas,:) = new_outloc(1,:) toutloc(ibas,:) = new_outloc(1,:)
tlat(ibas) = new_lat(1) tlat(ibas) = new_lat(1)
...@@ -1810,14 +1862,26 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz, ...@@ -1810,14 +1862,26 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
! Trung: Should come back here soon (small bug can come from ! Trung: Should come back here soon (small bug can come from
! Trung: this point) ! Trung: this point)
! !
!IF ( new_outbas(ip) .EQ. 1 .AND. ip .NE. mp ) THEN IF ( new_outbas(ip) .EQ. 1 .AND. ip .NE. mp ) THEN
! toutbas(nbb+ip-1) = ibas toutbas(nbb+ip-1) = ibas
!ELSE IF ( ip .NE. mp ) THEN ELSE IF ( ip .NE. mp ) THEN
! toutbas(nbb+ip-1) = new_outbas(ip)+nbb-1 toutbas(nbb+ip-1) = new_outbas(ip)+nbb-1
!ELSE 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) ! toutbas(nbb+ip-1) = new_outbas(ip)
!ENDIF !END IF
toutbas(nbb+ip-1) = new_outbas(ip)
! !
toutloc(nbb+ip-1,:) = new_outloc(ip,:) toutloc(nbb+ip-1,:) = new_outloc(ip,:)
tlat(nbb+ip-1) = new_lat(ip) tlat(nbb+ip-1) = new_lat(ip)
...@@ -1858,7 +1922,7 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz, ...@@ -1858,7 +1922,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, & 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) & tpts, trip, basin, fac, lon, lat, orog, rlen, rdz, topoindex_stream, distance, dorog)
! !
IMPLICIT NONE IMPLICIT NONE
! !
...@@ -1883,6 +1947,7 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz, ...@@ -1883,6 +1947,7 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
! !
!! OUTPUT VARIABLES !! OUTPUT VARIABLES
REAL(r_std), INTENT(out) :: topoindex_stream !! REAL(r_std), INTENT(out) :: topoindex_stream !!
REAL(r_std), INTENT(out) :: distance, dorog !!
! !
!! LOCAL VARIABLES !! LOCAL VARIABLES
REAL(r_std), PARAMETER :: flag=-9999. !! REAL(r_std), PARAMETER :: flag=-9999. !!
...@@ -1894,7 +1959,6 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz, ...@@ -1894,7 +1959,6 @@ SUBROUTINE routing_reg_divbas_divide(nb_htu, nbv,nbi,nbj, ijdimmax, tbname, tsz,
CHARACTER(LEN=13) :: afmtr !! CHARACTER(LEN=13) :: afmtr !!
! !
REAL(r_std), DIMENSION(nbi,nbj) :: factmp_glo !! REAL(r_std), DIMENSION(nbi,nbj) :: factmp_glo !!
REAL(r_std) :: distance, dorog !!
INTEGER(i_std), DIMENSION(nbi,nbj) :: triptmp, triptemp !! INTEGER(i_std), DIMENSION(nbi,nbj) :: triptmp, triptemp !!
! !
INTEGER(i_std) :: il, jl, ill, jll, jp !! INTEGER(i_std) :: il, jl, ill, jll, jp !!
...@@ -3547,8 +3611,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b ...@@ -3547,8 +3611,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) & 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_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_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 !basin_rdz(ib, totakeover) = (basin_rdz(ib, totakeover)+basin_rdz(ib, tokill))*0.5
! !
ELSE IF ( kill_option .EQ. 2 ) THEN ELSE IF ( kill_option .EQ. 2 ) THEN
IF (outflow_basin(ib,tokill) == totakeover .OR. outflow_basin(ib,totakeover) == tokill ) THEN IF (outflow_basin(ib,tokill) == totakeover .OR. outflow_basin(ib,totakeover) == tokill ) THEN
...@@ -3556,23 +3620,30 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b ...@@ -3556,23 +3620,30 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
IF ( ABS(basin_orog_mean(ib,tokill)-basin_orog_mean(ib,totakeover)) > meandz ) THEN IF ( ABS(basin_orog_mean(ib,tokill)-basin_orog_mean(ib,totakeover)) > meandz ) THEN
! HTUs flow into each other ! HTUs flow into each other
basin_topoind(ib, totakeover) = basin_topoind(ib, totakeover)+basin_topoind(ib, tokill) 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_rlen(ib, totakeover) = basin_rlen(ib, totakeover)+basin_rlen(ib, tokill)
basin_rdz(ib, totakeover) = basin_rdz(ib, totakeover)+basin_rdz(ib, tokill) !basin_rdz(ib, totakeover) = basin_rdz(ib, totakeover)+basin_rdz(ib, tokill)
ELSE ELSE
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) & 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_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_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 !basin_rdz(ib, totakeover) = (basin_rdz(ib, totakeover)+basin_rdz(ib, tokill))*0.5
ENDIF ENDIF
ELSE ELSE
basin_topoind(ib, totakeover) = (basin_topoind(ib, totakeover)*basin_area(ib, totakeover) & 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_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_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 !basin_rdz(ib, totakeover) = (basin_rdz(ib, totakeover)+basin_rdz(ib, tokill))*0.5
ENDIF ENDIF
ELSE ELSE
CALL ipslerr_p(3,'routing_reg_killbas','The selected kill_option does not exist.','','') CALL ipslerr_p(3,'routing_reg_killbas','The selected kill_option does not exist.','','')
ENDIF ENDIF
! Keep the stream rlen / rdz
! See if we keep one only rlen / rdz corresponding to the stream