Skip to content

GitLab

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

The river segment length and river segment elevation change per HTU is now... 

The river segment length and river segment elevation change per HTU is now carried until routing_reg_globalize. This allows to explore other methods of computing the topoindex.
parent 3b997f43
Hide whitespace changes
Inline Side-by-side
Showing with 187 additions and 88 deletions
F90subroutines/Makefile
View file @ c8b9a77d
......@@ -11,7 +11,7 @@ FPY = f2py
#
EXT_SUFFIX := $(shell python3-config --extension-suffix)
#
all : routing_interface.so diagst.so
all : routing_interface$(EXT_SUFFIX) diagst$(EXT_SUFFIX)
routing_interface : routing_interface$(EXT_SUFFIX)
diagst : diagst$(EXT_SUFFIX)
......
F90subroutines/routing_interface.f90
View file @ c8b9a77d
SUBROUTINE initatmgrid(rank, nbcore, nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
SUBROUTINE initatmgrid(rankmpi, nbcore, nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_in)
!
USE grid
USE ioipsl
......@@ -7,7 +7,7 @@ SUBROUTINE initatmgrid(rank, nbcore, nbpt, nbseg, polygons_in, centre_in, area_i
IMPLICIT NONE
!
!! INPUT VARIABLES
INTEGER, INTENT(in) :: rank !! rank of processor
INTEGER, INTENT(in) :: rankmpi !! rank of processor
INTEGER, INTENT(in) :: nbcore
INTEGER, INTENT(in) :: nbpt !! Atmospheric Domain size (unitless)
INTEGER, INTENT(in) :: nbseg !! Number of segments in the polygone
......@@ -29,8 +29,8 @@ SUBROUTINE initatmgrid(rank, nbcore, nbpt, nbseg, polygons_in, centre_in, area_i
!
CALL grid_init(nbpt, nbseg, polygons_in, centre_in, area_in, contfrac_in, neighbours_loc)
!
numout=100+rank
WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",rank,"_from_",nbcore,".txt"
numout=100+rankmpi
WRITE(outfname,'(A15,I4.4,A6,I4.4,A4)') "Out_RoutingReg_",rankmpi,"_from_",nbcore,".txt"
OPEN(unit=numout, file=outfname)
!
END SUBROUTINE initatmgrid
......@@ -44,8 +44,9 @@ SUBROUTINE closeinterface()
END SUBROUTINE closeinterface
SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_index, sub_area, max_basins, &
& min_topoind, lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy,orog,floodp, nbi, nbj, area_bx, &
& trip_bx, basin_bx, topoind_bx, fac_bx, hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
& min_topoind, meanrlen, meanrdz, 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, &
& hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
!
USE ioipsl
USE grid
......@@ -60,7 +61,8 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
INTEGER, INTENT(in) :: sub_pts(nbpt) !! Number of high resolution points on this 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(in) :: min_topoind !! The current minimum of topographic index (m)
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) :: sub_area(nbpt,nbvmax_in) !! Area on the fine grid (m^2)
!
REAL, INTENT(in) :: lon_rel(nbpt,nbvmax_in) !!
......@@ -69,6 +71,8 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
REAL, INTENT(inout) :: trip(nbpt,nbvmax_in) !! The trip field (unitless)
REAL, INTENT(inout) :: basins(nbpt,nbvmax_in) !! data on the fine grid
REAL, INTENT(inout) :: topoindex(nbpt,nbvmax_in) !! Topographic index of the residence time (m)
REAL, INTENT(inout) :: rlen(nbpt,nbvmax_in) !!
REAL, INTENT(inout) :: rdz(nbpt,nbvmax_in) !!
REAL, INTENT(inout) :: hierarchy(nbpt,nbvmax_in) !! data on the fine grid
REAL, INTENT(inout) :: fac(nbpt,nbvmax_in) !! data on the fine grid
REAL, INTENT(inout) :: orog(nbpt,nbvmax_in) !! data on the fine grid
......@@ -83,6 +87,9 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
REAL, INTENT(out) :: lat_bx(nbpt,ijdimmax,ijdimmax) !!
REAL, INTENT(out) :: lshead_bx(nbpt,ijdimmax,ijdimmax) !! Large scale heading for outflow points.
REAL, INTENT(out) :: topoind_bx(nbpt,ijdimmax,ijdimmax) !! Topographic index of the residence time for each of the smaller boxes (m)
REAL, INTENT(out) :: rlen_bx(nbpt,ijdimmax,ijdimmax) !!
REAL, INTENT(out) :: rdz_bx(nbpt,ijdimmax,ijdimmax) !!
REAL, INTENT(out) :: rweight_bx(nbpt,ijdimmax,ijdimmax) !!
INTEGER, INTENT(out) :: trip_bx(nbpt,ijdimmax,ijdimmax) !! The trip field for each of the smaller boxes (unitless)
INTEGER, INTENT(out) :: basin_bx(nbpt,ijdimmax,ijdimmax) !!
REAL, INTENT(out) :: orog_bx(nbpt,ijdimmax,ijdimmax) !! Orography (m)
......@@ -107,18 +114,21 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, maxpercent_in, sub_pts, sub_i
DO ib=1,nbpt
CALL routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, fac, hierarchy, &
& orog, floodp, nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), basin_bx(ib,:,:), &
& topoind_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:),orog_bx(ib,:,:),floodp_bx(ib,:,:),&
& lon_bx(ib,:,:), lat_bx(ib,:,:), lshead_bx(ib,:,:))
& meanrlen, meanrdz, 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,:,:), &
& 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,:,:), &
& lshead_bx(ib,:,:))
ENDDO
!
rweight_bx(:,:,:) = zero
!
END SUBROUTINE gethydrogrid
SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_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)
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)
!
USE ioipsl
USE grid
......@@ -135,11 +145,14 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
REAL, INTENT(in) :: fac_bx(nbpt,ijdimmax,ijdimmax) !! Flow accumulation
REAL, INTENT(in) :: hierarchy_bx(nbpt,ijdimmax,ijdimmax) !! Level in the basin of the point
REAL, INTENT(inout) :: topoind_bx(nbpt,ijdimmax,ijdimmax) !! Topographic index of the residence time for each of the smaller boxes (m)
REAL, INTENT(inout) :: rlen_bx(nbpt,ijdimmax,ijdimmax) !!
REAL, INTENT(inout) :: rdz_bx(nbpt,ijdimmax,ijdimmax) !!
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.
!
!
! OUTPUT VARIABLES
INTEGER, INTENT(out) :: nb_basin(nbpt) !! Number of sub-basins (unitless)
INTEGER, INTENT(out) :: nb_basin(nbpt) !! Number of sub-basins (unitless)
INTEGER, INTENT(out) :: basin_inbxid(nbpt,nb_htu) !!
REAL, INTENT(out) :: basin_outlet(nbpt,nb_htu,2) !!
REAL, INTENT(out) :: basin_outtp(nbpt,nb_htu) !!
......@@ -166,20 +179,20 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, ijdimmax, nbi, nbj, trip_bx, basin_bx,
DO ib=1,nbpt
CALL routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi(ib), nbj(ib), trip_bx(ib,:,:), &
& basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), &
& topoind_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_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), topoind_bx(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_bbout(ib,:), basin_pts(ib,:,:,:), basin_lshead(ib,:), coast_pts(ib,:), lontmp(ib,:,:), lattmp(ib,:,:))
ENDDO
END SUBROUTINE findbasins
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, &
& 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, &
& basin_bbout, lshead, coast_pts, nwbas, 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_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
& 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_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
!
USE ioipsl
......@@ -199,10 +212,13 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: hierarchy_bx !! Level in the basin of the point
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: orog_bx !! Orography (m)
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: floodp_bx !! Floodplains area (m^2)
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: fac_bx !! Level in the basin of the point
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: topoind_bx !! Topographic index of the residence time for each of the smaller boxes (m)
REAL(r_std), INTENT (in) :: min_topoind !! The current minimum of topographic index (m)
INTEGER(i_std), INTENT (in), DIMENSION(nbpt) :: nb_basin !! Number of sub-basins (unitless)
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: fac_bx !! Level in the basin of the point
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: topoind_bx !! Topographic index of the residence time for each of the smaller boxes (km)
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: rlen_bx !!
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: rdz_bx !!
REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: rweight_bx !!
REAL(r_std), INTENT (in) :: min_topoind !! The current minimum of topographic index (km)
INTEGER(i_std), INTENT (in), DIMENSION(nbpt) :: nb_basin !! Number of sub-basins (unitless)
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu,2) :: basin_outlet !! Outlet coordinate of each subgrid basin
REAL(r_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: basin_outtp !!
......@@ -227,7 +243,9 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_orog_max !!
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_floodp !!
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_fac !!
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (km)
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_rlen
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_rdz
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_lshead !! Large scale heading out of the grid box.
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_beta_fp !! Beta parameter for floodplains
INTEGER(i_std), INTENT (out), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
......@@ -241,12 +259,13 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, ijdimmax, area_bx, lon_bx, lat_bx, trip_
CALL routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, area_bx(ib,:,:),&
& lon_bx(ib,:,:), lat_bx(ib,:,:), trip_bx(ib,:,:), &
& hierarchy_bx(ib,:,:), orog_bx(ib,:,:), floodp_bx(ib,:,:), fac_bx(ib,:,:), &
& topoind_bx(ib,:,:), min_topoind, nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), &
& basin_sz(ib,:), basin_pts(ib,:,:,:), basin_bxout(ib,:), basin_bbout(ib,:), lshead(ib,:), coast_pts(ib,:), nwbas, &
& basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy,&
& topoind_bx(ib,:,:), rlen_bx(ib,:,:), rdz_bx(ib,:,:), rweight_bx(ib,:,:), min_topoind, &
& nb_basin(ib), basin_inbxid(ib,:), basin_outlet(ib,:,:), basin_outtp(ib,:), basin_sz(ib,:), &
& basin_pts(ib,:,:,:), basin_bxout(ib,:), basin_bbout(ib,:), lshead(ib,:), &
& coast_pts(ib,:), nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy,&
& basin_orog_mean, basin_orog_min, basin_orog_max, basin_floodp, basin_fac,&
& basin_topoind, basin_id, basin_coor, &
& basin_type, basin_flowdir, basin_lshead, basin_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
& 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)
ENDDO
!
END SUBROUTINE globalize
......
F90subroutines/routing_reg.f90
View file @ c8b9a77d
......@@ -13,6 +13,8 @@ MODULE routing_reg
INTEGER(i_std), SAVE :: nbpt_save
INTEGER(i_std), SAVE :: nbasmax_save
REAL(r_std), SAVE :: topoindexmin = 0.0
REAL(r_std), SAVE :: hydro_meanlen = 0.0
REAL(r_std), SAVE :: hydro_meandz = 0.0
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)
......@@ -88,9 +90,9 @@ CONTAINS
!_ ================================================================================================================================
SUBROUTINE routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, fac, hierarchy, &
& orog, floodp,&
& nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, fac_bx, hierarchy_bx,orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
& meanrlen, meanrdz, 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, &
& hierarchy_bx, orog_bx, floodp_bx, lon_bx, lat_bx, lshead_bx)
!
IMPLICIT NONE
!
......@@ -102,6 +104,7 @@ CONTAINS
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(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) :: sub_area(nbpt,nbvmax) !! Area on the fine grid (m^2)
!
REAL(r_std), INTENT(in) :: lon_rel(nbpt,nbvmax) !!
......@@ -114,6 +117,8 @@ 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) :: 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
......@@ -128,6 +133,8 @@ CONTAINS
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) !!
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) !!
......@@ -164,12 +171,16 @@ CONTAINS
inc(8,2) = -1
!
topoindexmin = min_topoind
hydro_meanlen = meanrlen
hydro_meandz = meanrdz
!
! Set everything to undef to locate easily empty points
!
trip_bx(:,:) = undef_int
basin_bx(:,:) = undef_int
topoind_bx(:,:) = undef_sechiba
rlen_bx(:,:) = undef_sechiba
rdz_bx(:,:) = undef_sechiba
area_bx(:,:) = undef_sechiba
hierarchy_bx(:,:) = undef_sechiba
fac_bx(:,:) = undef_sechiba
......@@ -221,6 +232,8 @@ CONTAINS
basin_bx(iloc, jloc) = NINT(basins(ib, ip))
area_bx(iloc, jloc) = sub_area(ib, ip)
topoind_bx(iloc, jloc) = topoindex(ib, ip)
rlen_bx(iloc, jloc) = rlen(ib, ip)
rdz_bx(iloc, jloc) = rdz(ib, ip)
hierarchy_bx(iloc, jloc) = hierarchy(ib, ip)
fac_bx(iloc, jloc) = fac(ib, ip)
orog_bx(iloc, jloc) = orog(ib, ip)
......@@ -244,6 +257,8 @@ CONTAINS
max_basins = max_basins + 1
area_bx(iloc, jloc) = resolution(ib,1)*resolution(ib,2)*contfrac(ib)
topoind_bx(iloc, jloc) = min_topoind
rlen_bx(iloc, jloc) = hydro_meanlen
rdz_bx(iloc, jloc) = hydro_meandz
! Be careful here: We deal with basin_hierarchy(ib,ij) = min_topoind later!
hierarchy_bx(iloc, jloc) = min_topoind
fac_bx(iloc, jloc) = zero
......@@ -384,9 +399,9 @@ CONTAINS
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_findbasins(nb_htu, nbv, ib, ijdimmax, nbi, nbj, trip, basin, fac, hierarchy, topoind, lshead, diaglalo, &
& nb_basin, basin_inbxid, basin_outlet, basin_outtp, basin_sz, basin_bxout, basin_bbout, basin_pts, &
& basin_lshead, coast_pts, lontmp, lattmp)
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)
!
IMPLICIT NONE
!
......@@ -405,7 +420,10 @@ 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 (m)
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(out) :: rweight(ijdimmax,ijdimmax) !!
!
! For debug and get coordinate of river outlet
!
......@@ -454,7 +472,6 @@ CONTAINS
INTEGER(i_std) :: sortind(nbvmax) !!
!
INTEGER(i_std) :: rivpas(ijdimmax,ijdimmax) !! Number of passage through grid when following rivers
INTEGER(i_std) :: rivlen(ijdimmax,ijdimmax) !! number of grids along the river
INTEGER(i_std) :: color(ijdimmax,ijdimmax)
!
INTEGER(i_std) :: itrans !!
......@@ -988,10 +1005,12 @@ CONTAINS
ENDDO
ENDDO
!
! Add current point to the weight calculation.
rweight(:,:) = 1.0
!
! Work only for topo_option number 3 or 4
!
IF ( topo_option .EQ. 3 ) THEN
rivlen(:,:) = 1
rivpas(:,:) = 0
! Loop for all sub-basins and each point belongs to that sub-basin
DO ij=1,nb_basin
......@@ -1013,9 +1032,11 @@ CONTAINS
! this point still belongs to current sub-basin. Here I
! assume that the sub-routine routing_reg_divbas worked well.
! So we don't need to check.
topoind(ip,jp)=topoind(ip,jp)+topoind(ill,jll)
topoind(ip,jp) = topoind(ip,jp)+topoind(ill,jll)
rlen(ip,jp) = rlen(ip,jp)+rlen(ill,jll)
rdz(ip,jp) = rdz(ip,jp)+rdz(ill,jll)
rivpas(ill,jll) = rivpas(ill,jll)+1
rivlen(ip,jp) = rivlen(ip,jp)+1
rweight(ip,jp) = rweight(ip,jp)+1.0
!
il = ill ; jl = jll
p = trip(il,jl)
......@@ -1030,11 +1051,12 @@ CONTAINS
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
IF ( rivpas(ip,jp) == 0 ) THEN
! Found a river head
topoind(ip,jp) = topoind(ip,jp) * rivlen(ip,jp)
cnt = cnt + rivlen(ip,jp)
cnt = cnt + rweight(ip,jp)
ELSE
topoind(ip,jp) = 0.0
topoind(ip,jp) = zero
rlen(ip,jp) = zero
rdz(ip,jp) = zero
rweight(ip,jp) = zero
ENDIF
ENDDO
! Finish the weight computation of the rivers followed
......@@ -1042,7 +1064,7 @@ CONTAINS
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
topoind(ip,jp) = topoind(ip,jp)/cnt
rweight(ip,jp) = rweight(ip,jp)/cnt
ENDDO
ENDIF
ENDDO
......@@ -1557,11 +1579,11 @@ CONTAINS
!_ ================================================================================================================================
SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, area_bx, lon_bx, lat_bx, trip_bx, &
& hierarchy_bx, orog_bx, floodp_bx, fac_bx, topoind_bx, &
& hierarchy_bx, orog_bx, floodp_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, &
& basin_bbout, lshead, coast_pts, nwbas, 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_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
& 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_beta_fp, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
!
IMPLICIT NONE
......@@ -1582,16 +1604,19 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
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) :: 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
REAL(r_std), DIMENSION(nb_htu,2) :: basin_outlet !! Outlet coordinate of each subgrid basin (lat,lon)
REAL(r_std), DIMENSION(nbvmax) :: basin_outtp !!
INTEGER(i_std), DIMENSION(nb_htu,nbv,2) :: basin_pts !! Points in each basin
INTEGER(i_std), DIMENSION(nb_htu) :: basin_bxout !! outflow direction
INTEGER(i_std), DIMENSION(nb_htu) :: basin_bbout !! outflow sub-basin
REAL(r_std), DIMENSION(nb_htu) :: lshead !! Large scale heading of outflow.
INTEGER(i_std) :: coast_pts(nb_htu) !! The coastal flow points (unitless)
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) :: 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
REAL(r_std), DIMENSION(nb_htu,2) :: basin_outlet !! Outlet coordinate of each subgrid basin (lat,lon)
REAL(r_std), DIMENSION(nbvmax) :: basin_outtp !!
INTEGER(i_std), DIMENSION(nb_htu,nbv,2) :: basin_pts !! Points in each basin
INTEGER(i_std), DIMENSION(nb_htu) :: basin_bxout !! outflow direction
INTEGER(i_std), DIMENSION(nb_htu) :: basin_bbout !! outflow sub-basin
REAL(r_std), DIMENSION(nb_htu) :: lshead !! Large scale heading of outflow.
INTEGER(i_std) :: coast_pts(nb_htu) !! The coastal flow points (unitless)
! global maps
INTEGER(i_std) :: nwbas !!
INTEGER(i_std), DIMENSION(nbpt) :: basin_count !!
......@@ -1608,7 +1633,9 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
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 (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_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 !!
......@@ -1618,6 +1645,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: coastal_basin !!
!
INTEGER(i_std) :: ij, iz !! Indices (unitless)
INTEGER(i_std) :: ip, jp, cnt
CHARACTER(LEN=4) :: hierar_method = 'OUTP' !!
!
! To calculate topoindex in new method
......@@ -1645,6 +1673,8 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
ENDIF
!
basin_topoind(ib,:) = undef_sechiba
basin_rlen(ib,:) = undef_sechiba
basin_rdz(ib,:) = undef_sechiba
!
DO ij=1, nb_basin
!
......@@ -1696,6 +1726,8 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!
END SELECT
basin_topoind(ib,ij) = zero
basin_rlen(ib,ij) = zero
basin_rdz(ib,ij) = zero
!
DO iz=1,basin_sz(ij)
!
......@@ -1786,7 +1818,6 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
!
! Calculating topoindex depends on options
!
basin_topoind(ib,ij) = zero
! Constant value (1000.)
IF ( topo_option .EQ. 1 ) THEN
basin_topoind(ib,ij) = 1000.
......@@ -1798,13 +1829,37 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
ENDDO
basin_topoind(ib,ij) = basin_topoind(ib,ij)/REAL(basin_sz(ij),r_std)
ENDIF
! All topoindex with HTU have already been weighted and those not at the root of the river set to zero.
! Weight averaged topoindex of rivers within HTU. Those not at the root of the river set to zero.
IF ( topo_option .EQ. 3 ) THEN
DO iz=1,basin_sz(ij)
basin_topoind(ib,ij) = basin_topoind(ib,ij) + topoind_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
basin_topoind(ib,ij) = basin_topoind(ib,ij) + rweight_bx(ip,jp)*topoind_bx(ip,jp)
ENDDO
ENDIF
!
! Compute mean length of rivers in HTU
!
cnt = 0
DO iz=1,basin_sz(ij)
ip = basin_pts(ij,iz,1)
jp = basin_pts(ij,iz,2)
! Both rlen and rdz are zero outside of the first point in a river of the HTU
IF ( rlen_bx(ip,jp) .GT. 0 ) THEN
cnt = cnt + 1
basin_rlen(ib,ij) = basin_rlen(ib,ij)+rlen_bx(ip,jp)
basin_rdz(ib,ij) = basin_rdz(ib,ij)+rdz_bx(ip,jp)
ENDIF
ENDDO
IF ( cnt > 0 ) THEN
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.
basin_rlen(ib,ij) = hydro_meanlen
basin_rdz(ib,ij) = hydro_meandz
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
......
HydroGrid.py
View file @ c8b9a77d
......@@ -109,6 +109,10 @@ def calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) :
land = np.where(nf.variables["orog"][jstr:jend,istr:iend] < missing)
topoind=np.zeros((jend-jstr,iend-istr))
topoind[:,:] = np.nan
rlen=np.zeros((jend-jstr,iend-istr))
rlen[:,:] = np.nan
rdz=np.zeros((jend-jstr,iend-istr))
rdz[:,:] = np.nan
#
for j,i in zip(land[:][0],land[:][1]) :
trip = int(nf.variables["trip"][jstr+j,istr+i]-1)
......@@ -132,6 +136,8 @@ def calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) :
rx,ry = hydrogrid.resolution(i,j)
dl = (rx+ry)*0.5*0.5
topoind[j,i]=np.sqrt(dl**3./(dz*1.0e6))
rlen[j,i]=dl
rdz[j,i]=dz
#
mintopo = np.nanmin(topoind)
maxtopo = np.nanmax(topoind)
......@@ -139,7 +145,11 @@ def calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid) :
ERROR("Topoind close to zero encoutered.")
sys.exit()
#
return gather(topoind, index, default = 10), "Topographic index which serves for the residence time", mintopo, maxtopo
meanrlen = np.nanmean(rlen)
meanrdz = np.nanmean(rdz)
#
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
#
class HydroGrid :
def __init__(self, lolacorners, wfile) :
......@@ -254,9 +264,12 @@ class HydroData :
self.topoinddesc = "desc"
mintopo = 1
else:
self.topoind, self.topoinddesc, mintopo, maxtopo = calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid)
self.topoind, self.rlen, self.rdz, self.topoinddesc, mintopo, maxtopo, mrlen, mrdz = \
calctopoindex(nf, istr, iend, jstr, jend, index, hydrogrid)
self.topoindmin = part.domainmin(mintopo)
self.topoindmax = part.domainmax(maxtopo)
self.meanlength = part.domainmean(mrlen)
self.meandz = part.domainmean(mrdz)
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.topoind=gather(nf.variables["topoind"][jstr:jend,istr:iend].astype(np.float32), index, default = 10)
......
Interface.py
View file @ c8b9a77d
......@@ -91,7 +91,8 @@ class HydroOverlap :
#
ijdim=[]
for ib in range(nbpt) :
ijdim.append(max(np.max(sub_index[ib,:sub_pts[ib],0])-np.min(sub_index[ib,:sub_pts[ib],0])+