Improving and organizing diagnostic tools

Diags/Routing_Test_and_Correction.py → 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,:]

Diags/hydrosuper.f90 → 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