diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index 6540154cb18838d0f002d235d50c972ac7f74582..4c7b94445b23ff82c23891f4aeff41980e291c4b 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2030,7 +2030,6 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
DO sp=1,nbpt
!
DO sb=1,basin_count(sp)
- !
!
! We only work on this point if it does not flow into the ocean
! or flow to another sub-basin in the same grid box.
@@ -2087,6 +2086,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
ENDIF
!
IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
+ IF ( outflow_grid(sp,sb) == 0 ) WRITE(numout,*) "Linkup 1.0 -- Flow out of Halo zone"
WRITE(numout,*) "Linkup 1.0 - In routing_reg_linkup working on sp & sb:",sp, sb
WRITE(numout,*) "Linkup 1.0 - Coords : ", lalo_g(sp,2), lalo_g(sp,1)
WRITE(numout,*) "Linkup 1.0 - outflow_flowdir = ", basin_flowdir(sp,sb)
@@ -2388,8 +2388,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
WRITE(numout,*) "Linkup 4.0 - Tested outflow basin = ", gridbasin(1:nbtotest)
WRITE(numout,*) "Linkup 4.0 - outflow_grid =", outflow_grid(sp,sb)
WRITE(numout,*) "Linkup 4.0 - outflow_basin = ", outflow_basin(sp,sb)
- WRITE(numout,*) "Linkup 4.0 - Number of ocean options : ", nbocean
- IF (nbocean > 0 ) WRITE(numout,*) "Linkup 4.0 - ocean direction = ", wocean/nbocean
+ WRITE(numout,*) "Linkup 4.0 - Number of ocean options : ", nbocean, nbtotest
IF ( basin_id(sp,sb) < invented_basins ) THEN
WRITE(numout,*) "Linkup 4.0 - Hmin for basin :", hatoutflow(basin_id(sp,sb))
WRITE(numout,*) "Linkup 4.0 - Hmin location :", minhloc(basin_id(sp,sb),1), minhloc(basin_id(sp,sb),2)
@@ -2423,12 +2422,17 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
! In this phase we should also consider the size of the basin. If it small then we can
! neglect its and put it into the ocean.
!
- IF ( wocean/nbocean < 90.0 .OR. (basin_area(sp,sb)/area_g(sp)*100 < 0.5) ) THEN
+ IF ( wocean/nbocean < 90.0 .OR. (basin_area(sp,sb)/area_g(sp)*100 < 1.0) ) THEN
outflow_grid(sp,sb) = -2
outflow_basin(sp,sb) = undef_int
solved(sp,4) = solved(sp,4)+1
IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 4.0 - ", sp, sb, "becomes a coastal flow basin."
+ WRITE(numout,*) "Linkup 4.1 - In routing_reg_linkup working on sp & sb:", sp, sb
+ WRITE(numout,*) "Linkup 4.1 - basin_id = ", basin_id(sp,sb)
+ WRITE(numout,*) "Linkup 4.1 - hierarchy = ", basin_hierarchy(sp,sb)
+ WRITE(numout,*) "Linkup 4.1 - crit =", crit
+ WRITE(numout,*) "Linkup 4.1 - ocean direction = ", wocean/nbocean
+ WRITE(numout,*) "Linkup 4.1 - other criteria : ", basin_area(sp,sb)/area_g(sp)*100
ENDIF
ENDIF
ENDIF
@@ -2441,12 +2445,12 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
!
! These ocean points should not be too far away from the large scale direction.
!
- IF ( wocean/nbocean < 90.0 .OR. (basin_area(sp,sb)/area_g(sp)*100 < 0.5) ) THEN
+ IF ( wocean/nbocean < 180.0 .OR. (basin_area(sp,sb)/area_g(sp)*100 < 25.0) ) THEN
outflow_grid(sp,sb) = -2
outflow_basin(sp,sb) = undef_int
solved(sp,4) = solved(sp,4)+1
IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
- WRITE(numout,*) "Linkup 4.0 - ", sp, sb, "becomes a coastal flow basin. Relaxed condition"
+ WRITE(numout,*) "Linkup 4.1 - ", sp, sb, "becomes a coastal flow basin. Relaxed condition"
ENDIF
ENDIF
ENDIF
@@ -2463,6 +2467,12 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
! We probably have not yet solved all points. So we go through all points and basins again
! to see if we cannot find some local solution as a last resort.
!
+ IF ( debug .AND. routing_diagbox_g(sp,diaglalo) ) THEN
+ WRITE(numout,*) "Linkup 5.0 - In routing_reg_linkup working on sp & sb:", sp, sb
+ WRITE(numout,*) "Linkup 5.0 - outflow_basin :", outflow_basin(sp,sb)
+ WRITE(numout,*) "Linkup 5.0 - outflow_grid :", outflow_grid(sp,sb)
+ WRITE(numout,*) "Linkup 5.0 - basin_flowdir :", basin_flowdir(sp,sb)
+ ENDIF
!
DO sp=1,nbpt
!
@@ -2646,7 +2656,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, basin_count, basin_area,
!
inp = outflow_grid(sp,sb)
sbl = outflow_basin(sp,sb)
- IF ( inp .GE. 0 ) THEN
+ IF ( inp .GT. 0 ) THEN
IF ( basin_count(inp) .LT. sbl ) THEN
WRITE(numout,*) 'Point :', sp, ' Basin :', sb
WRITE(numout,*) 'Flows into point :', inp, ' basin :', sbl
@@ -3367,12 +3377,14 @@ SUBROUTINE routing_reg_truncate(nbpt, nbasmax, gridarea, contfrac, nwbas, num_la
DO ij=1,basin_count(ib)
ibf = route_togrid_glo(ib,ij)
ijf = route_tobasin_glo(ib,ij)
- IF ( ijf .GT. basin_count(ibf) .AND. ijf .LE. nbasmax) THEN
- WRITE(numout,*) 'Second check'
- WRITE(numout,*) 'point :', ib, ' basin :', ij
- WRITE(numout,*) 'Flows into point :', ibf, ' basin :', ijf
- WRITE(numout,*) 'But it flows into now here as basin count = ', basin_count(ibf)
- CALL ipslerr_p(3,'routing_reg_truncate','Problem with routing..','','')
+ IF ( ibf .GT. 0 ) THEN
+ IF ( ijf .GT. basin_count(ibf) .AND. ijf .LE. nbasmax) THEN
+ WRITE(numout,*) 'Second check'
+ WRITE(numout,*) 'point :', ib, ' basin :', ij
+ WRITE(numout,*) 'Flows into point :', ibf, ' basin :', ijf
+ WRITE(numout,*) 'But it flows into now here as basin count = ', basin_count(ibf)
+ CALL ipslerr_p(3,'routing_reg_truncate','Problem with routing..','','')
+ ENDIF
ENDIF
ENDDO
ENDDO
diff --git a/ModelGrid.py b/ModelGrid.py
index a253481756013228f478c7e5f6274ff910a3c489..e525ea6d527d5670a27dd5382c4025a267898278 100644
--- a/ModelGrid.py
+++ b/ModelGrid.py
@@ -35,7 +35,7 @@ def mindist(coord,lon,lat) :
#
# Function to gather all land points but while keeping also the neighbour information.
#
-def gatherland(lon,lat,land,indP,indFi,indFj) :
+def gatherland(lon, lat, land, indP, indFi, indFj) :
nj,ni=lon.shape
coord=[]
indF_land=[]
@@ -55,6 +55,12 @@ def gatherland(lon,lat,land,indP,indFi,indFj) :
for j in range(nj) :
if (land[j,i] > 0 ) :
ntmp=[]
+ #
+ # Indices of neighbouring points are in C and thus +1 will be performed
+ # For FORTRAN.
+ # -1 will become 0 and indicate point is outside of domain.
+ # -2 will become -1 and indicate ocean.
+ #
for r in rose :
nnj=j+r[0]
nni=i+r[1]
@@ -64,7 +70,7 @@ def gatherland(lon,lat,land,indP,indFi,indFj) :
else :
ntmp.append(-2)
else :
- ntmp.append(-2)
+ ntmp.append(-1)
neighbours.append(ntmp)
return nbland, coord,neighbours,indP_land,indF_land
@@ -256,7 +262,8 @@ class ModelGrid :
# Gather all land points
#
self.nbland, self.coordll,self.neighbours,self.indP,indF = gatherland(self.lon_full,self.lat_full,self.land,ind,\
- np.reshape(indFi[:],self.lon_full.shape),np.reshape(indFj[:],self.lon_full.shape))
+ np.reshape(indFi[:],self.lon_full.shape),\
+ np.reshape(indFj[:],self.lon_full.shape))
INFO("Shape of region :"+str(ni)+" x "+str(nj)+" with nbland="+str(self.nbland))
if self.nbland < 1 or self.nbland > self.nj*self.ni :
INFO("Found impossible number of land points : "+str(self.nbland))