diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index bba1e001e0af1eceea9b9a20be7e9f184070a643..4d4bd97699558d4cb97299517893ca91d80d8995 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -43,7 +43,7 @@ SUBROUTINE closeinterface()
!
END SUBROUTINE closeinterface
-SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
+SUBROUTINE gethydrogrid(nbpt, nbvmax_in, ijdimmax, 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)
!
@@ -54,7 +54,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
!
IMPLICIT NONE
!
- INTEGER, INTENT(in) :: nbvmax_in, nbxmax_in
+ INTEGER, INTENT(in) :: nbvmax_in, ijdimmax
!
INTEGER, INTENT(in) :: nbpt !! Domain size (unitless)
INTEGER, INTENT(in) :: sub_pts(nbpt) !! Number of high resolution points on this grid (unitless)
@@ -76,26 +76,28 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
!
!! OUTPUT VARIABLES
INTEGER, INTENT(out) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless)
- REAL, INTENT(out) :: area_bx(nbpt,nbxmax_in,nbxmax_in) !! Area of each small box in the grid box (m^2)
- REAL, INTENT(out) :: hierarchy_bx(nbpt,nbxmax_in,nbxmax_in) !! Level in the basin of the point
- REAL, INTENT(out) :: fac_bx(nbpt,nbxmax_in,nbxmax_in) !! Flow accumulation
- REAL, INTENT(out) :: lon_bx(nbpt,nbxmax_in,nbxmax_in) !!
- REAL, INTENT(out) :: lat_bx(nbpt,nbxmax_in,nbxmax_in) !!
- REAL, INTENT(out) :: lshead_bx(nbpt,nbxmax_in,nbxmax_in) !! Large scale heading for outflow points.
- 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)
+ REAL, INTENT(out) :: area_bx(nbpt,ijdimmax,ijdimmax) !! Area of each small box in the grid box (m^2)
+ REAL, INTENT(out) :: hierarchy_bx(nbpt,ijdimmax,ijdimmax) !! Level in the basin of the point
+ REAL, INTENT(out) :: fac_bx(nbpt,ijdimmax,ijdimmax) !! Flow accumulation
+ REAL, INTENT(out) :: lon_bx(nbpt,ijdimmax,ijdimmax) !!
+ 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)
+ 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)
+ REAL, INTENT(out) :: floodp_bx(nbpt,ijdimmax,ijdimmax) !! floodplains (m^2)
!
INTEGER :: ii, ib
REAL :: resolution(nbpt,2)
!
- ! These two parameters are part of the routing_reg module saved variables. They are transfered here.
+ ! nbvmax is still used to dimension wome variables in routing_reg.f90.
+ ! It is transfered here but should be argument to the various subroutines.
!
nbvmax = nbvmax_in
- nbxmax = nbxmax_in
+ !
+ WRITE(numout,*) "Memory Mgt getgrid : nbvmax, ijdimmax = ", nbvmax, ijdimmax
!
DO ii=1,nbpt
resolution(ii,1) = SQRT(area(ii))
@@ -103,7 +105,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
ENDDO
DO ib=1,nbpt
- CALL routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
+ CALL routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& 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,:,:),&
@@ -113,7 +115,7 @@ SUBROUTINE gethydrogrid(nbpt, nbvmax_in, nbxmax_in, sub_pts, sub_index, sub_area
END SUBROUTINE gethydrogrid
-SUBROUTINE findbasins(nbpt, nb_htu, nbv, nbxmax_in, nbi, nbj, trip_bx, basin_bx, fac_bx, hierarchy_bx, topoind_bx, lshead_bx, &
+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)
!
@@ -124,15 +126,15 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, nbxmax_in, nbi, nbj, trip_bx, basin_bx,
!
! Arguments
!
- INTEGER, INTENT(in) :: nb_htu, nbv, nbxmax_in
+ INTEGER, INTENT(in) :: nb_htu, nbv, ijdimmax
INTEGER, INTENT(in) :: nbpt !! Domain size (unitless)
INTEGER, INTENT(in) :: nbi(nbpt), nbj(nbpt) !! Number of point in x and y within the grid (unitless)
- INTEGER, INTENT(inout) :: trip_bx(nbpt,nbxmax_in,nbxmax_in) !! The trip field for each of the smaller boxes (unitless)
- INTEGER, INTENT(inout) :: basin_bx(nbpt,nbxmax_in,nbxmax_in) !!
- REAL, INTENT(in) :: fac_bx(nbpt,nbxmax_in,nbxmax_in) !! Flow accumulation
- REAL, INTENT(in) :: hierarchy_bx(nbpt,nbxmax_in,nbxmax_in) !! Level in the basin of the point
- REAL, INTENT(inout) :: topoind_bx(nbpt,nbxmax_in,nbxmax_in) !! Topographic index of the residence time for each of the smaller boxes (m)
- REAL, INTENT(in) :: lshead_bx(nbpt,nbxmax_in,nbxmax_in) !! Large scale heading for outflow points.
+ INTEGER, INTENT(inout) :: trip_bx(nbpt,ijdimmax,ijdimmax) !! The trip field for each of the smaller boxes (unitless)
+ INTEGER, INTENT(inout) :: basin_bx(nbpt,ijdimmax,ijdimmax) !!
+ 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(in) :: lshead_bx(nbpt,ijdimmax,ijdimmax) !! Large scale heading for outflow points.
!
!
! OUTPUT VARIABLES
@@ -149,8 +151,8 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, nbxmax_in, nbi, nbj, trip_bx, basin_bx,
!
! For debug and get coordinate of river outlet
!
- REAL, INTENT(in) :: lontmp(nbpt,nbxmax_in,nbxmax_in) !! Longitude
- REAL, INTENT(in) :: lattmp(nbpt,nbxmax_in,nbxmax_in) !! Latitude
+ REAL, INTENT(in) :: lontmp(nbpt,ijdimmax,ijdimmax) !! Longitude
+ REAL, INTENT(in) :: lattmp(nbpt,ijdimmax,ijdimmax) !! Latitude
!
! Local
!
@@ -158,11 +160,8 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, nbxmax_in, nbi, nbj, trip_bx, basin_bx,
!
!diaglalo(1,:) = (/ 39.6791, 2.6687 /)
!
- IF ( nbxmax_in .NE. nbxmax ) THEN
- WRITE(*,*) "nbxmax has changed !!"
- STOP
- ENDIF
-
+ WRITE(numout,*) "Memory Mgt findbasin : nbvmax, nb_htu, nbv = ", nbvmax, nb_htu, nbv
+
DO ib=1,nbpt
CALL routing_reg_findbasins(nb_htu, nbv, ib, nbi(ib), nbj(ib), trip_bx(ib,:,:), &
& basin_bx(ib,:,:), fac_bx(ib,:,:), hierarchy_bx(ib,:,:), &
@@ -170,10 +169,10 @@ SUBROUTINE findbasins(nbpt, nb_htu, nbv, nbxmax_in, nbi, nbj, trip_bx, basin_bx,
& 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, nbxmax_in, 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, &
& 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, &
@@ -190,16 +189,16 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, nbxmax_in, area_bx, lon_bx, lat_bx, trip
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
- INTEGER(i_std), INTENT(in) :: nbv,nb_htu, nbxmax_in
- REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: area_bx !! Area of each small box in the grid box (m^2)
- REAL(r_std), INTENT (in), DIMENSION(nbpt,nbxmax_in,nbxmax_in) :: lon_bx !! Longitude of each small box in the grid box
- 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)
+ INTEGER(i_std), INTENT(in) :: nbv, nb_htu, ijdimmax
+ REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: area_bx !! Area of each small box in the grid box (m^2)
+ REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: lon_bx !! Longitude of each small box in the grid box
+ REAL(r_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: lat_bx !! Latitude of each small box in the grid box
+ INTEGER(i_std), INTENT (in), DIMENSION(nbpt,ijdimmax,ijdimmax) :: trip_bx !! The trip field for each of the smaller boxes (unitless)
+ 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)
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,nb_htu) :: basin_inbxid, basin_sz !! ID of basin, number of points in the basin
@@ -233,13 +232,10 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, nbxmax_in, area_bx, lon_bx, lat_bx, trip
!!
INTEGER(i_std) :: ib
!!
- IF ( nbxmax_in .NE. nbxmax ) THEN
- WRITE(*,*) "GLOBALIZE : nbxmax has changed !!"
- STOP
- ENDIF
+ WRITE(numout,*) "Memory Mgt globalize : nbvmax, ijdimmax, nbv, nwbas, nb_htu = ", nbvmax, ijdimmax, nbv, nwbas, nb_htu
!!
DO ib=1,nbpt
- CALL routing_reg_globalize(nbpt,nb_htu, nbv, ib, neighbours, area_bx(ib,:,:),&
+ 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,:), &
@@ -251,7 +247,7 @@ SUBROUTINE globalize(nbpt, nb_htu, nbv, nbxmax_in, area_bx, lon_bx, lat_bx, trip
!
END SUBROUTINE globalize
-SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, basin_flowdir, &
+SUBROUTINE linkup(nbpt, ijdimmax, nwbas, inflowmax, basin_count, basin_area, basin_id, basin_flowdir, &
& basin_lshead, basin_hierarchy, basin_fac, outflow_grid, outflow_basin, inflow_number, inflow_grid, &
& inflow_basin, nbcoastal, coastal_basin, invented_basins)
!
@@ -264,7 +260,7 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
- INTEGER(i_std), INTENT (in) :: nbxmax_in
+ INTEGER(i_std), INTENT (in) :: ijdimmax, inflowmax
REAL(r_std), INTENT(in) :: invented_basins !!
!
INTEGER(i_std), INTENT (in) :: nwbas !!
@@ -278,20 +274,18 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_basin !!
!
- INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax_in*8) :: inflow_number !!
- INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax_in*8,nbxmax_in*8) :: inflow_basin !!
- INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax_in*8,nbxmax_in*8) :: inflow_grid !!
+ INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas) :: inflow_number !!
+ INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_basin !!
+ INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_grid !!
!
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt) :: nbcoastal !!
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: coastal_basin !!
!
!
- IF ( nbxmax_in .NE. nbxmax ) THEN
- WRITE(*,*) "LINKUP : nbxmax has changed !!"
- STOP
- ENDIF
- CALL routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, &
+ WRITE(numout,*) "Memory Mgt Linkup : nbvmax, ijdimmax, nwbas, inflowmax = ", nbvmax, ijdimmax, nwbas, inflowmax
+
+ CALL routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basin_count, basin_area, basin_id, basin_flowdir, &
& basin_lshead, basin_hierarchy, basin_fac, diaglalo, outflow_grid, outflow_basin, inflow_number, inflow_grid, &
& inflow_basin, nbcoastal, coastal_basin, invented_basins)
!
@@ -431,7 +425,7 @@ 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, &
+SUBROUTINE finish_truncate(nbpt, inflowmax, 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, &
@@ -447,7 +441,7 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
- INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax
+ INTEGER(i_std), INTENT (in) :: inflowmax, nbasmax
INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), INTENT (in) :: num_largest
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea
@@ -468,11 +462,11 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
!
- ! Changed nwbas to nbxmax*4 and *8 to save the memory
!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in), INTENT(inout) :: inflow_number !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_basin !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_grid !!
+ !
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: inflow_number !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_basin !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid !!
!
! Output variables
!
@@ -494,7 +488,7 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
!
!
- CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, num_largest, &
+ CALL routing_reg_end_truncate(nbpt, nbasmax, gridarea, cfrac, nwbas, inflowmax, 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)
@@ -519,7 +513,7 @@ SUBROUTINE finish_truncate(nbpt, nbxmax_in, nbasmax, nwbas, num_largest, gridare
END SUBROUTINE finish_truncate
-SUBROUTINE finish_inflows(nbpt,nbxmax_in, nbasmax, inf_max, basin_count, inflow_number, inflow_grid, inflow_basin,&
+SUBROUTINE finish_inflows(nbpt, nwbas, nbasmax, inf_max, basin_count, inflow_number, inflow_grid, inflow_basin,&
& route_innum, route_ingrid, route_inbasin)
!
USE ioipsl
@@ -530,15 +524,15 @@ SUBROUTINE finish_inflows(nbpt,nbxmax_in, nbasmax, inf_max, basin_count, inflow_
! Arguments
!
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
- INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax
+ INTEGER(i_std), INTENT (in) :: nwbas, nbasmax
INTEGER(i_std), INTENT (in) :: inf_max !!
!
- INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in), INTENT(in) :: inflow_number !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,inf_max), INTENT(in) :: inflow_basin !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,inf_max), INTENT(in) :: inflow_grid !!
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: inflow_number !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inf_max), INTENT(in) :: inflow_basin !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inf_max), INTENT(in) :: inflow_grid !!
!
- REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_innum
+ REAL(r_std), DIMENSION(nbpt,nbasmax), INTENT(out) :: route_innum
REAL(r_std), DIMENSION(nbpt,nbasmax,inf_max), INTENT(out) :: route_ingrid
REAL(r_std), DIMENSION(nbpt,nbasmax,inf_max), INTENT(out) :: route_inbasin
!
@@ -567,7 +561,7 @@ SUBROUTINE finish_inflows(nbpt,nbxmax_in, nbasmax, inf_max, basin_count, inflow_
END SUBROUTINE finish_inflows
-SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numops, basin_count, basin_area, &
+SUBROUTINE killbas(nbpt, inflowmax, 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)
@@ -580,7 +574,7 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
- INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax
+ INTEGER(i_std), INTENT (in) :: inflowmax, nbasmax
INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), INTENT (in) :: ops !!
INTEGER(i_std), DIMENSION(nbpt, ops), INTENT (in) :: tokill, totakeover
@@ -600,11 +594,11 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
!
- ! Changed nwbas to nbxmax*4 and *8 to save the memory
+ !
!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in), INTENT(inout) :: inflow_number !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_basin !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax_in,nbxmax_in), INTENT(inout) :: inflow_grid !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: inflow_number !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_basin !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid !!
!! LOCAL
INTEGER(i_std) :: ib, op, tok, totak, igrif, ibasf
@@ -633,7 +627,7 @@ SUBROUTINE killbas(nbpt, nbxmax_in, nbasmax, nwbas,ops, tokill, totakeover, numo
igrif = it
END IF
END DO
- CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, basin_count, basin_area, &
+ CALL routing_reg_killbas(nbpt, ib, tok, totak, nwbas, inflowmax, 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)
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 4c061942a47f2772df237228d43dd48f67fb4274..1ed64bc3377bba2084ec7b1609b01b34f42fb334 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -8,7 +8,6 @@ MODULE routing_reg
IMPLICIT NONE
INTEGER(i_std), SAVE :: nbvmax=440 !! The maximum number of basins we can handle at any time during the generation of the maps (unitless)
- INTEGER(i_std), SAVE :: nbxmax=63 !! The maximum number of points in one direction (NS or WE) using in routing_reg_getgrid (unitless)
INTEGER(i_std), SAVE :: maxpercent=2 !! The maximum area percentage of a sub-basin in a grib-box (default = 2%)
INTEGER(i_std), SAVE :: nbpt_save
@@ -73,7 +72,7 @@ CONTAINS
!! \n
!_ ================================================================================================================================
- SUBROUTINE routing_reg_getgrid(nbpt, ib, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
+ SUBROUTINE routing_reg_getgrid(nbpt, ib, ijdimmax, sub_pts, sub_index, sub_area, max_basins, min_topoind, &
& 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)
@@ -83,6 +82,7 @@ CONTAINS
!! INPUT VARIABLES
INTEGER(i_std), INTENT(in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT(in) :: ib !! Current basin (unitless)
+ INTEGER(i_std), INTENT(in) :: ijdimmax !! Maximum dimension in i or j of the underlaying hydro grid.
INTEGER(i_std), INTENT(in) :: sub_pts(nbpt) !! Number of high resolution points on this grid (unitless)
INTEGER(i_std), INTENT(in) :: sub_index(nbpt,nbvmax,2) !! Indices of the points we need on the fine grid (unitless)
REAL(r_std), INTENT(inout) :: max_basins !! The current maximum of basins
@@ -106,25 +106,25 @@ CONTAINS
!
!! OUTPUT VARIABLES
INTEGER(i_std), INTENT(out) :: nbi, nbj !! Number of point in x and y within the grid (unitless)
- REAL(r_std), INTENT(out) :: area_bx(nbxmax,nbxmax) !! Area of each small box in the grid box (m^2)
- REAL(r_std), INTENT(out) :: hierarchy_bx(nbxmax,nbxmax) !! Level in the basin of the point
- REAL(r_std), INTENT(out) :: fac_bx(nbxmax,nbxmax) !! Flow accumulation
- REAL(r_std), INTENT(out) :: lon_bx(nbxmax,nbxmax) !!
- REAL(r_std), INTENT(out) :: lat_bx(nbxmax,nbxmax) !!
- REAL(r_std), INTENT(out) :: lshead_bx(nbxmax,nbxmax) !! Large scale heading for outflow points.
- 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) !!
+ REAL(r_std), INTENT(out) :: area_bx(ijdimmax,ijdimmax) !! Area of each small box in the grid box (m^2)
+ REAL(r_std), INTENT(out) :: hierarchy_bx(ijdimmax,ijdimmax) !! Level in the basin of the point
+ REAL(r_std), INTENT(out) :: fac_bx(ijdimmax,ijdimmax) !! Flow accumulation
+ REAL(r_std), INTENT(out) :: lon_bx(ijdimmax,ijdimmax) !!
+ 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)
+ 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) !!
+ REAL(i_std), INTENT(out) :: floodp_bx(ijdimmax,ijdimmax) !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ip, jp, ll(1), iloc, jloc !! Indices (unitless)
INTEGER(i_std) :: ipp, jpp
INTEGER(i_std), DIMENSION(8,2) :: inc
REAL(r_std) :: cenlon, cenlat, dlon, dlat, deslon, deslat, facti, factj
- REAL(r_std) :: lonstr(nbxmax*nbxmax) !!
- REAL(r_std) :: latstr(nbxmax*nbxmax) !!
+ REAL(r_std) :: lonstr(ijdimmax*ijdimmax) !!
+ REAL(r_std) :: latstr(ijdimmax*ijdimmax) !!
!
!_ ================================================================================================================================
@@ -167,9 +167,9 @@ CONTAINS
IF ( sub_pts(ib) > 0 ) THEN
!
DO ip=1,sub_pts(ib)
- IF ( ip >nbxmax*nbxmax ) THEN
- CALL ipslerr_p(3,'routing_reg_getgrid','nbxmax too small when filling lonstr',&
- & 'Please change method to estimate nbxmax','')
+ IF ( ip >ijdimmax*ijdimmax ) THEN
+ CALL ipslerr_p(3,'routing_reg_getgrid','ijdimmax too small when filling lonstr',&
+ & 'Please change method to estimate ijdimmax','')
ENDIF
lonstr(ip) = lon_rel(ib, ip)
latstr(ip) = lat_rel(ib, ip)
@@ -182,9 +182,9 @@ CONTAINS
!
! Verify dimension and allocated space
!
- IF ( nbi > nbxmax .OR. nbj > nbxmax ) THEN
- WRITE(numout,*) "size of area : nbi=",nbi,"nbj=",nbj, "nbxmax=", nbxmax
- CALL ipslerr_p(3,'routing_reg_getgrid','nbxmax too small','Please change method to estimate nbxmax','')
+ IF ( nbi > ijdimmax .OR. nbj > ijdimmax ) THEN
+ WRITE(numout,*) "size of area : nbi=",nbi,"nbj=",nbj, "ijdimmax=", ijdimmax
+ CALL ipslerr_p(3,'routing_reg_getgrid','ijdimmax too small','Please change method to estimate ijdimmax','')
ENDIF
!
! Transfer the data in such a way that (1,1) is the North Western corner and
@@ -1649,8 +1649,8 @@ CONTAINS
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, neighbours, area_bx, lon_bx, lat_bx, trip_bx, hierarchy_bx, &
- & orog_bx, floodp_bx, fac_bx, topoind_bx, &
+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, &
& 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, &
@@ -1662,18 +1662,19 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, neighbours, area_bx, lon
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt, nb_htu, nbv !! Domain size (unitless)
INTEGER(i_std), INTENT (in) :: ib !! Current basin (unitless)
+ INTEGER(i_std), INTENT (in) :: ijdimmax !! Maximum dimension in i or j of the underlaying hydro grid.
INTEGER(i_std), INTENT(in) :: neighbours(nbpt,NbNeighb)!! Vector of neighbours for each grid point
!! (1=North and then clockwise)
!! LOCAL VARIABLES
- REAL(r_std), DIMENSION(nbxmax,nbxmax) :: area_bx !! Area of each small box in the grid box (m^2)
- REAL(r_std), DIMENSION(nbxmax,nbxmax) :: lon_bx !! Longitude of each small box in the grid box
- 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), DIMENSION(ijdimmax,ijdimmax) :: area_bx !! Area of each small box in the grid box (m^2)
+ REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: lon_bx !! Longitude of each small box in the grid box
+ REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: lat_bx !! Latitude of each small box in the grid box
+ INTEGER(i_std), DIMENSION(ijdimmax,ijdimmax) :: trip_bx !! The trip field for each of the smaller boxes (unitless)
+ REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: hierarchy_bx !! Level in the basin of the point
+ REAL(r_std), DIMENSION(ijdimmax,ijdimmax) :: orog_bx !! Level in the basin of the point
+ 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
@@ -1966,7 +1967,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, neighbours, area_bx, lon
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area, basin_id, basin_flowdir, &
+SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basin_count, basin_area, basin_id, basin_flowdir, &
& basin_lshead, basin_hierarchy, basin_fac, diaglalo, outflow_grid, outflow_basin, inflow_number, inflow_grid, &
& inflow_basin, nbcoastal, coastal_basin, invented_basins)
!
@@ -1975,6 +1976,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT (in), DIMENSION(nbpt,NbNeighb) :: neighbours !!
+ INTEGER(i_std), INTENT (in) :: ijdimmax, inflowmax
REAL(r_std), INTENT(in) :: invented_basins !!
!
INTEGER(i_std), INTENT (in) :: nwbas !!
@@ -1989,9 +1991,9 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: outflow_basin !!
!
- INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax*8) :: inflow_number !!
- INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_basin !!
- INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_grid !!
+ INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas) :: inflow_number !!
+ INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_basin !!
+ INTEGER(i_std), INTENT(out), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_grid !!
!
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt) :: nbcoastal !!
INTEGER(i_std), INTENT(inout), DIMENSION(nbpt,nwbas) :: coastal_basin !!
@@ -2078,7 +2080,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
IF ( outflow_grid(sp,sb) .EQ. -4 ) THEN
! Flow into a basin of the same grid
bop = outflow_basin(sp,sb)
- CALL routing_updateflow(sp, sb, sp, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ CALL routing_updateflow(sp, sb, sp, bop, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,1) = solved(sp,1) + 1
@@ -2115,7 +2117,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
!
IF ( bop .LT. undef_int ) THEN
!
- CALL routing_updateflow(sp, sb, inp, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ CALL routing_updateflow(sp, sb, inp, bop, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,1) = solved(sp,1) + 1
@@ -2233,7 +2235,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
IF ( dop > 0 .AND. dop < undef_int .AND. bop < undef_int ) THEN
!
- CALL routing_updateflow(sp, sb, dop, bop, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ CALL routing_updateflow(sp, sb, dop, bop, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid, inflow_basin)
!
IF ( outflow_basin(sp,sb) == bop ) THEN
@@ -2303,7 +2305,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
IF (( sbint < undef_int ) .AND. (sbint .GT. 0) .AND. (sbint .NE. sb)) THEN
basin_hierarchy(sp,sb) = basin_hierarchy(sp,sbint)
- CALL routing_updateflow(sp, sb, sp, sbint, nbpt, nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ CALL routing_updateflow(sp, sb, sp, sbint, nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, &
& inflow_number, inflow_grid,inflow_basin)
found = 1
WRITE (numout,*) "Lowest basin hierarchy in the grid file"
@@ -2321,7 +2323,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
bop = minhloc(basin_id(sp,sb),2)
IF ((dop < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0)) THEN
- CALL routing_updateflow(sp, sb, dop, bop, nbpt,nwbas, nbxmax*8, outflow_grid, outflow_basin, &
+ CALL routing_updateflow(sp, sb, dop, bop, nbpt,nwbas, inflowmax, outflow_grid, outflow_basin, &
& 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
@@ -2583,7 +2585,7 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, nbcore, corepts, b
!! \n
!_ ================================================================================================================================
-SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_largest, &
+SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, inflowmax, 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)
@@ -2601,7 +2603,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: gridarea !! The size of each grid box in X and Y (m)
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box (unitless;0-1)
!
- INTEGER(i_std), INTENT(in) :: nwbas !!
+ INTEGER(i_std), INTENT(in) :: nwbas, inflowmax !!
INTEGER(i_std), INTENT (in) :: num_largest
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
@@ -2619,11 +2621,10 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, nu
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_basin !!
!
- ! Changed nwbas to nbxmax*4 and *8 to save the memory
!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax*8), INTENT(inout) :: inflow_number !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8), INTENT(inout) :: inflow_basin !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8), INTENT(inout) :: inflow_grid !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: inflow_number !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_basin !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(inout) :: inflow_grid !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij, ic, ibf, ijf, igrif, ibasf, cnt, pold, bold, ff(2) !! Indices (unitless)
@@ -2789,7 +2790,7 @@ 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, &
+SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, 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)
@@ -2802,7 +2803,7 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
INTEGER(i_std) :: nbpt !! Domain size (unitless)
INTEGER(i_std) :: ib !! Current basin (unitless)
!
- INTEGER(i_std) :: nwbas !!
+ INTEGER(i_std) :: nwbas, inflowmax !!
INTEGER(i_std), DIMENSION(nbpt) :: basin_count !!
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: basin_id !!
REAL(r_std), DIMENSION(nbpt,nwbas,2) :: basin_coor !!
@@ -2817,11 +2818,10 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_grid !! Type of outflow on the grid box (unitless)
INTEGER(i_std), DIMENSION(nbpt,nwbas) :: outflow_basin !!
!
- ! Saving memory: changed nwbas to nbxmax*4 and *8
!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax*8) :: inflow_number !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_basin !!
- INTEGER(i_std), DIMENSION(nbpt,nbxmax*8,nbxmax*8) :: inflow_grid !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas) :: inflow_number !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_basin !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax) :: inflow_grid !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: inf, ibs, ing, inb, ibasf, igrif, it !! Indices (unitless)
@@ -2832,6 +2832,8 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
! Update the information needed in the basin "totakeover"
! For the moment only area
!
+ WRITE(numout, *) "XXXXX Arguments : inflow_number = ", inflow_number(ib,1:10)
+ !
basin_cg(ib, totakeover, 1) = (basin_cg(ib, totakeover, 1)*basin_area(ib, totakeover) &
& + basin_cg(ib, tokill, 1)*basin_area(ib, tokill))/(basin_area(ib, totakeover) + basin_area(ib, tokill))
!
@@ -2880,10 +2882,12 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
! Redirect the flows which went into the basin to be killed before we change everything
!
DO inf = 1, inflow_number(ib, tokill)
- outflow_basin(inflow_grid(ib, tokill, inf), inflow_basin(ib, tokill, inf)) = totakeover
- inflow_number(ib, totakeover) = inflow_number(ib, totakeover) + 1
- inflow_grid(ib, totakeover, inflow_number(ib, totakeover)) = inflow_grid(ib, tokill, inf)
- inflow_basin(ib, totakeover, inflow_number(ib, totakeover)) = inflow_basin(ib, tokill, inf)
+ IF ( inflow_grid(ib, tokill, inf) > 0 .AND. inflow_basin(ib, tokill, inf) > 0 ) THEN
+ outflow_basin(inflow_grid(ib, tokill, inf), inflow_basin(ib, tokill, inf)) = totakeover
+ inflow_number(ib, totakeover) = inflow_number(ib, totakeover) + 1
+ inflow_grid(ib, totakeover, inflow_number(ib, totakeover)) = inflow_grid(ib, tokill, inf)
+ inflow_basin(ib, totakeover, inflow_number(ib, totakeover)) = inflow_basin(ib, tokill, inf)
+ ENDIF
ENDDO
!
! Take out the basin to be killed from the list of inflow basins of the downstream basin
@@ -2957,7 +2961,17 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, basin_count,
!
DO ibs=tokill+1, basin_count(ib)
DO inf = 1, inflow_number(ib, ibs)
- outflow_basin(inflow_grid(ib, ibs, inf), inflow_basin(ib, ibs, inf)) = ibs-1
+ IF (inflow_grid(ib, ibs, inf) < 1 .OR. inflow_grid(ib, ibs, inf) > nbpt) THEN
+ WRITE(numout,*) "XXXX ib, tokill+1, basin_count : ", ib, tokill+1, basin_count(ib)
+ WRITE(numout,*) "XXXX ibs, inf, _grid =", ibs, inf, inflow_grid(ib, ibs, inf)
+ ENDIF
+ IF (inflow_basin(ib, ibs, inf) < 1 .OR. inflow_basin(ib, ibs, inf) > nwbas ) THEN
+ WRITE(numout,*) "XXXX ib, tokill+1, basin_count : ", ib, tokill+1, basin_count(ib)
+ WRITE(numout,*) "XXXX ibs, inf, _basin =", ib, ibs, inf, inflow_basin(ib, ibs, inf)
+ ENDIF
+ IF (inflow_grid(ib, ibs, inf) > 0 .AND. inflow_grid(ib, ibs, inf) > 0) THEN
+ outflow_basin(inflow_grid(ib, ibs, inf), inflow_basin(ib, ibs, inf)) = ibs-1
+ ENDIF
ENDDO
ENDDO
!
diff --git a/Interface.py b/Interface.py
index 14ce51df7e0792dca92f9441ec18f92778e2554b..6055acb0aaafb4495801bbfb2670db5b50e68100 100644
--- a/Interface.py
+++ b/Interface.py
@@ -262,15 +262,23 @@ class HydroOverlap :
trip_tmp[np.isnan(trip_tmp)] = undef_int
basins_tmp[np.isnan(trip_tmp)] = undef_int
#
+ # Compute nbxmax
+ #
+ ijdim=[]
+ for ib in range(nbpt) :
+ ijdim.append(max(max(sub_index[ib,:,0])-min(sub_index[ib,:,0])+1,max(sub_index[ib,:,1])-min(sub_index[ib,:,1])+1))
+ ijdimmax = max(ijdim)
+ #
# Go to the call of the FORTRAN interface
#
print("GETHYDROGRID : nbpt = ", nbpt, nbvmax)
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, \
+ routing_interface.gethydrogrid(ijdimmax, sub_pts, sub_index, sub_area, \
hydrodata.basinsmax, hydrodata.topoindmin, sub_lon, sub_lat, trip_tmp, basins_tmp, topoind_tmp, fac_tmp,\
hierarchy_tmp, orog_tmp, floodp_tmp)
#
@@ -296,15 +304,14 @@ class HydroSuper :
#
# nb_htu can be adjusted with self.nwbas
# nb_htu can be lowered with a larger maxpercent (routing_reg.f90)
- nb_htu = 600
nb_htu = nbvmax
nbv = nbvmax
#
# Call findbasins
#
nb_basin, basin_inbxid, basin_outlet, basin_outtp, self.basin_sz, basin_bxout, basin_bbout, self.basin_pts, basin_lshead, coast_pts = \
- routing_interface.findbasins(nbpt = self.nbpt, nb_htu = nb_htu, nbv = nbv, nbi = hydrooverlap.nbi, nbj = hydrooverlap.nbj, \
- trip_bx = hydrooverlap.trip_bx, \
+ routing_interface.findbasins(nbpt = self.nbpt, nb_htu = self.nbhtuext, nbv = nbv, nbi = hydrooverlap.nbi, \
+ nbj = hydrooverlap.nbj, trip_bx = hydrooverlap.trip_bx, \
basin_bx = hydrooverlap.basin_bx, fac_bx = hydrooverlap.fac_bx, \
hierarchy_bx = hydrooverlap.hierarchy_bx, \
topoind_bx = hydrooverlap.topoind_bx, lshead_bx = hydrooverlap.lshead_bx, \
@@ -313,7 +320,13 @@ class HydroSuper :
# Adjust nwbas to the maximum found over the domain
#
self.nwbas = part.domainmax(np.max(nb_basin))
+ # Set the number of inflows per basin. For the moment twice the maximum number of basins.
+ self.inflowmax = self.nwbas*2
print("Maximum number of basin created : {0}".format(self.nwbas))
+ ijdim=[]
+ for i in range(self.nbpt) :
+ ijdim.append(max(hydrooverlap.area_bx[i,:,:].shape))
+ self.ijdimmax = max(ijdim)
#
# Call Globalize
#
@@ -321,12 +334,13 @@ 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_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(nbpt = self.nbpt, nb_htu = nb_htu,nbv = nbv, area_bx = hydrooverlap.area_bx, lon_bx = lon_bx_tmp, \
- lat_bx = lat_bx_tmp, trip_bx = hydrooverlap.trip_bx, \
+ routing_interface.globalize(nbpt = self.nbpt, nb_htu = self.nbhtuext, nbv = nbv, ijdimmax = self.ijdimmax, \
+ area_bx = hydrooverlap.area_bx, lon_bx = lon_bx_tmp, lat_bx = lat_bx_tmp, trip_bx = hydrooverlap.trip_bx, \
hierarchy_bx = hydrooverlap.hierarchy_bx, orog_bx = hydrooverlap.orog_bx, floodp_bx = hydrooverlap.floodp_bx,\
fac_bx = hydrooverlap.fac_bx, topoind_bx = hydrooverlap.topoind_bx, min_topoind = hydrodata.topoindmin, \
nb_basin = nb_basin, basin_inbxid = basin_inbxid, basin_outlet = basin_outlet, basin_outtp = basin_outtp, \
@@ -343,11 +357,15 @@ class HydroSuper :
# Call the linkup routine in routing_reg.
#
print("Invented basins =", hydrodata.basinsmax)
- self.inflow_number,self.inflow_grid,self.inflow_basin = routing_interface.linkup(nbxmax, self.basin_count, self.basin_area, self.basin_id, \
- self.basin_flowdir, self.basin_lshead, self.basin_hierarchy, \
- self.basin_fac, self.outflow_grid, self.outflow_basin, \
- self.nbcoastal, self.coastal_basin, float(hydrodata.basinsmax))
+ self.inflow_number,self.inflow_grid,self.inflow_basin = \
+ routing_interface.linkup(self.ijdimmax, self.inflowmax, self.basin_count, self.basin_area, \
+ self.basin_id, self.basin_flowdir, self.basin_lshead, self.basin_hierarchy, \
+ self.basin_fac, self.outflow_grid, self.outflow_basin, \
+ self.nbcoastal, self.coastal_basin, float(hydrodata.basinsmax))
self.nbxmax_in = self.inflow_number.shape[1]
+ #
+ #
+ #
return
#
@@ -415,8 +433,10 @@ class HydroSuper :
#
def killbas(self, tokill, totakeover, numops):
ops = tokill.shape[1]
-
- routing_interface.killbas(nbpt = self.nbpt, nbxmax_in = self.nbxmax_in, \
+ #
+ nbxmax_tmp = self.inflow_grid.shape[2]
+ #
+ routing_interface.killbas(nbpt = self.nbpt, inflowmax = nbxmax_tmp, \
nbasmax = self.nbasmax, nwbas = self.nwbas, ops = ops, tokill = tokill,\
totakeover = totakeover, numops = numops, basin_count = self.basin_count,\
basin_area = self.basin_area, basin_orog = self.basin_orog, basin_floodp = self.basin_floodp, \
@@ -551,12 +571,13 @@ class HydroGraph :
self.nbasmax = nbasmax
self.nbpt = hydrosuper.basin_count.shape[0]
nwbas = hydrosuper.basin_topoind.shape[1]
- nbxmax_in = hydrosuper.inflow_grid.shape[1]
+ nbxmax_in = hydrosuper.inflow_grid.shape[2]
+ #
#
self.routing_area, self.routing_orog, self.routing_floodp, self.routing_cg, self.topo_resid, self.route_nbbasin,\
self.route_togrid, self.route_tobasin, self.route_nbintobas, self.global_basinid, \
self.route_outlet, self.route_type, self.origin_nbintobas, self.routing_fetch = \
- routing_interface.finish_truncate(nbpt = self.nbpt, nbxmax_in = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
+ routing_interface.finish_truncate(nbpt = self.nbpt, inflowmax = nbxmax_in, nbasmax = nbasmax, nwbas = nwbas, \
num_largest = hydrosuper.num_largest, gridarea = modelgrid.area, \
cfrac = modelgrid.contfrac, basin_count = hydrosuper.basin_count, \
basin_notrun = hydrosuper.basin_notrun, basin_area = hydrosuper.basin_area, \
@@ -579,12 +600,10 @@ class HydroGraph :
# Inflows
self.max_inflow = part.domainmax(np.max(hydrosuper.inflow_number))
gingrid = part.l2glandindex( hydrosuper.inflow_grid[:,:,:self.max_inflow])
- self.route_innum, self.route_ingrid, self.route_inbasin = routing_interface.finish_inflows(nbpt = self.nbpt, nbxmax_in = nbxmax_in, \
- nbasmax = nbasmax, inf_max = self.max_inflow, \
- basin_count = hydrosuper.basin_count, \
- inflow_number = hydrosuper.inflow_number, \
- inflow_grid = gingrid, \
- inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
+ self.route_innum, self.route_ingrid, self.route_inbasin = \
+ routing_interface.finish_inflows(nbpt = self.nbpt, nwbas = nwbas, nbasmax = nbasmax, inf_max = self.max_inflow, \
+ basin_count = hydrosuper.basin_count, inflow_number = hydrosuper.inflow_number, \
+ inflow_grid = gingrid, inflow_basin = hydrosuper.inflow_basin[:,:,:self.max_inflow])
return
#
diff --git a/tests/Iberia_n48/BuildHTUs_IP.pbs b/tests/Iberia_n48/BuildHTUs_IP.pbs
index 11dbcddc5954f459a62a920625ea7b647cf89bca..697f13fca5610b820d127308ed10ab5ee3d1ebc6 100644
--- a/tests/Iberia_n48/BuildHTUs_IP.pbs
+++ b/tests/Iberia_n48/BuildHTUs_IP.pbs
@@ -19,6 +19,19 @@ source ../../Environment
#
/bin/rm -f DocumentationInterface *graph.nc *.txt
#
+# If weights file does not exist, then compute it quickly
+#
+if [ ! -e Weights.nc ] ; then
+ mpirun -n ${NSLOTS} python ../../WeightsOnly.py
+ if [ $? -gt 0 ] ; then
+ exit
+ else
+ echo "============================================="
+ echo "= Weight calculations successful ="
+ echo "============================================="
+ fi
+fi
+#
# Run the Python code to generate the HTUs and write them into a netCDF file.
#
mpirun -n ${NSLOTS} python ../../RoutingPreProc.py
diff --git a/tests/Iberia_n48_regular/BuildHTUs_IP.pbs b/tests/Iberia_n48_regular/BuildHTUs_IP.pbs
index 8e9e843b684a7bbe6648e3b17c22da11d008c938..888578004c015425a6fdc5b2038f7f2e5b515b01 100644
--- a/tests/Iberia_n48_regular/BuildHTUs_IP.pbs
+++ b/tests/Iberia_n48_regular/BuildHTUs_IP.pbs
@@ -3,10 +3,10 @@
#PBS -N BuildHTUs_IP
#
#PBS -j oe
-#PBS -l nodes=1:ppn=48
+#PBS -l nodes=1:ppn=64
#PBS -l walltime=12:00:00
-#PBS -l mem=160gb
-#PBS -l vmem=160gb
+#PBS -l mem=240gb
+#PBS -l vmem=240gb
#
cd ${PBS_O_WORKDIR}
export NSLOTS=$(($PBS_NUM_NODES*$PBS_NUM_PPN))
@@ -20,6 +20,7 @@ source ../../Environment
/bin/rm -f DocumentationInterface *graph.nc *.txt
#
# Run the Python code to generate the HTUs and write them into a netCDF file.
+# It will also generate the weights if needed.
#
mpirun -n ${NSLOTS} python ../../RoutingPreProc.py
if [ $? -gt 0 ] ; then