diff --git a/Diags/HTUs_test_and_corr.pbs b/Diags/Test_and_correction/HTUs_test_and_corr.pbs
similarity index 100%
rename from Diags/HTUs_test_and_corr.pbs
rename to Diags/Test_and_correction/HTUs_test_and_corr.pbs
diff --git a/Diags/Makefile b/Diags/Test_and_correction/Makefile
similarity index 100%
rename from Diags/Makefile
rename to Diags/Test_and_correction/Makefile
diff --git a/Diags/Routing_Test_and_Correction.py b/Diags/Test_and_correction/Routing_Test_and_Correction.py
similarity index 95%
rename from Diags/Routing_Test_and_Correction.py
rename to Diags/Test_and_correction/Routing_Test_and_Correction.py
index 29c527ea2d9050ece016f7befcd85cc1fddf89d8..b24cd108363212c5251f709b3f2904d25435d6ec 100644
--- a/Diags/Routing_Test_and_Correction.py
+++ b/Diags/Test_and_correction/Routing_Test_and_Correction.py
@@ -3,6 +3,7 @@
 import sys
 sys.path.append("/home/anthony/Documents/LOCAL/RoutingPP/ROUTING_Analysis_tools/Programs")
 import numpy as np
+import numpy.ma as ma
 from netCDF4 import Dataset as NetCDFFile
 from hydrosuper import check_fetch_order, check_inflows_in_outflow, check_outflow_in_inflows, test_outflow, compute_fetch, compute_inflows
 #
@@ -34,7 +35,8 @@ class routing_test:
       self.nbpt_glo = self.nc.variables["nbpt_glo"][:]
 
       self.num_pt = int(np.nanmax(self.nbpt_glo))
-      self.conv_land_ji = [[np.where(self.nbpt_glo == i)[0][0], np.where(self.nbpt_glo == i)[1][0]] for i in range(1, self.num_pt+1)]
+      a = np.unravel_index(ma.argsort(self.nbpt_glo, axis = None)[:self.num_pt], self.nbpt_glo.shape)
+      self.conv_land_ji = [[a[0][i], a[1][i]] for i in range(a[0].shape[0])]
 
       self.basin_count_land = self.var1d("routenbintobas")
       self.max_bas_count = int(np.max(self.basin_count_land))
@@ -51,17 +53,17 @@ class routing_test:
 
    def var1d(self, varn):
        var = self.nc.variables[varn][:]
-       output = np.array([var[j,i] for j,i in self.conv_land_ji])
+       output = np.array([var[j,i] for j,i in self.conv_land_ji], order = 'F')
        return output
 
    def var2d(self, varn):
        var = self.nc.variables[varn][:]
-       output = np.array([var[:self.max_bas_count,j,i] for j,i in self.conv_land_ji])
+       output = np.array([var[:self.max_bas_count,j,i] for j,i in self.conv_land_ji], order = 'F')
        return output
 
    def var3d(self, varn):
        var = self.nc.variables[varn][:]
-       output = np.zeros((self.num_pt,self.max_bas_count, self.max_innum))
+       output = np.zeros((self.num_pt,self.max_bas_count, self.max_innum), order = 'F')
        for l, ji in enumerate(self.conv_land_ji):
            j,i = ji
            output[l,:,:] = np.transpose(var[:self.max_innum,:self.max_bas_count,j,i])
@@ -103,7 +105,7 @@ class routing_test:
    ####
    def get_var_corr(self, X):
        if X.ndim == 2:
-           A = np.zeros(self.nc.variables["fetch"][:].shape)
+           A = np.zeros(self.nc.variables["fetch"][:].shape, order = 'F')
            for k,ji in enumerate(self.conv_land_ji):
                j,i = ji 
                A[:self.max_bas_count,j,i] = X[k,:]
diff --git a/Diags/hydrosuper.f90 b/Diags/Test_and_correction/hydrosuper.f90
similarity index 96%
rename from Diags/hydrosuper.f90
rename to Diags/Test_and_correction/hydrosuper.f90
index 79453d7d4bcdf6862e167fafa0f7bfdd1d420abf..6db2615d0e81cdad03ba27b686a157cd7fa1a2e5 100644
--- a/Diags/hydrosuper.f90
+++ b/Diags/Test_and_correction/hydrosuper.f90
@@ -17,7 +17,7 @@ SUBROUTINE check_fetch_order(nbpt, nbbas, fetch, outgrid, outbas, basin_count)
   DO ig=1,nbpt
     ibmx = basin_count(ig)
     DO ib=1,ibmx
-      IF (outbas(ig,ib) .LE. 35) THEN
+      IF (outbas(ig,ib) .LE. nbbas) THEN
         f1 = fetch(ig,ib)
         f2 = fetch(outgrid(ig,ib), outbas(ig,ib))
         IF (f1 .GT. f2) THEN 
@@ -54,7 +54,7 @@ SUBROUTINE check_inflows_in_outflow(nbpt, nbbas, nbin, outgrid, outbas, basin_co
     DO ib=1,ibmx
       og = outgrid(ig,ib)
       ob = outbas(ig,ib)
-      IF ( ob .LE. 35) THEN
+      IF ( ob .LE. nbbas) THEN
         iinmx = innum(og,ob)
         test = 0
         ! Look through the inflows
@@ -150,12 +150,12 @@ SUBROUTINE test_outflow(nbpt, nbbas, outgrid, outbas, basin_count)
       ind = 0
       og = outgrid(ig,ib)
       ob = outbas(ig,ib)
-      DO WHILE (ob .LE. 35)
+      DO WHILE (ob .LE. nbbas)
         it = outgrid(og,ob)
         ob = outbas(og,ob)
         og = it
         ind = ind+1
-        IF (ind .GE. 5000) THEN
+        IF (ind .GE. nbbas*nbpt) THEN
           og = 0
           WRITE(*,*) "ERROR", ig, ib
         END IF
@@ -202,7 +202,7 @@ SUBROUTINE compute_fetch(nbpt, nbbas, outgrid, outbas, basin_count, basin_area,
 
       og = outgrid(ig,ib)
       ob = outbas(ig,ib)
-      DO WHILE (ob .LE. 35)
+      DO WHILE (ob .LE. nbbas)
         fetch_new(og,ob) = fetch_new(og,ob) + area
         it = outgrid(og,ob)
         ob = outbas(og,ob)
@@ -243,7 +243,7 @@ SUBROUTINE compute_inflows(nbpt, nbbas, nbin, outgrid, outbas, basin_count, ingr
     DO ib=1,ibmx
       og = outgrid(ig,ib)
       ob = outbas(ig,ib)
-      IF ( ob .LE. 35) THEN
+      IF ( ob .LE. nbbas) THEN
         innum_new(og,ob) = innum_new(og,ob) +1
         iin = innum_new(og,ob)
         ingrid_new(og,ob,iin) = ig