Commit 47b7a8a4 authored by Anthony Schrapffer's avatar Anthony Schrapffer
Browse files

Merge branch 'master' of gitlab.in2p3.fr:jan.polcher/routingpp

parents 0c98f5c9 8317f5fb
......@@ -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
......
......@@ -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))
......
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