From e1d15ce4ada41c745fab8878e2a5cb0b23e45f1a Mon Sep 17 00:00:00 2001
From: Anthony <anthony.schrapffer@polytechnique.fr>
Date: Tue, 5 May 2020 16:03:05 +0200
Subject: [PATCH] Subroutine to check the inflows and light correction

---
 F90subroutines/routing_interface.f90 | 54 ++++++++++++++++++++++++++++
 Interface.py                         | 10 ++++--
 2 files changed, 61 insertions(+), 3 deletions(-)

diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 4d4bd97..e741f22 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -715,6 +715,60 @@ SUBROUTINE checkrouting(nbpt, nwbas, outflow_grid, outflow_basin, basin_count)
 END SUBROUTINE checkrouting
 
 
+SUBROUTINE check_inflows(nbpt, nwbas, inflowmax, outflow_grid, outflow_basin, basin_count, inflow_number, inflow_grid, inflow_basin)
+
+  !
+  USE ioipsl
+  USE grid
+  USE routing_tools
+  USE routing_reg
+  !
+  !! INPUT VARIABLES
+  INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
+  INTEGER(i_std), INTENT (in) :: nwbas !!
+  INTEGER(i_std), INTENT (in) :: inflowmax
+
+  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in)       :: outflow_grid !! Type of outflow on the grid box (unitless)
+  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in)       :: outflow_basin !!
+  INTEGER(i_std), DIMENSION(nbpt),       INTENT(in)       :: basin_count !!
+  INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in)          :: inflow_number
+  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(in) :: inflow_basin
+  INTEGER(i_std), DIMENSION(nbpt,nwbas,inflowmax), INTENT(in) :: inflow_grid
+
+  ! LOCAL
+  INTEGER(i_std) :: ig, nbas, ib, og, ob, inf, found
+
+  WRITE(numout,*) "Checking if the HTUs are in the inflows of their outflow"  
+
+  DO ig=1,nbpt
+    nbas = basin_count(ig)
+    DO ib=1,nbas
+      og = outflow_grid(ig,ib)
+      ob = outflow_basin(ig,ib)
+      if (og .GT. 0) THEN
+        IF (inflow_number(og,ob) .EQ. 0) THEN
+          WRITE(numout,*) ig, ib, "Error : outflow has no inflow"
+          WRITE(numout,*) og, ob
+        ELSE
+          found = 0
+          DO inf = 1,inflow_number(og,ob)
+            IF ((inflow_grid(og,ob,inf) .EQ. ig) .AND. (inflow_basin(og,ob,inf) .EQ. ib)) THEN
+              found = 1
+            END IF
+          END DO
+          IF (found .EQ. 0) THEN
+            WRITE(numout,*) ig,ib, "Error, not found the inflows of its outflow"
+            WRITE(numout,*) og, ob
+          END IF
+        END IF
+      END IF 
+    END DO
+  END DO
+
+
+
+END SUBROUTINE check_inflows
+
 SUBROUTINE checkfetch(nbpt, nwbas, fetch_basin, outflow_grid, outflow_basin, basin_count)
   !
   USE ioipsl
diff --git a/Interface.py b/Interface.py
index 6055acb..f8f0e67 100644
--- a/Interface.py
+++ b/Interface.py
@@ -227,7 +227,7 @@ class HydroOverlap :
         #
         # Reshape stuff so that it fits into arrays
         #
-        sub_index = np.zeros((nbpt,nbvmax,2), dtype=np.int8, order='F')
+        sub_index = np.zeros((nbpt,nbvmax,2), dtype=np.int32, order='F')
         sub_area = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
         sub_lon = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
         sub_lat = np.zeros((nbpt,nbvmax), dtype=np.float32, order='F')
@@ -265,8 +265,8 @@ class HydroOverlap :
         # Compute nbxmax
         #
         ijdim=[]
-        for ib in range(nbpt) :
-            ijdim.append(max(max(sub_index[ib,:,0])-min(sub_index[ib,:,0])+1,max(sub_index[ib,:,1])-min(sub_index[ib,:,1])+1))
+        for ib in range(nbpt) :   
+            ijdim.append(max(np.max(sub_index[ib,:,0])-np.min(sub_index[ib,:,0])+1,np.max(sub_index[ib,:,1])-np.min(sub_index[ib,:,1])+1))
         ijdimmax = max(ijdim)
         #
         # Go to the call of the FORTRAN interface
@@ -429,6 +429,10 @@ class HydroSuper :
                                        basin_count = self.basin_count)
 
         return
+
+    def check_inflows(self):
+        nbxmax_tmp = self.inflow_grid.shape[2]
+        routing_interface.check_inflows(nbpt = self.nbpt, nwbas = self.nwbas, inflowmax = nbxmax_tmp, outflow_grid = self.outflow_grid, outflow_basin = self.outflow_basin, basin_count = self.basin_count, inflow_number = self.inflow_number, inflow_grid = self.inflow_grid, inflow_basin = self.inflow_basin)
     #
     #
     def killbas(self, tokill, totakeover, numops):
-- 
GitLab