Merge gitlab.in2p3.fr:ipsl/lmd/intro/routingpp

F90subroutines/routing_interface.f90

@@ -885,15 +885,13 @@ SUBROUTINE get_floodcri(nbpt, nwbas,inflowmax, inflow_number, inflow_basin, infl
   INTEGER :: ig, ib, ug, ub, i_inf, nbas
   REAL    :: diff, orog, d
 
-  floodcri(:,:) = 0
+  floodcri(:,:) = 2
 
   DO ig=1,nbpt
     nbas = basin_count(ig)
     DO ib=1,nbas
       IF (basin_floodp(ig,ib) .GT. 0) THEN
-        IF (inflow_number(ig,ib) .EQ. 0) THEN
-          diff = 2
-        ELSE
+        IF (inflow_number(ig,ib) .GE. 0) THEN
           diff = 999
           orog = basin_orog_min(ig,ib)
           DO i_inf=1,inflow_number(ig,ib)

F90subroutines/routing_reg.f90

@@ -1670,7 +1670,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
        basin_orog_max(ib,ij) = zero
        basin_floodp(ib,ij) = zero
        basin_orog_std(ib,ij) = 0
-       basin_beta_fp(ib,ij) = 0
+       basin_beta_fp(ib,ij) = 2
        !
        SELECT CASE (hierar_method)
           !
@@ -1749,7 +1749,7 @@ SUBROUTINE routing_reg_globalize(nbpt, nb_htu, nbv, ib, ijdimmax, neighbours, ar
          DO iz=1,basin_sz(ij)
            temp =  temp + (orog_bx(basin_pts(ij,iz,1),basin_pts(ij,iz,2)) -basin_orog_mean(ib,ij))**2
          END DO
-         IF (basin_sz(ij) .GT. 0) basin_orog_std(ib,ij) = SQRT(temp)/(basin_sz(ij))
+         basin_orog_std(ib,ij) = SQRT(temp)/(basin_sz(ij))
          IF (basin_orog_std(ib,ij) .LT. lowlim_std) basin_orog_std(ib,ij) = lowlim_std
          IF (basin_orog_std(ib,ij) .GT. uplim_std) basin_orog_std(ib,ij) = uplim_std
        END IF
@@ -2768,7 +2768,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
           routing_orog_min_glo(ib,ij) = basin_orog_min(ib,ij)
           routing_orog_max_glo(ib,ij) = basin_orog_max(ib,ij)
           ! Now the final floodplains is converted into fraction
-          routing_floodp_glo(ib,ij) = basin_floodp(ib,ij)/basin_area(ib,ij)
+          routing_floodp_glo(ib,ij) = MIN(basin_floodp(ib,ij)/basin_area(ib,ij),1.)
           routing_beta_glo(ib,ij) = basin_beta_fp(ib,ij)
           routing_cg_glo(ib,ij,:) = basin_cg(ib,ij,:)
           !
@@ -2782,7 +2782,7 @@ SUBROUTINE routing_reg_end_truncate(nbpt, nbasmax, gridarea, contfrac, gridcente
           !
           route_togrid_glo(ib,ij) = outflow_grid(ib,ij)
           route_tobasin_glo(ib,ij) = outflow_basin(ib,ij)
-          routing_floodcri_glo(ib,ij) = floodcri(ib,ij)*2000
+          routing_floodcri_glo(ib,ij) = floodcri(ib,ij)*1000
           !
        ENDDO
     ENDDO
@@ -2946,6 +2946,14 @@ SUBROUTINE routing_reg_killbas(nbpt, ib, tokill, totakeover, nwbas, inflowmax, b
     !basin_orog_min(ib, totakeover) =
     basin_orog_max(ib, totakeover) = MAX(basin_orog_mean(ib, totakeover), basin_orog_mean(ib, tokill))
     !
+    IF ((basin_beta_fp(ib, totakeover)>0) .AND. (basin_beta_fp(ib,tokill)> 0)) THEN
+        basin_beta_fp(ib, totakeover) = (basin_beta_fp(ib,totakeover)*basin_area(ib, totakeover) &
+            & + basin_beta_fp(ib,tokill) * basin_area(ib, tokill))/(basin_area(ib,totakeover) + basin_area(ib, tokill))
+    ELSE IF (basin_beta_fp(ib,tokill)> 0) THEN
+        ! IF beta_fp from totakeover is 0 we keep beta_fp from tokill
+        basin_beta_fp(ib, totakeover) = basin_beta_fp(ib,tokill) 
+    END IF
+    !
     ! Add the fetch of the basin will kill to the one which gets the water
     !
     fetch_basin(ib, totakeover) = fetch_basin(ib, totakeover) + fetch_basin(ib, tokill)

Locations.py

@@ -156,7 +156,7 @@ class Locations :
     #
     def addtonetcdf(self, nbpt, nbasmax, outnf, procgrid, part, spacecoord, loccoord, vtyp) :
         #
-        locname = "locations"
+        locname = "st_locations"
         loctitle = "Basic information of located stations or structures"
         locunits = "-"
         monname = "HTUmonitor"

Tools/HydroDATA_shapefile/lib/GRID.py

@@ -32,11 +32,11 @@ class GlobalGrid:
    def __init__(self, fdir, input_type):
       # Get the coordinates
       nc = Dataset(fdir, "r")
-      griddesc, self.lon, self.lat, self.res_lon, self.res_lat = MG.getcoordinates(nc, 0, -1, 0,  -1)
-      self.nj,self.ni = self.lon.shape
       #
-      self.land = np.full(self.lon.shape, 1); self.nbland = np.sum(self.land)
       if input_type == "geogrid":
+          griddesc, self.lon, self.lat, self.res_lon, self.res_lat = MG.getcoordinates(nc, 0, -1, 0,  -1)
+          self.nj,self.ni = self.lon.shape
+          
           # Get the projection
           if griddesc['type'] == "lambert_conformal_conic"  :
               self.proj = PS.LambertC(griddesc['dx'], griddesc['known_lon'], griddesc['known_lat'], griddesc['truelat1'], \
@@ -49,7 +49,12 @@ class GlobalGrid:
           else :
               ERROR("Unknown grid type")
               sys.exit()
-      elif input_grid == "hydro":
+      elif input_type == "hydro":
+         self.lon=np.copy(nc.variables["nav_lon"][:,:])
+         self.lat=np.copy(nc.variables["nav_lat"][:,:])
+         self.nj,self.ni = self.lon.shape
+         self.res_lon = np.mean(np.abs(np.diff(nc.variables["nav_lon"][0,:])))
+         self.res_lat = np.mean(np.abs(np.diff(nc.variables["nav_lat"][:,0])))
          try:
              self.proj=PS.RegLonLat(self.res_lon, self.res_lat, self.lon[0], self.lat[0]) 
          except:

Tools/HydroDATA_shapefile/main_multiprocess.py

@@ -49,7 +49,7 @@ if multiprocessing.cpu_count() < nbcore:
 
 ##############
 
-gg = GlobalGrid(dinput, "geogrid")
+gg = GlobalGrid(dinput, input_type)
 
 landP= landPolygon(shpdir)
 

Tools/HydroDATA_shapefile/run.def

@@ -16,12 +16,12 @@ Debug = False
 
 # [DEF] input_type : geogrid, hydro
 #
-input_type = 
-dinput = 
+input_type = hydro
+dinput = /bdd/ORCHIDEE_Forcing/Routing/Hydro4ORCH/MERIT_Global.nc 
 
 # output file direction
 doutput = ./output.nc
 
 # Shapefile direction
-shpdir = 
+shpdir = ./shapefile/ne_10m_land/ne_10m_land.shp