Skip to content

GitLab

Commit 8b2082f5 authored by Anthony's avatar Anthony
Browse files

Add the orography and the floodplains

parent 5b2959e8
Hide whitespace changes
Inline Side-by-side
Showing with 180 additions and 85 deletions
F90subroutines/routing_interface.f90
View file @ 8b2082f5
......@@ -44,8 +44,8 @@ SUBROUTINE closeinterface()
END SUBROUTINE closeinterface
SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& lon_rel, lat_rel, trip, basins, topoindex, fac, hierarchy, nbi, nbj, area_bx, trip_bx, basin_bx, &
& topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx)
& 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)
!
USE ioipsl
USE grid
......@@ -71,6 +71,8 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
REAL, INTENT(inout) :: topoindex(nbpt,nbvmax_in) !! Topographic index of the residence time (m)
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
REAL, INTENT(inout) :: floodp(nbpt,nbvmax_in) !! data on the fine grid
!
!! OUTPUT VARIABLES
INTEGER, INTENT(out) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless)
......@@ -83,6 +85,8 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
REAL, INTENT(out) :: topoind_bx(nbpt,nbxmax_in,nbxmax_in) !! Topographic index of the residence time for each of the smaller boxes (m)
INTEGER, INTENT(out) :: trip_bx(nbpt,nbxmax_in,nbxmax_in) !! The trip field for each of the smaller boxes (unitless)
INTEGER, INTENT(out) :: basin_bx(nbpt,nbxmax_in,nbxmax_in) !!
REAL, INTENT(out) :: orog_bx(nbpt,nbxmax_in,nbxmax_in) !! Orography (m)
REAL, INTENT(out) :: floodp_bx(nbpt,nbxmax_in,nbxmax_in) !! floodplains (m^2)
!
INTEGER :: ii, ib
REAL :: resolution(nbpt,2)
......@@ -101,8 +105,9 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
DO ib=1,nbpt
CALL routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, fac, hierarchy, &
& nbi(ib), nbj(ib), area_bx(ib,:,:), trip_bx(ib,:,:), basin_bx(ib,:,:), topoind_bx(ib,:,:), &
& fac_bx(ib,:,:), hierarchy_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), lshead_bx(ib,:,:))
& 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,:,:))
ENDDO
!
END SUBROUTINE gethydrogrid
......@@ -167,9 +172,11 @@ SUBROUTINE findbasins(nbpt, nbvmax_in, nbxmax_in, nbi, nbj, trip_bx, basin_bx, f
END SUBROUTINE findbasins
SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, fac_bx, topoind_bx, &
SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, orog_bx, floodp_bx, &
& fac_bx, topoind_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, basin_floodp, &
& basin_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
& outflow_grid, outflow_basin, nbcoastal, coastal_basin)
!
......@@ -188,6 +195,8 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: lat_bx !! Latitude of each small box in the grid box
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: trip_bx !! The trip field for each of the smaller boxes (unitless)
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: hierarchy_bx !! Level in the basin of the point
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: orog_bx !! Orography (m)
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: floodp_bx !! Floodplains area (m^2)
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: fac_bx !! Level in the basin of the point
REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: 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)
......@@ -211,6 +220,8 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_area !!
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas,2) :: basin_cg !! Centre of gravity of the HTU
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_hierarchy !!
REAL(r_std), INTENT (out), DIMENSION(nbpt,nwbas) :: basin_orog !!
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_lshead !! Large scale heading out of the grid box.
......@@ -228,10 +239,11 @@ SUBROUTINE globalize(nbpt, nbvmax_in, nbxmax_in, area_bx, lon_bx, lat_bx, trip_b
!!
DO ib=1,nbpt
CALL routing_reg_globalize(nbpt, ib, neighbours, area_bx(ib,:,:), lon_bx(ib,:,:), lat_bx(ib,:,:), trip_bx(ib,:,:), &
& hierarchy_bx(ib,:,:), fac_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, basin_fac, basin_topoind, basin_id, basin_coor, &
& basin_count, basin_notrun, basin_area, basin_cg, basin_hierarchy,&
& basin_orog, basin_floodp, basin_fac, basin_topoind, basin_id, basin_coor, &
& basin_type, basin_flowdir, basin_lshead, outflow_grid, outflow_basin, nbcoastal, coastal_basin)
ENDDO
!
......@@ -417,10 +429,12 @@ SUBROUTINE rivclassification(nbpt, nwbas, nbcore, corepts, basin_count, outflow_
END SUBROUTINE rivclassification
SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, basin_notrun, basin_area, &
& basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin, routing_area, routing_cg, topo_resid, route_nbbasin, route_togrid, &
& route_tobasin, route_nbintobas, global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridarea, cfrac, basin_count, &
& basin_notrun, basin_area, basin_orog, basin_floodp, basin_cg, basin_topoind, fetch_basin, &
& basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin, routing_area, routing_orog, routing_floodp, &
& routing_cg, topo_resid, route_nbbasin, route_togrid, route_tobasin, route_nbintobas, &
& global_basinid, route_outlet, route_type, origin_nbintobas, routing_fetch)
!
USE ioipsl
USE grid
......@@ -444,6 +458,8 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_type !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_flowdir !! Water flow directions in the basin (unitless)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_orog !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_cg
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
......@@ -459,6 +475,8 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
! Output variables
!
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_area !! Surface of basin (m^2)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_orog !! Mean Orography (m)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_floodp !! Surface of Floodplains (m^2)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: routing_fetch !! Upstream are flowing into HTU (m^2)
REAL(r_std), DIMENSION(nbpt,nbasmax,2), INTENT(out) :: routing_cg !! Centre of gravity of HTU (Lat, Lon)
REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: topo_resid !! Topographic index of the retention time (m)
......@@ -474,11 +492,14 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
!
!
CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, basin_count, &
& basin_notrun, basin_area, basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, &
& basin_count, basin_notrun, basin_area, basin_orog, basin_floodp, basin_cg, &
& basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
routing_area(:,:) = routing_area_glo(:,:)
routing_orog(:,:) = routing_orog_glo(:,:)
routing_floodp(:,:) = routing_floodp_glo(:,:)
routing_cg(:,:,:) = routing_cg_glo(:,:,:)
topo_resid(:,:) = topo_resid_glo(:,:)
route_nbbasin(:) = route_count_glo(:)
......@@ -545,6 +566,7 @@ END SUBROUTINE finish_inflows
SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numops, basin_count, basin_area, &
& basin_orog, basin_floodp, &
& basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
!
......@@ -568,6 +590,8 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_type !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_flowdir !! Water flow directions in the basin (unitless)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_orog !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_cg
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
......@@ -589,7 +613,7 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
IF (basin_count(ib) > nbasmax) THEN
tok = tokill(ib,op)
totak = totakeover(ib,op)
IF (tok .GT. 0) THEN
IF (tok .GT. 0) THEN
WRITE(numout,*) "nbpt", ib, "tokill", tok, "totakover", totak
! Test if tokill is downstream of totakeover (avoid loop)
igrif = outflow_grid(ib,totak)
......@@ -603,11 +627,12 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
tok = it
ELSE
it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END IF
END DO
CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, basin_area, basin_cg, basin_topoind,&
CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, basin_area, &
& basin_orog, basin_floodp, basin_cg, basin_topoind,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
& inflow_grid, inflow_basin)
END IF
......@@ -615,8 +640,6 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
END DO
END DO
! REPRENDRE LA FIN DE TRUNCATE
END SUBROUTINE killbas
......@@ -633,7 +656,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
! Local
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: flag !!
......@@ -641,13 +664,13 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
flag(:,:) = 0
WRITE(numout,*) "Checking routing"
test = 0
DO ig=1,nbpt
basnum = basin_count(ig)
DO ib=1,basnum
DO ib=1,basnum
IF (flag(ig,ib) .EQ. 0) THEN
flag(ig,ib) = -1
igrif = outflow_grid(ig,ib)
......@@ -655,7 +678,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
!
DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -1) .AND. (flag(igrif,ibasf) .NE. -99))
flag(igrif, ibasf) = -1
it = outflow_grid(igrif,ibasf)
it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
......@@ -668,7 +691,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
DO WHILE ((igrif .GT. 0) .AND. (flag(igrif,ibasf) .NE. -99))
flag(igrif, ibasf) = -99
it = outflow_grid(igrif,ibasf)
it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
......@@ -681,7 +704,7 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
ibasf = outflow_basin(ig,ib)
DO WHILE (flag(igrif,ibasf) .EQ. -1)
flag(igrif, ibasf) = -99
it = outflow_grid(igrif,ibasf)
it = outflow_grid(igrif,ibasf)
ibasf = outflow_basin(igrif,ibasf)
igrif = it
END DO
......@@ -720,11 +743,11 @@ SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, outflow_grid, outflow_basin, bas
DO ig=1,nbpt
basnum = basin_count(ig)
DO ib=1,basnum
DO ib=1,basnum
igrif = outflow_grid(ig,ib)
ibasf = outflow_basin(ig,ib)
IF (igrif .GT. 0) THEN
IF (igrif .GT. 0) THEN
IF (fetch_basin(ig,ib) .GT. fetch_basin(igrif,ibasf)) THEN
WRITE(numout,*) "****ERROR Fetch, nbpt:",ig
WRITE(numout,*) igrif, ibasf
......
F90subroutines/routing_reg.f90
View file @ 8b2082f5
......@@ -15,6 +15,8 @@ MODULE routing_reg
INTEGER(i_std), SAVE :: nbasmax_save
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_area_glo !! Surface of basin (m^2)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_orog_glo !! Mean topography (m)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: routing_floodp_glo !! Surface of floodplains (m^2)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: routing_cg_glo !! Centre of gravity of HTU (Lat, Lon)
REAL(r_std), SAVE, ALLOCATABLE, DIMENSION(:,:) :: topo_resid_glo !! Topographic index of the retention time (m)
INTEGER(i_std), SAVE, ALLOCATABLE, DIMENSION(:) :: route_count_glo !! Number of basins per grid
......@@ -73,7 +75,8 @@ CONTAINS
SUBROUTINE routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& lon_rel, lat_rel, lalo, resolution, contfrac, trip, basins, topoindex, fac, hierarchy, &
& nbi, nbj, area_bx, trip_bx, basin_bx, topoind_bx, fac_bx, hierarchy_bx, lon_bx, lat_bx, lshead_bx)
& 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)
!
IMPLICIT NONE
!
......@@ -98,6 +101,8 @@ CONTAINS
REAL(r_std), INTENT(inout) :: topoindex(nbpt,nbvmax) !! Topographic index of the residence time (m)
REAL(r_std), INTENT(inout) :: hierarchy(nbpt,nbvmax) !! data on the fine grid
REAL(r_std), INTENT(inout) :: fac(nbpt,nbvmax) !! data on the fine grid
REAL(r_std), INTENT(inout) :: orog(nbpt,nbvmax) !! data on the fine grid
REAL(r_std), INTENT(inout) :: floodp(nbpt,nbvmax) !! data on the fine grid
!
!! OUTPUT VARIABLES
INTEGER(i_std), INTENT(out) :: nbi, nbj !! Number of point in x and y within the grid (unitless)
......@@ -110,6 +115,8 @@ CONTAINS
REAL(r_std), INTENT(out) :: topoind_bx(nbxmax,nbxmax) !! Topographic index of the residence time for each of the smaller boxes (m)
INTEGER(i_std), INTENT(out) :: trip_bx(nbxmax,nbxmax) !! The trip field for each of the smaller boxes (unitless)
INTEGER(i_std), INTENT(out) :: basin_bx(nbxmax,nbxmax) !!
REAL(i_std), INTENT(out) :: orog_bx(nbxmax,nbxmax) !!
REAL(i_std), INTENT(out) :: floodp_bx(nbxmax,nbxmax) !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc !! Indices (unitless)
......@@ -152,6 +159,8 @@ CONTAINS
lon_bx(:,:) = undef_sechiba
lat_bx(:,:) = undef_sechiba
lshead_bx(:,:) = undef_sechiba
orog_bx(:,:) = undef_sechiba
floodp_bx(:,:) = undef_sechiba
cenlon = zero
cenlat = zero
!
......@@ -193,6 +202,8 @@ CONTAINS
topoind_bx(iloc, jloc) = topoindex(ib, ip)
hierarchy_bx(iloc, jloc) = hierarchy(ib, ip)
fac_bx(iloc, jloc) = fac(ib, ip)
orog_bx(iloc, jloc) = orog(ib, ip)
floodp_bx(iloc, jloc) = floodp(ib, ip) * sub_area(ib, ip)
lon_bx(iloc, jloc) = lon_rel(ib, ip)
lat_bx(iloc, jloc) = lat_rel(ib, ip)
cenlon = cenlon + lon_rel(ib, ip)/sub_pts(ib)
......@@ -1635,9 +1646,11 @@ CONTAINS
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, fac_bx, topoind_bx, &
SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, &
& orog_bx, floodp_bx, fac_bx, topoind_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, basin_floodp, &
& basin_fac, basin_topoind, basin_id, basin_coor, basin_type, basin_flowdir, basin_lshead, &
& outflow_grid, outflow_basin, nbcoastal, coastal_basin)
!
......@@ -1654,6 +1667,8 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: lat_bx !! Latitude of each small box in the grid box
INTEGER(i_std), DIMENSION(nbxmax,nbxmax) :: trip_bx !! The trip field for each of the smaller boxes (unitless)
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: hierarchy_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: orog_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: floodp_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: fac_bx !! Level in the basin of the point
REAL(r_std), DIMENSION(nbxmax,nbxmax) :: 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)
......@@ -1677,6 +1692,8 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area !!
REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_cg !! Centre of gravity of the HTU in latitude, longitude
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_hierarchy !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog !!
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_lshead !! Large scale heading out of the grid box.
......@@ -1755,6 +1772,8 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
basin_cg(ib,ij,:) = zero
basin_hierarchy(ib,ij) = zero
basin_fac(ib,ij) = zero
basin_orog(ib,ij) = zero
basin_floodp(ib,ij) = zero
!
SELECT CASE (hierar_method)
!
......@@ -1767,6 +1786,10 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
DO iz=1,basin_sz(ij)
!
basin_area(ib,ij) = basin_area(ib,ij) + area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
! Numerator for mean orography, then need to be divided by the HTUs area
basin_orog(ib,ij) = basin_orog(ib,ij) + orog_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) &
& * area_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
basin_floodp(ib,ij) = basin_floodp(ib,ij) + floodp_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2))
!
! Compute centre of gravity
!
......@@ -1813,6 +1836,9 @@ SUBROUTINE routing_reg_globalize(nbpt, ib, neighbours, area_bx, lon_bx, lat_bx,
!
ENDDO
!
! Finish the calculation of mean orography
basin_orog(ib,ij) = basin_orog(ib,ij)/basin_area(ib,ij)
!
SELECT CASE (hierar_method)
!
CASE("MEAN")
......@@ -2054,7 +2080,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,1) = solved(sp,1) + 1
ENDIF
WRITE(numout,*) sp,sb,"flow in the same grid, basin:",bop
WRITE(numout,*) sp,sb,"flow in the same grid, basin:",bop
!
ELSE IF ( outflow_grid(sp,sb) .EQ. -3 ) THEN
! Return flow
......@@ -2296,20 +2322,20 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
& inflow_number, inflow_grid, inflow_basin)
! It is possible that the lowest hierarchy is in two grid cells
! In this case, if we have an error for this routing_updateflow
! We go to next step and make it a coastal flow
! We go to next step and make it a coastal flow
IF (outflow_basin(sp, sb) < undef_int) THEN
found = 1
WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points"
ELSE
WRITE(numout,*) "Lowest hierarchy may be in two different grid points"
END IF
WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb)
WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb)
WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop)
END IF
END IF
END DO
END IF
ENDIF
ENDIF
!
! Last option : coastal outflow
IF (found .EQ. 0) THEN
......@@ -2554,8 +2580,9 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, basin_count, &
& basin_notrun, basin_area, basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, &
& basin_count, basin_notrun, basin_area, basin_orog, basin_floodp, &
& basin_cg, basin_topoind, fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, &
& outflow_grid, outflow_basin, inflow_number, inflow_grid, inflow_basin)
!
IMPLICIT NONE
......@@ -2581,6 +2608,8 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_type !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_flowdir !! Water flow directions in the basin (unitless)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_orog !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2), INTENT(inout) :: basin_cg !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
......@@ -2647,6 +2676,8 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
DO ij=1,basin_count(ib)
!
routing_area_glo(ib,ij) = basin_area(ib,ij)
routing_orog_glo(ib,ij) = basin_orog(ib,ij)
routing_floodp_glo(ib,ij) = basin_floodp(ib,ij)
routing_cg_glo(ib,ij,:) = basin_cg(ib,ij,:)
topo_resid_glo(ib,ij) = basin_topoind(ib,ij)
global_basinid_glo(ib,ij) = basin_id(ib,ij)
......@@ -2751,7 +2782,8 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
!! \n
!_ ================================================================================================================================
SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, basin_cg, basin_topoind,&
SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count, basin_area, &
& basin_orog, basin_floodp, basin_cg, basin_topoind,&
& fetch_basin, basin_id, basin_coor, basin_type, basin_flowdir, outflow_grid, outflow_basin, inflow_number,&
& inflow_grid, inflow_basin)
!
......@@ -2770,6 +2802,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_type !!
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_flowdir !! Water flow directions in the basin (unitless)
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_area !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_orog !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_floodp !!
REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_cg !!
REAL(r_std), DIMENSION(nbpt,nwbas) :: basin_topoind !! Topographic index of the residence time for a basin (m)
REAL(r_std), DIMENSION(nbpt,nwbas) :: fetch_basin !!
......@@ -2802,6 +2836,11 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
!
basin_area(ib, totakeover) = basin_area(ib, totakeover) + basin_area(ib,tokill)
!
basin_floodp(ib, totakeover) = basin_floodp(ib, totakeover) + basin_floodp(ib,tokill)
!
basin_orog(ib, totakeover) = (basin_orog(ib, totakeover)*basin_area(ib, totakeover) &
& + basin_orog(ib, tokill)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
! Add the fetch of the basin will kill to the one which gets the water
!
fetch_basin(ib, totakeover) = fetch_basin(ib, totakeover) + fetch_basin(ib, tokill)
......@@ -2937,6 +2976,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
nbasmax_save = nbasmax
ALLOCATE (routing_area_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (routing_orog_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (routing_floodp_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (route_fetch_glo(nbpt,nbasmax), stat=ier)
ALLOCATE (routing_cg_glo(nbpt,nbasmax,2), stat=ier)
ALLOCATE (topo_resid_glo(nbpt,nbasmax), stat=ier)
......
HydroGrid.py
View file @ 8b2082f5
......@@ -49,10 +49,10 @@ def corners(lon, lat) :
#
centersll.append([lon[j,i], lat[j,i]])
radiusll.append(RPP.maxradius([lon[j,i], lat[j,i]], lons, lats))
#
#
index.append([j,i])
cornersll.append(boxll)
return cornersll, cornerspoly, centersll, radiusll, index
#
def gather(x, index) :
......@@ -130,3 +130,15 @@ class HydroData :
#
self.fac=gather(nf.variables["fac"][jstr:jend,istr:iend].astype(np.float32), index)
self.facdesc=nf.variables["fac"].long_name
#
if "orog" in nf.variables.keys():
self.orog = gather(nf.variables["orog"][jstr:jend,istr:iend].astype(np.float32), index)
self.orogdesc=nf.variables["orog"].long_name
else:
self.orog = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
#
if "floodplains" in nf.variables.keys():
self.floodplains = gather(nf.variables["floodplains"][jstr:jend,istr:iend].astype(np.float32), index)
self.floodplainsdesc=nf.variables["floodplains"].long_name
else:
self.floodplains = gather(np.zeros((jend-jstr,iend-istr)).astype(np.float32), index)
Interface.py
View file @ 8b2082f5
......@@ -38,7 +38,7 @@ prec = 100.0
#
# Print the documentation for the FORTRAN interface
#
if gendoc.lower() == "true" :
if gendoc.lower() == "true" :
docwrapper = open('DocumentationInterface', 'w')
docwrapper.write(routing_interface.initatmgrid.__doc__)
docwrapper.write("====================================================================\n")
......@@ -212,9 +212,9 @@ def finalfetch(part, routing_area, basin_count, route_togrid, route_tobasin, fet
#
fetch_error = np.sum(np.abs(fetch_out[part.landcorelist,:]-fetch_in[part.landcorelist,:]), axis=1)\
/np.sum(routing_area[part.landcorelist,:], axis=1)
if np.max(fetch_error) > prec :
if np.max(fetch_error) > prec :
print("Rank :"+str(part.rank)+" Too large fetch error (fraction of greid area) : ", fetch_error)
print("Total fetch error in fraction of grid box : ", np.sum(fetch_error))
#
return fetch_out
......@@ -245,6 +245,8 @@ class HydroOverlap :
topoind_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
fac_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
hierarchy_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
orog_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
floodp_tmp = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
#
trip_tmp[:,:] = np.nan
basins_tmp[:,:] = np.nan
......@@ -254,6 +256,8 @@ class HydroOverlap :
topoind_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.topoind[ib][:])
fac_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.fac[ib][:])
hierarchy_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.disto[ib][:])
orog_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.orog[ib][:])
floodp_tmp[ib,0:sub_pts[ib]] = np.asarray(hydrodata.floodplains[ib][:])
#
trip_tmp[np.isnan(trip_tmp)] = undef_int
basins_tmp[np.isnan(trip_tmp)] = undef_int
......@@ -264,9 +268,11 @@ class HydroOverlap :
print("GETHYDROGRID : nbvmax = ", nbvmax)
print("GETHYDROGRID : nbxmax = ", nbxmax)
self.nbi, self.nbj, self.area_bx, self.trip_bx, self.basin_bx, self.topoind_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(nbxmax, sub_pts, sub_index, sub_area, \
hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp, hierarchy_tmp)
hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp,\
hierarchy_tmp, orog_tmp, floodp_tmp)
#
# Plot some diagnostics for the hydrology grid within the atmospheric meshes.
#
......@@ -302,16 +308,18 @@ class HydroSuper :
lon_bx_tmp[np.isnan(lon_bx_tmp)] = undef_int
lat_bx_tmp = hydrooverlap.lat_bx
lat_bx_tmp[np.isnan(lat_bx_tmp)] = undef_int
self.basin_count, self.basin_notrun, self.basin_area, self.basin_cg, self.basin_hierarchy, self.basin_fac, self.basin_topoind, \
self.basin_count, self.basin_notrun, self.basin_area, self.basin_cg, self.basin_hierarchy, \
self.basin_orog, self.basin_floodp, self.basin_fac, self.basin_topoind, \
self.basin_id, self.basin_outcoor, self.basin_type, self.basin_flowdir, \
self.basin_lshead, self.outflow_grid, self.outflow_basin, self.nbcoastal, self.coastal_basin = \
routing_interface.globalize(hydrooverlap.area_bx, lon_bx_tmp, lat_bx_tmp, hydrooverlap.trip_bx, \
hydrooverlap.hierarchy_bx, hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \
hydrooverlap.hierarchy_bx, hydrooverlap.orog_bx, hydrooverlap.floodp_bx,\
hydrooverlap.fac_bx, hydrooverlap.topoind_bx, hydrodata.topoindmin, \