diff --git a/VERSION b/VERSION
index ecedc98d..84ed6d2a 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-4.3.1
\ No newline at end of file
+4.3.1 MCT
\ No newline at end of file
diff --git a/docs/3_step3_preparing-setting-file/index.md b/docs/3_step3_preparing-setting-file/index.md
index d9a3ec9d..43c63a8a 100644
--- a/docs/3_step3_preparing-setting-file/index.md
+++ b/docs/3_step3_preparing-setting-file/index.md
@@ -378,6 +378,41 @@ These parameters are all related to the [routing of water in the channels](https
+### Diffusive wave routing parameters
+
+The following parameters are related to the [diffusive wave routing](https://ec-jrc.github.io/lisflood-model/3_14_optLISFLOOD_diffusive-wave/) in river channels. The multiplier *CalChanMan3* can be used to fine-tune the diffusive wave propagation when using the Muskingum-Cunge-Todini (MCT) routing, and it can be defined as either a single value or a map. The map *ChannelsMCT* is a Bolean map with the mask of rivers where MCT wave routing must be used. The parameter *ChanGradMaxMCT* defines the maximum riverbed slope for river grid cells using the MCT wave routing. The parameter is provided as a single number and it is recommended to set it to values < 0.001 and > *ChanGradMin*
+```xml
+
+ **************************************************************
+ MUSKINGUM-CUNGE-TODINI ROUTING PARAMETERS
+ **************************************************************
+
+
+
+ default: 3.0 [-]
+ Multiplier [-] applied to Channel Manning's n for MCT routing
+
+
+
+
+ Boolean map with value 1 at channel pixels where MCT is
+ used, and 0 at all other pixels
+
+
+
+
+ Maximum channel gradient for channels using MCT routing [-] (for MCT wave: slope cannot be 0)
+ Default: 0.001
+
+
+```
+
+- **CalChanMan3** is a multiplier that is applied to the Manning’s roughness map of the [channel system](https://ec-jrc.github.io/lisflood-model/2_16_stdLISFLOOD_channel-routing/) [-] for the grid cells where MCT routing is used
+- **ChannelsMCT** is Bolean mask including the rivers grid cells using the MCT wave routing [-]
+- **ChanGradMaxMCT** is a upper limit for the channel gradient used in the calculation of the MCT wave routing [m/m]
+
+
+
### Parameters related to numerics
This category only contains one parameter at the moment, which can only be a single value. We strongly recommend keeping this parameter at its default value.
@@ -530,8 +565,8 @@ Note that a couple of variables can be [initialized internally](https://ec-jrc.g
-**************************************************************
- The following variables can also be initialized in the model internally. if you want this to happen set them to bogus value of -9999
+ **************************************************************
+ The following variables can also be initialized in the model internally. If you want this to happen set them to bogus value of -9999
**************************************************************
@@ -571,7 +606,19 @@ Note that a couple of variables can be [initialized internally](https://ec-jrc.g
only needed for lakes reservoirs and transmission loss
-9999: use discharge of half bankfull
-
+
+
+
+ Courant number at previous step for MCT routing
+ Cold start: -9999: use 1
+
+
+
+
+ Reynolds number at previous step for MCT routing
+ Cold start: -9999: use 0
+
+
```
- **WaterDepthInitValue** is the initial amount of water on the soil surface $[mm]$
@@ -604,6 +651,10 @@ Note that a couple of variables can be [initialized internally](https://ec-jrc.g
- **PrevDischarge** is the initial discharge from previous run $[\frac{m^3} {s}]$ used for lakes, reservoirs and transmission loss (only needed if option is on for lakes or reservoirs or transmission loss). Note that PrevDischarge is discharge as an average over the time step (a flux) . A value of **-9999** sets the initial amount of discharge to equivalent of half bankfull.
+- **PrevCmMCTInitValue** is the Courant number at the end of the previous step and it is only used for MCT wave routing [-]. A value of **-9999 ** sets the initial value to 1.
+
+- **PrevDmMCTInitValue** is the Reynols number at the end of the previous step and it is only used for MCT wave routing [-]. A value of **-9999 ** sets the initial value to 0.
+
```xml
**************************************************************
diff --git a/docs/3_step4_preparing-input-files/index.md b/docs/3_step4_preparing-input-files/index.md
index cc7cf9bf..b0d545ad 100644
--- a/docs/3_step4_preparing-input-files/index.md
+++ b/docs/3_step4_preparing-input-files/index.md
@@ -31,7 +31,7 @@ The section [Static Maps](../4_Static-Maps-introduction) provides detailed guide
+ [land use maps](../4_Static-Maps_land-use/): fraction of forest; fraction of irrigated crops; fraction of rice crops; fraction of inland water; fraction of sealed surfaces; fraction of other land uses.
+ [land use depending](../4_Static-Maps_land-use-depending/):crop coefficient; crop group number; Manning/s's surface roughness; soil depth.
+ [soil hydraulic properties](../4_Static-Maps_soil-hydraulic-properties/): saturated hydraulic conductivity; soil water content at saturation; residual soil water content; parameters alpha and lambda of Van Genuchten's equations.
-+ [channel geometry](../4_Static-Maps_channel-geometry/): channels mask; channels side slope; channels length; channels gradient; Manning's rougheness coefficient of the channels; channels bottom width; floodplain width; bankfull channels depth.
++ [channel geometry](../4_Static-Maps_channel-geometry/): channels mask; channels side slope; channels length; channels gradient; Manning's rougheness coefficient of the channels; channels bottom width; floodplain width; bankfull channels depth; MCT diffusive wave routing channels.
+ [leaf area index](../4_Static-Maps_leaf-area-index/): evolution of vegetation over time (leaf area index) for land covers forest, irrigated areas, others.
+ [reservoirs and lakes](../4_Static-Maps_reservoirs-lakes/): lake mask; lakes ID points; reservoirs ID points. These maps are required only upon activation of the [lakes module](https://ec-jrc.github.io/lisflood-model/3_02_optLISFLOOD_lakes/) and/or of the [reservoirs module](https://ec-jrc.github.io/lisflood-model/3_03_optLISFLOOD_reservoirs/).
+ [rice calendar](../4_Static-Maps_rice-calendar/): rice calendar for planting and harvesting seasons. These maps are required only when activating the [rice irrigation module](https://ec-jrc.github.io/lisflood-model/2_17_stdLISFLOOD_irrigation/)
@@ -39,9 +39,9 @@ The section [Static Maps](../4_Static-Maps-introduction) provides detailed guide
+ water demand maps: domestic, energetic, livestock, industrial water use. These maps represent the time series of spatially distributed values of water demand for domestic, energetic, livestock, and industrial water use. These maps are required only when activating the [water use module](https://ec-jrc.github.io/lisflood-model/2_18_stdLISFLOOD_water-use/)
+ outlet points: locations and IDs of the points for which LISFLOOD provides the time series of discharge values.
-#### Role of "mask" and "channels" maps
+#### Role of "mask", "channels" ans "channelsMCT" maps
-The mask map (i.e. "area.map") defines the model domain. In order to avoid unexpected results, **it is vital that all maps that are related to topography, land use and soil are defined** (i.e. don't contain a missing value) for each pixel that is "true" (has a Boolean 1 value) on the mask map. The same applies for all meteorological input and the Leaf Area Index maps. Similarly, all pixels that are "true" on the channels map must have some valid (non-missing) value on each of the channel parameter maps. Undefined pixels can lead to unexpected behaviour of the model, output that is full of missing values, loss of mass balance and possibly even model crashes. Some maps needs to have values in a defined range e.g. the gradient map has to be greater than 0.
+The mask map (i.e. "area.map") defines the model domain. In order to avoid unexpected results, **it is vital that all maps that are related to topography, land use and soil are defined** (i.e. don't contain a missing value) for each pixel that is "true" (has a Boolean 1 value) on the mask map. The same applies for all meteorological input and the Leaf Area Index maps. Similarly, all pixels that are "true" on the channels map must have some valid (non-missing) value on each of the channel parameter maps. At the same time, all pixels that have value "true" in the MCT rivers mask must also belong to the "channels" map. Undefined pixels can lead to unexpected behaviour of the model, output that is full of missing values, loss of mass balance and possibly even model crashes. Some maps needs to have values in a defined range e.g. the gradient map has to be greater than 0.
#### Geometrical properties of the computational grid cell
diff --git a/docs/4_Static-Maps_channel-geometry/index.md b/docs/4_Static-Maps_channel-geometry/index.md
index 638e5b0a..f8588a2c 100644
--- a/docs/4_Static-Maps_channel-geometry/index.md
+++ b/docs/4_Static-Maps_channel-geometry/index.md
@@ -2,6 +2,7 @@
In the LISFLOOD model flow through the channel is simulated using the kinematic wave equations. Channel maps describe the sub grid information of the channel geometry, i.e. the length, slope, width and depth of the main channel inside a grid-cell.
+ The **channel mask** map is the Boolean field with '1' for all grid-cells with channels and NoData for all grid-cells with no channels.
++ The **MCT channel mask** map is the Boolean field with '1' for all river grid-cells using MCT wave routing and '0' for other channel grid cells.
+ The **channel side slope** map (referred as 's' in Figure 41) defines the slope of the channel banks.
+ The **channel length** map is the length of the river in each grid-cell, and the value can exceed grid-size to account for meandering rivers.
+ The **channel gradient** (or channel slope) map is the average gradient of the main river inside a cell.
@@ -23,6 +24,7 @@ Channel characteristics, explained above, are shown in the Figure 41 below.
Type: Boolean | Units: -;
Range: NoData or 1 |Boolean map that identifies the channel grid-cells |
+|MCT Channel mask | chanmct.nc;
Type: Boolean | Units: -;
Range: [0-1] |Boolean map that identifies the channel grid-cells using the MCT diffusive river routing |
|Side slope |chans.nc;
Type: Float32 | Units: m;
Range>0 |Channel side slope|
|Channel length |chanlength.nc;
Type: Float32 | Units: m;
Range>0 |Channel length (value can exceed grid size, to account for meandering rivers)|
|Channel gradient |changrad.nc;
Type: Float32 |Units: m/m;
Range: [0-1] |Channel longitudinal gradient|
@@ -45,6 +47,9 @@ Channel characteristics, explained above, are shown in the Figure 41 below.
R.: map>0 | Saturated conductivity layer 2 |
| **CHANNEL GEOMETRY** | | | |
| Channels | chan.map | U.: [-]
R.: 0 or 1 | Map with Boolean 1 for all channel pixels, and Boolean 0 for all other pixels on MaskMap |
+| MCT Channels | chanmct.map | U.: [-]
R.: 0 or 1 | Map with Boolean 1 for channel pixels using MCT diffusive wave routing, and Boolean 0 for all other pixels on MaskMap |
| ChanGrad | changrad.map | U.: $\frac{m} {m}$
R.: map > 0
!!! | Channel gradient |
| ChanMan | chanman.map | U.: $m^{-1/3} s$
R.: map > 0 | Manning's roughness coefficient for channels |
| ChanLength | chanleng.map | U.: $m$
R.: map > 0 | Channel length (can exceed grid size, to account for meandering rivers) |
diff --git a/docs/4_annex_output-files/index.md b/docs/4_annex_output-files/index.md
index 2811aa3d..6843e6e5 100644
--- a/docs/4_annex_output-files/index.md
+++ b/docs/4_annex_output-files/index.md
@@ -9,7 +9,8 @@ LISFLOOD can generate a wide variety of outputs. Output files can be time series
| Description | Units | File name |
| ------------------------------------------------------------ | ---------------- | --------------------- |
| **RATE VARIABLES AT GAUGES** | | |
-| $^{1,2}$ channel discharge | $\frac{m^3} {s}$ | dis.tss |
+| $^{1,2}$ average river discharge | $\frac{m^3} {s}$ | dis.tss |
+| $^{1,2}$ instantaneous rivers discharge | $\frac{m^3} {s}$ | chanq.tss |
| **NUMERICAL CHECKS** | | |
| $^2$ cumulative mass balance error | $m^3$ | mbError.tss |
| $^2$ cumulative mass balance error, expressed as mm water slice (average over catchment) | $mm$ | mbErrorMm.tss |
diff --git a/docs/4_annex_settings_and_options/index.md b/docs/4_annex_settings_and_options/index.md
index 06199206..182f24a4 100644
--- a/docs/4_annex_settings_and_options/index.md
+++ b/docs/4_annex_settings_and_options/index.md
@@ -20,6 +20,7 @@ The Tables below presents a comprehensive list of setting options, inputs, and o
| lfoptions | | WATER ABSTRACTION | | useWaterDemandAveYear | | | switch | input | Use "average" year for water demand and loop it over years |
| lfoptions | | TRANSMISSION LOSS | | TransLoss | | | switch | input | Activate transmission loss |
| lfoptions | | DOUBLE KINEMATIC WAVE | | SplitRouting | | | switch | input | Activate double kinematic wave routing |
+| lfoptions | | MCT DIFFUSIVE WAVE | | MCTRouting | | | switch | input | Activate MCT diffusive wave routing |
| lfoptions | | WATER ABSTRACTION | | wateruse | | | switch | input | Activate water use computation |
| lfoptions | | GROUNDWATER | | groundwaterSmooth | | | switch | input | Activate smoothing for groundwater to correct error by using windowtotal, based on groundwater bodies and catchments |
| lfoptions | | WATER ABSTRACTION | | wateruseRegion | | | switch | input | Use water regions in water use module |
@@ -127,6 +128,7 @@ The Tables below presents a comprehensive list of setting options, inputs, and o
| lfuser | ROUTING | KINEMATIC WAVE | n | CalChanMan | $(PathParams)/params_CalChanMan1 | $(PathParams)/params_CalChanMan1 | value/map | input calib par | It is a multiplier that is applied to the Manning's roughness map of the channel system default: 2.0 $(PathParams)/params_CalChanMan1.nc |
| lfuser | SNOW AND FROST | SNOW | Cm | SnowMeltCoef | $(PathParams)/params_SnowMeltCoef | $(PathParams)/params_SnowMeltCoef | value/map | input calib par | Snowmelt coefficient [mm/deg C /day]. It is the degree-day factor that controls the rate of snowmelt default: 4.0 $(PathParams)/params_SnowMeltCoef.nc SRM: 0.45 cm/C/day ( = 4.50 mm/C/day), Kwadijk: 18 mm/C/month (= 0.59 mm/C/day) See also Martinec et al., 1998. |
| lfuser | ROUTING | DOUBLE KINEMATIC WAVE | n? | CalChanMan2 | $(PathParams)/params_CalChanMan2 | $(PathParams)/params_CalChanMan2 | value/map | input calib par | PBchange Multiplier applied to Channel Manning's n for second routing line default: 3.0 $(PathParams)/params_CalChanMan2.nc |
+| lfuser | ROUTING | MCT DIFFUSIVE WAVE | n? | CalChanMan3 | $(PathParams)/params_CalChanMan3 | $(PathParams)/params_CalChanMan3 | value/map | input calib par | Multiplier [-] applied to Channel Manning's n for MCT diffusive wave routing default: 3.0 $(PathParams)/params_CalChanMan3.nc |
| lfuser | CALIBRATION | LAKES | | LakeMultiplier | $(PathParams)/params_LakeMultiplier | $(PathParams)/params_LakeMultiplier | value/map | input calib par | default: 1.0 $(PathParams)/params_LakeMultiplier.nc Multiplier applied to the lake parameter A |
| lfuser | CALIBRATION | RESERVOIRS | | adjust_Normal_Flood | $(PathParams)/params_adjust_Normal_Flood | $(PathParams)/params_adjust_Normal_Flood | value/map | input calib par | default: 0.8 $(PathParams)/params_adjust_Normal_Flood.nc adjusting the balance between normal and flood storage big values (= closer to 1.0 = close to flood) results in keeping the outflow longer at rnormq Range 0.01 - 0.99 |
| lfuser | CALIBRATION | RESERVOIRS | | ReservoirRnormqMult | $(PathParams)/params_ReservoirRnormqMult | $(PathParams)/params_ReservoirRnormqMult | value/map | input calib par | default: 1.0 $(PathParams)/params_ReservoirRnormqMult.nc Reservoir rnormq (normal outflow) fraction of the value in table rnormq.txt range: 0.5 - 2 |
@@ -164,6 +166,8 @@ The Tables below presents a comprehensive list of setting options, inputs, and o
| lfuser | INITIAL CONDITIONS | SOIL | | ThetaInit3Value | -9999 | $(PathInit)/thc | value/map | input initial/internal | initial soil moisture content layer 2 -9999: use field capacity values It is the initial moisture content [mm3 mm-3] of the lower soil layer. |
| lfuser | INITIAL CONDITIONS | DOUBLE KINEMATIC WAVE | | CrossSection2AreaInitValue | -9999 | $(PathInit)/ch2cr | value/map | input initial/internal | initial channel cross-sectional area [m2] of the water in the river channels for 2nd routing channel -9999: use 0 |
| lfuser | INITIAL CONDITIONS | DOUBLE KINEMATIC WAVE | | PrevSideflowInitValue | -9999 | $(PathInit)/chside | value/map | input initial/internal | initial inflow from each pixel to the channel [mm]. A value of -9999 sets the initial amount of sideflow to thye channel to 0 |
+| lfuser | INITIAL CONDITIONS | MCT DIFFUSIVE WAVE | | PrevCmMCTInitValue | -9999 | $(PathInit)/prevcm | value/map | input initial/internal | Courant number at previous step for MCT diffusive routing [-]. A value of -9999 sets the initial value to 0 |
+| lfuser | INITIAL CONDITIONS | MCT DIFFUSIVE WAVE | | PrevDmMCTInitValue | -9999 | $(PathInit)/prevdm | value/map | input initial/internal | Reynolds number at previous step for MCT diffusive routing [-]. A value of -9999 sets the initial value to 1 |
| lfuser | INITIAL CONDITIONS | LAKES | | LakeInitialLevelValue | -9999 | $(PathInit)/lakeh | value/map | input initial/internal | Initial lake level [m] -9999 sets initial value to steady-state level |
| lfuser | INITIAL CONDITIONS | KINEMATIC WAVE | | PrevDischarge | -9999 | $(PathInit)/chanq | value/map | input initial/internal | initial discharge from previous run for lakes, reservoirs and transmission loss only needed for lakes reservoirs and transmission loss -9999: use 0 It is the initial discharge from previous run [m3s-1] used for lakes, reservoirs and transmission loss (only needed if option is on for lakes or reservoirs or transmission loss). Note that PrevDischarge is discharge as an average over the time step (a flux) . |
| lfuser | INITIAL CONDITION FOREST | INTERCEPTION | | CumIntForestInitValue | 0 | $(PathInit)/cumf | value/map | input initial/internal | cumulative interception forest [mm] |
@@ -241,6 +245,8 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| lfuser | ROUTING | KINEMATIC WAVE | | OFDepRef | 5 | 5 | value | input par | It is a reference flow depth from which the flow velocity of the surface runoff is calculated [mm] Reference depth of overland flow [mm], used to compute overland flow Alpha for kin. wave |
| lfuser | ROUTING | KINEMATIC WAVE | | GradMin | 0.001 | 0.001 | value | input par | Minimum slope gradient of the surface (for kin. wave: slope cannot be 0) It is a lower limit for the slope gradient used in the calculation of the surface runoff flow velocity [m m-1] |
| lfuser | ROUTING | KINEMATIC WAVE | | ChanGradMin | 0.0001 | 0.0001 | value | input par | Minimum channel gradient (for kin. wave: slope cannot be 0) It is a lower limit for the channel gradient used in the calculation of the channel flow velocity [m m-1] |
+| lfuser | ROUTING | MCT DIFFUSIVE WAVE | | ChannelsMCT | $(PathRoot)/maps/chanmct | $(PathRoot)/maps/chanmct | map | input | Boolean map with value 1 at channel pixels where MCT is used, and 0 at all other pixels |
+| lfuser | ROUTING | MCT DIFFUSIVE WAVE | | ChanGradMaxMCT | 0.001 | 0.001 | value | input par | Maximum channel gradient for channels using MCT routing [-] (for MCT wave: slope cannot be 0) [m m-1] |
| lfuser | ROUTING | DOUBLE KINEMATIC WAVE | | QSplitMult | 2 | 2 | value | input par | PBchange Multiplier applied to average Q to split into a second line of routing |
| lfuser | NUMERICS | SOIL | | CourantCrit | 0.4 | 0.4 | value | input par | Minimum value for Courant condition in soil moisture routine. Always less than or equal to 1. Small values result in improved numerical accuracy, at the expense of increased computing time (more sub-steps needed). If reported time series of soil moisture contain large jumps, lowering CourantCrit should fix this |
| lfuser | OPTION RESERVOIR | RESERVOIRS | | DtSecReservoirs | 3600 | 3600 | value | input | Sub time step used for reservoir simulation [s]. Within the model, the smallest out of DtSecReservoirs and DtSec is used. |
@@ -295,11 +301,14 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| lfbinding | TRANSMISSION LOSS | TRANSMISSION | | UpAreaTrans | $(UpAreaTrans) | $(UpAreaTrans) | 0 | input | upstream area for transmission loss |
| lfbinding | ROUTING | KINEMATIC WAVE | | CalChanMan | $(CalChanMan) | $(CalChanMan) | 0 | input | It is a multiplier that is applied to the Manning's roughness map of the channel system default: 2.0 $(PathParams)/params_CalChanMan1.nc |
| lfbinding | ROUTING | DOUBLE KINEMATIC WAVE | | CalChanMan2 | $(CalChanMan2) | $(CalChanMan2) | value/map | input | PBchange Multiplier applied to Channel Manning's n for second routing line default: 3.0 $(PathParams)/params_CalChanMan2.nc |
+| lfbinding | ROUTING | MCT DIFFUSIVE WAVE | | CalChanMan3 | $(CalChanMan3) | $(CalChanMan3) | value/map | input | Multiplier [-] applied to Channel Manning's n for MCT routing default: 3.0 $(PathParams)/params_CalChanMan3.nc |
| lfbinding | ROUTING | DOUBLE KINEMATIC WAVE | | QSplitMult | $(QSplitMult) | $(QSplitMult) | value | input | PBchange Multiplier applied to average Q to split into a second line of routing |
| lfbinding | ROUTING | KINEMATIC WAVE | | beta | $(beta) | $(beta) | 0 | input | It is the routing coefficient in Manning's equation (2/3). kinematic wave parameter: 0.6 is for broad sheet flow |
| lfbinding | ROUTING | KINEMATIC WAVE | | OFDepRef | $(OFDepRef) | $(OFDepRef) | 0 | input | It is a reference flow depth from which the flow velocity of the surface runoff is calculated [mm] Reference depth of overland flow [mm], used to compute overland flow Alpha for kin. wave |
| lfbinding | ROUTING | KINEMATIC WAVE | | GradMin | $(GradMin) | $(GradMin) | 0 | input | Minimum slope gradient of the surface (for kin. wave: slope cannot be 0) It is a lower limit for the slope gradient used in the calculation of the surface runoff flow velocity [m m-1] |
| lfbinding | ROUTING | KINEMATIC WAVE | | ChanGradMin | $(ChanGradMin) | $(ChanGradMin) | 0 | input | Minimum channel gradient (for kin. wave: slope cannot be 0) It is a lower limit for the channel gradient used in the calculation of the channel flow velocity [m m-1] |
+| lbinding | ROUTING | MCT DIFFUSIVE WAVE | | ChannelsMCT | $(ChannelsMCT) | $(ChannelsMCT) | map | input | Boolean map with value 1 at channel pixels where MCT is used, and 0 at all other pixels |
+| lbinding | ROUTING | MCT DIFFUSIVE WAVE | | ChanGradMaxMCT | $(ChanGradMaxMCT) | $(ChanGradMaxMCT) | map | input | Maximum channel gradient for channels using MCT routing [-] (for MCT wave: slope cannot be 0) |
| lfbinding | NUMERICS | SOIL | | CourantCrit | $(CourantCrit) | $(CourantCrit) | value | input | Minimum value for Courant condition in soil moisture routine. Always less than or equal to 1. Small values result in improved numerical accuracy, at the expense of increased computing time (more sub-steps needed). If reported time series of soil moisture contain large jumps, lowering CourantCrit should fix this |
| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | GROUNDWATER | | LZAvInflowMap | $(PathMaps)/lzavin.map | $(PathMaps)/lzavin.map | map | input | $(PathInit)/lzavin.map Reported map of average percolation rate from upper to lower groundwater zone (reported for end of simulation) |
| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | KINEMATIC WAVE | | AvgDis | $(PathMaps)/avgdis.map | $(PathMaps)/avgdis.map | map | input | $(PathInit)/avgdis.map CHANNEL split routing in two lines Average discharge map [m3/s] |
@@ -322,6 +331,8 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | SOIL | | ThetaInit3Value | $(ThetaInit3Value) | $(ThetaInit3Value) | 0 | input | initial soil moisture content layer 2 -9999: use field capacity values |
| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | DOUBLE KINEMATIC WAVE | | CrossSection2AreaInitValue | $(CrossSection2AreaInitValue) | $(CrossSection2AreaInitValue) | value/map | input | initial channel crosssection for 2nd routing channel -9999: use 0 |
| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | DOUBLE KINEMATIC WAVE | | PrevSideflowInitValue | $(PrevSideflowInitValue) | $(PrevSideflowInitValue) | value/map | input | initial inflow from each pixel to the channel [mm]. A value of -9999 sets the initial amount of sideflow to thye channel to 0 |
+| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | MCT DIFFUSIVE WAVE | | PrevCmMCTInitValue | $(PrevCmMCTInitValue) | $(PrevCmMCTInitValue) | value/map | input initial/internal | Courant number at previous step for MCT diffusive routing [-]. A value of -9999 sets the initial value to 0 |
+| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | MCT DIFFUSIVE WAVE | | PrevDmMCTInitValue | $(PrevDmMCTInitValue) | $(PrevDmMCTInitValue) | value/map | input initial/internal | Reynolds number at previous step for MCT diffusive routing [-]. A value of -9999 sets the initial value to 1 |
| lfbinding | INITIAL CONDITIONS WATER BALANCE MODEL | KINEMATIC WAVE | | PrevDischarge | $(PrevDischarge) | $(PrevDischarge) | 0 | input | initial discharge from previous run for lakes, reservoirs and transmission loss only needed for lakes reservoirs and transmission loss -9999: use 0 |
| lfbinding | INITIAL CONDITION FOREST | INTERCEPTION | | CumIntForestInitValue | $(CumIntForestInitValue) | $(CumIntForestInitValue) | 0 | input | cumulative interception forest [mm] |
| lfbinding | INITIAL CONDITION FOREST | GROUNDWATER | | UZForestInitValue | $(UZForestInitValue) | $(UZForestInitValue) | 0 | input | Initial water storage water in upper groundwater zone for forest [mm] |
@@ -403,6 +414,8 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | DOUBLE KINEMATIC WAVE | | CrossSection2State | $(PathOut)/ch2cr | $(PathOut)/ch2cr | map | output/state | Cross section area for split routing [m2] |
| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | DOUBLE KINEMATIC WAVE | | ChSideState | $(PathOut)/chside | $(PathOut)/chside | map | output/state | Reported sideflow to channel for first line of routing [m3/s] |
| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | ROUTING | | ChanQState | $(PathOut)/chanq | $(PathOut)/chanq | map | output/state | Reported istantaneous discharge at end of computation step [cu m/s] ChanQ |
+| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | ROUTING | | PrevCmMCTState | $(PathOut)/prevcm | $(PathOut)/prevcm | map | output/state | Reported Courant number at previous step for MCT routing |
+| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | ROUTING | | PrevDmMCTState | $(PathOut)/prevdm | $(PathOut)/prevdm | map | output/state | Reported Reynolds number at previous step for MCT routing |
| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | SOIL | | DSLRForestState | $(PathOut)/dslf | $(PathOut)/dslf | map | output/state | Reported days since last rain for forest |
| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | INTERCEPTION | | CumInterceptionForestState | $(PathOut)/cumf | $(PathOut)/cumf | map | output/state | Reported interception storage for forest |
| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | SOIL | | Theta1ForestState | $(PathOut)/thfa | $(PathOut)/thfa | map | output/state | theta for soil layer 1a forest fraction |
@@ -438,7 +451,10 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| lfbinding | REPORTED OUTPUT MAPS (END) | DOUBLE KINEMATIC WAVE | | CrossSection2End | $(PathOut)/ch2cr.end | $(PathOut)/ch2cr.end | map | output/end | Cross section area for split routing [m2] |
| lfbinding | REPORTED OUTPUT MAPS (END) | DOUBLE KINEMATIC WAVE | | ChSideEnd | $(PathOut)/chside.end | $(PathOut)/chside.end | map | output/end | Reported channel side flow |
| lfbinding | REPORTED OUTPUT MAPS (END) | ROUTING | | ChanQEnd | $(PathOut)/chanq.end | $(PathOut)/chanq.end | map | output/end | Reported istantaneous discharge at end of computation step [cu m/s] ChanQ |
-| lfbinding | REPORTED OUTPUT MAPS (END) | KINEMATIC WAVE | | DischargeEnd | $(PathOut)/dis.end | $(PathOut)/dis.end | map | output/end | Reported discharge [m3/s] for lakes and reservoirs |
+| lfbinding | REPORTED OUTPUT MAPS (END) | KINEMATIC WAVE | | DischargeEnd | $(PathOut)/dis.end | $(PathOut)/dis.end | map | output/end | Reported average discharge [m3/s] for lakes and reservoirs |
+| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | ROUTING | | PrevCmMCTEnd | $(PathOut)/prevcm.end | $(PathOut)/prevcm.end | map | output/end | Reported Courant number at previous step for MCT routing |
+| lfbinding | REPORTED OUTPUT MAPS (STATE VARIABLES AT SELECTED TIME STEPS) | ROUTING | | PrevDmMCTEnd | $(PathOut)/prevdm.end | $(PathOut)/prevdm.end
+| Reported Reynolds number at previous step for MCT routing |
| lfbinding | REPORTED OUTPUT MAPS (END) | SOIL | | DSLRForestEnd | $(PathOut)/dslf.end | $(PathOut)/dslf.end | map | output/end | Reported days since last rain for forest |
| lfbinding | REPORTED OUTPUT MAPS (END) | INTERCEPTION | | CumInterceptionForestEnd | $(PathOut)/cumf.end | $(PathOut)/cumf.end | map | output/end | Reported interception storage for forest |
| lfbinding | REPORTED OUTPUT MAPS (END) | SOIL | | Theta1ForestEnd | $(PathOut)/thfa.end | $(PathOut)/thfa.end | map | output/end | Reported volumetric soil moisture content for soil layer 1a for forest [V/V] |
@@ -487,6 +503,7 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| lfoptions | | | | useWaterDemandAveYear | 0 | WATER ABSTRACTION | switch | input | 0 | | | | Use "average" year for water demand and loop it over years |
| lfoptions | | | | TransLoss | 0 | TRANSMISSION LOSS | switch | input | 0 | | | | Activate transmission loss |
| lfoptions | | | | SplitRouting | 0 | KINEMATIC WAVE | switch | input | 0 | | | | Activate double kinematic wave routing |
+| lfoptions | | | | MCTRouting | 0 | MCT DIFFUSIVE WAVE| switch | input | 0 | | | | Activate MCT diffusive wave routing |
| lfoptions | | | | wateruse | 0 | WATER ABSTRACTION | switch | input | 0 | | | | Activate water use computation |
| lfoptions | | | | groundwaterSmooth | 0 | GROUNDWATER | switch | input | 0 | | | | Activate smoothing for groundwater to correct error by using windowtotal, based on groundwater bodies and catchments |
| lfoptions | | | | wateruseRegion | 0 | WATER ABSTRACTION | switch | input | 0 | | | | Use water regions in water use module |
@@ -579,7 +596,7 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| key | outputVar | where | repoption | restrictoption | operation | Description |
| -------------------- | ------------------------------------------------------------ | -------------- | ---------------------- | --------------------------- | --------------- | ------------------------------------------------------------ |
-| DisTS | ChanQAvg | Gauges | repDischargeTs | | | Reported discharge at gauges [m3/s] |
+| DisTS | ChanQAvg | Gauges | repDischargeTs | | | Reported average discharge at gauges [m3/s] |
| WaterLevelTS | WaterLevel | Gauges | repWaterLevelTs | nonInit,simulateWaterLevels | | Reported water level [m] |
| StepsSoilTS | NoSubSteps | 1 | repInternalCom | | mapmaximum | Number of sub-steps needed for soil moisture routine |
| StepsWaterUseTS | NoWaterUseExe | 1 | repInternalCom | nonInit,wateruse | | number of loops needed for water use routine |
@@ -705,6 +722,8 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| CumIntSealedState | CumInterSealed | mm | repStateMaps | | | | nonInit | Reported cumulative depressions storage [mm] |
| DischargeMaps | ChanQAvg | m3/s | repStateMaps | | | | nonInit | Reported discharge [cu m/s] (average over the timesteps) |
| ChanQState | ChanQ | m3/s | repStateMaps | | | | nonInit | Reported istantaneous discarge at end of time step |
+| PrevCmMCTState | PrevCm0 | - | repStateMaps | | | | nonInit | Reported Courant number at previous step for MCT routing |
+| PrevDmMCTState | PrevDm0 | - | repStateMaps | | | | nonInit | Reported Reynolds number at previous step for MCT routing |
| LakeLevelState | LakeLevel | m | repStateMaps | | | | nonInit,simulateLakes | Output map(s) with lake level [m] |
| ReservoirFillState | ReservoirFill | - | repStateMaps | | | | nonInit,simulateReservoirs | Output map(s) with Reservoir Filling [-] |
| CrossSection2State | CrossSection2Area | m | repStateMaps | | | | nonInit,SplitRouting | Cross section area for split routing [m2] |
@@ -741,6 +760,8 @@ param| lfuser | IRRIGATION AND WATER ABSTRACTION
| CumIntSealedEnd | CumInterSealed | mm | | | repEndMaps | | nonInit | Reported cumulative depressions storage [mm] |
| DischargeEnd | ChanQAvg | m3/s | | | repEndMaps | | nonInit | Reported discharge [cu m/s] (average over the timesteps) |
| ChanQEnd | ChanQ | m3/s | | | repEndMaps | | nonInit | Reported istantaneous discarge at end of time step |
+| PrevCmMCTEnd | PrevCm0 | - | | | repEndMaps | | nonInit | Reported Courant number at previous step for MCT routing |
+| PrevDmMCTEnd | PrevDm0 | - | | | repEndMaps | | nonInit | Reported Reynolds number at previous step for MCT routing |
| LakeLevelEnd | LakeLevel | m | | | repEndMaps | | nonInit,simulateLakes | Reported lake level |
| ReservoirFillEnd | ReservoirFill | m | | | repEndMaps | | nonInit,simulateReservoirs | Reported reservoir filling |
| CrossSection2End | CrossSection2Area | m | | | repEndMaps | | nonInit,SplitRouting | Cross section area for split routing [m2] |
diff --git a/docs/4_annex_tests/index.md b/docs/4_annex_tests/index.md
index 04343343..a0aaed0d 100644
--- a/docs/4_annex_tests/index.md
+++ b/docs/4_annex_tests/index.md
@@ -785,6 +785,27 @@ Verifies chunking files using NetCDFTimeChunks with values `1`, `10`, `auto` and
| test_chunking_24h | 86400 |30/07/2016 06:00 - 01/09/2016 06:00 | compares output chunking files with reference files |
| test_chunking_6h | 21600 |30/07/2016 06:00 - 01/09/2016 06:00 | compares output chunking files with reference files |
+
+### Testing MCT diffusive wave routing
+
+A set of tests is available to verify the correct functioning of the MCT diffusive wave routing. In these tests all OSLisflood options are switched off and the model runs with core components and river routing. Two sets of tests are available: 1) 'fast' tests that can be run routinely and only take a few minutes; 2) 'slow' tests that need to run before merging new model developments.
+
+#### Testing warm start for MCT diffusive routing
+
+Test ensures that a long cold run is equivalent to an initial cold start + repeated warm starts, when MCT diffusive wave routing is used. The test runs OSLisflood continuously for a long period of time and saves all outputs in a reference folder. It then runs OSLisflood on the same period but restarting the model at every step (start and stop). Outputs of the two runs are compared and maps are expected to be identical.
+During this test, both 'slow' and 'fast' tests generate the reference dataset dynamically and are designed to work when model developments change the model outputs.
+
+#### Testing inflow for MCT diffusive routing
+
+This test verifies the correct functioning of the option inflow, when using MCT diffusive routing. The test runs OSLisflood continuously for a long period of time and saves river discharge time series for a point along the main river (point A) and the point immediatly downstream (point B). It then runs OSLisflood again on the sub-catchment downstream of point B and using the generated river discharge at point A as an inflow at point B. Discharge time series from the two runs at point B are compared, differences should be <0.001 m3/s.
+During this test, both 'slow' and 'fast' tests generate the reference dataset dynamically and are designed to work when model developments change the model outputs.
+
+#### Testing results for MCT diffusive routing
+
+These tests perform short and long (max 1 year) execution of OSlisflood on the Po test catchments clipped from the full setup (EFAS) and assert that results are within defined tolerances if compared with test oracle. All tests have all switches turned off and run with core OSLisflood components and river routing. Note that the tests will fail if model development changes outputs from any of the core model componets.
+
+
+
## Release test for EFAS and GloFAS
At each release, in addition to pass all tests described above, OSLisflood is tested also with full domains of EFAS and GLOFAS.
diff --git a/requirements.txt b/requirements.txt
index 297cc0c0..4e0395cb 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -4,6 +4,7 @@ cftime>=1.0.4.2
coverage>=6.0
dask>=2021.10.0
future>=0.18.0
+#GDAL==3.5.3
GDAL>=3.2.0
importlib-metadata<5.0.0
iniconfig>=1.1.1
diff --git a/src/lisflood/Lisflood_dynamic.py b/src/lisflood/Lisflood_dynamic.py
index 70f4c188..d547ce81 100644
--- a/src/lisflood/Lisflood_dynamic.py
+++ b/src/lisflood/Lisflood_dynamic.py
@@ -24,6 +24,7 @@
import gc
from pcraster.framework import DynamicModel
+
import numpy as np
from lisflood.global_modules.errors import LisfloodWarning
@@ -31,6 +32,10 @@
from .global_modules.settings import CDFFlags, LisSettings, MaskInfo
+from pcraster import upstream
+from lisflood.global_modules.add1 import decompress, compressArray
+
+
class LisfloodModel_dyn(DynamicModel):
# =========== DYNAMIC ====================================================
@@ -175,6 +180,7 @@ def splitlanduse(array1, array2=None, array3=None):
# ************************************************************
maskinfo = MaskInfo.instance()
self.sumDisDay = maskinfo.in_zero()
+
# sums up discharge of the sub steps
for NoRoutingExecuted in range(self.NoRoutSteps):
self.routing_module.dynamic(NoRoutingExecuted)
@@ -191,27 +197,28 @@ def splitlanduse(array1, array2=None, array3=None):
# if option['simulatePolders']:
# ChannelToPolderM3=ChannelToPolderM3Old;
- if option['InitLisflood'] or (not(option['SplitRouting'])):
- # kinematic routing
- self.ChanM3 = self.ChanM3Kin.copy()
- # Total channel storage [cu m], equal to ChanM3Kin
- else:
- # split routing
- self.ChanM3 = self.ChanM3Kin + self.Chan2M3Kin - self.Chan2M3Start
+ # # moved to routing.py
+ # # Calculate total water storage in river channel
+ # if option['InitLisflood'] or (not(option['SplitRouting'])):
+ # # kinematic routing
+ # self.ChanM3 = self.ChanM3Kin.copy()
+ # # Total channel storage [m3], equal to ChanM3Kin fir kinematic routing only at t+dt
+ # else:
+ # # split routing
+ # self.ChanM3 = self.ChanM3Kin + self.Chan2M3Kin - self.Chan2M3Start
+ # # Total channel storage [m3] = Volume in main channel (ChanM3Kin) + volume above bankfull (Chan2M3Kin - Chan2M3Start)
+ # # at t+dt
- # Total channel storage [cu m], equal to ChanM3Kin
- # sum of both lines
- # CrossSection2Area = pcraster.max(scalar(0.0), (self.Chan2M3Kin - self.Chan2M3Start) / self.ChanLength)
self.TotalCrossSectionArea = self.ChanM3 * self.InvChanLength
+ # Total river channel cross-section area at t+dt
self.sumDis += self.sumDisDay
+ # Accumulated average discharge over the entire simulation period
+
self.ChanQAvg = self.sumDisDay/self.NoRoutSteps
+ # Average channel outflow over the model computation step
- # Total volume of water in channel per inv channel length
- # New cross section area (kinematic wave)
- # This is the value after the kinematic wave, so we use ChanM3Kin here
- # (NOT ChanQKin, which is average discharge over whole step, we need state at the end of all iterations!)
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if not(option['dynamicWave']):
@@ -222,12 +229,16 @@ def splitlanduse(array1, array2=None, array3=None):
# Set water level dynamic wave to dummy value (needed
if option['InitLisflood'] or option['repAverageDis']:
- self.CumQ += self.ChanQ
+ # self.CumQ += self.ChanQ
+ self.CumQ += self.ChanQAvg
+ #cmcheck - we should use ChanQAvg here not ChanQ
self.avgdis = self.CumQ/self.TimeSinceStart
- # to calculate average discharge
+ # to calculate average discharge over the entire simulation
- self.DischargeM3Out += np.where(self.AtLastPointC ,self.ChanQ * self.DtSec,0)
- # Cumulative outflow out of map
+ #self.DischargeM3Out += np.where(self.AtLastPointC ,self.ChanQ * self.DtSec,0)
+ self.DischargeM3Out += np.where(self.AtLastPointC, self.ChanQAvg * self.DtSec, 0)
+ # Cumulative outflow out of map
+ # cmcheck - we should use ChanQAvg here not ChanQ
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Calculate water level
diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py
index 8a84816a..346c1dff 100644
--- a/src/lisflood/Lisflood_initial.py
+++ b/src/lisflood/Lisflood_initial.py
@@ -217,6 +217,10 @@ def __init__(self):
# CHANNEL INITIAL SPLIT UP IN SECOND CHANNEL
self.surface_routing_module.initialSecond()
+ if option.get('MCTRouting'):
+ self.routing_module.initialMCT()
+ # initialising Muskingum-Cunge-Todini routing for channel
+
self.evapowater_module.initial()
self.riceirrigation_module.initial()
self.waterabstraction_module.initial()
diff --git a/src/lisflood/global_modules/default_options.py b/src/lisflood/global_modules/default_options.py
index 84c88181..a7a715d2 100644
--- a/src/lisflood/global_modules/default_options.py
+++ b/src/lisflood/global_modules/default_options.py
@@ -8,6 +8,7 @@
'InitLisfloodwithoutSplit': False,
'MonteCarlo': False,
'SplitRouting': False,
+ 'MCTRouting': False,
'TemperatureInKelvin': False,
'TransLoss': False,
'TransientLandUseChange': False,
@@ -677,8 +678,21 @@
output_var='SeepTopToSubB[2]', unit='mm',
end=[], steps=[], all=['repPercolationMaps'],
restrictoption=['nonInit'], monthly=False,
- yearly=False),
-
+ yearly=False),
+ 'PrevCmMCTEnd': ReportedMap(name='PrevCmMCTEnd', output_var='PrevCm0', unit='-',
+ end=['repEndMaps'], steps=[], all=[],
+ restrictoption=[], monthly=False, yearly=False),
+ 'PrevCmMCTState': ReportedMap(name='PrevCmMCTState', output_var='PrevCm0', unit='-',
+ end=[], steps=['repStateMaps'], all=[],
+ restrictoption=['nonInit'], monthly=False,
+ yearly=False),
+ 'PrevDmMCTEnd': ReportedMap(name='PrevDmMCTEnd', output_var='PrevDm0', unit='-',
+ end=['repEndMaps'], steps=[], all=[],
+ restrictoption=[], monthly=False, yearly=False),
+ 'PrevDmMCTState': ReportedMap(name='PrevDmMCTState', output_var='PrevDm0', unit='-',
+ end=[], steps=['repStateMaps'], all=[],
+ restrictoption=['nonInit'], monthly=False,
+ yearly=False),
'PolderLevelEnd': ReportedMap(name='PolderLevelEnd', output_var='PolderLevel',
unit='m', end=['repEndMaps'], steps=[],
all=[],
diff --git a/src/lisflood/global_modules/settings.py b/src/lisflood/global_modules/settings.py
index a4b44afe..d77cf7ba 100755
--- a/src/lisflood/global_modules/settings.py
+++ b/src/lisflood/global_modules/settings.py
@@ -242,11 +242,15 @@ def __init__(self, mask, maskmap):
self.mask_all.mask = self.flat
self.info = self.Info(mask, mask.shape, self.flat, self.flat.shape, self.mask_compressed.shape, self.mask_all)
self._in_zero = np.zeros(self.mask_compressed.shape)
+ self._in_one = np.ones(self.mask_compressed.shape)
self.maskmap = maskmap
def in_zero(self):
return self._in_zero.copy()
+ def in_one(self):
+ return self._in_one.copy()
+
def __iter__(self):
return iter(self.info)
diff --git a/src/lisflood/hydrological_modules/evapowater.py b/src/lisflood/hydrological_modules/evapowater.py
index 85757549..4f917d83 100755
--- a/src/lisflood/hydrological_modules/evapowater.py
+++ b/src/lisflood/hydrological_modules/evapowater.py
@@ -60,7 +60,8 @@ def initial(self):
maskinfo = MaskInfo.instance()
if option['openwaterevapo']:
LakeMask = loadmap('LakeMask', pcr=True)
- lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresKinematic, 5)
+ # lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresKinematic, 5)
+ lmask = ifthenelse(LakeMask != 0, self.var.LddStructuresChan, 5)
LddEva = lddrepair(lmask)
lddC = compressArray(LddEva)
inAr = decompress(np.arange(maskinfo.info.mapC[0], dtype="int32"))
@@ -132,13 +133,15 @@ def dynamic(self):
UpstreamEva = self.var.EWRef * self.var.MMtoM3 * self.var.WaterFraction
# evaporation for loop is amount of water per timestep [cu m]
# Volume of potential evaporation from water surface per time step (conversion to [m3])
- ChanMIter = self.var.ChanM3Kin.copy()
+ # ChanMIter = self.var.ChanM3Kin.copy()
+ ChanMIter = self.var.ChanM3.copy()
# for Iteration loop: First value is amount of water in the channel
- # amount of water in bankful (first line of routing)
+
ChanLeft = ChanMIter * 0.1
- # 10% of the discharge must stay in the river
+ # 10% of the volume must stay in the river
+
self.var.EvaAddM3 = MaskInfo.instance().in_zero()
- # real water consumption is set to 0
+ # real water consumption is set to 0
for NoEvaExe in range(self.var.maxNoEva):
ChanHelp = np.maximum(ChanMIter - UpstreamEva, ChanLeft)
diff --git a/src/lisflood/hydrological_modules/indicatorcalc.py b/src/lisflood/hydrological_modules/indicatorcalc.py
index b70aa0ed..a6b6b500 100644
--- a/src/lisflood/hydrological_modules/indicatorcalc.py
+++ b/src/lisflood/hydrological_modules/indicatorcalc.py
@@ -132,7 +132,8 @@ def dynamic(self):
# EXTERNAL FLOW
wreg = decompress(self.var.WUseRegionC)
# to pcraster map because of the following expression!
- self.var.RegionMonthExternalInflowM3 = compressArray(areatotal(cover(ifthen(self.var.WaterRegionInflowPoints != 0, upstream(self.var.LddStructuresKinematic,decompress(self.var.MonthDisM3))),0),wreg))
+ # self.var.RegionMonthExternalInflowM3 = compressArray(areatotal(cover(ifthen(self.var.WaterRegionInflowPoints != 0, upstream(self.var.LddStructuresKinematic,decompress(self.var.MonthDisM3))),0),wreg))
+ self.var.RegionMonthExternalInflowM3 = compressArray(areatotal(cover(ifthen(self.var.WaterRegionInflowPoints != 0, upstream(self.var.LddStructuresChan,decompress(self.var.MonthDisM3))),0),wreg))
self.var.RegionMonthAbstractionRequiredAllSourcesM3 = np.take(np.bincount(self.var.WUseRegionC,weights=self.var.MonthAbstractionRequiredAllSourcesM3),self.var.WUseRegionC)
self.var.RegionMonthAbstractionRequiredSurfaceGroundWaterM3 = np.take(np.bincount(self.var.WUseRegionC,weights=self.var.MonthAbstractionRequiredSurfaceGroundWaterM3),self.var.WUseRegionC)
self.var.RegionMonthAbstractionRequiredSurfaceWaterM3 = np.take(np.bincount(self.var.WUseRegionC,weights=self.var.MonthAbstractionRequiredSurfaceWaterM3),self.var.WUseRegionC)
diff --git a/src/lisflood/hydrological_modules/inflow.py b/src/lisflood/hydrological_modules/inflow.py
index 1678339f..5a71e557 100755
--- a/src/lisflood/hydrological_modules/inflow.py
+++ b/src/lisflood/hydrological_modules/inflow.py
@@ -47,7 +47,7 @@ def __init__(self, inflow_variable):
# --------------------------------------------------------------------------
# --------------------------------------------------------------------------
def initial(self):
- """ initial part of the inflow module
+ """ Initialization part of the inflow module.
"""
# ************************************************************
# ***** INFLOW INIT
@@ -55,14 +55,17 @@ def initial(self):
settings = LisSettings.instance()
option = settings.options
if option['inflow']:
- self.var.InflowPoints = loadmap('InflowPoints')
+ self.var.InflowPoints = loadmap('InflowPoints') #1D array size is pixels belonging to basin mask
+
self.var.QInM3Old = np.where(self.var.InflowPoints > 0, self.var.ChanQ * self.var.DtSec, 0)
+ # inflow volume for model step
# read inflow map
inflowmapprc = loadmap('InflowPoints', pcr=True)
- inflowmapnp = pcr2numpy(inflowmapprc, -9999)
+ inflowmapnp = pcr2numpy(inflowmapprc, -9999) #2D array
inflowmapnp = np.where(inflowmapnp > 0, inflowmapnp, 0)
+
# get outlets ids from outlets map
inflow_ids = np.unique(inflowmapnp)
# drop negative values (= missing data in pcraster map)
@@ -84,19 +87,20 @@ def initial(self):
warnings.warn(LisfloodWarning("Inflow point was removed ID: %d\n" % inf_id))
# substitute indexes to id in map
- self.var.InflowPointsMap = np.copy(inflowmapnp)
+ self.var.InflowPointsMap = np.copy(inflowmapnp) #np
for k, v in iteritems(id_dict):
- self.var.InflowPointsMap[inflowmapnp == k] = v
+ self.var.InflowPointsMap[inflowmapnp == k] = v #np
# convert map to pcraster format
- self.var.InflowPointsMap = numpy2pcr(Nominal, self.var.InflowPointsMap, -9999)
+ self.var.InflowPointsMap = numpy2pcr(Nominal, self.var.InflowPointsMap, -9999) #pcr
+
# Initialising cumulative output variables
# These are all needed to compute the cumulative mass balance error
- # inflow substep amount
+
def dynamic_init(self):
- """ dynamic part of the inflow module
+ """ Initialization of the dynamic part of the inflow module
init inflow before sub step routing
"""
@@ -121,6 +125,7 @@ def dynamic(self):
# Get inflow hydrograph at each inflow point [m3/s]
QIn = compressArray(QIn)
QIn[np.isnan(QIn)] = 0
+ #cmcheck - inflow
self.var.QInM3 = QIn * self.var.DtSec
# Convert to [m3] per time step
self.var.TotalQInM3 += self.var.QInM3
@@ -140,7 +145,7 @@ def dynamic_inloop(self, NoRoutingExecuted):
settings = LisSettings.instance()
option = settings.options
- if option['inflow']:
+ if option['inflow']: #cmcheck - inflow
self.var.QInDt = (self.var.QInM3Old + (NoRoutingExecuted + 1) * self.var.QDelta) * self.var.InvNoRoutSteps
# flow from inlets per sub step
diff --git a/src/lisflood/hydrological_modules/kinematic_wave_parallel.py b/src/lisflood/hydrological_modules/kinematic_wave_parallel.py
index 48f51cfd..1f4d5fd2 100755
--- a/src/lisflood/hydrological_modules/kinematic_wave_parallel.py
+++ b/src/lisflood/hydrological_modules/kinematic_wave_parallel.py
@@ -83,6 +83,8 @@ def streamLookups(flow_dir, land_mask):
'''
flow_dir[~land_mask] = 8 # exceeds number of rows of IX_ADDS
num_pixs = land_mask.sum()
+ # Create 2D array of indices of land pixels (each index is unique)
+ # The IDs are numbered from 0 to num_pixs-1
land_points = -np.ones(land_mask.shape, int)
land_points[land_mask] = np.arange(num_pixs, dtype=int)
downstream_lookup, upstream_lookup = kwpt.upDownLookups(flow_dir, np.ascontiguousarray(land_mask).astype(np.uint8), land_points, num_pixs, IX_ADDS)
@@ -120,16 +122,21 @@ def __init__(self, compressed_encoded_ldd, land_mask, alpha_channel, beta, space
self.kinematic_wave_warning_printed=False
self.flagnancheck=flagnancheck
# Parameters for the solution of the discretised Kinematic wave continuity equation
+ self.alpha_channel = alpha_channel
self.space_delta = space_delta
+ self.inv_time_delta = 1 / time_delta
self.beta = beta
self.inv_beta = 1 / beta
self.b_minus_1 = beta - 1
self.a_dx_div_dt_channel = alpha_channel * space_delta / time_delta
self.b_a_dx_div_dt_channel = beta * self.a_dx_div_dt_channel
+ self.a_dx = alpha_channel * space_delta
# If split-routing (floodplains)
if alpha_floodplains is not None:
+ self.alpha_floodplains = alpha_floodplains
self.a_dx_div_dt_floodplains = alpha_floodplains * space_delta / time_delta
self.b_a_dx_div_dt_floodplains = beta * self.a_dx_div_dt_floodplains
+ self.a_dx_floodplains = alpha_floodplains * space_delta
# Process flow direction matrix: downstream and upstream lookups, and routing orders
flow_dir = decodeFlowMatrix(rebuildFlowMatrix(compressed_encoded_ldd, land_mask))
self.downstream_lookup, self.upstream_lookup = streamLookups(flow_dir, land_mask)
@@ -139,7 +146,7 @@ def __init__(self, compressed_encoded_ldd, land_mask, alpha_channel, beta, space
def _setRoutingOrders(self):
"""Compute the kinematic wave routing order. Pixels are grouped in sets with the same order.
- Pixels in the same se are independent and can be routed in parallel. Sets must be processed in series, starting from order 0.
+ Pixels in the same set are independent and can be routed in parallel. Sets must be processed in series, starting from order 0.
Pixels are ordered topologically starting from the outlets, as in:
Liu et al. (2014), A layered approach to parallel computing for spatially distributed hydrological modeling,
Environmental Modelling & Software 51, 221-227.
@@ -152,12 +159,19 @@ def _setRoutingOrders(self):
self.pixels_ordered = self.pixels_ordered.sort_values(["order", "pixels"]).set_index("order").squeeze()
except: # FOR COMPATIBILITY WITH OLDER PANDAS VERSIONS
self.pixels_ordered = self.pixels_ordered.sort(["order", "pixels"]).set_index("order").squeeze()
+ # Ensure output is always a Series, even if only one element
+ if not isinstance(self.pixels_ordered, pd.DataFrame):
+ self.pixels_ordered = pd.Series(self.pixels_ordered)
+
order_counts = self.pixels_ordered.groupby(self.pixels_ordered.index).count()
+ # calc number of sets (orders)
stop = order_counts.cumsum()
- self.order_start_stop = np.column_stack((np.append(0, stop[:-1]), stop)).astype(int)
- self.pixels_ordered = self.pixels_ordered.values.astype(int)
+ self.order_start_stop = np.column_stack((np.append(0, stop[:-1]), stop)).astype(int)
+ # find start and end pixel in each order
+ self.pixels_ordered = self.pixels_ordered.values.astype(int)
+ # array containing all pixels organised by computation order
- def kinematicWaveRouting(self, discharge, specific_lateral_inflow, section="main_channel"):
+ def kinematicWaveRouting(self, discharge_avg, discharge, specific_lateral_inflow, section="main_channel"):
"""Kinematic wave routing: wrapper around kinematic_wave_parallel_tools.kinematicWave"""
# Lateral inflow (m3 s-1)
lateral_inflow = nx.evaluate("q * dx", local_dict={"q": specific_lateral_inflow, "dx": self.space_delta})
@@ -165,18 +179,20 @@ def kinematicWaveRouting(self, discharge, specific_lateral_inflow, section="main
if section == "main_channel":
a_dx_div_dt = self.a_dx_div_dt_channel
b_a_dx_div_dt = self.b_a_dx_div_dt_channel
+ a_dx = self.a_dx
elif section == "floodplains":
a_dx_div_dt = self.a_dx_div_dt_floodplains
b_a_dx_div_dt = self.b_a_dx_div_dt_floodplains
+ a_dx = self.a_dx_floodplains
else:
raise Exception("The section parameter must be either 'main_channel' or 'floodplain'!")
# Constant term in f(x) evaluation for Newton-Raphson method
local_dict = {"a_dx_div_dt": a_dx_div_dt, "Qold": discharge, "b": self.beta, "lateral_inflow": lateral_inflow}
constant = nx.evaluate("a_dx_div_dt * Qold ** b + lateral_inflow", local_dict=local_dict)
# Solve the Kinematic wave equation
- kwpt.kinematicRouting(discharge, lateral_inflow, constant, self.upstream_lookup,\
+ kwpt.kinematicRouting(discharge_avg, discharge, lateral_inflow, constant, self.upstream_lookup,\
self.num_upstream_pixels, self.pixels_ordered, self.order_start_stop,\
- self.beta, self.inv_beta, self.b_minus_1, a_dx_div_dt, b_a_dx_div_dt)
+ self.inv_time_delta,self.beta, self.inv_beta, self.b_minus_1, a_dx_div_dt, b_a_dx_div_dt, a_dx)
if self.flagnancheck:
if self.kinematic_wave_warning_printed==False:
if np.all(np.isfinite(discharge))==False:
diff --git a/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py b/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py
index 8646cefa..bfd95d8e 100755
--- a/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py
+++ b/src/lisflood/hydrological_modules/kinematic_wave_parallel_tools.py
@@ -13,7 +13,7 @@
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the Licence for the specific language governing permissions and limitations under the Licence.
-Collection of functions called by the Python module kinematic_wave_parallel.py.
+Collection of functions called by the Python module kinematic_wave_parallel.py.
"""
from math import fabs
@@ -32,65 +32,141 @@
# ROUTING FUNCTIONS
# -------------------------------------------------------------------------------------------------
@njit(parallel=True, fastmath=False, cache=True)
-def kinematicRouting(discharge, lateral_inflow, constant, upstream_lookup,\
- num_upstream_pixels, ordered_pixels, start_stop, beta, inv_beta,\
- b_minus_1, a_dx_div_dt, b_a_dx_div_dt):
- """"""
+def kinematicRouting(discharge_avg, discharge, lateral_inflow, constant, upstream_lookup,\
+ num_upstream_pixels, ordered_pixels, start_stop, inv_time_delta, beta, inv_beta,\
+ b_minus_1, a_dx_div_dt, b_a_dx_div_dt, a_dx):
+ """
+ This function performs the kinematic wave routing algorithm to simulate the movement of water through a network of interconnected channels.
+ :param discharge_avg:
+ :param discharge:
+ :param lateral_inflow:
+ :param constant:
+ :param upstream_lookup:
+ :param num_upstream_pixels:
+ :param ordered_pixels:
+ :param start_stop:
+ :param inv_time_delta: 1/DtRouting
+ :param beta:
+ :param inv_beta:
+ :param b_minus_1:
+ :param a_dx_div_dt:
+ :param b_a_dx_div_dt:
+ :param a_dx: ChannelAlpha * ChanLength
+ :return:
+ """
num_orders = start_stop.shape[0]
- # Iterate through routing orders (sets of pixels for which the kinemativc wave can be solved independently and thus in parallel)
+ # Iterate through each routing order (sets of pixels for which the kinemativc wave can be solved independently and thus in parallel)
for order in range(num_orders):
first = start_stop[order,0]
last = start_stop[order,1]
- for index in prange(first, last):
- solve1Pixel(ordered_pixels[index], discharge, lateral_inflow, constant, upstream_lookup,\
- num_upstream_pixels, a_dx_div_dt, b_a_dx_div_dt, beta, inv_beta, b_minus_1)
+ # Iterate through each pixel in the current order in parallel
+ for index in prange(first, last): # this is a parallel loop
+ # Solve the kinematic wave for the current pixel
+ pix = ordered_pixels[index]
+ discharge_before = discharge[pix]
+ if solve1Pixel(pix, discharge, constant[pix], upstream_lookup,\
+ num_upstream_pixels, a_dx_div_dt[pix], b_a_dx_div_dt[pix], beta, inv_beta, b_minus_1):
+ solve1PixelAvg(pix, discharge_avg, discharge_before, discharge[pix], lateral_inflow[pix], upstream_lookup,\
+ num_upstream_pixels, inv_time_delta, beta, a_dx[pix])
+
@njit(nogil=True, fastmath=False, cache=True)
-def solve1Pixel(pix, discharge, lateral_inflow, constant,\
- upstream_lookup, num_upstream_pixels, a_dx_div_dt,\
- b_a_dx_div_dt, beta, inv_beta, b_minus_1):
- """"""
- count = 0
- previous_estimate = -1.0
- upstream_inflow = 0.0
+def solve1Pixel(
+ pix, discharge, constant,
+ upstream_lookup, num_upstream_pixels, a_dx_div_dt,
+ b_a_dx_div_dt, beta, inv_beta, b_minus_1):
+ """
+ Te Chow et al. 1988 - Applied Hydrology - Chapter 9.6
+ :param pix:
+ :param discharge: instantaneous outflow discharge
+ :param constant:
+ :param upstream_lookup:
+ :param num_upstream_pixels:
+ :param a_dx_div_dt:
+ :param b_a_dx_div_dt:
+ :param beta:
+ :param inv_beta:
+ :param b_minus_1:
+ :return:
+ """
+
# Inflow from upstream pixels
+ upstream_inflow = constant # upstream_inflow + alpha*dx/dt*Qold**beta + dx*specific_lateral_inflow
for ups_ix in range(num_upstream_pixels[pix]):
upstream_inflow += discharge[upstream_lookup[pix,ups_ix]]
- const_plus_ups_infl = upstream_inflow + constant[pix] # upstream_inflow + alpha*dx/dt*Qold**beta + dx*specific_lateral_inflow
+
# If old discharge, upstream inflow and lateral inflow are below accuracy: set discharge to 0 and exit
- if const_plus_ups_infl <= NEWTON_TOL:
+ if upstream_inflow <= NEWTON_TOL:
discharge[pix] = 0
- return
+ return False
+
# Initial discharge guess using analytically derived boundary values
- a_cpui_pow_b_m_1 = b_a_dx_div_dt[pix] * const_plus_ups_infl**b_minus_1
+ a_cpui_pow_b_m_1 = b_a_dx_div_dt * upstream_inflow**b_minus_1
if a_cpui_pow_b_m_1 <= 1:
- secant_bound = const_plus_ups_infl / (1 + a_cpui_pow_b_m_1)
+ secant_bound = upstream_inflow / (1 + a_cpui_pow_b_m_1)
else:
- secant_bound = const_plus_ups_infl / (1 + a_cpui_pow_b_m_1**inv_beta)
- other_bound = ((const_plus_ups_infl - secant_bound) / a_dx_div_dt[pix])**inv_beta
- discharge[pix] = (secant_bound + other_bound) / 2
- error = closureError(discharge[pix], const_plus_ups_infl, a_dx_div_dt[pix], beta)
+ secant_bound = upstream_inflow / (1 + a_cpui_pow_b_m_1**inv_beta)
+ other_bound = ((upstream_inflow - secant_bound) / a_dx_div_dt)**inv_beta
+ dis = (secant_bound + other_bound) / 2
+ error = dis + a_dx_div_dt * dis**beta - upstream_inflow
+
# Iterations
- while fabs(error) > NEWTON_TOL and discharge[pix] != previous_estimate and count < MAX_ITERS: # is previous_estimate useful?
- previous_estimate = discharge[pix]
- discharge[pix] -= error / (1 + b_a_dx_div_dt[pix] * discharge[pix]**b_minus_1)
- discharge[pix] = max(discharge[pix], NEWTON_TOL)
- error = closureError(discharge[pix], const_plus_ups_infl, a_dx_div_dt[pix], beta)
+ count = 0
+ previous_estimate = -1.0
+ while count < MAX_ITERS:
+ if fabs(error) <= NEWTON_TOL or dis == previous_estimate:
+ break
+ previous_estimate = dis
+ factor = (1 + b_a_dx_div_dt * dis**b_minus_1)
+ dis = max(dis - error / factor, NEWTON_TOL)
+ error = dis + a_dx_div_dt * dis**beta - upstream_inflow
count += 1
+
# If iterations converge to NEWTON_TOL, set value to 0
- if discharge[pix] == NEWTON_TOL:
+ if dis == NEWTON_TOL:
discharge[pix] = 0
- # to simulate inf or nan: discharge[pix] = 1.0/0.0
- # with gil:
- # got_valid_value = np.isfinite(discharge[pix])
- # if got_valid_value==False:
- # print(str(discharge[pix]) + ' found at pix:' + str(pix))
+ return False
+
+ discharge[pix] = dis
+ return True
+
@njit(nogil=True, fastmath=False, cache=True)
-def closureError(discharge, upper_bound, a_dx_div_dt, beta):
- """"""
- return discharge + a_dx_div_dt * discharge**beta - upper_bound
+def solve1PixelAvg(
+ pix, discharge_avg, discharge_before, discharge_after, lateral_inflow,
+ upstream_lookup, num_upstream_pixels,
+ inv_time_delta, beta, a_dx):
+ """
+ :param pix:
+ :param discharge_avg: average outflow discharge
+ :param discharge_before: instantaneous outflow discharge before solve1Pixel computation
+ :param discharge_after: instantaneous outflow discharge after solve1Pixel computation
+ :param lateral_inflow:
+ :param upstream_lookup:
+ :param num_upstream_pixels:
+ :param inv_time_delta: 1/DtRouting
+ :param beta:
+ :param a_dx: ChannelAlpha * ChanLength
+ :return:
+ """
+ upstream_inflow_avg = 0.0
+
+ # Inflow from upstream pixels
+ for ups_ix in range(num_upstream_pixels[pix]):
+ upstream_inflow_avg += discharge_avg[upstream_lookup[pix, ups_ix]]
+
+ # volume of water in channel at beginning of computation step
+ channel_volume_start = discharge_before**beta
+
+ # volume of water in channel at end of computation step
+ channel_volume_end = discharge_after**beta
+ # integration on control volume to calc average outflow (channel water mass balance)
+ discharge_avg[pix] = upstream_inflow_avg + lateral_inflow + a_dx*(channel_volume_start - channel_volume_end) * inv_time_delta
+
+ # avoid negative average discharge
+ if discharge_avg[pix] < 0:
+ discharge_avg[pix] = 0
# -------------------------------------------------------------------------------------------------
@@ -99,36 +175,113 @@ def closureError(discharge, upper_bound, a_dx_div_dt, beta):
@njit(nogil=True, fastmath=False, cache=True)
def updateUpstreamRouted(pixels_routed, downstream_lookup,\
upstream_lookup, upstream_routed):
- """Called by kinematic_wave_parallel.orderKinematicRouting to implement greedy pixel ordering starting from headwaters."""
+ """
+ Called by kinematic_wave_parallel.orderKinematicRouting to implement greedy pixel ordering starting from headwaters.
+ Update the boolean array upstream_routed to mark the upstream pixels that have already been routed.
+ The function loops over each pixel in pixels_routed, identifies its downstream pixel, and sets the corresponding
+ element of upstream_routed to True. This is achieved by finding the index of the current pixel in the upstream_lookup
+ array of the downstream pixel and setting the corresponding element of upstream_routed to True.
+
+ Parameters:
+ -----------
+ pixels_routed : array-like
+ 1D array of integers representing pixels that have already been routed.
+ downstream_lookup : array-like
+ 1D array of integers representing the downstream pixel index for each pixel.
+ upstream_lookup : array-like
+ 2D array of integers representing the upstream pixel index for each downstream pixel and each of its upstream pixels.
+ upstream_routed : array-like
+ 2D boolean array representing whether each upstream pixel has already been routed.
+
+ Returns:
+ --------
+ None
+ :param pixels_routed:
+ :param downstream_lookup:
+ :param upstream_lookup:
+ :param upstream_routed:
+ :return:
+ """
for pix in pixels_routed:
+ # Find the downstream pixel index for the current pixel
downs_pix = downstream_lookup[pix]
if downs_pix != -1:
+ # Find the index of the current pixel among its upstream pixels
ups_index = 0
while upstream_lookup[downs_pix,ups_index] != pix:
ups_index += 1
+ # Mark the upstream pixel as routed
upstream_routed[downs_pix,ups_index] = True
+
@njit(nogil=True, fastmath=False, cache=True)
def upDownLookups(flow_d8, land_mask, land_points, num_pixs, ix_adds):
- '''Called by catchment.streamLookups'''
+ """
+ Called by catchment.streamLookups
+ Create lookups of pixels upstream and downstream of each pixel in a catchment.
+ These arrays are used to determine the flow paths of water from a source pixel to all the pixels downstream and upstream of it, respectively.
+ The function loops through each pixel in the dataset and checks whether it has a downstream neighbor in the direction of its flow_d8 value.
+ If it does, it calculates the ID of the downstream pixel using the ix_adds array and checks whether the downstream pixel is also a land pixel
+ (i.e. it is not water or outside the dataset). If both pixels are land pixels, the function updates the downstream_lookup and upstream_lookup
+ arrays with the IDs of the upstream and downstream pixels, respectively. The function keeps track of the number of upstream neighbors for
+ each downstream pixel using the ups_count array.
+
+ Parameters:
+ -----------
+ flow_d8 : long[:,::1]
+ 2D array of D8 flow direction codes
+ land_mask : unsigned char[:,::1]
+ 2D array of 0/1 indicating whether each pixel is land or not
+ land_points : long[:,::1]
+ 2D array of indices of land pixels (each index is unique).
+ The IDs are numbered from 0 to num_pixs - 1.
+ num_pixs : long
+ Number of pixels in the catchment
+ ix_adds : long[:,::1]
+ 2D array of pixel coordinates to add to a source pixel to get destination pixel
+ in a given direction, given as (row, col) pairs
+
+ Returns:
+ --------
+ downstream_lookup : long[::1]
+ 1D array of indices of downstream pixels for each pixel in the catchment
+ upstream_lookup : long[:,::1]
+ 2D array of indices of upstream pixels for each pixel in the catchment
+ :param flow_d8:
+ :param land_mask:
+ :param land_points:
+ :param num_pixs:
+ :param ix_adds:
+ :return:
+ """
downstream_lookup = -np.ones(num_pixs)
upstream_lookup = -np.ones((num_pixs, 8))
ups_count = np.zeros(num_pixs, np.uintc)
num_row = flow_d8.shape[0]
num_col = flow_d8.shape[1]
+ # Loop over every pixel in the catchment
for src_row in range(num_row):
for src_col in range(num_col):
+ # If the flow direction is less than 8 (i.e. the pixel is part of the catchment)
if flow_d8[src_row,src_col] < 8:
+ # Compute the destination pixel coordinates based on flow direction
dst_row = src_row + ix_adds[flow_d8[src_row,src_col],0]
dst_col = src_col + ix_adds[flow_d8[src_row,src_col],1]
+ # Check if the destination pixel is within the boundaries of the land_mask
+ # and is itself a land pixel
if dst_row != -1 and dst_col != -1 and dst_row != num_row and dst_col != num_col and land_mask[dst_row,dst_col]:
+ # Get the indices of the source and destination pixels in the land_points array
ups_p = land_points[src_row,src_col]
downs_p = land_points[dst_row,dst_col]
+ # Assign the downstream pixel index to the upstream_lookup array for the source pixel
downstream_lookup[ups_p] = downs_p
+ # Assign the upstream pixel index to the upstream_lookup array for the destination pixel
+ # and increment the counter for the number of upstream pixels
upstream_lookup[downs_p,ups_count[downs_p]] = ups_p
ups_count[downs_p] += 1
- return np.asarray(downstream_lookup), np.asarray(upstream_lookup)
+ # Return the downstream and upstream lookup arrays
+ return np.asarray(downstream_lookup), np.asarray(upstream_lookup)
# -------------------------------------------------------------------------------------------------
@@ -136,11 +289,37 @@ def upDownLookups(flow_d8, land_mask, land_points, num_pixs, ix_adds):
# -------------------------------------------------------------------------------------------------
@njit(nogil=True, fastmath=False, cache=True)
def immediateUpstreamInflow(discharge, upstream_lookup, num_upstream_pixels):
- """"""
+ """Calculate the immediate upstream inflow for each pixel in a river network.
+ The implementation of the function is as follows:
+ A new 1D NumPy array of type double called inflow is created with the same length as discharge, and is initialized to all zeros.
+ A loop over all the pixels in the river network is started, using the range() function and the num_pixels variable (which is set to the length of discharge).
+ Inside the loop, another loop over the upstream pixels for the current pixel is started, using the range() function and the num_upstream_pixels[pix] variable.
+ Inside the inner loop, the inflow for the current pixel is updated by adding the discharge of the upstream pixel, which is retrieved using the upstream_lookup array.
+ Finally, the inflow array is returned as a NumPy array.
+
+ Args:
+ discharge: A 1D NumPy array of type double that represents the discharge (flow rate) for each pixel in the river network.
+ upstream_lookup: A 2D NumPy array of type long that represents the upstream routing of each pixel in the river network.
+ num_upstream_pixels: A 1D NumPy array of type long that represents the number of upstream pixels for each pixel in the river network.
+
+ Returns:
+ A 1D NumPy array of type double that represents the immediate upstream inflow for each pixel in the river network.
+
+ :param discharge:
+ :param upstream_lookup:
+ :param num_upstream_pixels:
+ :return:
+ """
num_pixels = discharge.size
inflow = np.zeros(num_pixels)
+ # Loop over all the pixels in the river network
for pix in range(num_pixels):
+
+ # Loop over the upstream pixels for the current pixel
for ups_ix in range(num_upstream_pixels[pix]):
+
+ # Update the inflow for the current pixel by adding the discharge of the upstream pixel
inflow[pix] += discharge[upstream_lookup[pix,ups_ix]]
- return np.asarray(inflow)
+ # Return the inflow array as a NumPy array
+ return np.asarray(inflow)
diff --git a/src/lisflood/hydrological_modules/lakes.py b/src/lisflood/hydrological_modules/lakes.py
index fc4a4b7a..1b3fe103 100755
--- a/src/lisflood/hydrological_modules/lakes.py
+++ b/src/lisflood/hydrological_modules/lakes.py
@@ -80,6 +80,9 @@ def initial(self):
self.var.IsStructureKinematic = np.where(LakeSitesC > 0, np.bool8(1), self.var.IsStructureKinematic)
# Add lake locations to structures map (used to modify LddKinematic
# and to calculate LddStructuresKinematic)
+ self.var.IsStructureChan = np.where(LakeSitesC > 0, np.bool8(1), self.var.IsStructureChan)
+ # Add lake locations to structures map (used to modify LddChan
+ # and to calculate LddStructuresChan)
# PCRaster part
# -----------------------
@@ -87,7 +90,8 @@ def initial(self):
LakeSitePcr = pcraster.ifthen((pcraster.defined(LakeSitePcr) & pcraster.boolean(decompress(self.var.IsChannel))), LakeSitePcr)
IsStructureLake = pcraster.boolean(LakeSitePcr)
# additional structure map only for lakes to calculate water balance
- self.var.IsUpsOfStructureLake = pcraster.downstream(self.var.LddKinematic, pcraster.cover(IsStructureLake, 0))
+ # self.var.IsUpsOfStructureLake = pcraster.downstream(self.var.LddKinematic, pcraster.cover(IsStructureLake, 0))
+ self.var.IsUpsOfStructureLake = pcraster.downstream(self.var.LddChan, pcraster.cover(IsStructureLake, 0))
# Get all pixels just upstream of lakes
# -----------------------
diff --git a/src/lisflood/hydrological_modules/mct.py b/src/lisflood/hydrological_modules/mct.py
new file mode 100644
index 00000000..56c90484
--- /dev/null
+++ b/src/lisflood/hydrological_modules/mct.py
@@ -0,0 +1,549 @@
+import numpy as np
+import pandas as pd
+from numba import njit, prange
+
+from .kinematic_wave_parallel import rebuildFlowMatrix, decodeFlowMatrix, streamLookups, topoDistFromSea
+
+
+class MCTWave:
+ """Build pixels loop for MCT channels"""
+
+ def __init__(
+ self,
+ compressed_encoded_ldd, # MCT Ldd
+ land_mask, # MCT mask
+ ChanLength, # Channel length
+ ChanGrad, # Channel riverbed slope
+ ChanBottomWidth, # Riverbed bottom width
+ ChanManMCT, # MCT pixels Mannings' coefficient
+ ChanSdXdY, # Riverbed side slope
+ dt, # computation time step for routing [s]
+ river_router, # class
+ mapping_mct, # MCT pixels mapping
+ ):
+
+ # Process flow direction matrix: downstream and upstream lookups, and routing orders
+ flow_dir = decodeFlowMatrix(rebuildFlowMatrix(compressed_encoded_ldd, land_mask))
+ self.downstream_lookup, self.upstream_lookup = streamLookups(flow_dir, land_mask)
+ self.num_upstream_pixels = (self.upstream_lookup != -1).sum(1).astype(int) # astype for cython import in windows (to avoid 'long long' buffer dtype mismatch)
+ # Routing order: decompose domain into batches; within each batch, pixels can be routed in parallel
+ self._setMCTRoutingOrders()
+
+ self.ChanLength = ChanLength
+ self.ChanGrad = ChanGrad
+ self.ChanBottomWidth = ChanBottomWidth
+ self.ChanManMCT = ChanManMCT
+ self.ChanSdXdY = ChanSdXdY
+ self.dt = dt
+ self.river_router = river_router
+ self.mapping_mct = mapping_mct
+
+
+ def _setMCTRoutingOrders(self):
+ """Compute the MCT wave routing order. Pixels are grouped in sets with the same order.
+ Pixels in the same set are independent and can be routed in parallel. Sets must be processed in series, starting from order 0.
+ Sets contain MCT pixels only, but pixels can receive contribution from upstream kinematic pixels.
+ Pixels are ordered topologically starting from the outlets, as in:
+ Liu et al. (2014), A layered approach to parallel computing for spatially distributed hydrological modeling,
+ Environmental Modelling & Software 51, 221-227.
+ Order MAX is given to pixels with no downstream relations (outlets); order MAX-1 is given to
+ pixels whose downstream pixels are all of order MAX; and so on."""
+
+ ocean_topo_distance = topoDistFromSea(self.downstream_lookup, self.upstream_lookup)
+ routing_order = ocean_topo_distance.max() - ocean_topo_distance
+ self.pixels_ordered = pd.DataFrame({"pixels": np.arange(routing_order.size), "order": routing_order})
+ try:
+ self.pixels_ordered = self.pixels_ordered.sort_values(["order", "pixels"]).set_index("order").squeeze()
+ except: # FOR COMPATIBILITY WITH OLDER PANDAS VERSIONS
+ self.pixels_ordered = self.pixels_ordered.sort(["order", "pixels"]).set_index("order").squeeze()
+ # Ensure output is always a Series, even if only one element
+ if not isinstance(self.pixels_ordered, pd.Series):
+ self.pixels_ordered = pd.Series(self.pixels_ordered)
+
+ order_counts = self.pixels_ordered.groupby(self.pixels_ordered.index).count()
+ # calc number of sets (orders)
+ stop = order_counts.cumsum()
+ self.order_start_stop = np.column_stack((np.append(0, stop[:-1]), stop)).astype(int) # astype for cython import in windows (see above)
+ # find start and end pixel in each order
+ self.pixels_ordered = self.pixels_ordered.values.astype(int) # astype for cython import in windows (see above)
+ # array containing all pixels organised by computation order
+
+ def routing(
+ self,
+ ChanQ_0, # -> used to calc q00
+ ChanM3_0, # V00
+ SideflowChanMCT, # sideflow for MCT pixels [m3/s]
+ # THESE ARE USED AS INPUTS AND OUTPUTS
+ ChanQ, # -> in input: used to calc q01; in output: q11 outflow (x+dx) at time t+dt (instant)
+ ChanQAvgDt, # -> in input: used to calculate q0m; in output: q1m outflow (x+dx) at time t+dt (average)
+ PrevCm0, # Courant number at the end of previous routing step t
+ PrevDm0, # Reynolds number at the end of previous routing step t
+ ChanM3, # Channel storage volume. In input: at time t V00; in output: at time t+dt V11
+ ):
+
+
+ mct_routing(
+ self.ChanLength, # Channel length
+ self.ChanGrad, # Channel riverbed slope
+ self.ChanBottomWidth, # Riverbed bottom width
+ self.ChanManMCT, # MCT pixels Mannings' coefficient
+ self.ChanSdXdY, # Riverbed side slope
+ self.dt, # computation time step for routing [s]
+ self.order_start_stop, # index of first and last pixel in the set
+ self.pixels_ordered, # list of pixel indexes in order of routing
+ self.river_router.upstream_lookup, # indexes of upstream contributing pixels
+ self.river_router.num_upstream_pixels, # number of upstream contributing pixels
+ self.mapping_mct, # MCT pixels mapping
+ ChanQ_0, # -> used to calc q00
+ ChanM3_0, # V00
+ SideflowChanMCT, # sideflow for MCT pixels [m3/s]
+ # THESE ARE USED AS INPUTS AND OUTPUTS
+ ChanQ, # -> in input: used to calc q01; in output: q11 outflow (x+dx) at time t+dt (instant)
+ ChanQAvgDt, # -> in input: used to calculate q0m; in output: q1m outflow (x+dx) at time t+dt (average)
+ PrevCm0, # Courant number in input: at time t; in output: at time t+dt
+ PrevDm0, # Reynolds number in input: at time t; in output: at time t+dt
+ ChanM3, # V11 as output
+ )
+
+
+@njit(parallel=True, fastmath=False, cache=True)
+def mct_routing(
+ # static inputs (not changing between time steps)
+ ChanLength, # Channel length
+ ChanGrad, # Channel riverbed slope
+ ChanBottomWidth, # Riverbed bottom width
+ ChanManMCT, # MCT pixels Mannings' coefficient
+ ChanSdXdY, # Riverbed side slope
+ dt, # computation time step for channel [s]
+ mct_order_start_stop, # index of first and last pixel in the set
+ mct_pixels_ordered, # list of pixel indexes in order of routing
+ upstream_lookup, # indexes of upstream contributing pixels
+ num_upstream_pixels, # number of upstream contributing pixels
+ mapping_mct, # MCT pixels mapping
+ # dynamic inputs
+ ChanQ_0, # -> used to calc q00
+ ChanM3_0, # V00
+ SideflowChanMCT, # sideflow for MCT pixels [m3/s]
+ # THESE ARE USED AS INPUTS AND OUTPUTS
+ ChanQ, # -> in input: used to calc q01; in output: q11 outflow (x+dx) at time t+dt (instant)
+ ChanQAvgDt, # -> in input: used to calculate q0m; in output: q1m outflow (x+dx) at time t+dt (average)
+ PrevCm0, # Courant number in input: at time t; in output: at time t+dt
+ PrevDm0, # Reynolds number in input: at time t; in output: at time t+dt
+ ChanM3, # V11 as output
+):
+ """This function implements Muskingum-Cunge-Todini routing method
+ MCT routing is calculated on MCT pixels only but gets inflow from both Kinematic/Split and MCT upstream pixels.
+ Function compress_mct is used to compress arrays with all river channel pixels to arrays containing MCT pixels only.
+ References:
+ Todini, E. (2007). A mass conservative and water storage consistent variable parameter Muskingum-Cunge approach. Hydrol. Earth Syst. Sci.
+ (Chapter 5)
+ Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030
+ (Appendix B)
+
+ :return:
+ ChanQ, # q11 outflow (x+dx) at time t+dt (instant)
+ ChanQAvgDt, # q1m outflow (x+dx) at time t+dt (average)
+ PrevCm0, # Courant number at time t+dt
+ PrevDm0, # Reynolds number at time t+dt
+ ChanM3, # V11 channel storage volume at t+dt (instant)
+ """
+
+ num_orders = mct_order_start_stop.shape[0]
+
+ # loop on orders
+ for order in range(num_orders):
+ first = mct_order_start_stop[order, 0]
+ last = mct_order_start_stop[order, 1]
+ # loop on pixels in the order
+ for index in prange(first, last): # this is a parallel loop
+ # get MCT pixel id in the MCT LDD
+ mctpix = mct_pixels_ordered[index]
+ # Find the corresponding pixel id in the full LDD
+ kinpix = mapping_mct[mctpix]
+ # Find id of upstream contributing pixels (from full LDD)
+ upstream_pixels = upstream_lookup[kinpix]
+
+ # get inflow from upstream pixels for current (t+dt) and previous steps (t) for the MCT pixel
+ q00 = 0.0
+ q0m = 0.0
+ q01 = 0.0
+ for ups_ix in range(num_upstream_pixels[kinpix]):
+ ups_pix = upstream_pixels[ups_ix] # upstream pixel id
+ q00 += ChanQ_0[ups_pix] # Inflow (x) to the pixel at previous step t (instant)
+ q0m += ChanQAvgDt[ups_pix] # Average inflow (x) to the pixel at previous step t (average)
+ q01 += ChanQ[ups_pix] # Inflow (x) at current step t+dt (instant)
+
+ # get outflow from the pixel at previous step t
+ q10 = ChanQ_0[kinpix] # Outflow (x+dx) from the pixel at previous step t (instant)
+
+ V00 = ChanM3_0[kinpix] # Channel storage volume at the end of previous step t (instant)
+
+ Cm0 = PrevCm0[kinpix] # Courant number at the end of previous step t
+ Dm0 = PrevDm0[kinpix] # Reynolds number at the end of previous step t
+
+ ql = SideflowChanMCT[kinpix] # Sideflow during step dt
+
+ # static data
+ xpix = ChanLength[kinpix] # Channel length
+ s0 = ChanGrad[kinpix] # Channel riverbed slope
+ Balv = ChanBottomWidth[kinpix] # Riverbed bottom width
+ Nalv = ChanManMCT[kinpix] # MCT pixels Mannings' coefficient
+ ANalv = np.arctan(1 / ChanSdXdY[kinpix]) # angle of the riverbed side [rad]
+
+ # calling MCT function for single cell idpx
+ q11, q1m, V11, Cm1, Dm1 = MCTRouting_single(
+ V00, q10, q01, q00, ql, q0m, Cm0, Dm0,
+ dt, xpix, s0, Balv, ANalv, Nalv)
+
+ ChanQ[kinpix] = q11 # Outflow (x+dx) at the end of routing step t+dt (instant)
+ ChanQAvgDt[kinpix] = q1m # Average outflow (x+dx) at the end of routing step t+dt (average)
+ ChanM3[kinpix] = V11 # Channel storage at the end of routing step t+dt (instant)
+ PrevCm0[kinpix] = Cm1 # Courant number at the end of routing step t+dt (instant)
+ PrevDm0[kinpix] = Dm1 # Reynolds number at the end of routing step t+dt (instant)
+
+
+@njit(nogil=True, fastmath=False, cache=True)
+def MCTRouting_single(
+ V00, q10, q01, q00, ql, q0mm, Cm0, Dm0, dt, xpix, s0, Balv, ANalv, Nalv
+):
+ """
+ This function implements Muskingum-Cunge-Todini routing method for a single pixel.
+ References:
+ Todini, E. (2007). A mass conservative and water storage consistent variable parameter Muskingum-Cunge approach. Hydrol. Earth Syst. Sci.
+ (Chapter 5)
+ Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030
+ (Appendix B)
+
+ :param V00: channel storage volume at t (beginning of computation step)
+ :param q10: O(t) - outflow (x+dx) at time t
+ :param q01: I(t+dt) - inflow (x) at time t+dt (end of computation step)
+ :param q00: I(t) - inflow (x) at time t
+ :param ql: lateral flow over time dt [m3/s]
+ :param q0mm: I - average inflow during step dt
+ :param Cm0: Courant number at time t
+ :param Dm0: Reynolds number at time t
+ :param dt: time interval step
+ :param xpix: channel length
+ :param s0: channel slope
+ :param Balv: channel bankfull width
+ :param ANalv: angle of the riverbed side [rad]
+ :param Nalv: channel Manning roughness coefficient
+ :return:
+ q11: Outflow (x+dx) at O(t+dt)
+ V11: channel storage volume at t+dt
+ Cm1: Courant number at t+1 for state file
+ Dm1: Reynolds number at t+1 for state file
+ """
+
+ eps = 1e-06
+
+ # Calc O' first guess for the outflow at time t+dt
+ # O'(t+dt)=O(t)+(I(t+dt)-I(t))
+ q11 = q10 + (q01 - q00)
+
+ # check for negative and zero discharge values
+ # zero outflow is not allowed
+ if q11 < eps: # cmcheck <=0
+ q11 = eps
+
+ # calc reference discharge at time t
+ # qm0 = (I(t)+O(t))/2
+ # qm0 = (q00 + q10) / 2.
+
+ # Calc O(t+dt)=q11 at time t+dt using MCT equations
+ for i in range(2): # repeat 2 times for accuracy
+
+ # reference I discharge at x=0
+ qmx0 = (q00 + q01) / 2.0
+ if qmx0 < eps: # cmcheck ==0
+ qmx0 = eps
+ hmx0 = hoq(qmx0, s0, Balv, ANalv, Nalv)
+
+ # reference O discharge at x=1
+ qmx1 = (q10 + q11) / 2.0
+ if qmx1 < eps: # cmcheck ==0
+ qmx1 = eps
+ hmx1 = hoq(qmx1, s0, Balv, ANalv, Nalv)
+
+ # Calc riverbed slope correction factor
+ cor = 1 - (1 / s0 * (hmx1 - hmx0) / xpix)
+ sfx = s0 * cor
+ if sfx < (0.8 * s0):
+ sfx = 0.8 * s0 # In case of instability raise from 0.5 to 0.8
+
+ # Calc reference discharge time t+dt
+ # Q(t+dt)=(I(t+dt)+O'(t+dt))/2
+ qm1 = (q01 + q11) / 2.0
+ # cm
+ if qm1 < eps: # cmcheck ==0
+ qm1 = eps
+ # cm
+ hm1 = hoq(qm1, s0, Balv, ANalv, Nalv)
+ dummy, Ax1, Bx1, Px1, ck1 = qoh(hm1, s0, Balv, ANalv, Nalv)
+ if ck1 <= eps:
+ ck1 = eps
+
+ # Calc correcting factor Beta at time t+dt
+ Beta1 = ck1 / (qm1 / Ax1)
+ # calc corrected cell Reynolds number at time t+dt
+ Dm1 = qm1 / (sfx * ck1 * Bx1 * xpix) / Beta1
+ # corrected Courant number at time t+dt
+ Cm1 = ck1 * dt / xpix / Beta1
+
+ # Calc MCT parameters
+ den = 1 + Cm1 + Dm1
+ c1 = (-1 + Cm1 + Dm1) / den
+ c2 = (1 + Cm0 - Dm0) / den * (Cm1 / Cm0)
+ c3 = (1 - Cm0 + Dm0) / den * (Cm1 / Cm0)
+ c4 = (2 * Cm1) / den
+
+ # cmcheck
+ # Calc outflow q11 at time t+1
+ # Mass balance equation without lateral flow
+ # q11 = c1 * q01 + c2 * q00 + c3 * q10
+ # Mass balance equation that takes into consideration the lateral flow
+ q11 = c1 * q01 + c2 * q00 + c3 * q10 + c4 * ql
+
+ if q11 < eps: # cmcheck <=0
+ q11 = eps
+
+ #### end of for loop
+
+ # # cmcheck
+ # calc_t = xpix / ck1
+ # if calc_t < dt:
+ # print('xpix/ck1 < dt')
+
+ # k1 = dt / Cm1
+ # x1 = (1. - Dm1) / 2.
+
+ # Calc the corrected mass-conservative expression for the reach segment storage at time t+dt
+ # The lateral inflow ql is only explicitly accounted for in the mass balance equation, while it is not in the equation expressing
+ # the storage as a weighted average of inflow and outflow.The rationale of this approach lies in the fact that the outflow
+ # of the reach implicitly takes the effect of the lateral inflow into account.
+ if q11 == 0:
+ V11 = V00 + (q00 + q01 - q10 - q11) * dt / 2
+ else:
+ V11 = (1 - Dm1) * dt / (2 * Cm1) * q01 + (1 + Dm1) * dt / (2 * Cm1) * q11
+ # V11 = k1 * (x1 * q01 + (1. - x1) * q11) # MUST be the same as above!
+
+ ### calc integration on the control volume (pixel)
+ # calc average discharge outflow q1m for MCT channels during routing sub step dt
+ # Calculate average outflow using water balance for MCT channel grid cell over sub-routing step
+ q1mm = q0mm + ql + (V00 - V11) / dt
+
+ # cmcheck
+ # q1m cannot be smaller than eps or it will cause instability
+ if q1mm < eps:
+ q1mm = eps
+ V11 = V00 + (q0mm + ql - q1mm) * dt
+
+ # q11 Outflow at O(t+dt)
+ # q1m average outflow in time dt
+ # V11 water volume at t+dt
+ # Cm1 Courant number at t+dt for state files
+ # Dm1 Reynolds numbers at t+dt for state files
+ return q11, q1mm, V11, Cm1, Dm1
+
+
+@njit(nogil=True, fastmath=False, cache=True)
+def hoq(q, s0, Balv, ANalv, Nalv):
+ """Water depth h from discharge q.
+ Given a generic cross-section (rectangular, triangular or trapezoidal) and a steady-state discharge q=Q*, it computes
+ water depth (y), wet contour (Bx), wet area (Ax) and wave celerity (cel) using Newton-Raphson method.
+ Reference:
+ Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method.
+ Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030
+
+ :param q: steady-state river discharge [m3/s]
+ :param s0: river bed slope (tan B)
+ :param Balv : width of the riverbed [m]
+ :param ChanSdXdY : slope dx/dy of riverbed side
+ :param ANalv : angle of the riverbed side [rad]
+ :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3]
+
+ :return:
+ y: water depth referred to the bottom of the riverbed [m]
+ """
+
+ alpha = 5.0 / 3.0 # exponent (5/3)
+ eps = 1.0e-06
+ max_tries = 1000
+
+ rs0 = np.sqrt(s0)
+ usalpha = 1.0 / alpha
+
+ # cotangent(angle of the riverbed side - dXdY)
+ if ANalv < np.pi / 2:
+ # triangular or trapezoid cross-section
+ c = cotan(ANalv)
+ else:
+ # rectangular corss-section
+ c = 0.0
+
+ # sin(angle of the riverbed side - dXdY)
+ if ANalv < np.pi / 2:
+ # triangular or trapezoid cross-section
+ s = np.sin(ANalv)
+ else:
+ # rectangular cross-section
+ s = 1.0
+
+ # water depth first approximation y0 based on steady state q
+ if Balv == 0:
+ # triangular cross-section
+ y = (Nalv * q / rs0) ** (3.0 / 8.0) * (2 / s) ** 0.25 / c ** (5.0 / 8.0)
+ else:
+ # rectangular cross-section and first approx for trapezoidal cross-section
+ y = (Nalv * q / (rs0 * Balv)) ** usalpha
+
+ if (Balv != 0) and (ANalv < np.pi / 2):
+ # trapezoid cross-section
+ y = (Nalv * q / rs0) ** usalpha * (Balv + 2.0 * y / s) ** 0.4 / (Balv + c * y)
+
+ for tries in range(1, max_tries):
+ # calc Q(y) for the different tries of y
+ q0, Ax, Bx, Px, cel = qoh(y, s0, Balv, ANalv, Nalv)
+ # Ax: wet area[m2]
+ # Bx: cross-section width at water surface[m]
+ # Px: cross-section wet contour [m]
+ # cel: wave celerity[m/s]
+
+ # this is the function we want to find the 0 for f(y)=Q(y)-Q*
+ fy = q0 - q
+ # calc first derivative of f(y) f'(y)=Bx(y)*cel(y)
+ dfy = Bx * cel
+ # calc update for water depth y
+ dy = fy / dfy
+ # update yt+1=yt-f'(yt)/f(yt)
+ y = y - dy
+ # stop loop if correction becomes too small
+ if np.abs(dy) < eps:
+ break
+
+ return y
+
+
+@njit(nogil=True, fastmath=False, cache=True)
+def qoh(y, s0, Balv, ANalv, Nalv):
+ """Discharge q from water depth h.
+ Given a generic river cross-section (rectangular, triangular and trapezoidal)
+ and a water depth (y [m]) referred to the bottom of the riverbed, it uses Manning’s formula to calculate:
+ q: steady-state discharge river discharge [m3/s]
+ a: wet area [m2]
+ b: cross-section width at water surface [m]
+ p: cross-section wet contour [m]
+ cel: wave celerity [m/s]
+ Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030
+
+ :param y: river water depth [m]
+ :param s0: river bed slope (tan B)
+ :param Balv : width of the riverbed [m]
+ :param ANalv : angle of the riverbed side [rad]
+ :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3]
+
+ :return:
+ q,a,b,p,cel
+ """
+
+ alpha = 5.0 / 3.0 # exponent (5/3)
+
+ rs0 = np.sqrt(s0)
+ alpham = alpha - 1.0
+
+ # np.where(ANalv < np.pi/2, triang. or trapeiz., rectangular)
+ # cotangent(angle of the riverbed side - dXdY)
+ c = np.where(
+ ANalv < np.pi / 2,
+ # triangular or trapezoid cross-section
+ cotan(ANalv),
+ # rectangular cross-section
+ 0.0,
+ )
+ # sin(angle of the riverbed side - dXdY)
+ s = np.where(
+ ANalv < np.pi / 2,
+ # triangular or trapezoid cross-section
+ np.sin(ANalv),
+ # rectangular corss-section
+ 1.0,
+ )
+
+ a = (Balv + y * c) * y # wet area [m2]
+ b = Balv + 2.0 * y * c # cross-section width at water surface [m]
+ p = Balv + 2.0 * y / s # cross-section wet contour [m]
+ q = rs0 / Nalv * a**alpha / p**alpham # steady-state discharge [m3/s]
+ cel = (q / 3.0) * (5.0 / a - 4.0 / (p * b * s)) # wave celerity [m/s]
+
+ return q, a, b, p, cel
+
+
+@njit(nogil=True, fastmath=False, cache=True)
+def hoV(V, xpix, Balv, ANalv):
+ """Water depth h from volume V.
+ Given a generic river cross-section (rectangular, triangular and trapezoidal) and a river channel volume V,
+ it calculates the water depth referred to the bottom of the riverbed [m] (y).
+ Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030
+
+ :param V : volume of water in channel riverbed
+ :param xpix : dimension along the flow direction [m]
+ :param Balv : width of the riverbed [m]
+ :param ANalv : angle of the riverbed side [rad]
+
+ :return:
+ y : channel water depth [m]
+ """
+ eps = 1e-6
+
+ c = np.where(
+ ANalv < np.pi / 2, # angle of the riverbed side dXdY [rad]
+ cotan(ANalv), # triangular or trapezoidal cross-section
+ 0.0,
+ ) # rectangular cross-section
+
+ a = V / xpix # wet area [m2]
+
+ # np.where(c < 1.d-6, rectangular, triangular or trapezoidal)
+ y = np.where(
+ np.abs(c) < eps,
+ a / Balv, # rectangular cross-section
+ (-Balv + np.sqrt(Balv**2 + 4 * a * c)) / (2 * c),
+ ) # triangular or trapezoidal cross-section
+
+ return y
+
+
+@njit(nogil=True, fastmath=False, cache=True)
+def qoV(V, xpix, s0, Balv, ANalv, Nalv):
+ """Discharge q from river channel volume V.
+ Given a generic river cross-section (rectangular, triangular and trapezoidal)
+ and a water volume (V [m3]), it uses Manning’s formula to calculate the corresponding discharge (q [m3/s]).
+ Reference: Reggiani, P., Todini, E., & Meißner, D. (2016). On mass and momentum conservation in the variable-parameter Muskingum method. Journal of Hydrology, 543, 562–576. https://doi.org/10.1016/j.jhydrol.2016.10.030
+
+ :param V : volume of water in channel riverbed
+ :param xpix : dimension along the flow direction [m]
+ :param s0: river bed slope (tan B)
+ :param Balv : width of the riverbed [m]
+ :param ANalv : angle of the riverbed side [rad]
+ :param Nalv : channel mannings coefficient n for the riverbed [s/m1/3]
+
+ :return:
+ y : channel water depth [m]
+ """
+ y = hoV(V, xpix, Balv, ANalv)
+ q, a, b, p, cel = qoh(y, s0, Balv, ANalv, Nalv)
+ return q
+
+
+@njit(nogil=True, fastmath=False, cache=True)
+def cotan(x):
+ """There is no cotangent function in numpy"""
+ return np.cos(x) / np.sin(x)
+
+
+@njit(nogil=True, fastmath=False, cache=True)
+def rad_from_dxdy(dxdy):
+ """Calculate radians"""
+ rad = np.arctan(1 / dxdy)
+ # angle = np.rad2deg(rad)
+ return rad
diff --git a/src/lisflood/hydrological_modules/reservoir.py b/src/lisflood/hydrological_modules/reservoir.py
index c2e15cd5..0c4ec944 100755
--- a/src/lisflood/hydrological_modules/reservoir.py
+++ b/src/lisflood/hydrological_modules/reservoir.py
@@ -114,19 +114,24 @@ def initial(self):
# Add reservoir locations to structures map
# (used to modify LddKinematic and to calculate LddStructuresKinematic)
self.var.IsStructureKinematic = np.where(self.var.ReservoirSitesC > 0, np.bool8(1), self.var.IsStructureKinematic)
-
- # load reservoirs in PCRaster
- reservoirs_pcr = loadmap('ReservoirSites', pcr=True)
- self.var.ReservoirSites = reservoirs_pcr
- # filter out reservoirs that are not part of the channel network
+ # Add reservoir locations to structures map (used to modify LddKinematic
+ # and to calculate LddStructuresKinematic)
+ self.var.IsStructureChan = np.where(self.var.ReservoirSitesC > 0, np.bool8(1), self.var.IsStructureChan)
+ # Add reservoir locations to structures map (used to modify LddChan
+ # and to calculate LddStructuresChan)
+
+ ReservoirSitePcr = loadmap('ReservoirSites', pcr=True)
+ self.var.ReservoirSites = ReservoirSitePcr
+ ReservoirSitePcr = ifthen((defined(ReservoirSitePcr) & boolean(decompress(self.var.IsChannel))), ReservoirSitePcr)
+ # Get rid of any reservoirs that are not part of the channel network
# (following logic of 'old' code the inflow into these reservoirs is
# always zero, so either change this or leave them out!)
- reservoirs_pcr = ifthen((defined(reservoirs_pcr) & boolean(decompress(self.var.IsChannel))), reservoirs_pcr)
+ ReservoirSitePcr = ifthen((defined(ReservoirSitePcr) & boolean(decompress(self.var.IsChannel))), ReservoirSitePcr)
# RESERVOIR CHARACTERISTICS
# reservoir storage capacity [m3]
- total_storage = lookupscalar(str(binding['ReservoirTotalStorage']), reservoirs_pcr)
+ total_storage = lookupscalar(str(binding['ReservoirTotalStorage']), ReservoirSitePcr)
total_storage = compressArray(total_storage)
self.var.TotalReservoirStorageM3C = np.where(np.isnan(total_storage), 0, total_storage)
self.var.TotalReservoirStorageM3 = np.compress(self.var.ReservoirSitesC > 0, self.var.TotalReservoirStorageM3C)
@@ -140,7 +145,7 @@ def initial(self):
# flood storage limit (fraction of total storage [-])
if str(binding['ReservoirFloodStorage']).endswith('txt'):
- flood_storage = lookupscalar(str(binding['ReservoirFloodStorage']), reservoirs_pcr)
+ flood_storage = lookupscalar(str(binding['ReservoirFloodStorage']), ReservoirSitePcr)
flood_storage = compressArray(flood_storage)
else:
flood_storage = loadmap('ReservoirFloodStorage')
@@ -149,7 +154,7 @@ def initial(self):
# factor of the flood outflow
if str(binding['ReservoirFloodOutflowFactor']).endswith('txt'):
- factor_outflow = lookupscalar(str(binding['ReservoirFloodStorage']), reservoirs_pcr)
+ factor_outflow = lookupscalar(str(binding['ReservoirFloodStorage']), ReservoirSitePcr)
factor_outflow = compressArray(factor_outflow)
else:
factor_outflow = loadmap('ReservoirFloodOutflowFactor')
@@ -170,17 +175,17 @@ def initial(self):
# RELEASE ATTRIBUTES
# minimum reservoir outflow [m3/s]
- MinReservoirOutflow = lookupscalar(str(binding['ReservoirMinOutflow']), reservoirs_pcr)
+ MinReservoirOutflow = lookupscalar(str(binding['ReservoirMinOutflow']), ReservoirSitePcr)
MinReservoirOutflowC = compressArray(MinReservoirOutflow)
self.var.MinReservoirOutflow = np.compress(self.var.ReservoirSitesC > 0, MinReservoirOutflowC)
# normal outflow [m3/s]
- normal_outflow = lookupscalar(str(binding['ReservoirNormalOutflow']), reservoirs_pcr)
+ normal_outflow = lookupscalar(str(binding['ReservoirNormalOutflow']), ReservoirSitePcr)
normal_outflow = compressArray(normal_outflow)
self.var.NormalReservoirOutflow = np.compress(self.var.ReservoirSitesC > 0, normal_outflow)
# flood-control outflow [m3/s]
- flood_outflow = lookupscalar(str(binding['ReservoirFloodOutflow']), reservoirs_pcr)
+ flood_outflow = lookupscalar(str(binding['ReservoirFloodOutflow']), ReservoirSitePcr)
flood_outflow = compressArray(flood_outflow)
flood_outflow = np.compress(self.var.ReservoirSitesC > 0, flood_outflow)
self.var.FloodReservoirOutflow = factor_outflow * flood_outflow
diff --git a/src/lisflood/hydrological_modules/routing.py b/src/lisflood/hydrological_modules/routing.py
index fd250084..fd157052 100644
--- a/src/lisflood/hydrological_modules/routing.py
+++ b/src/lisflood/hydrological_modules/routing.py
@@ -17,7 +17,9 @@
from __future__ import print_function, absolute_import
from pcraster import lddmask, accuflux, boolean, downstream, pit, path, lddrepair, ifthenelse, cover, nominal, uniqueid, \
- catchment, upstream
+ catchment, upstream, pcr2numpy
+
+import warnings
import numpy as np
@@ -27,8 +29,10 @@
from .inflow import inflow
from .transmission import transmission
from .kinematic_wave_parallel import kinematicWave, kwpt
+from .mct import MCTWave
from ..global_modules.settings import LisSettings, MaskInfo
+from ..global_modules.errors import LisfloodWarning
from ..global_modules.add1 import loadmap, loadmap_base, compressArray, decompress
from . import HydroModule
@@ -42,9 +46,11 @@ class routing(HydroModule):
"""
input_files_keys = {'all': ['beta', 'ChanLength', 'Ldd', 'Channels', 'ChanGrad', 'ChanGradMin',
'CalChanMan', 'ChanMan', 'ChanBottomWidth', 'ChanDepthThreshold',
- 'ChanSdXdY', 'TotalCrossSectionAreaInitValue', 'PrevDischarge'],
+ 'ChanSdXdY', 'TotalCrossSectionAreaInitValue', 'PrevDischarge',
+ 'ChanBottomWMult', 'ChanDepthTMult', 'ChanSMult'],
'SplitRouting': ['CrossSection2AreaInitValue', 'PrevSideflowInitValue', 'CalChanMan2'],
- 'dynamicWave': ['ChannelsDynamic']}
+ 'dynamicWave': ['ChannelsDynamic'],
+ 'MCTRouting': ['ChannelsMCT', 'ChanGradMaxMCT', 'PrevCmMCTInitValue', 'PrevDmMCTInitValue', 'CalChanMan3']}
module_name = 'Routing'
def __init__(self, routing_variable):
@@ -61,115 +67,134 @@ def __init__(self, routing_variable):
def initial(self):
""" initial part of the routing module
"""
+ settings = LisSettings.instance()
+ option = settings.options
maskinfo = MaskInfo.instance()
- self.var.avgdis = maskinfo.in_zero()
- self.var.Beta = loadmap('beta')
- self.var.InvBeta = 1 / self.var.Beta
- # Inverse of beta for kinematic wave
- self.var.ChanLength = loadmap('ChanLength').astype(float)
- self.var.InvChanLength = 1 / self.var.ChanLength
- # Inverse of channel length [1/m]
+
+ # ************************************************************
+ # ***** NUMBER OF ROUTING STEPS *********************
+ # ************************************************************
self.var.NoRoutSteps = int(np.maximum(1, round(self.var.DtSec / self.var.DtSecChannel,0)))
- # Number of sub-steps based on value of DtSecChannel,
- # or 1 if DtSec is smaller than DtSecChannel
- settings = LisSettings.instance()
- option = settings.options
+ # Number of routing sub-steps based on value of DtSecChannel or 1 if DtSec is smaller than DtSecChannel
+
if option['InitLisflood']:
self.var.NoRoutSteps = 1
- # InitLisflood is used!
- # so channel routing step is the same as the general time step
+ # InitLisflood is used! so channel routing step is the same as the general time step
+
self.var.DtRouting = self.var.DtSec / self.var.NoRoutSteps
- # Corresponding sub-timestep (seconds)
+ # Corresponding routing sub-timestep (seconds)
self.var.InvDtRouting = 1 / self.var.DtRouting
self.var.InvNoRoutSteps = 1 / float(self.var.NoRoutSteps)
# inverse for faster calculation inside the dynamic section
-# -------------------------- LDD
+ # ************************************************************
+ # ***** DRAINAGE NETWORK GEOMETRY - LDD *********************
+ # ************************************************************
- self.var.Ldd = lddmask(loadmap('Ldd', pcr=True, lddflag=True), self.var.MaskMap)
+ self.var.Ldd = lddmask(loadmap('Ldd', pcr=True, lddflag=True), self.var.MaskMap) #pcr
# Cut ldd to size of MaskMap (NEW, 29/9/2004)
# Prevents 'unsound' ldd if MaskMap covers sub-area of ldd
- # Count (inverse of) upstream area for each pixel
- # Needed if we want to calculate average values of variables
- # upstream of gauge locations
+ self.var.MaskMap = boolean(self.var.Ldd) #pcr
+ self.var.MaskMapNp=compressArray(self.var.MaskMap) #np
+ # Use boolean version of Ldd as calculation mask
+ # (important for correct mass balance check any water generated outside of Ldd won't reach channel anyway)
- self.var.UpArea = accuflux(self.var.Ldd, self.var.PixelAreaPcr)
+ self.var.UpArea = accuflux(self.var.Ldd, self.var.PixelAreaPcr) #pcr
# Upstream contributing area for each pixel
- # Note that you might expext that values of UpArea would be identical to
- # those of variable CatchArea (see below) at the outflow points.
- # This is NOT actually the case, because outflow points are shifted 1
+ # Note that you might expect that values of UpArea would be identical to those of variable CatchArea (see below)
+ # at the outflow points. This is NOT actually the case, because outflow points are shifted 1.
# cell in upstream direction in the calculation of CatchArea!
- self.var.InvUpArea = 1 / self.var.UpArea
+ self.var.InvUpArea = 1 / self.var.UpArea #pcr
+ # Count (inverse of) upstream area for each pixel
+ # Needed if we want to calculate average values of variables upstream of gauge locations
# Calculate inverse, so we can multiply in dynamic (faster than divide)
- self.var.IsChannelPcr = boolean(loadmap('Channels', pcr=True))
- self.var.IsChannel = np.bool8(compressArray(self.var.IsChannelPcr))
- # Identify channel pixels
- self.var.IsChannelKinematic = self.var.IsChannel.copy()
- # Identify kinematic wave channel pixels
- # (identical to IsChannel, unless dynamic wave is used, see below)
- self.var.IsStructureKinematic = np.bool8(maskinfo.in_zero())
-
- # Map that identifies special inflow/outflow structures (reservoirs, lakes) within the
- # kinematic wave channel routing. Set to (dummy) value of zero modified in reservoir and lake
- # routines (if those are used)
- LddChan = lddmask(self.var.Ldd, self.var.IsChannelPcr)
- # ldd for Channel network
- self.var.MaskMap = boolean(self.var.Ldd)
- # Use boolean version of Ldd as calculation mask
- # (important for correct mass balance check
- # any water generated outside of Ldd won't reach
- # channel anyway)
- self.var.LddToChan = lddrepair(ifthenelse(self.var.IsChannelPcr, 5, self.var.Ldd))
- # Routing of runoff (incl. ground water)en
- AtOutflow = boolean(pit(self.var.Ldd))
- # find outlet points...
+ self.var.IsChannelPcr = boolean(loadmap('Channels', pcr=True)) #pcr
+ self.var.IsChannel = np.bool8(compressArray(self.var.IsChannelPcr)) #bool
+ self.var.IsChannelPcr = boolean(decompress(self.var.IsChannel)) #pcr
+ # Identify grid cells containing a river channel
- if option['dynamicWave']:
- IsChannelDynamic = boolean(loadmap('ChannelsDynamic', pcr=True))
- # Identify channel pixels where dynamic wave is used
- self.var.IsChannelKinematic = (self.var.IsChannelPcr == 1) & (IsChannelDynamic == 0)
- # Identify (update) channel pixels where kinematic wave is used
- self.var.LddKinematic = lddmask(self.var.Ldd, self.var.IsChannelKinematic)
- # Ldd for kinematic wave: ends (pit) just before dynamic stretch
-
- # Following statements produce an ldd network that connects the pits in
- # LddKinematic to the nearest downstream dynamic wave pixel
-
- self.var.AtLastPoint = (downstream(self.var.Ldd, AtOutflow) == 1) & (AtOutflow != 1) & self.var.IsChannelPcr
-
- # NEW 23-6-2005
- # Dynamic wave routine gives no outflow out of pits, so we calculate this
- # one cell upstream (WvD)
- # (implies that most downstream cell is not taken into account in mass balance
- # calculations, even if dyn wave is not used)
- # Only include points that are on a channel (otherwise some small 'micro-catchments'
- # are included, for which the mass balance cannot be calculated
- # properly)
+ self.var.IsChannelKinematic = self.var.IsChannel.copy() #bool
+ # Identify channel pixels using kinematic or split routing
+ # (identical to IsChannel, unless MCT wave is used, see below)
+
+ self.var.IsStructureChan = np.bool8(maskinfo.in_zero()) #bool
+ # Initialise map that identifies special inflow/outflow structures (reservoirs, lakes) within the
+ # channel routing. Set to (dummy) value of zero modified in reservoirs and lakes functions (if those are used)
+
+ self.var.IsStructureKinematic = np.bool8(maskinfo.in_zero()) #bool
+ # Initialise map that identifies special inflow/outflow structures (reservoirs, lakes) within the
+ # kinematic wave channel routing. Set to (dummy) value of zero modified in reservoirs and lakes functions
+ # (if those are used)
+ self.var.LddChan = lddmask(self.var.Ldd, self.var.IsChannelPcr) #pcr
+ self.var.LddChanNp=compressArray(self.var.LddChan) #np
+ # LDD for channel network (excluding grid cells that do not have a channel)
+ # It is the same as Ldd if there is a channel on every grid cell
+ # We need a channel on each grid cell now but this could change in the future
+
+ self.var.LddToChan = lddrepair(ifthenelse(self.var.IsChannelPcr, 5, self.var.Ldd)) #pcr
+ self.var.LddToChanNp=compressArray(self.var.LddToChan) #np
+ # Routing of runoff (incl. groundwater) to the river channel
+
+ if option['dynamicWave']:
+ pass
+ # # legacy code
+ # AtOutflow = boolean(pit(self.var.Ldd)) # pcr
+ # AtOutflowNp = compressArray(AtOutflow) # np
+ # # find outlet points...
+ # # Function 'pit' assignes a unique number starting from 1 to pit cells (ldd=5) in the Ldd
+ #
+ # IsChannelDynamic = boolean(loadmap('ChannelsDynamic', pcr=True))
+ # # Identify channel pixels where dynamic wave is used
+ # self.var.IsChannelKinematic = (self.var.IsChannelPcr == 1) & (IsChannelDynamic == 0)
+ # # Identify (update) channel pixels where kinematic wave is used
+ # self.var.LddKinematic = lddmask(self.var.Ldd, self.var.IsChannelKinematic)
+ # # Ldd for kinematic wave: ends (pit) just before dynamic stretch
+ #
+ # # Following statements produce an ldd network that connects the pits in
+ # # LddKinematic to the nearest downstream dynamic wave pixel
+ #
+ # self.var.AtLastPoint = (downstream(self.var.Ldd, AtOutflow) == 1) & (AtOutflow != 1) & self.var.IsChannelPcr
+ #
+ # # NEW 23-6-2005
+ # # Dynamic wave routine gives no outflow out of pits, so we calculate this
+ # # one cell upstream (WvD)
+ # # (implies that most downstream cell is not taken into account in mass balance
+ # # calculations, even if dyn wave is not used)
+ # # Only include points that are on a channel (otherwise some small 'micro-catchments'
+ # # are included, for which the mass balance cannot be calculated
+ # # properly)
else:
- self.var.LddKinematic = LddChan
- # No dynamic wave, so kinematic ldd equals channel ldd
- self.var.AtLastPoint = AtOutflow
- self.var.AtLastPointC = np.bool8(compressArray(self.var.AtLastPoint))
- # assign unique identifier to each of them
- maskinfo = MaskInfo.instance()
- lddC = compressArray(self.var.LddKinematic)
- inAr = decompress(np.arange(maskinfo.info.mapC[0], dtype="int32"))
- # giving a number to each non missing pixel as id
- self.var.downstruct = (compressArray(downstream(self.var.LddKinematic, inAr))).astype("int32")
+ # For all routing options (kinematic, split and MCT)
+ self.var.LddKinematic = self.var.LddChan
+ self.var.LddKinematicNp = compressArray(self.var.LddKinematic) # np
+
+ self.var.AtLastPoint = boolean(pit(self.var.Ldd)) #pcr
+ # Assign True to each of the grid cells where there are outlet points
+ # Function 'pit' assigns a unique number starting from 1 to pit cells (ldd=5) in the Ldd
+ self.var.AtLastPointC = np.bool8(compressArray(self.var.AtLastPoint)) #np
+
+ lddC = compressArray(self.var.LddChan) #np
+ inAr = decompress(np.arange(maskinfo.info.mapC[0], dtype="int32")) #pcr
+ # Assign a number to each non-missing pixel as cell id, starting from 0
+
+ self.var.downstruct = (compressArray(downstream(self.var.LddChan, inAr))).astype("int32") #np
# each upstream pixel gets the id of the downstream pixel
- self.var.downstruct[lddC == 5] = maskinfo.info.mapC[0]
- # all pits gets a high number
- # upstream function in numpy
-
- OutflowPoints = nominal(uniqueid(self.var.AtLastPoint))
- # and assign unique identifier to each of them
- self.var.Catchments = (compressArray(catchment(self.var.Ldd, OutflowPoints))).astype(np.int32)
- CatchArea = np.bincount(self.var.Catchments, weights=self.var.PixelArea)[self.var.Catchments]
+ self.var.downstruct[lddC == 5] = maskinfo.info.mapC[0] #np
+ # assigning to all pits an id number higher than any of the other cells
+
+ OutflowPoints = nominal(uniqueid(self.var.AtLastPoint)) #pcr
+ OutflowPointsNp = compressArray(OutflowPoints) #np
+ # assigning id to the outflow points starting from 1
+
+ self.var.Catchments = (compressArray(catchment(self.var.Ldd, OutflowPoints))).astype(np.int32) #np
+ # assign outlet id to all pixel in its catchment
# define catchment for each outflow point
+
+ CatchArea = np.bincount(self.var.Catchments, weights=self.var.PixelArea)[self.var.Catchments] #np
# Compute area of each catchment [m2]
# Note: in earlier versions this was calculated using the "areaarea" function,
# changed to "areatotal" in order to enable handling of grids with spatially
@@ -182,23 +207,44 @@ def initial(self):
# ************************************************************
self.var.ChanGrad = np.maximum(loadmap('ChanGrad'), loadmap('ChanGradMin'))
- # avoid calculation of Alpha using ChanGrad=0: this creates MV!
+ # River bed slope
+ # Avoid calculation of Alpha using ChanGrad=0: this creates MV!
+ # Set minimum channel slope in rivers (ChanGradMin)
+
self.var.CalChanMan = loadmap('CalChanMan')
self.var.ChanMan = self.var.CalChanMan * loadmap('ChanMan')
- # Manning's n is multiplied by ChanManCal
- # enables calibration for peak timing
- self.var.ChanBottomWidth = loadmap('ChanBottomWidth')
- ChanDepthThreshold = loadmap('ChanDepthThreshold')
- ChanSdXdY = loadmap('ChanSdXdY')
- self.var.ChanUpperWidth = self.var.ChanBottomWidth + 2 * ChanSdXdY * ChanDepthThreshold
+ # Manning's roughtness coefficient n is multiplied by ChanManCal for calibration
+
+ self.var.CalChanBottomWidth = loadmap('ChanBottomWMult')
+ self.var.ChanBottomWidth = self.var.CalChanBottomWidth * loadmap('ChanBottomWidth')
+ # Riverbed width [m]
+
+ self.var.CalChanDepthThreshold = loadmap('ChanDepthTMult')
+ ChanDepthThreshold = self.var.CalChanDepthThreshold * loadmap('ChanDepthThreshold')
+ # Bankfull river depth [m]
+
+ self.var.CalChanSdXdY= loadmap('ChanSMult')
+ self.var.ChanSdXdY = self.var.CalChanSdXdY * loadmap('ChanSdXdY')
+ # Riverbed sides slope
+
+ # self.var.ChanSdXdY = maskinfo.in_zero() # rectangular cross-section
+ # self.var.ChanBottomWidth = maskinfo.in_zero() # triangular cross-section
+
+ self.var.ChanUpperWidth = self.var.ChanBottomWidth + 2 * self.var.ChanSdXdY * ChanDepthThreshold
# Channel upper width [m]
self.var.TotalCrossSectionAreaBankFull = 0.5 * \
ChanDepthThreshold * (self.var.ChanUpperWidth + self.var.ChanBottomWidth)
- # Area (sq m) of bank full discharge cross section [m2]
- # (trapezoid area equation)
- TotalCrossSectionAreaHalfBankFull = 0.5 * self.var.TotalCrossSectionAreaBankFull
+ # Area (sq m) of bank full discharge cross-section [m2] (trapezoid area equation)
+
+ # cmcheck - TotalCrossSectionAreaHalfBankFull is not 1/2 TotalCrossSectionAreaBankFull it's trapezoid
+ # ChanUpperWidthHalfBankFull = self.var.ChanBottomWidth + 2 * self.var.ChanSdXdY * 0.5 * ChanDepthThreshold
+ # TotalCrossSectionAreaHalfBankFull = 0.5 * \
+ # 0.5 * ChanDepthThreshold * (ChanUpperWidthHalfBankFull + self.var.ChanBottomWidth)
# Cross-sectional area at half bankfull [m2]
# This can be used to initialise channel flow (see below)
+ BankFullPerc = loadmap('BankFullPerc')
+ TotalCrossSectionAreaHalfBankFull = BankFullPerc * self.var.TotalCrossSectionAreaBankFull
+ # set BankFullPerc to 0.5 for half bankfull
TotalCrossSectionAreaInitValue = loadmap('TotalCrossSectionAreaInitValue')
self.var.TotalCrossSectionArea = np.where(TotalCrossSectionAreaInitValue == -9999, TotalCrossSectionAreaHalfBankFull, TotalCrossSectionAreaInitValue)
@@ -206,11 +252,10 @@ def initial(self):
# otherwise TotalCrossSectionAreaInitValue (typically end map from previous simulation)
if option['SplitRouting']:
- # in_zero = maskinfo.in_zero()
CrossSection2AreaInitValue = loadmap('CrossSection2AreaInitValue')
self.var.CrossSection2Area = np.where(CrossSection2AreaInitValue == -9999, maskinfo.in_zero(), CrossSection2AreaInitValue)
- # cross-sectional area [m2] for 2nd line of routing: if initial value in binding equals -9999 the value is set to 0
- # otherwise CrossSection2AreaInitValue (typically end map from previous simulation)
+ # cross-sectional area [m2] for 2nd line of routing (over bankfull only): if initial value in binding equals
+ # -9999 the value is set to 0 otherwise CrossSection2AreaInitValue (typically end map from previous simulation)
PrevSideflowInitValue = loadmap('PrevSideflowInitValue')
self.var.Sideflow1Chan = np.where(PrevSideflowInitValue == -9999, maskinfo.in_zero(), PrevSideflowInitValue)
@@ -218,47 +263,61 @@ def initial(self):
# otherwise PrevSideflowInitValue (typically end map from previous simulation)
# ************************************************************
- # ***** CHANNEL ALPHA (KIN. WAVE)*****************************
+ # ***** CHANNEL ALPHA (KIN. WAVE) ****************************
# ************************************************************
- # Following calculations are needed to calculate Alpha parameter in kinematic
- # wave. Alpha currently fixed at half of bankful depth (this may change in
- # future versions!)
+ # Manning's steady state flow equations
+ # from Ven The Chow - Applied Hydrology - page 283
+ # https: // wecivilengineers.files.wordpress.com / 2017 / 10 / applied - hydrology - ven - te - chow.pdf
+ # A = Alpha * Q ** Beta
+ # Q = (A/Alpha) ** (1/Beta) = (invAlpha * A)**invBeta
+
+ self.var.Beta = loadmap('beta')
+ # This is 3/5. Exponent of Manning's equation A=alpha*Q^beta -> Q= 1/alpha * A^(1/beta)
+ self.var.InvBeta = 1 / self.var.Beta
+ # Inverse of beta for kinematic wave
+ self.var.ChanLength = loadmap('ChanLength').astype(float)
+ # Channel length [m]
+ self.var.InvChanLength = 1 / self.var.ChanLength
+ # Inverse of channel length [1/m]
+
+ # Following calculations are needed to calculate Alpha parameter in kinewave.
+ # Alpha currently fixed at half of bankfull depth (this may change in future versions!)
ChanWaterDepthAlpha = np.where(self.var.IsChannel, 0.5 * ChanDepthThreshold, 0.0)
# Reference water depth for calculation of Alpha: half of bankfull
self.var.ChanWettedPerimeterAlpha = self.var.ChanBottomWidth + 2 * \
- np.sqrt(np.square(ChanWaterDepthAlpha) + np.square(ChanWaterDepthAlpha * ChanSdXdY))
- # Channel wetted perimeter [m](Pythagoras)
- AlpTermChan = (self.var.ChanMan / (np.sqrt(self.var.ChanGrad))) ** self.var.Beta
- self.var.AlpPow = 2.0 / 3.0 * self.var.Beta
+ np.sqrt(np.square(ChanWaterDepthAlpha) + np.square(ChanWaterDepthAlpha * self.var.ChanSdXdY))
+ # Channel wetted perimeter half bankfull [m](Pythagoras)
+ AlpTermChan = (self.var.ChanMan / (np.sqrt(self.var.ChanGrad))) ** self.var.Beta
+ self.var.AlpPow = 2.0 / 3.0 * self.var.Beta # this is 2/5
self.var.ChannelAlpha = (AlpTermChan * (self.var.ChanWettedPerimeterAlpha ** self.var.AlpPow)).astype(float)
+
self.var.InvChannelAlpha = 1 / self.var.ChannelAlpha
# ChannelAlpha for kinematic wave
# ************************************************************
- # ***** CHANNEL INITIAL DISCHARGE ****************************
+ # ***** CHANNEL INITIAL DISCHARGE (KINEMATIC) ****************
# ************************************************************
- self.var.ChanM3 = self.var.TotalCrossSectionArea * self.var.ChanLength
- # channel water volume [m3]
- self.var.ChanIniM3 = self.var.ChanM3.copy()
- self.var.ChanM3Kin = self.var.ChanIniM3.copy().astype(float)
+ self.var.ChanM3 = self.var.TotalCrossSectionArea * self.var.ChanLength #np
+ # Total water volume in river channel [m3]
+ self.var.ChanIniM3 = self.var.ChanM3.copy() #np
+ self.var.ChanM3Kin = self.var.ChanIniM3.copy().astype(float) #np
# Initialise water volume in kinematic wave channels [m3]
+
self.var.ChanQKin = np.where(self.var.ChannelAlpha > 0, (self.var.TotalCrossSectionArea / self.var.ChannelAlpha) ** self.var.InvBeta, 0).astype(float)
- # Initialise discharge at kinematic wave pixels (note that InvBeta is
- # simply 1/beta, computational efficiency!)
+ # Initialise discharge at kinematic wave pixels (note that InvBeta is simply 1/beta, computational efficiency!)
self.var.CumQ = maskinfo.in_zero()
- # ininialise sum of discharge to calculate average
+ # Initialise sum of discharge to calculate average
-# ************************************************************
-# ***** CHANNEL INITIAL DYNAMIC WAVE *************************
-# ************************************************************
+ # ************************************************************
+ # ***** CHANNEL INITIAL DYNAMIC WAVE *************************
+ # ************************************************************
if option['dynamicWave']:
pass
- # TODO !!!!!!!!!!!!!!!!!!!!
-
+ # legacy code
# lookchan = lookupstate(TabCrossSections, ChanCrossSections, ChanBottomLevel, self.var.ChanLength,
# DynWaveConstantHeadBoundary + ChanBottomLevel)
# ChanIniM3 = ifthenelse(AtOutflow, lookchan, ChanIniM3)
@@ -318,117 +377,238 @@ def initial(self):
# IsChannelDynamic, ChanQDyn, self.var.ChanQKin)
# Channel discharge: combine results of kinematic and dynamic wave
else:
+ # For all routing options (kinematic, split and MCT)
+ PrevDischarge = loadmap('PrevDischarge')
+ # Outflow (x+dx) Q at the end of previous computation step for full cross-section (instant)
+ # Used to calculated Inflow (x) from upstream pixels at the beginning of the computation step
+ self.var.ChanQ = np.where(PrevDischarge == -9999, self.var.ChanQKin, PrevDischarge) #np
+ # Initialise instantaneous channel discharge: cold start: equal to ChanQKin [m3/s]
- # ***** NO DYNAMIC WAVE *************************
- # Dummy code if dynamic wave is not used, in which case ChanQ equals ChanQKin
- # (needed only for polder routine)
+ # initialising average outflow for all routing options (kinematic, split and MCT)
+ self.var.ChanQAvgDt = maskinfo.in_zero()
+ self.var.ChanQKinAvgDt = maskinfo.in_zero()
+ # Outflow (x+dx) Q over the computation step for full cross-section (average)
+
+ self.var.Chan2QKinAvgDt = maskinfo.in_zero()
+ # Outflow (x+dx) Q over the computation step for second line of split routing (average)
+
+ # We do not need a state file to initialise the average outflow discharge (ChanQAvgDt and ChanQKinAvgDt).
+ # Initialisation would be necessary for pixels in order 0 (aka head pixels), but there is no upstream contribution for pixels in order 0.
+ # For pixels in order 1 and beyond, upstream contribution is calculated during the calculation step.
- PrevDischarge = loadmap('PrevDischarge')
- self.var.ChanQ = np.where(PrevDischarge == -9999, self.var.ChanQKin, PrevDischarge)
- # initialise channel discharge: cold start: equal to ChanQKin
- # [m3/s]
+ # ************************************************************
+ # ***** CUMULATIVE OUTPUT VARIABLES *************************
+ # ************************************************************
# Initialising cumulative output variables
- # These are all needed to compute the cumulative mass balance error
+ self.var.avgdis = maskinfo.in_zero()
+ # average discharge over the simulation period
+ # These are all needed to compute the cumulative mass balance error
self.var.DischargeM3Out = maskinfo.in_zero()
- # cumulative discharge at outlet [m3]
+ # Cumulative outflow volume at outlet [m3]
self.var.TotalQInM3 = maskinfo.in_zero()
- # cumulative inflow from inflow hydrographs [m3]
+ # Cumulative inflow volume from inflow hydrographs [m3]
self.var.sumDis = maskinfo.in_zero()
+ self.var.sumIn = maskinfo.in_zero()
+ # cmcheck - non so se sostituita da self.var.sumInWB
self.var.sumInWB = maskinfo.in_zero()
def initialSecond(self):
- """ initial part of the second channel routing module
+ """ initial part of the second line channel routing module for split routing
"""
settings = LisSettings.instance()
option = settings.options
flags = settings.flags
+ self.var.ChannelAlpha2 = None
+ # Default value, if split-routing is not active and water is routed only in the main channel
- self.var.ChannelAlpha2 = None # default value, if split-routing is not active and only water is routed only in the main channel
# ************************************************************
# ***** CHANNEL INITIAL SPLIT UP IN SECOND CHANNEL************
# ************************************************************
if option['SplitRouting']:
+ # calculate kin wave parameters for the (virtual) channel of the second line of routing
+
ChanMan2 = (self.var.ChanMan / self.var.CalChanMan) * loadmap('CalChanMan2')
+ # Manning's roughtness coefficient n for second line of routing
AlpTermChan2 = (ChanMan2 / (np.sqrt(self.var.ChanGrad))) ** self.var.Beta
self.var.ChannelAlpha2 = (AlpTermChan2 * (self.var.ChanWettedPerimeterAlpha ** self.var.AlpPow)).astype(float)
+ #cmcheck -> using channel wetted perimeter of half bankfull ChanWettedPerimeterAlpha ?
self.var.InvChannelAlpha2 = 1 / self.var.ChannelAlpha2
# calculating second Alpha for second (virtual) channel
if not(option['InitLisflood']):
-
- # use loadmap_base function as we don't want to cache avgdis it in the calibration
- self.var.QLimit = loadmap_base('AvgDis') * loadmap('QSplitMult')
-
+ # ************************************************************
+ # ***** INITIALISE SECOND LINE OF ROUTING (KINEMATIC) ********
+ # ************************************************************
# Over bankful discharge starts at QLimit
# lower discharge limit for second line of routing
- # set to mutiple of average discharge (map from prerun)
+ # set to multiple of average discharge (map from prerun)
# QSplitMult =2 is around 90 to 95% of Q
+ self.var.QLimit = loadmap_base('AvgDis') * loadmap('QSplitMult')
+ # Using loadmap_base function as we don't want to cache avgdis in the calibration
self.var.M3Limit = self.var.ChannelAlpha * self.var.ChanLength * (self.var.QLimit ** self.var.Beta)
# Water volume in bankful when over bankful discharge starts
-
- # QLimit should NOT be dependent on the NoRoutSteps (number of routing steps)
- # self.var.QLimit = self.var.QLimit / self.var.NoRoutSteps #original
-
- ###############################################
- # CM mod
- # TEMPORARY WORKAROUND FOR EFAS XDOM!!!!!!!!!!
- # This must be removed
- # self.var.QLimit = self.var.QLimit / 24.0
- ###############################################
+ # Manning's equation
self.var.Chan2M3Start = self.var.ChannelAlpha2 * self.var.ChanLength * (self.var.QLimit ** self.var.Beta)
- # virtual amount of water in the channel through second line
-
+ # Virtual (note we use ChannelAlpha2 now) amount of water in the main channel at the activation of second line of routing 'floodplains' (=> start using increased Manning coeff 2)
self.var.Chan2QStart = self.var.QLimit - compressArray(upstream(self.var.LddKinematic, decompress(self.var.QLimit)))
- # because kinematic routing with a low amount of discharge leads to long travel time:
+ # Virtual outflow from main channel at the activation of second line of routing (=> start using increased Manning coeff 2)
+ # Because kinematic routing with a low amount of discharge leads to long travel time:
# Starting Q for second line is set to a higher value
self.var.Chan2M3Kin = self.var.CrossSection2Area * self.var.ChanLength + self.var.Chan2M3Start
+ # Total (virtual) volume of water in river channel when second routing line is active (= using increased Manning coeff 2)
self.var.ChanM3Kin = self.var.ChanM3 - self.var.Chan2M3Kin + self.var.Chan2M3Start
-
- self.var.ChanM3Kin = np.where((self.var.ChanM3Kin < 0.0) & (self.var.ChanM3Kin > -0.0000001),0.0,self.var.ChanM3Kin) # this line prevents the warmstart from failing in case of small numerical imprecisions when writing and reading the maps
+ # (Real) Volume of water in main channel when second line of routing is active (= using riverbed Manning coeff)
+
+ self.var.ChanM3Kin = np.where((self.var.ChanM3Kin < 0.0) & (self.var.ChanM3Kin > -0.0000001),0.0,self.var.ChanM3Kin)
+ # this line prevents the warm start from failing in case of small numerical imprecisions when writing and reading the maps
self.var.Chan2QKin = (self.var.Chan2M3Kin * self.var.InvChanLength * self.var.InvChannelAlpha2) ** (self.var.InvBeta)
+ # Total (virtual) outflow from river channel when second routing line is active (= using increased Manning coeff 2)
self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta)
+ # (Real) outflow from main channel when second line of routing is active (= using riverbed Manning coeff 2)
+
+
+ # ************************************************************
+ # ***** INITIALISE PARALLEL KINEMATIC WAVE ROUTER ************
+ # ************************************************************
# Initialise parallel kinematic wave router: main channel-only routing if self.var.ChannelAlpha2 is None; else split-routing(main channel + floodplains)
+ # Initialization includes LDD for kinematic routing
maskinfo = MaskInfo.instance()
self.river_router = kinematicWave(compressArray(self.var.LddKinematic), ~maskinfo.info.mask, self.var.ChannelAlpha,
- self.var.Beta, self.var.ChanLength, self.var.DtRouting,
+ self.var.Beta, self.var.ChanLength, self.var.DtRouting,
alpha_floodplains=self.var.ChannelAlpha2, flagnancheck=flags['nancheck'])
-
- if option['InitLisflood'] and option['repMBTs']:
- self.var.StorageStepINIT= self.var.ChanM3Kin
- self.var.DischargeM3StructuresIni = maskinfo.in_zero()
- if option['simulateReservoirs']:
- self.var.StorageStepINIT += self.var.ReservoirStorageIniM3
- if option['simulateLakes']:
- self.var.StorageStepINIT += self.var.LakeStorageIniM3
- self.var.StorageStepINIT = np.take(np.bincount(self.var.Catchments, weights=self.var.StorageStepINIT), self.var.Catchments)
-
- if not option['InitLisflood'] and option['repMBTs']:
- DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
- if not(option['SplitRouting']):
- self.var.StorageStepINIT = self.var.ChanM3Kin
- if option['simulateReservoirs']:
- self.var.StorageStepINIT += self.var.ReservoirStorageIniM3
- DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
- if option['simulateLakes']:
- self.var.StorageStepINIT += self.var.LakeStorageIniM3
- DisStructure += np.where(compressArray(self.var.IsUpsOfStructureLake), 0.5 * self.var.ChanQ * self.var.DtRouting, 0)
- self.var.DischargeM3StructuresIni = np.take(np.bincount(self.var.Catchments, weights=DisStructure), self.var.Catchments)
- else:
- self.var.StorageStepINIT= self.var.ChanM3Kin+self.var.Chan2M3Kin-self.var.Chan2M3Start
- if option['simulateReservoirs']:
+
+ # ************************************************************
+ # ***** WATER BALANCE ************
+ # ************************************************************
+ if option['InitLisflood'] and option['repMBTs']:
+ # Calculate initial water storage in rivers (no lakes no reservoirs)
+ # self.var.StorageStepINIT= self.var.ChanM3Kin
+ self.var.StorageStepINIT = self.var.ChanM3
+ # Initial water volume in river channels
+ self.var.DischargeM3StructuresIni = maskinfo.in_zero()
+ if option['simulateReservoirs']:
self.var.StorageStepINIT += self.var.ReservoirStorageIniM3
- if option['simulateLakes']:
+ if option['simulateLakes']:
self.var.StorageStepINIT += self.var.LakeStorageIniM3
- self.var.StorageStepINIT = np.take(np.bincount(self.var.Catchments, weights=self.var.StorageStepINIT), self.var.Catchments)
+ self.var.StorageStepINIT = np.take(np.bincount(self.var.Catchments, weights=self.var.StorageStepINIT), self.var.Catchments)
+
+ if not option['InitLisflood'] and option['repMBTs']:
+ self.var.StorageStepINIT = self.var.ChanM3
+ # DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+ DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
+ if not(option['SplitRouting']):
+ # self.var.StorageStepINIT = self.var.ChanM3Kin
+ # self.var.StorageStepINIT = self.var.ChanM3
+ if option['simulateReservoirs']:
+ self.var.StorageStepINIT += self.var.ReservoirStorageIniM3
+ # DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+ DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
+ if option['simulateLakes']:
+ self.var.StorageStepINIT += self.var.LakeStorageIniM3
+ DisStructure += np.where(compressArray(self.var.IsUpsOfStructureLake), 0.5 * self.var.ChanQ * self.var.DtRouting, 0)
+ self.var.DischargeM3StructuresIni = np.take(np.bincount(self.var.Catchments, weights=DisStructure), self.var.Catchments)
+ else:
+ # self.var.StorageStepINIT= self.var.ChanM3Kin+self.var.Chan2M3Kin-self.var.Chan2M3Start
+ # self.var.StorageStepINIT = self.var.ChanM3
+ if option['simulateReservoirs']:
+ self.var.StorageStepINIT += self.var.ReservoirStorageIniM3
+ if option['simulateLakes']:
+ self.var.StorageStepINIT += self.var.LakeStorageIniM3
+ self.var.StorageStepINIT = np.take(np.bincount(self.var.Catchments, weights=self.var.StorageStepINIT), self.var.Catchments)
+
+ def initialMCT(self):
+ """ initial part of the Muskingum-Cunge-Todini routing module
+ """
+ settings = LisSettings.instance()
+ option = settings.options
+
+ if not (option['SplitRouting']):
+ self.var.ChannelAlpha2 = None
+ # default value, if split-routing is not active and water is routed only in the riverbed channel
+
+ # ************************************************************
+ # ***** INITIALISATION FOR MCT ROUTING ************
+ # ************************************************************
+
+ # even if MCT is active, it should be deactivated if there is no MCT cell in the domain
+ if option['MCTRouting']:
+ self.var.IsChannelMCTPcr = boolean(loadmap('ChannelsMCT', pcr=True)) #pcr
+ self.var.IsChannelMCT = np.bool8(compressArray(self.var.IsChannelMCTPcr)) #bool
+ if self.var.IsChannelMCT.sum()==0:
+ warnings.warn(LisfloodWarning('There are no MCT grid cell. MCT routing is deactivated'))
+ option['MCTRouting'] = False
+
+ if option['MCTRouting']:
+ maskinfo = MaskInfo.instance()
+
+ self.var.IsChannelMCTPcr = boolean(decompress(self.var.IsChannelMCT)) # pcr
+ # Identify channel pixels where Muskingum-Cunge-Todini is used
+
+ self.var.mctmask = np.bool8(pcr2numpy(self.var.IsChannelMCTPcr,0))
+ # Create a mask with cells using MCT
+
+ self.var.IsChannelKinematicPcr = (self.var.IsChannelPcr == 1) & (self.var.IsChannelMCTPcr == 0) #pcr
+ self.var.IsChannelKinematic = np.bool8(compressArray(self.var.IsChannelKinematicPcr)) #np
+ # Identify channel pixels where Kinematic wave is used instead of MCT
+
+ self.var.LddMCT = lddmask(self.var.LddChan, self.var.IsChannelMCTPcr) #pcr
+ # Ldd for MCT routing
+
+ self.var.LddKinematic = lddmask(self.var.LddChan, self.var.IsChannelKinematicPcr) #pcr
+ # Ldd for kinematic routing
+
+ ChanGradMaxMCT = loadmap('ChanGradMaxMCT')
+ # Maximum riverbed slope for MCT rivers
+ # Check where IsChannelMCT is True and values in ChanGrad > ChanGradMaxMCT
+ MCT_slope_mask = np.logical_and(self.var.IsChannelMCT, self.var.ChanGrad > ChanGradMaxMCT)
+ # Update values in ChanGrad where the condition is met
+ self.var.ChanGrad[MCT_slope_mask] = ChanGradMaxMCT
+ # set max channel slope for MCT pixels
+
+ # cmcheck
+ # This could become a calibration parameter if we want to use MCT+SplitRouting
+ self.var.ChanManMCT = (self.var.ChanMan / self.var.CalChanMan) * loadmap('CalChanMan3')
+ # Mannings coefficient for MCT pixels (same as second line of split routing)
+
+ PrevCmMCT = loadmap('PrevCmMCTInitValue')
+ self.var.PrevCm0 = np.where(PrevCmMCT == -9999, maskinfo.in_one(), PrevCmMCT) #np
+ # Courant numnber (Cm) for MCT at previous time step t0
+ PrevDmMCT = loadmap('PrevDmMCTInitValue')
+ self.var.PrevDm0 = np.where(PrevDmMCT == -9999, maskinfo.in_zero(), PrevDmMCT) #np
+ # Reynolds number (Dm) for MCT at previous time step t0
+
+ # ************************************************************
+ # ***** INITIALISE MUSKINGUM-CUNGE-TODINI WAVE ROUTER ********
+ # ************************************************************
+ mct_ldd = self.compress_mct(compressArray(self.var.LddMCT))
+ # Compress LddMCT to array with MCT pixels only
+
+ mapping_mct = self.compress_mct(range(len(self.var.ChanLength)))
+ # create mapping from global domain pixels index to MCT pixels index
+
+ self.mct_river_router = MCTWave(
+ mct_ldd, # MCT Ldd
+ self.var.mctmask, # MCT mask
+ self.var.ChanLength, # Channel length
+ self.var.ChanGrad, # Channel riverbed slope
+ self.var.ChanBottomWidth, # Riverbed bottom width
+ self.var.ChanManMCT, # MCT pixels Mannings' coefficient
+ self.var.ChanSdXdY, # Riverbed side slope
+ self.var.DtRouting, # computation time step for routing [s]
+ self.river_router, # class
+ mapping_mct # MCT pixels mapping
+ )
+
+
# --------------------------------------------------------------------------
# --------------------------------------------------------------------------
@@ -443,42 +623,50 @@ def dynamic(self, NoRoutingExecuted):
self.reservoir_module.dynamic_inloop(NoRoutingExecuted)
self.polder_module.dynamic_inloop()
- # End only with no Lisflood (no reservoirs, lakes and polder with
- # initLisflood)
-
self.inflow_module.dynamic_inloop(NoRoutingExecuted)
self.transmission_module.dynamic_inloop()
- # ************************************************************
- # ***** CHANNEL FLOW ROUTING: KINEMATIC WAVE ****************
- # ************************************************************
-
if not(option['dynamicWave']):
# ************************************************************
- # ***** SIDEFLOW
+ # ***** SIDEFLOW ****************
# ************************************************************
+ # Calculating water contribution to channels
SideflowChanM3 = self.var.ToChanM3RunoffDt.copy()
+ # Runoff contribution from overland flow, upper and lower groundwater zones to the channel during routing sub step [m3]
+
if option['openwaterevapo']:
SideflowChanM3 -= self.var.EvaAddM3Dt
+ # Evaporation from open water (rivers) - withdrawal [m3]
if option['wateruse']:
self.var.WUseAddM3Dt = (self.var.withdrawal_CH_actual_M3_routStep - self.var.returnflow_GwAbs2Channel_M3_routStep)
- SideflowChanM3 -= self.var.WUseAddM3Dt
+ SideflowChanM3 -= self.var.WUseAddM3Dt
+ # Water use abstraction from rivers - withdrawal [m3]
if option['inflow']:
SideflowChanM3 += self.var.QInDt
+ # Flow volume from inlets per sub step [m3]
if option['TransLoss']:
SideflowChanM3 -= self.var.TransLossM3Dt
+ # Transmission flow volume - withdrawal [m3]
+
if not(option['InitLisflood']): # only with no InitLisflood
if option['simulateLakes']:
- SideflowChanM3 += self.var.QLakeOutM3Dt
+ SideflowChanM3 += self.var.QLakeOutM3Dt
+ # Lakes outflow volume per routing sub step [m3]
if option['simulateReservoirs']:
- SideflowChanM3 += self.var.QResOutM3Dt
+ SideflowChanM3 += self.var.QResOutM3Dt
+ # Reservoirs outflow volume per routing sub step [m3]
if option['simulatePolders']:
SideflowChanM3 -= self.var.ChannelToPolderM3Dt
-
-
-
+ # Polders outflow volume per routing sub step [m3]
+
+ # NOTE: polders currently don't work with kinematic wave, but nevertheless
+ # ChannelToPolderM3 is already included in sideflow term (so it's there in case
+ # the polder routine is ever modified to make it work with kin. wave)
+ # Because of wateruse Sideflow might get even smaller than 0,
+ # instead of inflow its outflow
+
### check mass balance within routing ###
if option['repMBTs']:
if (NoRoutingExecuted<1):
@@ -497,210 +685,392 @@ def dynamic(self, NoRoutingExecuted):
self.var.AddedTRUN -= np.take(np.bincount(self.var.Catchments, weights=self.var.EvaAddM3Dt.copy()),self.var.Catchments)
if option['wateruse']:
self.var.AddedTRUN -= np.take(np.bincount(self.var.Catchments, weights=self.var.WUseAddM3Dt.copy()),self.var.Catchments)
-
- # Runoff (surface runoff + flow out of Upper and Lower Zone), outflow from
- # reservoirs and lakes and inflow from external hydrographs are added to the channel
- # system (here in [cu m])
- #
- # NOTE: polders currently don't work with kinematic wave, but nevertheless
- # ChannelToPolderM3 is already included in sideflow term (so it's there in case
- # the polder routine is ever modified to make it work with kin. wave)
- # Because of wateruse Sideflow might get even smaller than 0,
- # instead of inflow its outflow
-
+ # Sideflow contribution to kinematic and split routing grid cells expressed in [cu m /s / m channel length]
SideflowChan = np.where(self.var.IsChannelKinematic, SideflowChanM3 * self.var.InvChanLength * self.var.InvDtRouting,0)
# ************************************************************
- # ***** KINEMATIC WAVE ****************
+ # ***** RIVER ROUTING ****************
# ************************************************************
+ if option['InitLisflood']:
+ self.var.IsChannelKinematic = self.var.IsChannel.copy()
+ # Use kinematic routing in all grid cells
- if option['InitLisflood'] or (not(option['SplitRouting'])):
- # if InitLisflood no split routing is use
+ # only run kinematic for InitLisflood option
+ # either kinematic or split routing always runs for all pixels when using MCT routing
- # ---- Single Routing ---------------
- # No split routing
- # side flow consists of runoff (incl. groundwater), inflow from reservoirs (optional) and external inflow hydrographs (optional)
- SideflowChan[np.isnan(SideflowChan)] = 0 # TEMPORARY FIX - SEE DEBUG ABOVE!
-
- self.river_router.kinematicWaveRouting(self.var.ChanQKin, SideflowChan, "main_channel")
- self.var.ChanM3Kin = self.var.ChanLength * self.var.ChannelAlpha * self.var.ChanQKin**self.var.Beta
- # Volume in channel at end of computation step
+ # KINEMATIC ROUTING
+ if option['InitLisflood'] or not (option['SplitRouting']):
+ # Kinematic routing uses the same variables for inputs and outputs. From upstream to downstream, input
+ # values are used to calculate inflow from upstream pixels, then routing is solved and discharge values
+ # are updated.
- self.var.ChanM3Kin = np.maximum(self.var.ChanM3Kin, 0.0)
- # Check for negative volumes at the end of computation step
- self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta)
- # Correct negative discharge at the end of computation step
+ self.KinRouting(SideflowChan)
+ # Takes:
+ # self.var.ChanQKin = outflow (at x+dx) from pixel channels at time t in input, at time t+dt in output [m3/sec] (instant)
+ # self.var.ChanQKinAvgDt = average outflow (at x+dx) from pixel channels at time t in input, at time t+dt in output [m3/sec] (average)
+ # self.var.ChanM3Kin = amount of water stored in the channel at time t in input, at time t+dt in output [m3/sec] (instant) [m3]
+ # SideflowChan = lateral inflow into the channel segment (cell) during dt [m3/channellength/sec]
- self.var.ChanQ = self.var.ChanQKin.copy()
- # at single kin. ChanQ is the same
-
- self.var.sumDisDay += self.var.ChanQ
- # Total channel storage [cu m], equal to ChanM3Kin
-
-
+ # Store results of kinematic routing
- else:
+ ChanQ = self.var.ChanQKin.copy()
+ # Outflow (x+dx) Q at the end of computation step t+dt for full section (instant)
- # ---- Double Routing ---------------
- # routing is split in two (virtual) channels)
+ ChanM3 = self.var.ChanM3Kin.copy()
+ # Channel storage V at the end of computation step t+dt for full section (instant)
- # Ad
- SideflowRatio = np.where((self.var.ChanM3Kin + self.var.Chan2M3Kin) > 0, self.var.ChanM3Kin/(self.var.ChanM3Kin+self.var.Chan2M3Kin), 0.0)
-
+ ChanQAvgDt = self.var.ChanQKinAvgDt.copy()
+ # Average Outflow (x+dx) Q during computation step dt for full section (instant)
-
- # CM ##################################
- # self.var.Sideflow1Chan = np.where(self.var.ChanM3Kin > self.var.M3Limit, SideflowRatio*SideflowChan, SideflowChan)
- # This is creating instability because ChanM3Kin can be < M3Limit between two routing sub-steps
- # TO BY REPLACED WITH THE FOLLOWING
- self.var.Sideflow1Chan = np.where((self.var.ChanM3Kin + self.var.Chan2M3Kin - self.var.Chan2M3Start) > self.var.M3Limit,
- SideflowRatio*SideflowChan, SideflowChan)
- #######################################
-
+ # SPLIT ROUTING - no InitLisflood
+ else:
+ # Split routing uses the same variables for inputs and outputs. From upstream to downstream, input
+ # values are used to calculate inflow from upstream pixels, then routing is solved and discharge values
+ # are updated.
+
+ self.SplitRouting(SideflowChan)
+ # self.var.ChanQKin: outflow (x+dx) from the main channel at time t in input, at time t+dt in output (instant) [m3/s]
+ # self.var.ChanM3Kin: amount of water stored in the main channel at time t in input, at time t+dt in output (instant) [m3]
+ # self.var.Chan2QKin: outflow (x+dx) from the second line channel (including QLimit) at time t in input, at time t+dt in output (instant) [m3/s]
+ # self.var.Chan2M3Kin: amount of water stored in the second line channel at time t in input, at time t+dt in output (including self.var.Chan2M3Start) (instant) [m3]
+ # self.var.CrossSection2Area: cross-section area above bankfull for second line of routing at time t in input, at time t+dt in output (instant) [m2]
+
+ # Combine the two lines of routing together
+ # Store results of split routing
+
+ ChanQ = np.maximum(self.var.ChanQKin + self.var.Chan2QKin - self.var.QLimit, 0)
+ # (real) total outflow (at x+dx) at time t+dt end of step for the full cross-section (instant)
+ # Main channel routing and above bankfull routing from second line of routing
+
+ ChanM3 = self.var.ChanM3Kin + self.var.Chan2M3Kin - self.var.Chan2M3Start
+ # Total channel storage [m3] = Volume in main channel (ChanM3Kin) + volume above bankfull in second line (Chan2M3Kin - Chan2M3Start)
+ # Total channel storage V at the end of computation step t+dt for full section (instant)
+
+ # cmcheck
+ ChanQAvgDt = np.maximum(self.var.ChanQKinAvgDt + self.var.Chan2QKinAvgDt - self.var.QLimit, 0)
+ # (real) total outflow (at x+dx) at time t+dt end of step for the full cross-section (instant)
+ # Main channel routing and above bankfull routing from second line of routing
+
+ # MUSKINGUM-CUNGE-TODINI ROUTING - no InitLisflood
+ if option['MCTRouting']:
+ # MCT routing
+ # This is calculated for MCT grid cell only but takes the output of kinematic or split routing.
+ # First, Kinematic/Split routing is solved on all pixels (including MCT pixels) then results are updated
+ # for the MCT pixels.
+ # Sideflow contribution to MCT grid cells expressed in [m3/s]
+ SideflowChanMCT = np.where(self.var.IsChannelMCT, SideflowChanM3 * self.var.InvDtRouting, 0) #Ql
+
+ # Grab outflow at the end of the previous routing step t for all pixels) - current state of the MCT pixel
+ ChanQ_0 = self.var.ChanQ.copy() # Outflow (x+dx) at time t (end of previous routing step) (instant) -> used to calc q00
+ ChanM3_0 = self.var.ChanM3.copy() # Channel storage at time t (end of previous routing step) (instant) V00
+
+ # Grab outflow at the end of step t+dt for Kinematic/Split (all pixels) -> this is used to calc contribution from kin/split pixels to MCT pixels
+ # Put results of Kinematic/Split routing at the end of the routing step dt into MCT points
+ self.var.ChanQ = ChanQ # -> used to calc q01
+ self.var.ChanM3 = ChanM3
+ self.var.ChanQAvgDt = ChanQAvgDt # -> used to calc q0m
+
+ # Solve MCT routing and update current state at MCT pixels
+ self.mct_river_router.routing(
+ ChanQ_0, # -> used to calc q00
+ ChanM3_0, # V00
+ SideflowChanMCT, # sideflow for MCT pixels [m3/s]
+ # THESE ARE BOTH INPUTS AND OUTPUTS
+ self.var.ChanQ, # -> in input: used to calc q01; in output: q11 outflow (x+dx) at time t+dt (instant)
+ self.var.ChanQAvgDt, # -> in input: used to calculate q0m; in output: q1m outflow (x+dx) at time t+dt (average)
+ self.var.PrevCm0, # Courant number in input: at time t; in output: at time t+dt
+ self.var.PrevDm0, # Reynolds number in input: at time t; in output: at time t+dt
+ self.var.ChanM3 # Channel storage volume. In input: at time t V00; in output: at time t+dt V11
+ )
- self.var.Sideflow1Chan = np.where(np.abs(SideflowChan) < 1e-7, SideflowChan, self.var.Sideflow1Chan)
- # too small values are avoided
- Sideflow2Chan = SideflowChan - self.var.Sideflow1Chan
-
- Sideflow2Chan = Sideflow2Chan + self.var.Chan2QStart * self.var.InvChanLength # originale
- # as kinematic wave gets slower with less water
- # a constant amount of water has to be added
- # -> add QLimit discharge
-
- # --- Main channel routing ---
- self.river_router.kinematicWaveRouting(self.var.ChanQKin, self.var.Sideflow1Chan, "main_channel")
- self.var.ChanM3Kin = self.var.ChanLength * self.var.ChannelAlpha * self.var.ChanQKin ** self.var.Beta
-
- self.var.ChanM3Kin = np.maximum(self.var.ChanM3Kin, 0.0)
- # Check for negative volumes at the end of computation step
- self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta)
- # Correct negative discharge at the end of computation step
-
-
- # --- Floodplains routing ---
- self.river_router.kinematicWaveRouting(self.var.Chan2QKin, Sideflow2Chan, "floodplains")
- self.var.Chan2M3Kin = self.var.ChanLength * self.var.ChannelAlpha2 * self.var.Chan2QKin ** self.var.Beta
+ else:
+ # Store results of kinematic/split routing in the general variables
+ self.var.ChanQ = ChanQ
+ self.var.ChanM3 = ChanM3
+ self.var.ChanQAvgDt = ChanQAvgDt
+
+ # sum of average outflow Q on model sub-steps to give accumulated total outflow volume
+ self.var.sumDisDay += self.var.ChanQAvgDt
- diffM3 = self.var.Chan2M3Kin - self.var.Chan2M3Start
- self.var.Chan2M3Kin = np.where(diffM3 < 0.0, self.var.Chan2M3Start, self.var.Chan2M3Kin)
- # Check for negative volume in second line of routing at the end of routing substep
+ ### end of river routing calculation
- self.var.CrossSection2Area = (self.var.Chan2M3Kin - self.var.Chan2M3Start) * self.var.InvChanLength
- # Compute cross-section for second line of routing
-
- self.var.Chan2QKin = (self.var.Chan2M3Kin * self.var.InvChanLength * self.var.InvChannelAlpha2) ** (self.var.InvBeta)
- # Correct negative discharge at the end of computation step in second line
-
-
- self.var.ChanQ = np.maximum(self.var.ChanQKin + self.var.Chan2QKin - self.var.QLimit, 0.0)
-
- # Superposition Kinematic
- # Main channel routing and floodplains routing
- #print(np.sum(self.var.sumDisDay))
- self.var.sumDisDay_NOTlast=self.var.sumDisDay.copy()
- self.var.sumDisDay += self.var.ChanQ
- # ----------End splitrouting-------------------------------------------------
-
- ###
-
# ---- Uncomment lines 603-635 in order to compute the mass balance error within the routing module for the options (i) initial run or (ii) split routing off ----
- '''
- if option['repMBTs']:
- if option['InitLisflood'] or (not(option['SplitRouting'])):
- if NoRoutingExecuted == (self.var.NoRoutSteps-1):
- StorageStep = self.var.ChanM3Kin.copy()
- ChanQAvgR = self.var.sumDisDay/self.var.NoRoutSteps
- sum1=ChanQAvgR.copy()
- sum1[self.var.AtLastPointC == 0] = 0
- OutStep = np.take(np.bincount(self.var.Catchments,weights=sum1 * self.var.DtSec),self.var.Catchments)
-
- maskinfo = MaskInfo.instance()
- DisStructureR = maskinfo.in_zero()
- DischargeM3StructuresR = maskinfo.in_zero()
- if not option['InitLisflood']:
- if option['simulateReservoirs']:
- sum1 =self.var.ChanQ.copy()
- StorageStep = StorageStep + self.var.ReservoirStorageM3.copy()
- DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
- DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureR), self.var.Catchments)
- DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
- if not option['InitLisflood']:
- if option['simulateLakes']:
- sum1 =self.var.ChanQ.copy()
- StorageStep = StorageStep + self.var.LakeStorageM3Balance.copy()
- DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
- DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureR), self.var.Catchments)
- maskinfo = MaskInfo.instance()
- DisLake = maskinfo.in_zero()
- np.put(DisLake, self.var.LakeIndex, 0.5 * self.var.LakeInflowCC * self.var.DtRouting)
- DischargeM3Lake = np.take(np.bincount(self.var.Catchments, weights=DisLake),self.var.Catchments)
- DischargeM3StructuresR += DischargeM3Lake
- DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
- # Total Mass Balance Error in m3 per catchment for Initial Run OR Split Routing OFF
- MB =- np.sum(StorageStep)+np.sum(self.var.StorageStepINIT) - OutStep[0] -DischargeM3StructuresR[0] +self.var.AddedTRUN
- self.var.StorageStepINIT= np.sum(StorageStep) + DischargeM3StructuresR[0]
- '''
- if option['repMBTs']:
- if (not(option['InitLisflood'])) and (option['SplitRouting']):
- # compute the mass balance at the last of the sub-routing steps in order to account for the contributions of lakes and reservoirs
- if NoRoutingExecuted == (self.var.NoRoutSteps-1):
- ChanQAvgSR = self.var.sumDisDay/self.var.NoRoutSteps
- sum1=ChanQAvgSR.copy()
- sum1[self.var.AtLastPointC == 0] = 0
- OutStep = np.take(np.bincount(self.var.Catchments,weights=sum1 * self.var.DtSec),self.var.Catchments)
-
- StorageStep=[]
- StorageStep= self.var.ChanM3Kin.copy()+self.var.Chan2M3Kin.copy()-self.var.Chan2M3Start.copy()
-
-
- maskinfo = MaskInfo.instance()
- DisStructureSR = maskinfo.in_zero()
- DischargeM3StructuresR = maskinfo.in_zero()
-
- if option['simulateReservoirs']:
- sum1=[]
- sum1 =self.var.ChanQ.copy()
- StorageStep = StorageStep + self.var.ReservoirStorageM3.copy()
- DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
- DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureSR), self.var.Catchments)
- DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
-
- if option['simulateLakes']:
- sum1 =self.var.ChanQ.copy()
- StorageStep = StorageStep + self.var.LakeStorageM3Balance.copy()
- DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
- DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureSR), self.var.Catchments)
- DisLake = maskinfo.in_zero()
- np.put(DisLake, self.var.LakeIndex, 0.5 * self.var.LakeInflowCC * self.var.DtRouting)
- DischargeM3Lake = np.take(np.bincount(self.var.Catchments, weights=DisLake),self.var.Catchments)
- DischargeM3StructuresR += DischargeM3Lake
-
- DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
-
- # Mass Balance Error due to the Split Routing module
- StorageStep1=np.take(np.bincount(self.var.Catchments, weights=StorageStep), self.var.Catchments)
-
- self.var.MBErrorSplitRoutingM3 = - StorageStep1 + self.var.StorageStepINIT - OutStep - DischargeM3StructuresR + self.var.AddedTRUN
- # Discharge error at the outlet pointt [m3/s]
- QoutCorrection=self.var.MBErrorSplitRoutingM3/self.var.DtRouting
- QoutCorrection[self.var.AtLastPointC == 0] = 0
- self.var.OutletDischargeErrorSplitRouting = np.take(np.bincount(self.var.Catchments,weights=QoutCorrection),self.var.Catchments)
-
- self.var.StorageStepINIT= StorageStep1.copy()+DischargeM3StructuresR
-
-
-
- TotalCrossSectionArea = np.maximum(self.var.ChanM3Kin*self.var.InvChanLength,0.01)
-
+ #'''
+ # option['repMBTs']=True
+ if option['repMBTs']:
+ if option['InitLisflood'] or (not(option['SplitRouting'])):
+ # Kinematic routing and MCT
+ if NoRoutingExecuted == (self.var.NoRoutSteps-1):
+ # StorageStep = self.var.ChanM3Kin.copy()
+ StorageStep = self.var.ChanM3.copy()
+ # Water storage at t+dt end of routing step: rivers channels
+ # Using ChanM3 so it's OK for both MCT and KIN
+
+ ChanQAvgR = self.var.sumDisDay/self.var.NoRoutSteps
+ # average (of instantaneous) outflow (x+dx) at t+dt end of routing step
+ sum1=ChanQAvgR.copy()
+ sum1[self.var.AtLastPointC == 0] = 0
+ OutStepM3 = np.take(np.bincount(self.var.Catchments,weights=sum1 * self.var.DtSec),self.var.Catchments)
+ # average outflow volume (x+dx) volume at t+dt
+
+ maskinfo = MaskInfo.instance()
+ DisStructureR = maskinfo.in_zero()
+ DischargeM3StructuresR = maskinfo.in_zero()
+
+ if not option['InitLisflood']:
+ if option['simulateReservoirs']:
+ sum1 =self.var.ChanQ.copy()
+ StorageStep = StorageStep + self.var.ReservoirStorageM3.copy()
+ # DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+ DisStructureR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
+ DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureR), self.var.Catchments)
+ DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
+
+ if not option['InitLisflood']:
+ if option['simulateLakes']:
+ sum1 =self.var.ChanQ.copy()
+ StorageStep = StorageStep + self.var.LakeStorageM3Balance.copy()
+ # DisStructureR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+ DisStructureR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
+ DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureR), self.var.Catchments)
+ maskinfo = MaskInfo.instance()
+ DisLake = maskinfo.in_zero()
+ np.put(DisLake, self.var.LakeIndex, 0.5 * self.var.LakeInflowCC * self.var.DtRouting)
+ DischargeM3Lake = np.take(np.bincount(self.var.Catchments, weights=DisLake),self.var.Catchments)
+ DischargeM3StructuresR += DischargeM3Lake
+ DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
+
+ # Total Mass Balance Error in m3 per catchment for Initial Run OR Kinematic routing (Split Routing OFF)
+ # MB =-np.sum(StorageStep)+np.sum(self.var.StorageStepINIT) - OutStepM3[0] -DischargeM3StructuresR[0] +self.var.AddedTRUN
+
+ # # cmcheck
+ # if MB.any() > 1.e-12:
+ # print('Mass balance error MB > 1.e-12')
+
+
+ self.var.StorageStepINIT= np.sum(StorageStep) + DischargeM3StructuresR[0]
+ #'''
+
+ # ---- Uncomment lines in order to compute the mass balance error within the routing module for the options split routing ----
+ #'''
+ if option['repMBTs']:
+ if (not(option['InitLisflood'])) and (option['SplitRouting']):
+ # SplitRouting
+ # compute the mass balance at the last of the sub-routing steps in order to account for the contributions of lakes and reservoirs
+ if NoRoutingExecuted == (self.var.NoRoutSteps-1):
+ ChanQAvgSR = self.var.sumDisDay/self.var.NoRoutSteps #self.var.ChanQ
+ sum1=ChanQAvgSR.copy()
+ sum1[self.var.AtLastPointC == 0] = 0
+ OutStep = np.take(np.bincount(self.var.Catchments,weights=sum1 * self.var.DtSec),self.var.Catchments)
+
+ StorageStep=[]
+ # StorageStep= self.var.ChanM3Kin.copy()+self.var.Chan2M3Kin.copy()-self.var.Chan2M3Start.copy()
+ StorageStep = self.var.ChanM3
+
+
+ maskinfo = MaskInfo.instance()
+ DisStructureSR = maskinfo.in_zero()
+ DischargeM3StructuresR = maskinfo.in_zero()
+
+ if option['simulateReservoirs']:
+ sum1=[]
+ sum1 =self.var.ChanQ.copy()
+ StorageStep = StorageStep + self.var.ReservoirStorageM3.copy()
+ # DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+ DisStructureSR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
+ DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureSR), self.var.Catchments)
+ DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
+
+ if option['simulateLakes']:
+ sum1 =self.var.ChanQ.copy()
+ StorageStep = StorageStep + self.var.LakeStorageM3Balance.copy()
+ # DisStructureSR = np.where(self.var.IsUpsOfStructureKinematicC, sum1 * self.var.DtRouting, 0)
+ DisStructureSR = np.where(self.var.IsUpsOfStructureChanC, sum1 * self.var.DtRouting, 0)
+ DischargeM3StructuresR = np.take(np.bincount(self.var.Catchments, weights=DisStructureSR), self.var.Catchments)
+ DisLake = maskinfo.in_zero()
+ np.put(DisLake, self.var.LakeIndex, 0.5 * self.var.LakeInflowCC * self.var.DtRouting)
+ DischargeM3Lake = np.take(np.bincount(self.var.Catchments, weights=DisLake),self.var.Catchments)
+ DischargeM3StructuresR += DischargeM3Lake
+
+ DischargeM3StructuresR -= self.var.DischargeM3StructuresIni
+
+ # Mass Balance Error due to the Split Routing module
+ StorageStep1=np.take(np.bincount(self.var.Catchments, weights=StorageStep), self.var.Catchments)
+
+ self.var.MBErrorSplitRoutingM3 = - StorageStep1 + self.var.StorageStepINIT - OutStep - DischargeM3StructuresR + self.var.AddedTRUN
+ # Discharge error at the outlet pointt [m3/s]
+ QoutCorrection=self.var.MBErrorSplitRoutingM3/self.var.DtRouting
+ QoutCorrection[self.var.AtLastPointC == 0] = 0
+ self.var.OutletDischargeErrorSplitRouting = np.take(np.bincount(self.var.Catchments,weights=QoutCorrection),self.var.Catchments)
+
+ self.var.StorageStepINIT= StorageStep1.copy()+DischargeM3StructuresR
+ #'''
+
+ # Legacy
+ # TotalCrossSectionArea = np.maximum(self.var.ChanM3Kin * self.var.InvChanLength, 0.01)
+ TotalCrossSectionArea = np.maximum(self.var.ChanM3 * self.var.InvChanLength, 0.01)
self.var.FlowVelocity = np.minimum(self.var.ChanQKin/TotalCrossSectionArea, 0.36*self.var.ChanQKin**0.24)
# Channel velocity (m/s); dividing Q (m3/s) by CrossSectionArea (m2)
# avoid extreme velocities by using the Wollheim 2006 equation
# assume 0.1 for upstream areas (outside ChanLdd)
- self.var.FlowVelocity *= np.minimum(np.sqrt(self.var.PixelArea)*self.var.InvChanLength,1);
+ self.var.FlowVelocity *= np.minimum(np.sqrt(self.var.PixelArea)*self.var.InvChanLength,1)
# reduction for sinuosity of channels
- self.var.TravelDistance=self.var.FlowVelocity*self.var.DtSec;
+ self.var.TravelDistance=self.var.FlowVelocity*self.var.DtSec
# if flow is fast, Traveltime=1, TravelDistance is high: Pixellength*DtSec
# if flow is slow, Traveltime=DtSec then TravelDistance=PixelLength
# maximum set to 30km/day for 5km cell, is at DtSec/Traveltime=6, is at Traveltime 0,
+ self.var.ChanM3Kin / (self.var.ChanM3Kin + self.var.Chan2M3Kin), 0.0)
+ # Split sideflow between the two lines of routing
+
+ # self.var.Sideflow1Chan = np.where(self.var.ChanM3Kin > self.var.M3Limit, SideflowRatio*SideflowChan, SideflowChan)
+ # This is creating instability because ChanM3Kin can be < M3Limit between two routing sub-steps
+
+ self.var.Sideflow1Chan = np.where(
+ (self.var.ChanM3Kin + self.var.Chan2M3Kin - self.var.Chan2M3Start) > self.var.M3Limit,
+ SideflowRatio * SideflowChan, SideflowChan)
+ # sideflow to the main channel
+
+ self.var.Sideflow1Chan = np.where(np.abs(SideflowChan) < 1e-7, SideflowChan, self.var.Sideflow1Chan)
+ # too small values are avoided
+ Sideflow2Chan = SideflowChan - self.var.Sideflow1Chan
+ # sideflow to the 'second line' channel
+
+ Sideflow2Chan = Sideflow2Chan + self.var.Chan2QStart * self.var.InvChanLength
+ # as kinematic wave gets slower with less water a constant amount of water has to be added
+ # -> add QLimit discharge
+
+ # --- Main channel routing ---
+
+ self.river_router.kinematicWaveRouting(self.var.ChanQKinAvgDt, self.var.ChanQKin, self.var.Sideflow1Chan, "main_channel")
+ # sef.var.ChanQKin is outflow (x+dx) from main channel at time t in input and at time t+dt in output (instant)
+ # sef.var.ChanQKinAvgDt is average outflow (x+dx) from main channel during time dt (average)
+ self.var.ChanM3Kin = self.var.ChanLength * self.var.ChannelAlpha * self.var.ChanQKin ** self.var.Beta
+ # Volume in main channel at end of computation step (at t+dt) (instant)
+
+ self.var.ChanM3Kin = np.maximum(self.var.ChanM3Kin, 0.0)
+ # Check for negative volumes at the end of computation step in main channel
+ # Volume in main channel at t+dt
+ self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta)
+ # Correct negative discharge in main channel at the end of computation step
+ # Outflow (x+dx) at t+dt (instant)
+
+ # --- Second line channel routing (same channel geometry as main channel but increased Manning coeff) ---
+
+ self.river_router.kinematicWaveRouting(self.var.Chan2QKinAvgDt, self.var.Chan2QKin, Sideflow2Chan, "floodplains")
+ # sef.var.Chan2QKin is (virtual) total outflow (x+dx) at time t in input and at time t+dt in output (instant)
+ # (same channel geometry but using increased Manninig coeff)
+ # sef.var.Chan2QKinAvgDt is (virtual) average total outflow (x+dx) from second line channel during time dt (average)
+ self.var.Chan2M3Kin = self.var.ChanLength * self.var.ChannelAlpha2 * self.var.Chan2QKin ** self.var.Beta
+ # Total (virtual) volume of water in river channel when second routing line is active
+ # (same channel geometry but using increased Manninig coeff)
+
+ FldpM3 = self.var.Chan2M3Kin - self.var.Chan2M3Start
+ # Water volume over bankfull
+ self.var.Chan2M3Kin = np.where(FldpM3 < 0.0, self.var.Chan2M3Start, self.var.Chan2M3Kin)
+ # Check for negative volume over bankfull at the end of routing step for second line of routing
+ # Total (virtual) volume for second line of routing cannot be smaller than the bankfull volume for the second line
+ # of routing aka the bankfull volume calculated with increased Manning coeff (Chan2M3Start)
+
+ self.var.CrossSection2Area = (self.var.Chan2M3Kin - self.var.Chan2M3Start) * self.var.InvChanLength
+ # Compute cross-section area above bankfull for second line of routing
+
+ self.var.Chan2QKin = (self.var.Chan2M3Kin * self.var.InvChanLength * self.var.InvChannelAlpha2) ** (self.var.InvBeta)
+ # (virtual) total outflow (at x + dx) at time t + dt (instant) for second line of routing (using increased Manninig coeff)
+ # Correct negative discharge at the end of computation step in second line
+
+ # FldpQKin = self.var.Chan2QKin - self.var.QLimit
+ # Outflow at t+dt from above bankfull only
+
+ def upstream(self, var):
+ var_pcr = decompress(var)
+ var_ups = compressArray(upstream(self.var.LddChan, var_pcr))
+ return var_ups
+
+ def compress_mct(self, var):
+ """Compress to mct array.
+ For any array (var) with all catchment pixels, it finds the MCT pixels (x) and
+ reduces the dimension of the array (y).
+ :param var: array (dimension is the number of cells in the basin)
+ :return:
+ y: same as input array (var) but with only MCT pixels
+ """
+ x = np.ma.masked_where(self.var.IsChannelKinematic, var) # mask kinematic pixels
+ y = x.compressed()
+
+ return y
diff --git a/src/lisflood/hydrological_modules/structures.py b/src/lisflood/hydrological_modules/structures.py
index 60ced355..44dbeaa8 100755
--- a/src/lisflood/hydrological_modules/structures.py
+++ b/src/lisflood/hydrological_modules/structures.py
@@ -43,19 +43,41 @@ def __init__(self, structures_variable):
def initial(self):
""" initial part of the structures module
"""
- self.var.LddStructuresKinematic = self.var.LddKinematic
+ self.var.LddStructuresKinematic = self.var.LddKinematic #pcr map
+ LddStructuresKinematicNp = compressArray(self.var.LddStructuresKinematic)
+ # Unmodified version of LddKinematic is needed to connect inflow and outflow points
+ # of each structure (called LddStructuresKinematic now)
+
+ self.var.LddStructuresChan = self.var.LddChan #pcr map
+ LddStructuresChanNp = compressArray(self.var.LddStructuresChan)
+ # Unmodified version of LddChan is needed to connect inflow and outflow points
+ # of each structure (called LddStructuresChan now)
+
settings = LisSettings.instance()
option = settings.options
if not option['InitLisflood']:
# not done in Init Lisflood
- IsUpsOfStructureKinematic = downstream(
+ IsUpsOfStructureKinematic = downstream( #pcr map
self.var.LddKinematic,
cover(boolean(decompress(self.var.IsStructureKinematic)), boolean(0))
)
- # Get all pixels just upstream of kinematic structure locations
- self.var.IsUpsOfStructureKinematicC = compressArray(IsUpsOfStructureKinematic)
- # Unmodified version of LddKinematic is needed to connect inflow and outflow points
- # of each structure (called LddStructuresKinematic now)
- self.var.LddKinematic = lddrepair(ifthenelse(IsUpsOfStructureKinematic, 5, self.var.LddKinematic))
- # Cells just upstream of each structure are treated as pits in the kinematic wave
- # channel routing
+ # Downstream assigns to result the expression value of the neighbouring downstream cell
+ # Over is used to cover missing values on an expression with values taken from one or more different expression(s)
+ # Decompress is numpy2pcr
+ # Find location of pixels immediately upstream of a structure on the LddKinematic
+
+ IsUpsOfStructureChan = downstream( #pcr map
+ self.var.LddChan,
+ cover(boolean(decompress(self.var.IsStructureChan)), boolean(0))
+ )
+ # Find location of pixels immediately upstream of a structure on the LddChan
+
+ self.var.IsUpsOfStructureKinematicC = compressArray(IsUpsOfStructureKinematic) #np compressed array
+ # Location of pixels immediately upstream of a structure on the LddKinematic
+ self.var.IsUpsOfStructureChanC = compressArray(IsUpsOfStructureChan) #np compressed array
+ # Location of pixels immediately upstream of a structure on the LddChan
+
+ self.var.LddKinematic = lddrepair(ifthenelse(IsUpsOfStructureKinematic, 5, self.var.LddKinematic)) #pcr map
+ # Update LddKinematic by adding a pit in the pixel immediately upstream of a structure
+ self.var.LddChan = lddrepair(ifthenelse(IsUpsOfStructureChan, 5, self.var.LddChan)) #pcr map
+ # Update LddChan by adding a pit in the pixel immediately upstream of a structure
diff --git a/src/lisflood/hydrological_modules/surface_routing.py b/src/lisflood/hydrological_modules/surface_routing.py
index 00d54faa..c736bee2 100644
--- a/src/lisflood/hydrological_modules/surface_routing.py
+++ b/src/lisflood/hydrological_modules/surface_routing.py
@@ -94,6 +94,13 @@ def initial(self):
self.var.OFQOther = ((self.var.OFM3Other * self.var.InvPixelLength * self.var.InvOFAlpha.values[self.var.dim_runoff[1].index('Other')])**(self.var.InvBeta)).astype(float)
self.var.OFQForest = ((self.var.OFM3Forest * self.var.InvPixelLength * self.var.InvOFAlpha.values[self.var.dim_runoff[1].index('Forest')])**(self.var.InvBeta)).astype(float)
+ # cmcheck
+ # do I need to initialise from self.var.OFQDirect?
+ # Initial average overland discharge [m3 s-1]
+ self.var.OFQDirectAvg = maskinfo.in_zero()
+ self.var.OFQOtherAvg = maskinfo.in_zero()
+ self.var.OFQForestAvg = maskinfo.in_zero()
+
def initialSecond(self):
""" 2nd initialisation part of the surface routing module:
to be called after all needed parameters are set (PixelLength, LddToChan, Alpha...)
@@ -148,9 +155,10 @@ def dynamic(self):
SideflowOther = np.sum(self.var.SurfaceRunSoil.values[ilusevalues],self.var.SurfaceRunSoil.dims.index("landuse")) * self.var.MMtoM3 * self.var.InvPixelLength * self.var.InvDtSec
SideflowForest = self.var.SurfaceRunSoil.values[self.var.epic_settings.soil_uses.index('Forest')] * self.var.MMtoM3 * self.var.InvPixelLength * self.var.InvDtSec
# All surface runoff that is generated during current time step added as side flow [m3/s/m pixel-length]
- self.direct_surface_router.kinematicWaveRouting(self.var.OFQDirect, SideflowDirect)
- self.other_surface_router.kinematicWaveRouting(self.var.OFQOther, SideflowOther)
- self.forest_surface_router.kinematicWaveRouting(self.var.OFQForest, SideflowForest)
+
+ self.direct_surface_router.kinematicWaveRouting(self.var.OFQDirectAvg,self.var.OFQDirect, SideflowDirect)
+ self.other_surface_router.kinematicWaveRouting(self.var.OFQOtherAvg, self.var.OFQOther, SideflowOther)
+ self.forest_surface_router.kinematicWaveRouting(self.var.OFQForestAvg, self.var.OFQForest, SideflowForest)
# to PCRASTER
@@ -191,15 +199,20 @@ def dynamic(self):
self.var.OFM3Direct = self.var.PixelLength * self.var.OFAlpha.values[self.var.dim_runoff[1].index('Direct')] * self.var.OFQDirect**self.var.Beta
self.var.OFM3Other = self.var.PixelLength * self.var.OFAlpha.values[self.var.dim_runoff[1].index('Other')] * self.var.OFQOther**self.var.Beta
self.var.OFM3Forest = self.var.PixelLength * self.var.OFAlpha.values[self.var.dim_runoff[1].index('Forest')] * self.var.OFQForest**self.var.Beta
+ # Overland flow storage for different components at the end of routing step [m3]
self.var.Qall = self.var.OFQDirect + self.var.OFQOther + self.var.OFQForest
+ self.var.QallAvg = self.var.OFQDirectAvg + self.var.OFQOtherAvg + self.var.OFQForestAvg
self.var.M3all = self.var.OFM3Direct + self.var.OFM3Other + self.var.OFM3Forest
- # Overland flow storage [m3]
-
- self.var.OFToChanM3 = np.where(self.var.IsChannel, self.var.Qall * self.var.DtSec, 0)
+ # Total overland flow storage [m3]
+ # cmcheck
+ # this should be calculated using the average flow
+ self.var.OFToChanM3 = np.where(self.var.IsChannel, self.var.QallAvg * self.var.DtSec, 0)
+ # self.var.OFToChanM3 = np.where(self.var.IsChannel, self.var.Qall * self.var.DtSec, 0)
# Overland flow in channel pixels (in [m3])is added to channel
+
self.var.WaterDepth = self.var.M3all * self.var.M3toMM
# Update water depth [mm]
diff --git a/src/lisflood/hydrological_modules/waterabstraction.py b/src/lisflood/hydrological_modules/waterabstraction.py
index dfd348df..b3636c6c 100644
--- a/src/lisflood/hydrological_modules/waterabstraction.py
+++ b/src/lisflood/hydrological_modules/waterabstraction.py
@@ -159,14 +159,16 @@ def initial(self):
pitWuse2 = ifthen(pitWuseMax == self.var.UpArea, WUseRegion)
# search outlets in the inland water regions by using the maximum upstream area as criterium
- pitWuse3 = downstream(self.var.LddStructuresKinematic, WUseRegion)
+ # pitWuse3 = downstream(self.var.LddStructuresKinematic, WUseRegion)
+ pitWuse3 = downstream(self.var.LddStructuresChan, WUseRegion)
pitWuse3b = ifthen(pitWuse3 != WUseRegion, WUseRegion)
# search point where ldd leaves a water region
pitWuse = cover(pitWuse1b, pitWuse2, pitWuse3b, nominal(0))
# join all sources of pits
- LddWaterRegion = lddrepair(ifthenelse(pitWuse == 0, self.var.LddStructuresKinematic, 5))
+ # LddWaterRegion = lddrepair(ifthenelse(pitWuse == 0, self.var.LddStructuresKinematic, 5))
+ LddWaterRegion = lddrepair(ifthenelse(pitWuse == 0, self.var.LddStructuresChan, 5))
# create a Ldd with pits at every water region outlet
# this results in a interrupted ldd, so water cannot be transfered to the next water region
lddC = compressArray(LddWaterRegion)
@@ -185,8 +187,10 @@ def initial(self):
# outflowpoints to calculate upstream inflow for balances and Water Exploitation Index
# both inland outflowpoints to downstream subbasin, and coastal outlets
+ # WaterRegionInflow1 = boolean(
+ # upstream(self.var.LddStructuresKinematic, cover(scalar(self.var.WaterRegionOutflowPoints), 0)))
WaterRegionInflow1 = boolean(
- upstream(self.var.LddStructuresKinematic, cover(scalar(self.var.WaterRegionOutflowPoints), 0)))
+ upstream(self.var.LddStructuresChan, cover(scalar(self.var.WaterRegionOutflowPoints), 0)))
self.var.WaterRegionInflowPoints = ifthen(WaterRegionInflow1, boolean(1))
# inflowpoints to calculate upstream inflow for balances and Water Exploitation Index
else:
@@ -478,8 +482,13 @@ def dynamic(self):
# 10.1 Withdrawal required from channels (CH)
areatotal_withdrawal_CH_required_M3 = np.maximum(areatotal_withdrawal_SW_required - self.var.areatotal_withdrawal_LakRes_actual_M3, 0.)
# 10.2 Max abstractable volumes from channels, accounting for e-flow constraint
+
+ # PixelAvailableWaterFromChannelsM3 = np.maximum(
+ # self.var.ChanM3Kin - self.var.EFlowThreshold * self.var.DtSec, maskinfo.in_zero()) ### QUESTION! # * (1 - self.var.WUsePercRemain) THIS BIT IS COMMENTED FOR CONSISTENCY WITH EPIC, UNCOMMENT BEFORE THE FINAL MERGE
PixelAvailableWaterFromChannelsM3 = np.maximum(
- self.var.ChanM3Kin - self.var.EFlowThreshold * self.var.DtSec, maskinfo.in_zero()) ### QUESTION! # * (1 - self.var.WUsePercRemain) THIS BIT IS COMMENTED FOR CONSISTENCY WITH EPIC, UNCOMMENT BEFORE THE FINAL MERGE
+ self.var.ChanM3 - self.var.EFlowThreshold * self.var.DtSec, maskinfo.in_zero()) ### QUESTION! # * (1 - self.var.WUsePercRemain) THIS BIT IS COMMENTED FOR CONSISTENCY WITH EPIC, UNCOMMENT BEFORE THE FINAL MERGE
+ # using total water storage in river channel
+
self.var.AreaTotalAvailableWaterFromChannelsM3 = np.maximum(
np.take(np.bincount(self.var.WUseRegionC, weights=PixelAvailableWaterFromChannelsM3),
self.var.WUseRegionC),maskinfo.in_zero())
diff --git a/src/lisflood/hydrological_modules/waterbalance.py b/src/lisflood/hydrological_modules/waterbalance.py
index 3762f1bd..65d39b71 100755
--- a/src/lisflood/hydrological_modules/waterbalance.py
+++ b/src/lisflood/hydrological_modules/waterbalance.py
@@ -57,6 +57,7 @@ def initial(self):
# not done in Init Lisflood
ChannelInitM3 = self.var.ChanIniM3
+ # Initial total water volume in river channels [m3]
if not(option['InitLisflood']): # only with no InitLisflood
if option['simulateLakes']:
ChannelInitM3 += self.var.LakeStorageIniM3
@@ -89,7 +90,8 @@ def initial(self):
#DisStructure = cover(ifthen(
# self.var.IsUpsOfStructureKinematic, self.var.ChanQ * self.var.DtRouting), scalar(0.0))
- DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+ # DisStructure = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+ DisStructure = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
# Discharge upstream of structure locations (coded as pits) in [m3/time step]
# Needed for mass balance error calculations (see comment to calculation of WaterInit below)
@@ -202,6 +204,8 @@ def dynamic(self):
WaterStored_nextstep += np.take(np.bincount(self.var.Catchments,weights=HillslopeStoredM3),self.var.Catchments)
# Total water stored [m3]
+
+
# This goes out:
HillslopeOutM3 = (self.var.TaWB + self.var.TaInterceptionWB + self.var.ESActWB + self.var.GwLossWB) * self.var.MMtoM3
# Water that goes out of the system at the hillslope level [m3]
@@ -209,8 +213,12 @@ def dynamic(self):
##sum1 = self.var.sumDis.copy()
sum1 = self.var.ChanQAvg.copy()
+ # averge outflow during model time step
sum1[self.var.AtLastPointC == 0] = 0
WaterOut = np.take(np.bincount(self.var.Catchments,weights=sum1 * self.var.DtSec),self.var.Catchments)
+ # Water that goes out of the system at the channels level [m3]
+ # (rivers flow)
+
WaterOut += np.take(np.bincount(self.var.Catchments,weights=HillslopeOutM3),self.var.Catchments)
#WaterOut = areatotal(cover(ifthen(
@@ -232,7 +240,8 @@ def dynamic(self):
# for this double counting)
# new 12.11.09 PB
# added cumulative transmission loss
- DisStru = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+ # DisStru = np.where(self.var.IsUpsOfStructureKinematicC, self.var.ChanQ * self.var.DtRouting, 0)
+ DisStru = np.where(self.var.IsUpsOfStructureChanC, self.var.ChanQ * self.var.DtRouting, 0)
DischargeM3Structures = np.take(np.bincount(self.var.Catchments, weights=DisStru), self.var.Catchments)
# on the last time step lakes and reservoirs calculated with the previous routing results
@@ -258,6 +267,8 @@ def dynamic(self):
# Discharge just upstream of structure locations (coded as pits) in [cu m / time step]
# Needed for mass balance error calculations, because of double counting of structure
# storage and water in the channel.
+
+
# Mass balance:
self.var.MBError = self.var.WaterInit + WaterIn - WaterStored - WaterOut - DischargeM3Structures
# Totl mass balance error per catchment [cu m]. Mass balance error is computed for each computational time step.
diff --git a/src/lisfloodSettings_reference.xml b/src/lisfloodSettings_reference.xml
index 37b95aff..0053bfcc 100644
--- a/src/lisfloodSettings_reference.xml
+++ b/src/lisfloodSettings_reference.xml
@@ -48,6 +48,7 @@ You can use builtin path variables in this template and reference to other paths
+
# use inflow data
@@ -72,8 +73,10 @@ You can use builtin path variables in this template and reference to other paths
# report gauges and sites
+ # sites (pixel average)
+ # gauges (catchment average)
@@ -82,6 +85,9 @@ You can use builtin path variables in this template and reference to other paths
+ # report mass balance
+
+
#-----------------------------------------------------------
# report maps
#-----------------------------------------------------------
@@ -91,33 +97,69 @@ You can use builtin path variables in this template and reference to other paths
# report end maps
- # report maps
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
+ # report output maps
+ # forcings
+
+
+
+
+
+
+ # dicharge
+
+
+
+
+
+ # variables
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+ #Surface runoff
+ #Surface runoff + UZ [mm]
+ #Surface runoff + UZ [mm] + LZ
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
-
@@ -188,57 +230,37 @@ You can use builtin path variables in this template and reference to other paths
-
-
- /perm/mo/mocm/proj/efas/hprot/staticMaps/areawithdeadsea.map
- /perm/mo/mocm/proj/efas/hprot/staticMaps/work/DrinaMask.map
- /sbsDrina/DrinaMask.map
- /perm/mo/mocm/proj/efas/xdom/data/staticData/lisflood/area
- /perm/mo/mocm/proj/efas/hprot/staticMaps/areawithdeadsea.map
- /perm/mo/mocm/proj/efas/hprot/sbs/pontelagoscuro_sbs.map
- 70 40 0.1 -8.6 43.3 # Douro
- 30 30 0.1 6.0 52.0
- Clone map
- col row cellsize xupleft yupleft
- 300 200 0.1 -15 25 -> Niger
- 3600 1500 0.1 -180 90 -> World
- $(PathRoot)\areamaps\areaEurope.map
- $(PathRoot)\areamaps\AreaPo.map
- or pcraster maps or netcdf maps
- e.g.
- Rhine areaLobith.map
- Niger areaLokoja.map
- Po areaPo.map
- Danube areaDanube.map
- Sava areaSava.map
-
+
+
+ Clone map
+ col row cellsize xupleft yupleft
+ or pcraster maps or netcdf maps
+
-
+
Nominal map with gauge locations (i.e cells for which simulated discharge
is written to file(1,2,3 etc)
lat lon (lat2 lon2 ...) or pcraster maps or netcdf maps
Examples:
"6.15 51.85"
- "$(PathMaps)/outlets.nc"
+ "$(PathMaps)/outlets"
-
+
netcdf template used to copy metadata information for writing netcdf
- $(PathEvapo)/$(PrefixE0)
-
-
- latitude map to be used for snow/ice modelling
-
+
+
+ latitude map to be used for snow/ice modelling
+
-
**************************************************************
TIME-RELATED CONSTANTS
@@ -264,24 +286,26 @@ You can use builtin path variables in this template and reference to other paths
-
- 17280 5 times subrouting
- 21600 4 times
- Sub time step used for kinematic wave channel routing [seconds]
- Within the model,the smallest out of DtSecChannel and DtSec is used
-
+
+ 17280 5 times subrouting
+ 21600 6 hours
+ Sub time step used for kinematic and MCT wave channel routing [seconds]
+ Within the model,the smallest out of DtSecChannel and DtSec is used
+
-
- Number of first time step in simulation
-
+
+ 01/01/YEAR_START 12:00
+ Number of first time step in simulation
+
-
- Number of last time step in simulation
-
+
+ 01/01/YEAR_END 06:00
+ Number of last time step in simulation
+
@@ -338,7 +362,7 @@ You can use builtin path variables in this template and reference to other paths
default: 10
- $(PathParams)/params_UpperZoneTimeConstant.nc
+ $(PathParams)/params_UpperZoneTimeConstant
Time constant for water in upper zone [days]
This is the average time a water 'particle' remains in the reservoir
if we had a stationary system (average inflow=average outflow)
@@ -348,7 +372,7 @@ You can use builtin path variables in this template and reference to other paths
default: 100
- $(PathParams)/params_LowerZoneTimeConstant.nc
+ $(PathParams)/params_LowerZoneTimeConstant
Time constant for water in lower zone [days]
This is the average time a water 'particle' remains in the reservoir
if we had a stationary system (average inflow=average outflow)
@@ -358,7 +382,7 @@ You can use builtin path variables in this template and reference to other paths
default: 0.5
- $(PathParams)/params_GwPercValue.nc
+ $(PathParams)/GwPercValue
Maximum rate of percolation going from the Upper to the Lower
response box [mm/day]
@@ -367,15 +391,17 @@ You can use builtin path variables in this template and reference to other paths
default: 0.0
- $(PathParams)/params_GwLoss.nc
+ $(PathParams)/GwLoss
Maximum percolation rate from the lower groundwater zone [mm/days]. GWLoss
it’s lost beyond the catchment boundaries or to deep groundwater systems.
A value of 0 (closed lower boundary) is recommended as a starting value.
-
+
+ default: 10
+ $(PathParams)/params_LZThreshold
threshold value [mm] of the water storage in the lower groundwater zone.
If the water storage in the lower groundwater zone decreases below LZThreshold,
the flow from the lower zone to the nearby rivers (base-flow) stops.
@@ -385,7 +411,7 @@ You can use builtin path variables in this template and reference to other paths
default: 0.7 [-]
- $(PathParams)/params_b_Xinanjiang.nc
+ $(PathParams)/b_Xinanjiang
Power in Xinanjiang distribution function [-]
@@ -393,7 +419,7 @@ You can use builtin path variables in this template and reference to other paths
default: 3.5 [-]
- $(PathParams)/params_PowerPrefFlow.nc
+ $(PathParams)/params_PowerPrefFlow
Power that controls increase of proportion of preferential
flow with increased soil moisture storage [-]
@@ -402,7 +428,7 @@ You can use builtin path variables in this template and reference to other paths
default: 2.0 [-]
- $(PathParams)/params_CalChanMan1.nc
+ $(PathParams)/params_CalChanMan1
Multiplier [-] applied to Channel Manning's n
@@ -410,7 +436,7 @@ You can use builtin path variables in this template and reference to other paths
default: 4.0 [mm C-1 day-1]
- $(PathParams)/params_SnowMeltCoef.nc
+ $(PathParams)/params_SnowMeltCoef
SRM: 0.45 cm C-1 day-1 ( = 4.50 mm C-1 day-1), Kwadijk: 18 mm C-1 month-1 (= 0.59 mm C-1 day-1)
See also Martinec et al., 1998.
@@ -419,15 +445,23 @@ You can use builtin path variables in this template and reference to other paths
default: 3.0 [-]
- $(PathParams)/params_CalChanMan2.nc
+ $(PathParams)/params_CalChanMan2
Multiplier [-] applied to Channel Manning's n for second line of routing
+
+
+ default: 3.0 [-]
+ $(PathParams)/params_CalChanMan3
+ Multiplier [-] applied to Channel Manning's n for MCT routing
+
+
+
default: 1.0 [-]
- $(PathParams)/params_LakeMultiplier.nc
+ $(PathParams)/params_LakeMultiplier
Multiplier applied to the lake parameter A
@@ -464,28 +498,33 @@ You can use builtin path variables in this template and reference to other paths
- PBchange
- Multiplier applied to average Q to split into a second line of routing
- THIS NEEDS TO BE AMENDED BY ADDING THE RIGHT PATH TO $(PathParams)
+ default: 2.0
+ $(PathParams)/params_QSplitMult
+ Multiplier applied to average Q to split into a second line of routing
- PBchange
- Multiplier applied to ChanBottomWidth
+ default: 1.0
+ $(PathParams)/params_ChanBottomWMult
+ Multiplier applied to ChanBottomWidth
- Multiplier [] applied to ChanDepthThreshold
+ default: 1.0
+ $(PathParams)/params_ChanDepthTMult
+ Multiplier [] applied to ChanDepthThreshold
- Multiplier [] applied to ChanSdXdY
+ default: 1.0
+$(PathParams)/params_ChanSMult
+Multiplier [] applied to ChanSdXdY
@@ -500,12 +539,14 @@ You can use builtin path variables in this template and reference to other paths
+ RUNDIR
Output path (org=$(PathRoot)/out)
+ INITDIR
Path of the initial value maps e.g. lzavin.map
@@ -561,12 +602,12 @@ You can use builtin path variables in this template and reference to other paths
Water use maps path
- $(PathMaps)/waterdemand19902014
+ FORCINGSDIR
Meteo path
@@ -582,12 +623,10 @@ You can use builtin path variables in this template and reference to other paths
variable water fraction maps path
-
-
**************************************************************
INITIAL CONDITIONS
@@ -605,61 +644,72 @@ You can use builtin path variables in this template and reference to other paths
+ Cold start: 0
Initial overland flow storage [m3], direct runoff fraction
- Initial overland flow storage [m3], other fraction
+ Cold start: 0
+ Initial overland flow storage [m3], other fraction
+ Cold start: 0
Initial overland flow storage [m3], forest fraction
- initial snow depth in snow zone A [mm]
+ Cold start: 0
+ initial snow depth in snow zone A [mm]
+ Cold start: 0
initial snow depth in snow zone B [mm]
- initial snow depth in snow zone C [mm]
+ Cold start: 0
+ initial snow depth in snow zone C [mm]
+ Cold start: 0
initial Frost Index value [degC/day]
- cumulative interception [mm]
+ Cold start: 0
+ cumulative interception [mm]
- water in upper store [mm]
+ Cold start: 0
+ PRERUNDIR/uz.end.nc
+ water in upper store [mm]
- days since last rainfall [days]
+ Cold start: 1
+ days since last rainfall [days]
@@ -676,82 +726,99 @@ You can use builtin path variables in this template and reference to other paths
- water in lower store [mm]
- -9999: use steady-state storage
+ PRERUNDIR/lz.end.nc
+ water in lower store [mm]
+ Cold start: -9999: use steady-state storage
- initial cross-sectional area of flow in channel[m2]
- -9999: use half bankfull
+ initial cross-sectional area of flow in channel[m2]
+ Cold start: -9999: use half bankfull
- initial soil moisture content layer 1a (superficial soil layer) [mm3/mm3]
- -9999: use field capacity values
+ PRERUNDIR/tha.end.nc
+ initial soil moisture content layer 1a (superficial soil layer) [mm3/mm3]
+ Cold start: -9999: use field capacity values
+ PRERUNDIR/thb.end.nc
initial soil moisture content layer 1b (upper soil layer) [mm3/mm3]
- -9999: use field capacity values
+ Cold start: -9999: use field capacity values
- initial soil moisture content layer 2 (lower soil layer) [mm3/mm3]
- -9999: use field capacity values
+ PRERUNDIR/thc.end.nc
+ initial soil moisture content layer 2 (lower soil layer) [mm3/mm3]
+ Cold start: -9999: use field capacity values
- initial channel cross-section area [m2] for 2nd routing channel
- -9999: use 0
+ initial channel cross-section area [m2] for 2nd routing channel
+ Cold start: -9999: use 0
- initial lateral inflow for 1st line of routing [m2/s]
- -9999: use 0
+ initial lateral inflow for 1st line of routing [m2/s]
+ Cold start: -9999: use 0
- Initial lake level [m]
- -9999 sets initial value to steady-state level
+ Initial lake level [m]
+ Cold start: -9999 sets initial value to steady-state level
- Lake inflow at previous routing sub-step (ChanQ(t-1)) [m3/s]
- -9999 sets initial value equal to PrevDischarge (ChanQ(t))
+ Lake inflow at previous routing sub-step (ChanQ(t-1)) [m3/s]
+ Cold start: -9999 sets initial value equal to PrevDischarge (ChanQ(t))
- Lake outflow at previous routing sub-step (ChanQ(t-1)) [m3/s]
- -9999 sets initial value
+ Lake outflow at previous routing sub-step (ChanQ(t-1)) [m3/s]
+ Cold start: -9999 sets initial value
- initial discharge [m3/s] from previous run for lakes, reservoirs and transmission loss ChanQ
- only needed for lakes reservoirs and transmission loss
- -9999: use 0
+ initial discharge [m3/s] from previous run for lakes, reservoirs and transmission loss ChanQ
+ only needed for lakes reservoirs and transmission loss
+ -9999: use 0
+
+
+
+
+
+ Courant number at previous step for MCT routing
+ Cold start: -9999 sets initial value to 1
+
+
+ Reynolds number at previous step for MCT routing
+ Cold start: -9999 sets initial value to 0
+
+
@@ -764,48 +831,54 @@ You can use builtin path variables in this template and reference to other paths
-
- cumulative interception [mm]
-
+
+ cumulative interception [mm]
+ Cold start: 0
+
-
- water in upper groundwater store [mm]
-
+
+ PRERUNDIR/uzf.end.nc
+ water in upper groundwater store [mm]
+ Cold start: 0
+
-
- days since last rainfall
-
+
+ days since last rainfall
+ Cold start: 1
+
-
- initial soil moisture content [mm3/mm3] layer 1a (superficial soil layer)
- -9999: use field capacity values
-
+
+ PRERUNDIR/thfa.end.nc
+ initial soil moisture content [mm3/mm3] layer 1a (superficial soil layer)
+ Cold start: -9999: use field capacity values
+
-
- initial soil moisture content [mm3/mm3] layer 1b (upper soil layer)
- -9999: use field capacity values
-
+
+ PRERUNDIR/thfb.end.nc
+ initial soil moisture content [mm3/mm3] layer 1b (upper soil layer)
+ Cold start: -9999: use field capacity values
+
-
- initial soil moisture content [mm3/mm3] layer 2 (lower soil layer)
- -9999: use field capacity values
-
+
+ PRERUNDIR/thfc.end.nc
+ initial soil moisture content [mm3/mm3] layer 2 (lower soil layer)
+ Cold start: -9999: use field capacity values
+
-
**************************************************************
INITIAL CONDITIONS IRRIGATION
@@ -814,44 +887,50 @@ You can use builtin path variables in this template and reference to other paths
-
- cumulative interception [mm]
-
+
+ cumulative interception [mm]
+ Cold start: 0
+
-
- water in groundwater upper store [mm]
-
+
+ PRERUNDIR/uzi.end.nc
+ water in groundwater upper store [mm]
+ Cold start: 0
+
-
- days since last rainfall
-
+
+ days since last rainfall
+ Cold start: 1
+
-
- initial soil moisture content [mm3/mm3] layer 1a (superficial soil layer)
- -9999: use field capacity values
-
+
+ PRERUNDIR/thia.end.nc
+ initial soil moisture content [mm3/mm3] layer 1a (superficial soil layer)
+ Cold start: -9999: use field capacity values
+
-
- initial soil moisture content [mm3/mm3] layer 1b (upper soil layer)
- -9999: use field capacity values
-
+
+ PRERUNDIR/thib.end.nc
+ initial soil moisture content [mm3/mm3] layer 1b (upper soil layer)
+ Cold start: -9999: use field capacity values
+
-
- initial soil moisture content [mm3/mm3] layer 2 (lower soil layer)
- -9999: use field capacity values
-
+
+ PRERUNDIR/thic.end.nc
+ initial soil moisture content [mm3/mm3] layer 2 (lower soil layer)
+ Cold start: -9999: use field capacity values
+
-
@@ -864,9 +943,10 @@ You can use builtin path variables in this template and reference to other paths
-
- cumulative interception [mm]
-
+
+ cumulative interception [mm]
+ Cold start: 0
+
@@ -952,14 +1032,10 @@ You can use builtin path variables in this template and reference to other paths
-
- prefix variable water fraction
-
+
+ prefix variable water fraction
+
-
-
-
-
@@ -1094,9 +1170,9 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
+
- 0.75
+ default: 0.75
Field application irrigation efficiency [-]: max 1, ~0.90 drip irrigation, ~0.75 sprinkling
@@ -1144,7 +1220,7 @@ You can use builtin path variables in this template and reference to other paths
-
+
# Consumptive Use (1-Recycling ratio) for energy water use (0-1) [-]
# Source: Torcellini et al. (2003) "Consumptive Use for US Power Production"
@@ -1162,11 +1238,11 @@ You can use builtin path variables in this template and reference to other paths
-
+
- 0.2
- $(PathMaps)/leakage.map
- Fraction of leakage of public water supply (0=no leakage, 1=100% leakage)
+ Default : 0.2
+ $(PathRoot)/leakage
+ Fraction of leakage of public water supply (0=no leakage, 1=100% leakage)
@@ -1196,13 +1272,13 @@ You can use builtin path variables in this template and reference to other paths
-
+
Irrigation crop coefficient [-]
-
+
Irrigation crop group number [-]
@@ -1220,7 +1296,7 @@ You can use builtin path variables in this template and reference to other paths
-
+
Land Use Mask
@@ -1259,7 +1335,7 @@ You can use builtin path variables in this template and reference to other paths
-
+
@@ -1283,58 +1359,42 @@ You can use builtin path variables in this template and reference to other paths
- calibrate on GwLossFraction if necessary
-
-
length of the window used to smooth the LZ zone [number of cell length]
-
+
- $(PathMaps)/gwbodiesNew.map
- map of aquifers (0/1), used to smoothen LZ near extraction areas
+ map of aquifers (0/1), used to smoothen LZ near extraction areas
-
-
**************************************************************
EVAPORATION FROM OPEN WATER
**************************************************************
-
-
+
Fraction of maximum extend of water
-
+
Mask with Lakes from GLWD database
-
-
-
-
**************************************************************
TRANSMISSION LOSS PARAMETERS
@@ -1343,12 +1403,12 @@ You can use builtin path variables in this template and reference to other paths
Suggested parameterisation strategy:
- Use TransSub as calibration constant leave all other parameters at default values
-
+
-
- Transmission loss function parameter
-
+
+ Transmission loss function parameter
+
@@ -1365,17 +1425,15 @@ You can use builtin path variables in this template and reference to other paths
-
+
upstream area
-
-
**************************************************************
ROUTING PARAMETERS
@@ -1385,7 +1443,6 @@ You can use builtin path variables in this template and reference to other paths
- Use CalChanMan as calibration constant to fine-tune timing of
channel routing, leave all other parameters at default values
-
@@ -1401,27 +1458,25 @@ You can use builtin path variables in this template and reference to other paths
-
-
- Minimum slope gradient [-] (for kin. wave: slope cannot be 0)
-
+
+
+ Percentage [-] used to compute the TotalCrossSectionAreaHalfBankFull in the initialisation of the
+ total cross-sectional area
+ default: 0.5
+
-
-
- Minimum channel gradient [-] (for kin. wave: slope cannot be 0)
-
+
+
+ Minimum slope gradient [-] (for kin. wave: slope cannot be 0)
+
-
-
@@ -1437,14 +1492,13 @@ You can use builtin path variables in this template and reference to other paths
-
- Minimum value for Courant condition in soil moisture routine. Always less than or equal to 1.
- Small values result in improved numerical accuracy, at the expense of increased
- computing time (more sub-steps needed). If reported time series of soil moisture contain
- large jumps, lowering CourantCrit should fix this
-
+
+ Minimum value for Courant condition in soil moisture routine. Always less than or equal to 1.
+ Small values result in improved numerical accuracy, at the expense of increased
+ computing time (more sub-steps needed). If reported time series of soil moisture contain
+ large jumps, lowering CourantCrit should fix this
+
-
@@ -1473,13 +1527,12 @@ You can use builtin path variables in this template and reference to other paths
-
- Initial reservoir fill fraction [-]
- -9999 sets initial fill to 80% of the flood limit
- if you're not using the reservoir option, enter some bogus value
-
+
+ Initial reservoir fill fraction [-]
+ Cold start: -9999 sets initial fill to 80% of the flood limit
+ if you're not using the reservoir option, enter some bogus value
+
-
@@ -1491,18 +1544,17 @@ You can use builtin path variables in this template and reference to other paths
-
- Estimate of average net inflow into lake (=inflow - evaporation) [cu m / s]
- Used to calculated steady-state lake level in case LakeInitialLevelValue
- is set to -9999
-
+
+ Cold start: -9999
+ Estimate of average net inflow into lake (=inflow - evaporation) [cu m / s]
+ Used to calculated steady-state lake level in case LakeInitialLevelValue
+ is set to -9999
+
-
-
**************************************************************
POLDER OPTION
@@ -1516,11 +1568,10 @@ You can use builtin path variables in this template and reference to other paths
-
- Initial water level in polder [m]
-
+
+ Initial water level in polder [m]
+
-
@@ -1541,28 +1592,49 @@ You can use builtin path variables in this template and reference to other paths
+
+ Constant head [m] at most downstream pixel (relative to altitude
+ at most downstream pixel)
+
+
+
+
+
+
- Constant head [m] at most downstream pixel (relative to altitude
- at most downstream pixel)
+ **************************************************************
+ MUSKINGUM-CUNGE-TODINI OPTION
+ **************************************************************
+
+
+
+ Boolean map with value 1 at channel pixels where MCT is
+ used, and 0 at all other pixels
+
+
+
+ Maximum channel gradient for channels using MCT routing [-] (for MCT wave: slope cannot be 0)
+ Default: 0.001
+
+
-
-
+
**************************************************************
INFLOW HYDROGRAPH (OPTIONAL)
**************************************************************
-
-
- OPTIONAL: nominal map with locations of (measured)
- inflow hydrographs [cu m / s]
-
+
+
+ OPTIONAL: nominal map with locations of (measured)
+ inflow hydrographs [cu m / s]
+
@@ -1573,7 +1645,6 @@ You can use builtin path variables in this template and reference to other paths
InflowPoints map (also optional)
-
@@ -1650,13 +1721,13 @@ You can use builtin path variables in this template and reference to other paths
-
-
-
+
+
+
-
+
-
This is necessary when using projected coordinates (x,y)
@@ -1697,6 +1767,7 @@ You can use builtin path variables in this template and reference to other paths
latitude map to be used for snow/ice modelling
+
@@ -1745,7 +1816,6 @@ You can use builtin path variables in this template and reference to other paths
Number of last time step in simulation
-
@@ -1801,12 +1871,10 @@ You can use builtin path variables in this template and reference to other paths
which is not immediatly causing surface runoff [mm]
-
-
**************************************************************
SNOW AND FROST RELATED PARAMETERS
@@ -1882,12 +1950,10 @@ You can use builtin path variables in this template and reference to other paths
Molnau and Bissel found a value 56-85 for NW USA.
-
-
**************************************************************
INFILTRATION PARAMETERS
@@ -1911,12 +1977,10 @@ You can use builtin path variables in this template and reference to other paths
flow with increased soil moisture storage [-]
-
-
**************************************************************
GROUNDWATER RELATED PARAMETERS
@@ -1958,12 +2022,9 @@ You can use builtin path variables in this template and reference to other paths
A value of 0 (closed lower boundary) is recommended as a starting value.
-
-
-
**************************************************************
@@ -1971,7 +2032,6 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
water use daily maps with a (in this case negative) volume of water [cu m/s]
@@ -2015,15 +2075,10 @@ You can use builtin path variables in this template and reference to other paths
Time series of upstream water evaporation from open water at gauging stations
-
-
-
-
-
**************************************************************
TRANSMISSION LOSS PARAMETERS
@@ -2061,7 +2116,6 @@ You can use builtin path variables in this template and reference to other paths
upstream area
-
@@ -2090,6 +2144,12 @@ You can use builtin path variables in this template and reference to other paths
+
+
+ Multiplier applied to Channel Manning's n for MCT routing
+
+
+
Multiplier applied to average Q to split into a second line of routing
@@ -2127,6 +2187,14 @@ You can use builtin path variables in this template and reference to other paths
+
+
+ Percentage [-] used to compute the TotalCrossSectionAreaHalfBankFull in the initialisation of the
+ total cross-sectional area
+ default: 0.5
+
+
+
Minimum slope gradient (for kin. wave: slope cannot be 0)
@@ -2160,7 +2228,6 @@ You can use builtin path variables in this template and reference to other paths
large jumps, lowering CourantCrit should fix this
-
@@ -2174,8 +2241,7 @@ You can use builtin path variables in this template and reference to other paths
- $(PathMaps)/lzavin.map
- $(PathInit)/lzavin.map
+ PRERUNDIR/lzavin
Reported map of average percolation rate from upper to
lower groundwater zone (reported for end of simulation) [mm/day]
@@ -2183,8 +2249,7 @@ You can use builtin path variables in this template and reference to other paths
- $(PathMaps)/avgdis.map
- $(PathInit)/avgdis.map
+ PRERUNDIR/avgdis
CHANNEL split routing in two lines
Average discharge map [m3/s]
@@ -2257,9 +2322,9 @@ You can use builtin path variables in this template and reference to other paths
days since last rainfall
-
+
**************************************************************
@@ -2320,17 +2385,28 @@ You can use builtin path variables in this template and reference to other paths
- initial discharge [m3/s] from previous run for lakes, reservoirs and transmission loss
- only needed for lakes reservoirs and transmission loss
+ Outflow discharge at the end of previous step (instantaneous) [m3/s]
-9999: use 0
+
+
+ Courant number at previous step for MCT routing
+ Cold start: -9999 sets initial value to 1
+
+
+
+
+
+ Reynolds number at previous step for MCT routing
+ Cold start: -9999 sets initial value to 0
+
+
-
**************************************************************
INITIAL CONDITIONS FOREST
@@ -2386,7 +2462,6 @@ You can use builtin path variables in this template and reference to other paths
-
**************************************************************
INITIAL CONDITIONS IRRIGATION
@@ -2432,12 +2507,9 @@ You can use builtin path variables in this template and reference to other paths
-9999: use field capacity values
-
-
-
**************************************************************
@@ -2492,12 +2564,11 @@ You can use builtin path variables in this template and reference to other paths
leaf area index [m2/m2]
-
+
-
**************************************************************
INPUT WATER USE MAPS AND PARAMETERS
**************************************************************
@@ -2589,7 +2660,6 @@ You can use builtin path variables in this template and reference to other paths
-
Field application irrigation efficiency max 1, ~0.90 drip irrigation, ~0.75 sprinkling
@@ -2694,7 +2764,6 @@ You can use builtin path variables in this template and reference to other paths
-
@@ -2705,7 +2774,6 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
water amount in mm per day
@@ -2719,30 +2787,29 @@ You can use builtin path variables in this template and reference to other paths
-
+
map with starting day of the year
-
+
map with starting day of the year
-
+
map with starting day of the year
-
+
map with starting day of the year
-
@@ -2788,6 +2855,8 @@ You can use builtin path variables in this template and reference to other paths
Reported water level [m] % False by default
+
+
@@ -2895,7 +2964,6 @@ You can use builtin path variables in this template and reference to other paths
-
Reported cross section area 2nd line of routing [m2]
@@ -2914,6 +2982,18 @@ You can use builtin path variables in this template and reference to other paths
+
+
+ Courant number at previous step for MCT routing
+
+
+
+
+
+ Reynolds number at previous step for MCT routing
+
+
+
Reported days since last rain
@@ -3019,7 +3099,6 @@ You can use builtin path variables in this template and reference to other paths
Reported transpiration maps, weighted sum over the fractions of each pixel [mm]
-
@@ -3150,13 +3229,12 @@ You can use builtin path variables in this template and reference to other paths
-
+
Reported lake storage [m3]
-
Reported reservoir filling [-]
@@ -3183,7 +3261,19 @@ You can use builtin path variables in this template and reference to other paths
- Reported discharge [cu m/s] (average over the model timesteps) (AvgDis)
+ Reported average discharge [cu m/s] (average over the model timesteps) (AvgDis)
+
+
+
+
+
+ Courant number at previous step for MCT routing
+
+
+
+
+
+ Reynolds number at previous step for MCT routing
@@ -3192,23 +3282,27 @@ You can use builtin path variables in this template and reference to other paths
Reported days since last rain
+
Reported interception storage [mm]
+
Reported volumetric soil moisture content for
soil layer 1a (superficial layer) [V/V] [mm3/mm3]
+
Reported volumetric soil moisture content for both
soil layer 1b (upper layer) [V/V] [mm3/mm3]
+
Reported volumetric soil moisture content for both
@@ -3228,23 +3322,25 @@ You can use builtin path variables in this template and reference to other paths
-
Reported days since last rain
+
Reported interception storage [mm]
+
Reported volumetric soil moisture content for
soil layer 1a (superficial layer) [V/V] [mm3/mm3]
+
Reported volumetric soil moisture content for both
@@ -3265,7 +3361,6 @@ You can use builtin path variables in this template and reference to other paths
-
@@ -3306,9 +3401,9 @@ You can use builtin path variables in this template and reference to other paths
Potential evaporation from open water surface [mm per time step]
-
+
INDIVIDUAL STATE VARIABLES (AT EVERY TIME STEP)
@@ -3690,7 +3785,6 @@ You can use builtin path variables in this template and reference to other paths
-
Reported evaporation of intercepted water [mm]
@@ -3801,7 +3895,6 @@ You can use builtin path variables in this template and reference to other paths
Reported storage in upper response box [mm]
-
@@ -3851,10 +3944,6 @@ You can use builtin path variables in this template and reference to other paths
-
-
-
-
**************************************************************
@@ -3864,13 +3953,13 @@ You can use builtin path variables in this template and reference to other paths
- Reported discharge [cu m/s]
+ Reported average discharge over the model time step [cu m/s]
- Reported discharge over last routing sub-step [cu m/s]
+ Reported instantaneous discharge at the end of routing sub-step [cu m/s]
@@ -3879,6 +3968,8 @@ You can use builtin path variables in this template and reference to other paths
Reported water level [m]
+
+
@@ -3968,10 +4059,9 @@ You can use builtin path variables in this template and reference to other paths
Days since last rain [days]
-
-
+
OUTPUT AT SITES: RATE VARIABLES
@@ -4084,9 +4174,9 @@ You can use builtin path variables in this template and reference to other paths
Reported loss from lower zone [mm]
-
+
AVERAGE VALUES OF METEOROLOGICAL FORCING VARIABLES
@@ -4124,8 +4214,14 @@ You can use builtin path variables in this template and reference to other paths
+
+
+ ESActPixel+self.var.TaPixel+self.var.TaInterceptionAll+self.var.EvaAddM3*self.var.M3toMM
+
+
+
AVERAGE VALUES OF MODEL STATE VARIABLES
@@ -4180,10 +4276,9 @@ You can use builtin path variables in this template and reference to other paths
Days since last rain [days]
-
-
+
AVERAGE VALUES OF MODEL RATE VARIABLES
@@ -4293,14 +4388,14 @@ You can use builtin path variables in this template and reference to other paths
-
+
Reported loss from lower zone [mm]
-
+
NUMERICS
@@ -4368,10 +4463,8 @@ You can use builtin path variables in this template and reference to other paths
routine
-
-
@@ -4387,20 +4480,20 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
+
$(PathMaps)/ldd.map
local drain direction map (1-9)
-
+
slope gradient [-] (0.50 -> tangent=0.5, angle=26.5 degrees)
-
+
Elevation standard deviation [m], i.e. altitude difference elevation within pixel.
Used for sub-pixel modelling of snow accumulation
@@ -4416,31 +4509,28 @@ You can use builtin path variables in this template and reference to other paths
-
+
map with monitoring sites (1,2,3 etc)
-
+
optional map with pixel length [m], only needed in case map attributes
are not in [m] (e.g. lat / lon systems)
-
+
optional map with pixel area [m2], only needed in case map attributes
are not in [m] (e.g. lat / lon systems)
-
-
-
**************************************************************
SOIL PARAMETERS
@@ -4452,7 +4542,7 @@ You can use builtin path variables in this template and reference to other paths
-
+
Soil depth for 1st layer: 1a (superficial layer)
Minimum Soil Depth [mm] for superficial soil layer
@@ -4461,21 +4551,21 @@ You can use builtin path variables in this template and reference to other paths
-
+
Soil depth for 2nd layer: 1b (upper layer)
Minimum Soil Depth [mm] for second soil layer
-
+
Soil depth for 2nd layer: 2 (lower layer)
Minimum Soil Depth [mm] for second soil layer
-
+
Soil depth for 1st layer: 1a
Minimum Soil Depth [mm] for upper soil layer
@@ -4484,20 +4574,19 @@ You can use builtin path variables in this template and reference to other paths
-
+
Soil depth for 2nd layer: 1b
Minimum Soil Depth [mm] for lower soil layer
-
+
Soil depth for 2nd layer: 2
Minimum Soil Depth [mm] for lower soil layer 300
-
@@ -4508,9 +4597,9 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
+
- $(PathMapsLanduse)/fracsealed.map
+ $(PathMapsLanduse)/fracsealed
urban area (0-1)
@@ -4521,9 +4610,9 @@ You can use builtin path variables in this template and reference to other paths
-
+
- $(PathMapsLanduse)/fracforest.map
+ $(PathMapsLanduse)/fracforest
forest fraction of a pixel (0-1)
@@ -4535,9 +4624,9 @@ You can use builtin path variables in this template and reference to other paths
-
+
- $(PathMapsLanduse)/fracwater.map
+ $(PathMapsLanduse)/fracwater
forest fraction of a pixel (0-1)
@@ -4549,9 +4638,9 @@ You can use builtin path variables in this template and reference to other paths
-
+
- $(PathMapsLanduse)/fracother.map
+ $(PathMapsLanduse)/fracother
forest fraction of a pixel (0-1)
@@ -4563,9 +4652,9 @@ You can use builtin path variables in this template and reference to other paths
-
+
- $(PathMapsLanduse)/fracirrigated.map
+ $(PathMapsLanduse)/fracirrigated
irrigated area fraction of a pixel (0-1)
@@ -4584,7 +4673,7 @@ You can use builtin path variables in this template and reference to other paths
-
+
rice fraction of a pixel (0-1)
@@ -4596,13 +4685,13 @@ You can use builtin path variables in this template and reference to other paths
-
+
groundwater used fraction of a pixel (0-1)
-
+
NonConventionalWaterUsed (0-1)
@@ -4614,58 +4703,56 @@ You can use builtin path variables in this template and reference to other paths
-
+
$(PathMaps)/eflow.map
EFlowThreshold is map with m3/s discharge, e.g. the 10th percentile discharge of the baseline run
-
+
table with crop coefficient for each land use; values ranging from 0.75 to 1.25;
source: Supit, p.83 CFET
-
+
table with crop group number; table 6.1 and 6.2 in WOFOST 6.0, Supit, 1994: p. 86
-
+
- $(PathMapsTables)/manngs_2000_250m.map
+ $(PathMapsTables)/manngs_2000_250m
table with Manning's roughness [m^(1/3) s^(-1)] for each CORINE landcover class
-
+
table with crop coefficient for each land use; values ranging from 0.75 to 1.25;
source: Supit, p.83 CFET
-
+
table with crop group number; table 6.1 and 6.2 in WOFOST 6.0, Supit, 1994: p. 86
-
+
- $(PathMapsTables)/manngs_2000_250m.map
+ $(PathMapsTables)/manngs_2000_250m
table with Manning's roughness [m^(1/3) s^(-1)] for each CORINE landcover class
-
-
**************************************************************
Variable Channel NoSubStepChannel
@@ -4690,7 +4777,7 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
+
Boolean map with value 1 at channel pixels and 0 at
all other pixels
@@ -4699,51 +4786,51 @@ You can use builtin path variables in this template and reference to other paths
-
+
Channel gradient (fraction, dy/dx)
-
+
Channel Manning's n [m^(1/3) s^(-1)]
-
+
Channel length [m]
-
+
Channel bottom width [m]
-
+
Channel side slope, expressed as gradient
!! expressed as dx/dy !!! (NOT dy/dx, which would perhaps be more logical)
-
+
Floodplain Width [m]
-
+
Bankfull channel depth [m]
-
+
**************************************************************
@@ -4753,94 +4840,94 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
+
Soil moisture content at saturation [V/V] [mm3/mm3] for soil layer 1a (superficial soil layer)
-
+
Soil moisture content at saturation [V/V] [mm3/mm3] for soil layer 1b (upper soil layer)
-
+
Soil moisture content at saturation [V/V] [mm3/mm3] for soil layer 2 (lower soil layer)
-
+
Residual soil moisture content [V/V] [mm3/mm3] for soil layer 1a (superficial soil layer)
-
+
Residual soil moisture content [V/V] [mm3/mm3] for soil layer 1b (upper soil layer)
-
+
Residual soil moisture content [V/V] [mm3/mm3] for soil layer 2 (lower soil layer)
-
+
Pore-size index (Van Genuchten m and n are calculated from this)
Lambda is acually 1-n [-]
-
+
Pore-size index (Van Genuchten m and n are calculated from this)
Lambda is acually 1-n [-]
-
+
Pore-size index (Van Genuchten m and n are calculated from this)
Lambda is acually 1-n [-]
-
+
Van Genuchten parameter Alpha [1/cm] for soil layer 1a (superficial soil layer)
-
+
Van Genuchten parameter Alpha [1/cm] for soil layer 1b (upper soil layer)
-
+
Van Genuchten parameter Alpha [1/cm] for soil layer 2 (lower soil layer)
-
+
Saturated conductivity [mm/day] for soil layer 1a (superficial soil layer)
-
+
Saturated conductivity [mm/day] for soil layer 1b (upper soil layer)
-
+
Saturated conductivity [mm/day] for soil layer 2 (lower soil layer)
@@ -4848,6 +4935,7 @@ You can use builtin path variables in this template and reference to other paths
+
**************************************************************
@@ -4857,63 +4945,63 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
+
Soil moisture content at saturation [V/V] [mm3/mm3] for soil layer 1a (superficial soil layer)
-
+
Soil moisture content at saturation [V/V] [mm3/mm3]
-
+
Residual soil moisture content [V/V] [mm3/mm3] for soil layer 1a (superficial soil layer)
-
+
Residual soil moisture content [V/V] [mm3/mm3] for soil layer 1b (upper soil layer)
-
+
Pore-size index (Van Genuchten m and n are calculated from this) for soil layer 1a (superficial soil layer)
Lambda is acually 1-n [-] for soil layer 1b (upper soil layer)
-
+
Pore-size index (Van Genuchten m and n are calculated from this)
Lambda is acually 1-n [-] for soil layer 1b (upper soil layer)
-
+
Van Genuchten parameter Alpha [1/cm] for soil layer 1a (superficial soil layer)
-
+
Van Genuchten parameter Alpha [1/cm] for soil layer 1b (upper soil layer)
-
+
Saturated conductivity [mm/day] for soil layer 1a (superficial soil layer)
-
+
Saturated conductivity [mm/day] for soil layer 1b (upper soil layer)
@@ -4937,7 +5025,7 @@ You can use builtin path variables in this template and reference to other paths
-
+
Map with location of reservoirs (number)
@@ -5030,7 +5118,6 @@ You can use builtin path variables in this template and reference to other paths
reservoir routine switches to flood control mode.
-
@@ -5042,7 +5129,7 @@ You can use builtin path variables in this template and reference to other paths
**************************************************************
-
+
Map with location of lakes
@@ -5140,12 +5227,10 @@ You can use builtin path variables in this template and reference to other paths
-9999 sets initial value
-
-
**************************************************************
POLDERS(if used)
@@ -5158,7 +5243,7 @@ You can use builtin path variables in this template and reference to other paths
-
+
Map with polder locations
@@ -5223,7 +5308,6 @@ You can use builtin path variables in this template and reference to other paths
name of output map(s) with polder level
-
@@ -5249,7 +5333,6 @@ You can use builtin path variables in this template and reference to other paths
InflowPoints map (also optional)
-
@@ -5274,20 +5357,20 @@ You can use builtin path variables in this template and reference to other paths
-
+
Boolean map with value 1 at channel pixels where dynamic wave is
used, and 0 at all other pixels
-
+
Channel bottom level [m]
-
+
Nominal map with channel cross section IDs
@@ -5298,9 +5381,31 @@ You can use builtin path variables in this template and reference to other paths
Table with cross section parameters (H,A,P)
+
+
+
+
+ **************************************************************
+ MUSKINGUM-CUNGE-TODINI OPTION
+ **************************************************************
+
+
+
+
+ Boolean map with value 1 at channel pixels where MCT is
+ used, and 0 at all other pixels
+
+
+
+
+
+ Maximum channel gradient for channels using MCT routing [-] (for MCT wave: slope cannot be >0.001)
+
+
+
**************************************************************
@@ -5392,24 +5497,24 @@ You can use builtin path variables in this template and reference to other paths
Reported pF upper soil layer [-]
+
Reported pF lower soil layer [-]
+
Reported pF upper soil layer [-]
+
Reported pF lower soil layer [-]
-
-
-
@@ -5437,6 +5542,7 @@ You can use builtin path variables in this template and reference to other paths
map of aquifers (0/1), used to smoothen LZ near extraction areas
+
@@ -5452,7 +5558,6 @@ You can use builtin path variables in this template and reference to other paths
-
TotalAbstractionFromGroundwater [mm]
@@ -5616,7 +5721,6 @@ You can use builtin path variables in this template and reference to other paths
-
Abstraction from groundwater in m3 per timestep
@@ -5743,7 +5847,6 @@ You can use builtin path variables in this template and reference to other paths
-