Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 36fae2a5 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

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

parents 125ee7f0 dad41b15
......@@ -9,15 +9,18 @@ FFLAG = -fPIC
FPY = f2py
#PDBG = --debug
#
all : routing_interface.so
all : routing_interface.so diagst.so
clean :
rm -f compilation.log *.o *.so *.mod routing_interface.so *~
rm -f compilation.log *.o *.so *.mod routing_interface.so diagst.so *~
routing_interface.so : $(obj) routing_interface.f90
rm -f compilation.log
$(FPY) -c routing_interface.f90 $(PDBG) --fcompiler=gfortran -m routing_interface $(obj) > compilation.log
diagst.so : diagst.f90
$(FPY) -c diagst.f90 $(PDBG) --fcompiler=gfortran -m diagst > compilation.log
defprec.o : defprec.f90
$(FC) $(FFLAG) $(FDBG) -c defprec.f90
......
MODULE diagst
CONTAINS
SUBROUTINE upstmask(nbpt, nbasmax, routetogrid, routetobasin,&
& routenbbasin,basin_area, area, g0,b0, mask)
IMPLICIT NONE
!! INPUT
INTEGER(4), INTENT(in) :: nbpt, nbasmax
!
INTEGER(4), DIMENSION(nbpt, nbasmax), INTENT(in) :: routetogrid
INTEGER(4), DIMENSION(nbpt, nbasmax), INTENT(in) :: routetobasin
INTEGER(4), DIMENSION(nbpt), INTENT(in) :: routenbbasin
!
REAL(8), DIMENSION(nbpt, nbasmax), INTENT(in) :: basin_area
REAL(8), DIMENSION(nbpt), INTENT(in) :: area
!
INTEGER(4), INTENT(in) :: g0,b0
!
!! LOCAL
INTEGER(4) :: ig,ib,jg,jb,jt,nbas
!
!! OUTPUT
REAL(8), DIMENSION(nbpt), INTENT(out) :: mask
!!
!!
mask(:) = 0
DO ig =1,nbpt
!WRITE(*,*) ig, "/",nbpt
nbas = routenbbasin(ig)
DO ib=1,nbas
jg = ig
jb = ib
DO WHILE (jb .LE. nbasmax)
! CHECK if mask > 0
! valable si j'utilise mask ig,ib
IF (((jg .EQ. g0) .AND. (jb .EQ. b0))) THEN
! add area as a variable
mask(ig) = mask(ig) + basin_area(ig,ib)/area(ig)
jb = nbasmax+1
ELSE
jt = routetogrid(jg,jb)
jb = routetobasin(jg,jb)
jg = jt
END IF
END DO
END DO
END DO
END SUBROUTINE upstmask
END MODULE
#
#
import numpy as np
import numpy.ma as ma
from netCDF4 import Dataset
from netCDF4 import stringtochar, chartostring, stringtoarr
import sys
import codecs
import os
from inspect import currentframe, getframeinfo
import sys
localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
"""
F90=localdir+'/F90subroutines'
err=os.system("cd "+F90+"; make all")
if err != 0 :
print("Compilation error in the FORTRAN interface")
sys.exit()
"""
from diagst import *
#
import getargs
config = getargs.SetupConfig()
......@@ -278,9 +293,66 @@ class PlacedLocations :
self.FMONindex = locinfo[4,:].astype(int).tolist()
self.lonrange = [np.min(nc.variables["lon"][:,:]), np.max(nc.variables["lon"][:,:])]
self.latrange = [np.min(nc.variables["lat"][:,:]), np.max(nc.variables["lat"][:,:])]
# Parameters
basid = nc.variables["basinid"][:,:,:]
self.nbasmax = basid.shape[0]
land = nc.variables["land"][:,:]
self.nbpt = int(np.sum(land))
#
self.lnglat = {varn: nc.variables[varn][:] for varn in ["lon", "lat"]}
#
# Grid points
self.nbpt_glo = nc.variables["nbpt_glo"][:]
a = np.unravel_index(ma.argsort(self.nbpt_glo, axis = None)[:self.nbpt], self.nbpt_glo.shape)
self.conv_land2ij = [[int(a[0][i]), int(a[1][i])] for i in range(a[0].shape[0])]
# Variables
self.routetogrid = self.conv_2land(nc,"routetogrid")
self.routetobasin = self.conv_2land(nc,"routetobasin")
self.routenbintobas = self.conv_2land(nc,"routenbintobas")
self.basin_area = self.conv_2land(nc,"basin_area")
self.area = self.conv_2land(nc,"area")
nc.close()
return
def mask(self, i) :
def conv_2land(self, nc, varname):
"""
Convert the variable from nhtu/lat/lon to ig/ib format.
"""
var = nc.variables[varname][:]
if len(var.shape) == 3:
return np.array([var[:,j,i] for j,i in self.conv_land2ij], order = "F")
elif len(var.shape) == 2:
return np.array([var[j,i] for j,i in self.conv_land2ij], order = "F")
#
def upstmask(self, i, j, k):
"""
Get the upstream area mask from the ig/ib coordinates.
"""
g0 = self.nbpt_glo[int(j)-1,int(i)-1] # -1 Because F -> C
b0 = int(k)
mask = diagst.upstmask(
nbpt = self.nbpt,
nbasmax = self.nbasmax,
routetogrid = self.routetogrid,
routetobasin = self.routetobasin,
routenbbasin = self.routenbintobas,
basin_area = self.basin_area,
area = self.area,
g0 = g0,
b0 = b0,
)
return mask
def mask(self, index):
# Function to compute mask of area upstream of function with information in the GraphFile
mask=np.zeros(self.gridshape)
return mask
i = self.Fiindex[index]
j = self.Fjindex[index]
k = self.FHTUindex[index]
#
mask_raw = self.upstmask(i,j,k)
mask = np.zeros(self.gridshape)
for g in range(self.nbpt):
jj,ii = self.conv_land2ij[g]
mask[jj,ii] = min(mask_raw[g],1)
return mask
......@@ -88,6 +88,16 @@ nc.createDimension('x', L.gridshape[1])
nc.createDimension('y', L.gridshape[0])
nc.setncattr("RoutingGraph_File", GraphFile)
nc.setncattr("GRDCObservation_File", GRDCFile)
# Lon and Lat variables
D = {"lon": ["degrees east", "Longitude", "X"],
"lat": ["degrees north","Latitude", "Y"]}
NCFillValue = 1.0e20
for vname in ["lon", "lat"]:
var = nc.createVariable(vname, float, ('y','x'), fill_value=NCFillValue)
var[:] = L.lnglat[vname][:]
var.units="grid box centre "+D[vname][0]
var.title=D[vname][1]
var.axis=D[vname][2]
# Create variables only if the GRDC station has been placed on the graph
ncvar=[]
for i,loc in enumerate(Lind) :
......
......@@ -23,7 +23,7 @@ numop = 100
#
# Output
#
GraphFile = AmSud_A_graph.nc
GraphFile = AmSud_A_graph_HydroSHEDS_nbasmax35.nc
#
WeightFile = Weights_AmSud_HydroSHEDS.nc
#
......@@ -34,3 +34,5 @@ MaxDistErr = 25.0
# Maximum error in the upstream area in %
MaxUpstrErr = 25.0
GRDCFile = /homedata/aschrapffer/Observations/GRDC_Monthly_Jan20.nc
FloodplainsFile = /homedata/aschrapffer/FLOOPLAINS_INPUT/HYDROSHEDS_floodplains.nc
......@@ -4,7 +4,8 @@
#
EarthRadius = 6370000.
#
ModelGridFile = /home/lfita/estudios/RegIPSL/AmSud/geo_em_CentAmSud.nc
# ModelGridFile = /home/lfita/estudios/RegIPSL/AmSud/geo_em_CentAmSud.nc
ModelGridFile = /home/lfita/sandbox/copy/geo_em.d01.nc
# WEST_EAST = -90.75, -30.25
# SOUTH_NORTH = -55.25, 25.5
HydroFile = /bdd/ORCHIDEE_Forcing/Routing/Hydro4ORCH/HydroSHEDS_qGlobal.nc
......@@ -19,7 +20,7 @@ nbasmax = 3
#
# Number of operation of simplification performed together
#
numop = 20
numop = 50
#
# Output
#
......@@ -39,4 +40,5 @@ MaxDistErr = 25.0
# Maximum error in the upstream area in %
MaxUpstrErr = 25.0
GRDCFile = /homedata/aschrapffer/Observations/GRDC_Monthly_Jan20.nc
#
FloodplainsFile = /homedata/aschrapffer/FLOOPLAINS_INPUT/HYDROSHEDS_floodplains.nc
......@@ -23,7 +23,7 @@ numop =100
#
# Output
#
GraphFile = AmSud_A_graph_nbasmax15_t.nc
GraphFile = AmSud_A_graph_nbasmax15.nc
#
WeightFile = Weights_MERIT.nc
#
......
......@@ -4,7 +4,7 @@
#
#PBS -j oe
#PBS -l nodes=2:ppn=32
#PBS -l walltime=6:00:00
#PBS -l walltime=12:00:00
#PBS -l vmem=250g,mem=250g
#
cd ${PBS_O_WORKDIR}
......
......@@ -4,7 +4,8 @@
#
EarthRadius = 6370000.
#
ModelGridFile = /home/lfita/estudios/RegIPSL/AmSud/geo_em_CentAmSud.nc
# ModelGridFile = /home/lfita/estudios/RegIPSL/AmSud/geo_em_CentAmSud.nc
ModelGridFile = /home/lfita/sandbox/copy/geo_em.d01.nc
# WEST_EAST = -90.75, -30.25
# SOUTH_NORTH = -55.25, 25.5
HydroFile = /bdd/ORCHIDEE_Forcing/Routing/Hydro4ORCH/MERIT_Global.nc
......@@ -19,7 +20,10 @@ nbasmax = 3
#
# Number of operation of simplification performed together
#
numop = 10
numop = 20
maxpercent = 40
frac_lim = 0.5
limit_step5 = 40
#
# Output
#
......@@ -27,11 +31,13 @@ GraphFile = AmSud_A_CentAmSud_graph.nc
#
WeightFile = Weights_MERIT_CentAmSud.nc
#
FloodplainsFile = /homedata/aschrapffer/FLOOPLAINS_INPUT/MERIT_floodplains.nc
#
#
# File containing infrastructures to be placed.
# Maximum error in the distance of the station in km^2
MaxDistErr = 25.0
# Maximum error in the upstream area in %
MaxUpstrErr = 25.0
GRDCFile = /homedata/aschrapffer/Observations/GRDC_Monthly_Jan20.nc
GRDCFile = /homedata/aschrapffer/Observations/GRDC_Monthly_Jan20.nc
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment