Some corrections and clean up

F90subroutines/routing_reg.f90

@@ -2076,7 +2076,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
           found = 0
           IF (  outflow_grid(sp,sb) == 0 ) THEN
               found = 1
-              WRITE(numout,*) "Linkup 1.0 -- Flow out of Halo zone"
+              WRITE(numout,*) sp, sb, "Linkup 1.0 -- Flow out of Halo zone"
           ELSE
               found = 0
           END IF
@@ -2088,10 +2088,11 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                 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
+                   WRITE(numout,*) sp, sb, "flows in the same grid !"
                    found = 1
                    solved(sp,1) = solved(sp,1) + 1
                 ELSE
-                   WRITE(*,*) sp, sb, "flows in the same grid but has an error !"
+                   WRITE(numout,*) sp, sb, "flows in the same grid but has an error !"
                 ENDIF
                 !
              ELSE IF ( outflow_grid(sp,sb) .EQ. -3 ) THEN
@@ -2132,7 +2133,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                 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", sb, bop
+                   WRITE(numout,*) sp, sb, "Solution found in the original outflow_grid", sb, bop
                 ENDIF
              !
              ENDIF
@@ -2251,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. sp,sb=", sp, sb, " dop,bop=", dop, bop 
+                            WRITE (numout,*) sp, sb,"->Sol. in neighbours, output found at level:",nb,"dop,bop=",dop,bop 
                          ELSE
                             nb = nb+1
                          ENDIF
@@ -2321,7 +2322,7 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                         &                  inflow_number, inflow_grid,inflow_basin)
                    IF (outflow_basin(sp,sb) == sbint) THEN
                       found = 1
-                      WRITE (numout,*) "Lowest basin hierarchy in the grid file"
+                      WRITE (numout,*) sp, sb, "Lowest basin hierarchy in the grid file"
                    END IF
                 ENDIF
              ENDIF
@@ -2344,9 +2345,9 @@ SUBROUTINE routing_reg_linkup(nbpt, neighbours, nwbas, ijdimmax, inflowmax, basi
                             ! We go to next step and make it a coastal flow
                             IF (outflow_basin(sp, sb) .EQ. bop) THEN
                                found = 1
-                               WRITE(numout,*) "Lowest basin hierarchy in the neighbours grid points"
+                               WRITE(numout,*) sp, sb, "Lowest basin hierarchy in the neighbours grid points"
                             ELSE
-                               WRITE(numout,*) "Lowest hierarchy may be in two different grid points"
+                               WRITE(numout,*) sp, sb, "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)

Interface.py

@@ -450,15 +450,20 @@ class HydroSuper :
         fhalo = np.array([pt+1 for pt in range(self.nbpt) if pt not in part.landcorelist], order = "F")
         #
         # Outflow grid in global index and send to halo
+        hg = np.copy(self.outflow_grid)
         hg = part.l2glandindex(self.outflow_grid)
         part.landsendtohalo(hg, order='F')
         # Convert to local index
         outflows = np.unique(hg)
         outflows_out = [a for a in outflows if (a not in nbpt_glo and a>0)]
+        # Work in a copy to avoid error
+        hg1 = np.copy(hg)
         for a in outflows_out:
-          hg[hg == a] = 0
-        for i, b in enumerate(nbpt_glo):
-          hg[hg == b] = nbpt_loc[i]
+          hg1[hg == a] = 0
+        for loc, glo in zip(nbpt_loc,nbpt_glo):
+          hg1[hg == glo[0]] = loc[0]
+        hg = hg1
+        del hg1
         # Send Outflow basin to the halo and adapt it
         hb = np.copy(self.outflow_basin)
         part.landsendtohalo(hb, order='F')

Partition.py

@@ -113,8 +113,8 @@ def partitiondom(land, nbcore, rank) :
     test = np.zeros((njg,nig), dtype=np.int32)
     test[:,:] = land[:,:]
     n=0
-    for i in range(nig) :
-        for j in range(njg) :
+    for j in range(njg) :
+        for i in range(nig) :
             if (land[j,i] > 0 ) :
                 gindland[j,i] = n
                 n += 1
@@ -179,8 +179,8 @@ def landpartition(land, nbcore, rank) :
     glandcore=np.zeros((njg,nig), dtype=np.int32)
     glandcore[:,:]=-1
     n=0
-    for i in range(nig) :
-        for j in range(njg) :
+    for j in range(njg) :
+        for i in range(nig) :
             if (land[j,i] > 0 ) :
                 gindland[j,i] = n
                 cc = np.argmax(distppc > 0)