Skip to content

Commit 551ef40

Browse files
StefaniaGrimaldiStefaniaGrimaldi
committed
initializion of third soil layer
1 parent 06ce914 commit 551ef40

File tree

241 files changed

+34351
-39086
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

241 files changed

+34351
-39086
lines changed

‎src/lisflood/global_modules/decorators.py

+1-1
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,7 @@ def __call__(self, *args, **kwargs):
7878
return_data = self.cache[key]
7979
else:
8080
return_data = self.cache[key]
81-
self.found[self.name] += 1
81+
self.found[self.name] += 1
8282
return return_data
8383

8484
@classmethod

‎src/lisflood/global_modules/default_options.py

+27-12
Original file line numberDiff line numberDiff line change
@@ -848,6 +848,21 @@
848848
end=[], steps=[], all=['repSeepSubToGWMaps'],
849849
restrictoption=['nonInit'], monthly=False,
850850
yearly=False),
851+
'SeepTopToSubBAverageOtherMap': ReportedMap(name='SeepTopToSubBAverageOtherMap',
852+
output_var='SeepTopToSubBAv[0]', unit='mm/timestep',
853+
end=['InitLisflood'], steps=[], all=[],
854+
restrictoption=[], monthly=False,
855+
yearly=False),
856+
'SeepTopToSubBAverageForestMap': ReportedMap(name='SeepTopToSubBAverageForestMap',
857+
output_var='SeepTopToSubBAv[1]', unit='mm/timestep',
858+
end=['InitLisflood'], steps=[], all=[],
859+
restrictoption=[], monthly=False,
860+
yearly=False),
861+
'SeepTopToSubBAverageIrrigationMap': ReportedMap(name='SeepTopToSubBAverageIrrigationMap',
862+
output_var='SeepTopToSubBAv[2]', unit='mm/timestep',
863+
end=['InitLisflood'], steps=[], all=[],
864+
restrictoption=[], monthly=False,
865+
yearly=False),
851866
'SnowCoverAEnd': ReportedMap(name='SnowCoverAEnd', output_var='SnowCoverS[0]',
852867
unit='mm', end=['repEndMaps'], steps=[],
853868
all=[], restrictoption=[], monthly=False,
@@ -911,10 +926,10 @@
911926
end=[], steps=[], all=['repTavgMaps'],
912927
restrictoption=[], monthly=False, yearly=False),
913928
'Theta1End': ReportedMap(name='Theta1End', output_var='Theta1a[0]', unit='',
914-
end=['repEndMaps'], steps=[], all=[],
929+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
915930
restrictoption=[], monthly=False, yearly=False),
916931
'Theta1ForestEnd': ReportedMap(name='Theta1ForestEnd', output_var='Theta1a[1]',
917-
unit='', end=['repEndMaps'], steps=[],
932+
unit='', end=['repEndMaps','InitLisflood'], steps=[],
918933
all=[], restrictoption=[],
919934
monthly=False, yearly=False),
920935
'Theta1ForestMaps': ReportedMap(name='Theta1ForestMaps', output_var='Theta1a[1]',
@@ -929,7 +944,7 @@
929944
yearly=False),
930945
'Theta1IrrigationEnd': ReportedMap(name='Theta1IrrigationEnd',
931946
output_var='Theta1a[2]', unit='',
932-
end=['repEndMaps'], steps=[], all=[],
947+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
933948
restrictoption=[], monthly=False,
934949
yearly=False),
935950
'Theta1IrrigationMaps': ReportedMap(name='Theta1IrrigationMaps',
@@ -952,10 +967,10 @@
952967
all=[], restrictoption=['nonInit'], monthly=False,
953968
yearly=False),
954969
'Theta2End': ReportedMap(name='Theta2End', output_var='Theta1b[0]', unit='',
955-
end=['repEndMaps'], steps=[], all=[],
970+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
956971
restrictoption=[], monthly=False, yearly=False),
957972
'Theta2ForestEnd': ReportedMap(name='Theta2ForestEnd', output_var='Theta1b[1]',
958-
unit='', end=['repEndMaps'], steps=[],
973+
unit='', end=['repEndMaps','InitLisflood'], steps=[],
959974
all=[], restrictoption=[],
960975
monthly=False, yearly=False),
961976
'Theta2ForestMaps': ReportedMap(name='Theta2ForestMaps', output_var='Theta1b[1]',
@@ -970,7 +985,7 @@
970985
yearly=False),
971986
'Theta2IrrigationEnd': ReportedMap(name='Theta2IrrigationEnd',
972987
output_var='Theta1b[2]', unit='',
973-
end=['repEndMaps'], steps=[], all=[],
988+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
974989
restrictoption=[], monthly=False,
975990
yearly=False),
976991
'Theta2IrrigationMaps': ReportedMap(name='Theta2IrrigationMaps',
@@ -993,10 +1008,10 @@
9931008
all=[], restrictoption=['nonInit'], monthly=False,
9941009
yearly=False),
9951010
'Theta3End': ReportedMap(name='Theta3End', output_var='Theta2[0]', unit='',
996-
end=['repEndMaps'], steps=[], all=[],
1011+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
9971012
restrictoption=[], monthly=False, yearly=False),
9981013
'Theta3ForestEnd': ReportedMap(name='Theta3ForestEnd', output_var='Theta2[1]',
999-
unit='', end=['repEndMaps'], steps=[],
1014+
unit='', end=['repEndMaps','InitLisflood'], steps=[],
10001015
all=[], restrictoption=[],
10011016
monthly=False, yearly=False),
10021017
'Theta3ForestMaps': ReportedMap(name='Theta3ForestMaps', output_var='Theta2[1]',
@@ -1011,7 +1026,7 @@
10111026
yearly=False),
10121027
'Theta3IrrigationEnd': ReportedMap(name='Theta3IrrigationEnd',
10131028
output_var='Theta2[2]', unit='',
1014-
end=['repEndMaps'], steps=[], all=[],
1029+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
10151030
restrictoption=[], monthly=False,
10161031
yearly=False),
10171032
'Theta3IrrigationMaps': ReportedMap(name='Theta3IrrigationMaps',
@@ -1094,10 +1109,10 @@
10941109
all=[], restrictoption=['nonInit', 'wateruse'],
10951110
monthly=True, yearly=False),
10961111
'UZEnd': ReportedMap(name='UZEnd', output_var='UZ[0]', unit='mm',
1097-
end=['repEndMaps'], steps=[], all=[],
1112+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
10981113
restrictoption=[], monthly=False, yearly=False),
10991114
'UZForestEnd': ReportedMap(name='UZForestEnd', output_var='UZ[1]', unit='mm',
1100-
end=['repEndMaps'], steps=[], all=[],
1115+
end=['repEndMaps','InitLisflood'], steps=[], all=[],
11011116
restrictoption=[], monthly=False,
11021117
yearly=False),
11031118
'UZForestMaps': ReportedMap(name='UZForestMaps', output_var='UZ[1]', unit='mm',
@@ -1109,7 +1124,7 @@
11091124
all=[], restrictoption=['nonInit'], monthly=False,
11101125
yearly=False),
11111126
'UZIrrigationEnd': ReportedMap(name='UZIrrigationEnd', output_var='UZ[2]',
1112-
unit='mm', end=['repEndMaps'], steps=[],
1127+
unit='mm', end=['repEndMaps','InitLisflood'], steps=[],
11131128
all=[], restrictoption=[],
11141129
monthly=False, yearly=False),
11151130
'UZIrrigationMaps': ReportedMap(name='UZIrrigationMaps', output_var='UZ[2]',

‎src/lisflood/hydrological_modules/soil.py

+42-4
Original file line numberDiff line numberDiff line change
@@ -25,7 +25,13 @@
2525
from collections import OrderedDict
2626
#from global_modules.add1 import *
2727
import xarray as xr
28+
import numpy as np ##### 2024
29+
import scipy ##### 2024
30+
from scipy.optimize import least_squares ##### 2024
2831

32+
def function_estimate(satdegrthetastart,*data): ##### 2024
33+
SeepTopToSubBPixelAv, KSat2,GenuInvM2, GenuM2 = data
34+
return SeepTopToSubBPixelAv - KSat2 * np.sqrt(satdegrthetastart) * (1. - (1. - satdegrthetastart ** GenuInvM2) ** GenuM2) ** 2
2935

3036
def pressure2SoilMoistureFun(residual_sm, sat_sm, GenuA, GenuN, GenuM):
3137
"""Generate a function to compute soil moisture values corresponding to characteristic pressure head levels [cm].
@@ -71,6 +77,10 @@ def __init__(self, soil_variable):
7177
def initial(self):
7278
""" initial part of the soil module
7379
"""
80+
maskinfo = MaskInfo.instance()
81+
self.var.cumSeepTopToSubBAv = self.var.allocateVariableAllVegetation() #######2024
82+
self.var.SeepTopToSubBAv = self.var.allocateVariableAllVegetation() #######2024
83+
7484
def splitlanduse(array1, array2=None, array3=None):
7585
""" splits maps into the 3 different land use types - other , forest, irrigation
7686
"""
@@ -88,7 +98,6 @@ def splitlanduse(array1, array2=None, array3=None):
8898
settings = LisSettings.instance()
8999
option = settings.options
90100
binding = settings.binding
91-
maskinfo = MaskInfo.instance()
92101
if not option["cropsEPIC"]: # If EPIC is active, the rice fraction initialisation is handled by EPIC (setSoilFractions in EPIC_main.py)
93102
self.var.SoilFraction.values[self.var.vegetation.index('Rainfed_prescribed')] += self.var.RiceFraction
94103

@@ -264,18 +273,39 @@ def splitlanduse(array1, array2=None, array3=None):
264273
# Set to zero if soil depth is zero.
265274
# IMPORTANT: WInit1 and WInit2 represent the soil moisture in the *permeable* fraction of the pixel *only*
266275
# (since all soil moisture-related calculations are done for permeable fraction only!).
276+
if not option['InitLisflood']:
277+
self.var.SeepTopToSubBAv[0] = loadmap('SeepTopToSubBAverageOtherMap')
278+
self.var.SeepTopToSubBAv[1] = loadmap('SeepTopToSubBAverageForestMap')
279+
self.var.SeepTopToSubBAv[2] = loadmap('SeepTopToSubBAverageIrrigationMap')
267280

268281
for veg, luse in self.var.VEGETATION_LANDUSE.items():
269282
iveg = self.var.vegetation.index(veg)
270283
iluse = self.var.SOIL_USES.index(luse)
271284
ini_1a = np.where(ThetaInit1aValue[iveg] == -9999, self.var.WFC1a[iluse], ThetaInit1aValue[iveg] * self.var.SoilDepth1a[iluse])
272285
self.var.W1a[iveg] = np.where(self.var.PoreSpaceNotZero1a[iluse], ini_1a, 0)
273286
ini_1b = np.where(ThetaInit1bValue[iveg] == -9999, self.var.WFC1b[iluse], ThetaInit1bValue[iveg] * self.var.SoilDepth1b[iluse])
274-
self.var.W1b[iveg] = np.where(self.var.PoreSpaceNotZero1b[iluse], ini_1b, 0)
287+
self.var.W1b[iveg] = np.where(self.var.PoreSpaceNotZero1b[iluse], ini_1b, 0)
288+
275289
ini_2 = np.where(ThetaInit2Value[iveg] == -9999, self.var.WFC2[iluse], ThetaInit2Value[iveg] * self.var.SoilDepth2[iluse])
276290
self.var.W2[iveg] = np.where(self.var.PoreSpaceNotZero2[iluse], ini_2, 0)
277-
self.var.W1 = self.var.W1a + self.var.W1b
278-
291+
292+
293+
if not option['InitLisflood']: ## 2024
294+
295+
flag_coldstart = np.min( self.var.SeepTopToSubBAv[0] + self.var.SeepTopToSubBAv[1] + self.var.SeepTopToSubBAv[2])
296+
297+
if flag_coldstart > -9999.0:
298+
SOLVED = ((ThetaInit2Value[iveg] * self.var.SoilDepth2[iluse])-self.var.WRes2[iluse])/(self.var.WS2[iluse]-self.var.WRes2[iluse]) #### for extra safety - stef
299+
for ii in np.arange(len(ThetaInit2Value[iveg])):
300+
data = []
301+
analyticalcheckzero = []
302+
data=(self.var.SeepTopToSubBAv[iluse][ii], self.var.KSat2[iluse][ii], self.var.GenuInvM2[iluse][ii], self.var.GenuM2[iluse][ii])
303+
analyticalcheckzero = least_squares(function_estimate,((ThetaInit2Value[iveg][ii] * self.var.SoilDepth2[iluse][ii])-self.var.WRes2[iluse][ii])/(self.var.WS2[iluse][ii]-self.var.WRes2[iluse][ii]),bounds=(self.var.WFC2[iluse][ii]*0,self.var.WFC2[iluse][ii]*0+1.0), args=data)
304+
SOLVED[ii]=analyticalcheckzero.x
305+
ini_2 = SOLVED*(self.var.WS2[iluse]-self.var.WRes2[iluse])+self.var.WRes2[iluse]
306+
self.var.W2[iveg] = np.where(self.var.PoreSpaceNotZero2[iluse], ini_2, 0)
307+
308+
self.var.W1 = self.var.W1a + self.var.W1b
279309

280310
self.var.Sat1a = self.var.allocateVariableAllVegetation()
281311
self.var.Sat1b = self.var.allocateVariableAllVegetation()
@@ -508,6 +538,14 @@ def dynamic_perpixel(self):
508538
# Pixel-average seepage values in [mm] per timestep
509539
# (no seepage from direct runoff fraction)
510540

541+
self.var.cumSeepTopToSubBAv[0] += self.var.SeepTopToSubB[0] #### 2024
542+
self.var.SeepTopToSubBAv[0] = (self.var.cumSeepTopToSubBAv[0] * self.var.InvDtDay) / (self.var.TimeSinceStart) #### 2024
543+
self.var.cumSeepTopToSubBAv[1] += self.var.SeepTopToSubB[1] #### 2024
544+
self.var.SeepTopToSubBAv[1] = (self.var.cumSeepTopToSubBAv[1] * self.var.InvDtDay) / (self.var.TimeSinceStart) #### 2024
545+
self.var.cumSeepTopToSubBAv[2] += self.var.SeepTopToSubB[2] #### 2024
546+
self.var.SeepTopToSubBAv[2] = (self.var.cumSeepTopToSubBAv[2] * self.var.InvDtDay) / (self.var.TimeSinceStart) #### 2024
547+
548+
511549
# the variables below were added to report catchment-averaged soil moisture profiles
512550
self.var.Theta1aPixel = self.var.deffraction(self.var.Theta1a)
513551
self.var.Theta1bPixel = self.var.deffraction(self.var.Theta1b)

‎src/lisfloodSettings_reference.xml

+18
Original file line numberDiff line numberDiff line change
@@ -3667,6 +3667,24 @@ Reported seepage to groundwater(from soil layer 2 to groundwater), weighted sum
36673667
</comment>
36683668
</textvar>
36693669

3670+
<textvar name="SeepTopToSubBAverageOtherMap" value= "$(PathOut)/SeepTopToSubBAverageOtherMap">
3671+
<comment>
3672+
Reported infiltration from the soil layer 2 to soil layer 3, other land cover fraction, average flux over the simulation period [mm]
3673+
</comment>
3674+
</textvar>
3675+
3676+
<textvar name="SeepTopToSubBAverageForestMap" value= "$(PathOut)/SeepTopToSubBAverageForestMap">
3677+
<comment>
3678+
Reported infiltration from the soil layer 2 to soil layer 3, forest land cover fraction, average flux over the simulation period [mm]
3679+
</comment>
3680+
</textvar>
3681+
3682+
<textvar name="SeepTopToSubBAverageIrrigationMap" value= "$(PathOut)/SeepTopToSubBAverageIrrigationMap">
3683+
<comment>
3684+
Reported infiltration from the soil layer 2 to soil layer 3, irrigated land cover fraction, average flux over the simulation period [mm]
3685+
</comment>
3686+
</textvar>
3687+
36703688
<textvar name="AvailGWMaps" value="$(PathOut)/agw">
36713689
<comment>
36723690
Reported available groundwater LZ+ GWLoss [mm]

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/SeepTopToSubBAverageForestMap.nc

36.9 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/SeepTopToSubBAverageIrrigationMap.nc

36.9 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/SeepTopToSubBAverageOtherMap.nc

37 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/avgdis.nc

45 Bytes
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/lzavin.nc

45 Bytes
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/tha.end.nc

35.7 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thb.end.nc

36.1 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thc.end.nc

35.5 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thfa.end.nc

35.5 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thfb.end.nc

36 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thfc.end.nc

35.5 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thia.end.nc

35.5 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thib.end.nc

35.9 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/thic.end.nc

35.5 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/uz.end.nc

20.9 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/uzf.end.nc

20.8 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init/uzi.end.nc

21.3 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init_LF431/avgdis.nc

37.5 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/maps/safe_init_LF431/lzavin.nc

34.6 KB
Binary file not shown.

‎tests/data/LF_ETRS89_UseCase/reference/init_daily/chanqWin.tss ‎tests/data/LF_ETRS89_UseCase/out/chanqWin.tss

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
timeseries valuescale.scalar settingsfile: /BGFS/DISASTER/grimast/LF411/lisflood-code/tests/data/LF_ETRS89_UseCase/settings/prerun_testinit_createreference.xml date: Thu Jul 28 14:37:27 2022
1+
timeseries valuescale.scalar settingsfile: /BGFS/DISASTER/grimast/LF_DEVELOP/lisflood-code/tests/data/LF_ETRS89_UseCase/settings/prerun.xml_1c02bbb9-e5d4-4134-a294-99f6feae7078.xml date: Fri May 31 10:42:40 2024
22
31
33
timestep
44
334

‎tests/data/LF_ETRS89_UseCase/reference/init_daily/dis.tss ‎tests/data/LF_ETRS89_UseCase/out/dis.tss

+1-1
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
timeseries valuescale.scalar settingsfile: /BGFS/DISASTER/grimast/LF411/lisflood-code/tests/data/LF_ETRS89_UseCase/settings/prerun_testinit_createreference.xml date: Thu Jul 28 14:37:27 2022
1+
timeseries valuescale.scalar settingsfile: /BGFS/DISASTER/grimast/LF_DEVELOP/lisflood-code/tests/data/LF_ETRS89_UseCase/settings/prerun.xml_1c02bbb9-e5d4-4134-a294-99f6feae7078.xml date: Fri May 31 10:42:40 2024
22
31
33
timestep
44
334

‎tests/data/LF_ETRS89_UseCase/out/hLake.tss

+10
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
timeseries valuescale.scalar settingsfile: /BGFS/DISASTER/grimast/LF_DEVELOP/lisflood-code/tests/data/LF_ETRS89_UseCase/settings/full.xml_d06a5d5e-4c30-4bfa-858c-26d41e5297a4.xml date: Fri May 31 10:44:55 2024
2+
6
3+
timestep
4+
641
5+
643
6+
645
7+
646
8+
647
9+
6057 0.505187 1.11651 0.888341 0.650731 0.523296
10+
6058 0.522721 1.09129 0.876231 0.621849 0.517775

‎tests/data/LF_ETRS89_UseCase/out/qLakeIn.tss

+10
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
timeseries valuescale.scalar settingsfile: /BGFS/DISASTER/grimast/LF_DEVELOP/lisflood-code/tests/data/LF_ETRS89_UseCase/settings/full.xml_d06a5d5e-4c30-4bfa-858c-26d41e5297a4.xml date: Fri May 31 10:44:55 2024
2+
6
3+
timestep
4+
641
5+
643
6+
645
7+
646
8+
647
9+
6057 8.22618 41.5103 160.367 7.42939 5.52424
10+
6058 29.3635 58.4609 178.556 5.06521 1.19854

‎tests/data/LF_ETRS89_UseCase/out/qLakeOut.tss

+10
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
timeseries valuescale.scalar settingsfile: /BGFS/DISASTER/grimast/LF_DEVELOP/lisflood-code/tests/data/LF_ETRS89_UseCase/settings/full.xml_d06a5d5e-4c30-4bfa-858c-26d41e5297a4.xml date: Fri May 31 10:44:55 2024
2+
6
3+
timestep
4+
641
5+
643
6+
645
7+
646
8+
647
9+
6057 20.3273 105.669 214.421 26.3394 22.5795
10+
6058 19.566 99.1421 203.524 24.0527 22.127

0 commit comments

Comments
 (0)