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

Solved some issue in routing_reg_linkup. The HTU were flowing into themselves.

parent 8cb55bad
......@@ -2073,10 +2073,6 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
! At this point any of the outflows is designated by a negative values in
! outflow_grid
!
WRITE(numout,*) "*****"
WRITE(numout,*) "Linkup 0 - sp, sb = ", sp, sb
WRITE(numout,*) "Linkup 0 - outflow_grid = ", outflow_grid(sp,sb)
!
found = 0
IF ( outflow_grid(sp,sb) == 0 ) THEN
found = 1
......@@ -2130,14 +2126,14 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
& nbpt, nwbas, inp, basin_count, basin_id, basin_hierarchy, basin_fac, &
& basin_flowdir, nbcoastal, coastal_basin, bop, bopqual)
!
IF ( bop .LT. undef_int ) THEN
IF ( bop .LT. undef_int .AND. bop .NE. sb) THEN
!
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
found = 1
WRITE(numout,*) "Solution found in the original outflow_grid"
WRITE(numout,*) "Solution found in the original outflow_grid", sb, bop
ENDIF
!
ENDIF
......@@ -2190,7 +2186,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
! we know it's wrong but it will serve to make decision later
angp1 = routing_anglediff_g(sp, dp1, basin_flowdir(sp,sb))
! Check if grid point
IF (dp1 .GT. 0) THEN
IF (dp1 .GT. 0 .AND. dp1 .NE. sp) THEN
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), &
& basin_flowdir(sp,sb), invented_basins, &
& nbpt, nwbas, dp1, basin_count, basin_id, basin_hierarchy, basin_fac, &
......@@ -2205,7 +2201,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
dm1 = neighbours_g(sp,order_ref(nb*2+1))
! we know it's wrong but it will serve to make decision later
angm1 = routing_anglediff_g(sp, dm1, basin_flowdir(sp,sb))
IF (dm1 .GT. 0) THEN
IF (dm1 .GT. 0 .AND. dm1 .NE. sp) THEN
CALL routing_reg_bestsubbasin(sp, sb, basin_id(sp,sb), basin_hierarchy(sp,sb), basin_fac(sp,sb), &
& basin_flowdir(sp,sb), invented_basins, &
& nbpt, nwbas, dm1, basin_count, basin_id, basin_hierarchy, basin_fac, &
......@@ -2256,7 +2252,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
IF ( outflow_basin(sp,sb) == bop ) THEN
solved(sp,2) = solved(sp,2) + 1
found = 1
WRITE (numout,*) "Neighbours, output found at:" , nb, "level"
WRITE (numout,*) "Neighbours, output found at:" , nb, "level. sp,sb=", sp, sb, " dop,bop=", dop, bop
ELSE
nb = nb+1
ENDIF
......@@ -2340,7 +2336,8 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
IF (( neighbours_g(sp,order_ref(nb)) == minhloc(basin_id(sp,sb),1)) .AND. (found == 0)) THEN
dop = minhloc(basin_id(sp,sb),1)
bop = minhloc(basin_id(sp,sb),2)
IF ((dop < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0)) THEN
IF ((dop < undef_int ) .AND. (dop .GT. 0) .AND. (bop < undef_int) .AND. (bop .GT. 0) &
& .AND. (dop .NE. sp) .AND. (bop .NE. sb) ) THEN
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
......@@ -2352,6 +2349,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
ELSE
WRITE(numout,*) "Lowest hierarchy may be in two different grid points"
END IF
WRITE(numout,*) "sp,sb = ", sp,sb, " dop,bop = ", dop,bop
WRITE(numout,*) "hierarch of (sp,sb)", basin_hierarchy(sp,sb)
WRITE(numout,*) "Lowest basinid hierarch", basin_hierarchy(dop,bop)
END IF
......
......@@ -183,6 +183,8 @@ def getland (geo, ist, ni, jst, nj) :
vn=list(v.name for v in geo.variables.values())
if "LANDMASK" in vn :
land=geo.variables["LANDMASK"][0,jst:jst+nj,ist:ist+ni]
elif "Contfrac" in vn :
land=geo.variables["Contfrac"][jst:jst+nj,ist:ist+ni]
elif "elevation" in vn :
land=geo.variables["elevation"][jst:jst+nj,ist:ist+ni]
if "missing_value" in geo.variables["elevation"].ncattrs() :
......
......@@ -52,7 +52,7 @@ def fit_partition(partin, land):
for i, dom in enumerate(partin):
l = land[dom["jstart"]:dom["jstart"]+dom["nbj"],dom["istart"]:dom["istart"]+dom["nbi"]]
l1 = np.sum(l, axis = 0)
i0 = np.where(l1>0)[0][0]
i1 = np.where(l1>0)[0][-1]
......
......@@ -29,6 +29,7 @@ if [ $? -gt 0 ] ; then
echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
echo "X Run MECORDEX on 2 Proc failed X"
echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
exit
else
echo "======================================"
echo "= Run MEDCORDEX on 2 Proc successful ="
......@@ -36,9 +37,27 @@ else
ls -lt
fi
#
#
# 4 Proc
#
/bin/rm -f run.def *.txt
cp run_E2OFD.def run.def
mpirun -n 2 python ../../RoutingPreProc.py
if [ $? -gt 0 ] ; then
echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
echo "X Run E2OFD on 4 Proc failed X"
echo "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
exit
else
echo "=================================="
echo "= Run E2OFD on 4 Proc successful ="
echo "=================================="
ls -lt
fi
#
# More points on 32 proc
#
/bin/rm -f run.def
/bin/rm -f run.def *.txt
cp run_EuroSW.def run.def
mpirun -n 32 python ../../RoutingPreProc.py
if [ $? -gt 0 ] ; then
......
[OverAll]
#
#
EarthRadius = 6370000.
#
ModelGridFile = /bdd/MEDI/workspaces/polcher/NewRouting/MED_E2OFD_Contfrac.nc
HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
# Mallorca
WEST_EAST = 2.3, 3.75
SOUTH_NORTH = 39.00, 40.1
#
# FORTRAN interface parameters
#
Documentation = true
#
# Configuration for the graph to be generated
#
nbasmax = 35
#
# Number of operation of simplification performed together
#
numop = 200
#
# Output
#
GraphFile = E2OFD_test_graph.nc
#
......@@ -10,7 +10,7 @@ HydroFile = /bdd/MEDI/workspaces/polcher/NewRouting/routing_MED.nc
WEST_EAST = 2.2, 3.6
SOUTH_NORTH = 39.20, 40.10
#
##WeightFile = EuroSW_Weights.nc
WeightFile = EuroSW_Weights.nc
#
# FORTRAN interface parameters
#
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment