From ddb4ff9e7a7cab7ecb35e24e98fb2f9d9c91ef63 Mon Sep 17 00:00:00 2001 From: Helen Baron Date: Wed, 26 Aug 2026 14:22:42 +0100 Subject: [PATCH 1/4] Add reservoir variables --- doc/source/namelists/ancillaries.nml.rst | 48 +++- .../namelists/initial_conditions.nml.rst | 24 +- doc/source/namelists/jules_rivers.nml.rst | 13 + doc/source/output-variables.rst | 38 ++- .../jules-standalone/HEAD/rose-meta.conf | 22 +- rose-meta/jules-standalone/versions.py | 9 +- src/control/shared/jules_rivers_mod.F90 | 93 ++++++- .../allocate_river_arrays_mod.F90 | 10 +- .../shared/allocate_jules_arrays.F90 | 30 ++- .../init_rivers_process_data_mod.F90 | 81 +++++- .../ancillaries/init_rivers_props_mod.F90 | 26 +- .../ancillaries/jules_rivers_props_mod.F90 | 255 +++++++++++++++--- .../standalone/init_output_mod.F90 | 74 +++-- .../initial_conditions/init_ic_mod.F90 | 4 + src/io/dump/get_dim_info_mod.F90 | 3 +- src/io/dump/read_dump_mod.F90 | 3 +- src/io/dump/read_dump_var_rivers_mod.F90 | 3 + src/io/dump/required_vars_for_rivers_mod.F90 | 11 +- src/io/dump/write_dump_mod.F90 | 6 +- src/io/dump/write_dump_var_rivers_mod.F90 | 3 + src/io/model_interface/extract_var.inc | 42 +++ .../model_interface/model_interface_mod.F90 | 4 +- src/io/model_interface/populate_var.inc | 72 ++++- src/io/model_interface/variable_metadata.inc | 254 ++++++++++++++++- 24 files changed, 1010 insertions(+), 118 deletions(-) diff --git a/doc/source/namelists/ancillaries.nml.rst b/doc/source/namelists/ancillaries.nml.rst index 2b9bd9b2..bc4beb28 100644 --- a/doc/source/namelists/ancillaries.nml.rst +++ b/doc/source/namelists/ancillaries.nml.rst @@ -1093,7 +1093,7 @@ This namelist specifies the options available for initialising irrigated fractio .. nml:namelist:: JULES_RIVERS_PROPS -This namelist specifies how spatially varying river routing properties including overbank inundation should be set. +This namelist specifies how spatially varying river routing properties (including overbank inundation and reservoirs) should be set. .. note:: ``read_from_dump`` is not currently implemented for this namelist, meaning that river ancillary variables cannot be read from a dump file. Initial values of river prognostic variables can however be read from a dump file (see :nml:lst:`JULES_INITIAL`). @@ -1455,14 +1455,14 @@ The following table summarises river routing properties required to run RFM or T | ``longitude_2d`` | The longitude of each river grid point must be specified. This field is required only if the model | | | coordinates are latitude and longitude, i.e. if :nml:mem:`JULES_LATLON::l_coord_latlon` = FALSE. | +----------------------------+-----------------------------------------------------------------------------------------------------------+ -| ``logn_mean`` | Mean of ln(elevation-elev_min) for each grid cell (in units ln(m)) | +| ``logn_mean`` | Mean of ln(elevation-elev_min) for each grid cell (in units ln(m)). | | | | | | This is only used if :nml:mem:`JULES_OVERBANK::overbank_model` = 3. | | | | | | Note that elev_min is DEM minimum, not river/lake bed level (therefore large values close to water | | | bodies can occur in floodplain gridcells). | +----------------------------+-----------------------------------------------------------------------------------------------------------+ -| ``logn_stdev`` | Standard deviation of ln(elevation-elev_min) for each grid cell (in units ln(m)) | +| ``logn_stdev`` | Standard deviation of ln(elevation-elev_min) for each grid cell (in units ln(m)). | | | | | | This is only used if :nml:mem:`JULES_OVERBANK::overbank_model` = 3. | | | | @@ -1479,7 +1479,7 @@ The following table summarises river routing properties required to run RFM or T | | via OASIS the river outflow is distributed over the corresponding river outflow points on the ocean grid. | | | This is to ensure that water is conserved and rivers discharge into the correct ocean grid points. | +----------------------------+-----------------------------------------------------------------------------------------------------------+ -| ``rivers_storage`` | Water storage (kg) on the Rivers grid. | +| ``rivers_storage`` | Water storage (kg) on the rivers grid. | | | | | | This is only used if :nml:mem:`JULES_RIVERS::i_river_vn` = 3. | | | | @@ -1487,6 +1487,46 @@ The following table summarises river routing properties required to run RFM or T | | climatology, rather than the dump. If Rivers storage is requested from an ancillary then | | | ``rivers_sto_rp`` should not be requested via :nml:lst:`JULES_INITIAL`. | +----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_capacity_grid`` | Storage capacity of reservoirs (kg). | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_catch_grid`` | Catchment area of reservoirs (m2). | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_year_grid`` | Initialisation year of reservoirs. | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_critical_grid`` | Critical storage of reservoirs (kg). | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_flood_grid`` | Flood storage of reservoirs (kg). | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_emergency_grid`` | Emergency storage of reservoirs (kg). | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_normal_release_grid``| Normal release rate of reservoirs (kg s-1). | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ +| ``res_flood_release_grid`` | Flood release rate of reservoirs (kg s-1). | +| | | +| | This is only used if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | +| | | ++----------------------------+-----------------------------------------------------------------------------------------------------------+ Example of how to set up the river grid diff --git a/doc/source/namelists/initial_conditions.nml.rst b/doc/source/namelists/initial_conditions.nml.rst index 1fdffedb..79b2e709 100644 --- a/doc/source/namelists/initial_conditions.nml.rst +++ b/doc/source/namelists/initial_conditions.nml.rst @@ -402,18 +402,18 @@ The required variables for a particular configuration, along with their 'type' a | Required if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE, :nml:mem:`JULES_RIVERS::i_river_vn` = '2' and | | :nml:mem:`JULES_INITIAL::dump_file` = TRUE | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ -| ``rfm_surfstore_rp`` | Surface water storage on river routing points (m3) | none | +| ``rfm_surfstore_rp`` | Surface water storage on river routing points (m3). | none | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ -| ``rfm_substore_rp`` | Sub-surface water storage on river routing points (m3) | none | +| ``rfm_substore_rp`` | Sub-surface water storage on river routing points (m3). | none | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ -| ``rfm_flowin_rp`` | Surface flow into a grid box on river routing points (m3) | none | +| ``rfm_flowin_rp`` | Surface flow into a grid box on river routing points (m3). | none | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ -| ``rfm_bflowin_rp`` | Sub-surface flow into a grid box on river routing points (m3) | none | +| ``rfm_bflowin_rp`` | Sub-surface flow into a grid box on river routing points (m3). | none | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ | Required if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE, :nml:mem:`JULES_RIVERS::i_river_vn` = '1,3' and | | :nml:mem:`JULES_INITIAL::dump_file` = TRUE | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ -| ``rivers_sto_rp`` | Water storage (kg) | none | +| ``rivers_sto_rp`` | Water storage (kg). | none | | | | | | | **Alternatively**, can be initialised from an ancillary file via | | | | :nml:lst:`JULES_RIVERS_PROPS`. It cannot be requested via :nml:lst:`JULES_INITIAL` | | @@ -423,7 +423,12 @@ The required variables for a particular configuration, along with their 'type' a | :nml:mem:`JULES_INITIAL::dump_file` = TRUE and :nml:mem:`OASIS_RIVERS::send_fields` or | | :nml:mem:`JULES_OUTPUT_PROFILE::var` contains ``outflow_per_river``. | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ -| ``rivers_outflow_rp`` | River outflow on river routing points (kg s\ :sup:`-1`) | none | +| ``rivers_outflow_rp`` | River outflow on river routing points (kg s\ :sup:`-1`). | none | ++----------------------------------+-----------------------------------------------------------------------------------------+---------+ +| Required if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE, :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE and | +| :nml:mem:`JULES_INITIAL::dump_file` = TRUE | ++----------------------------------+-----------------------------------------------------------------------------------------+---------+ +| ``res_storage`` | Water stored in reservoirs, on river routing points (kg). | none | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ | Required if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE, :nml:mem:`JULES_RIVERS::i_river_vn` = '3', | | :nml:mem:`JULES_INITIAL::dump_file` = TRUE and :nml:mem:`OASIS_RIVERS::send_fields` contains ``inland_outflow``. | @@ -432,10 +437,9 @@ The required variables for a particular configuration, along with their 'type' a +----------------------------------+-----------------------------------------------------------------------------------------+---------+ | Required if :nml:mem:`JULES_VEGETATION::photo_acclim_model` = 2 or 3 | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ -| ``t_growth_gb`` | Running mean air temperature (K) | none | +| ``t_growth_gb`` | Running mean air temperature (K). | none | +----------------------------------+-----------------------------------------------------------------------------------------+---------+ - .. warning:: if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE, :nml:mem:`JULES_RIVERS::i_river_vn` = '2' and :nml:mem:`JULES_INITIAL::dump_file` = FALSE, ``rfm_surfstore_rp``, ``rfm_substore_rp``, ``rfm_flowin_rp`` and ``rfm_bflowin`` are initialised to zero. @@ -444,6 +448,10 @@ The required variables for a particular configuration, along with their 'type' a if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE, :nml:mem:`JULES_RIVERS::i_river_vn` = '1,3' and :nml:mem:`JULES_INITIAL::dump_file` = FALSE, ``rivers_sto_rp`` is initialised to zero **unless** it is initialised from an ancillary via :nml:lst:`JULES_RIVERS_PROPS`. +.. warning:: + if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE, :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE and :nml:mem:`JULES_INITIAL::dump_file` = FALSE, + ``res_storage`` is initialised to zero. + .. |mu| unicode:: μ .. u diff --git a/doc/source/namelists/jules_rivers.nml.rst b/doc/source/namelists/jules_rivers.nml.rst index f5c5cc0f..b7f6d97f 100644 --- a/doc/source/namelists/jules_rivers.nml.rst +++ b/doc/source/namelists/jules_rivers.nml.rst @@ -44,6 +44,19 @@ River routing introduces two more grids to a JULES run: the river routing input ``3`` Use a standalone JULES implementation of the TRIP model (see Oki et al. 1999). +.. nml:member:: l_reservoirs + + :type: logical + :default: F + + Switch for enabling routing through reservoirs. Only used if :nml:mem:`JULES_RIVERS::l_rivers` is TRUE. + + TRUE + River flow is routed through reservoirs. + + FALSE + No reservoirs are modelled. + .. nml:member:: l_riv_overbank :type: logical diff --git a/doc/source/output-variables.rst b/doc/source/output-variables.rst index 1e337f3c..78bb6822 100644 --- a/doc/source/output-variables.rst +++ b/doc/source/output-variables.rst @@ -471,13 +471,49 @@ Rivers | | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | | | :nml:mem:`JULES_RIVERS::i_river_vn` = 2. | | +-------------------------------+-----------------------------------------------------------------------------------------------+------------+ -| ``rivers_sto_rp`` | River routing gridbox river storage (kg) | np_rivers | +| ``rivers_sto_rp`` | River routing gridbox river storage (kg). | np_rivers | | | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | | | :nml:mem:`JULES_RIVERS::i_river_vn` = 3. | | +-------------------------------+-----------------------------------------------------------------------------------------------+------------+ | ``frac_fplain_rp`` | Overbank inundation area as a fraction of river routing gridcell area. | np_rivers | | | Only available if :nml:mem:`JULES_RIVERS::l_riv_overbank` = TRUE. | | +-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_capacity`` | Storage capacity of reservoirs (kg). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_catch`` | Catchment area of reservoirs (m\ :sup:`2`). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_year`` | Initialisation year of reservoirs. | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_critical`` | Critical storage of reservoirs (kg). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_flood`` | Flood storage of reservoirs (kg). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_emergency`` | Emergency storage of reservoirs (kg). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_normal_release`` | Normal release rate of reservoirs (kg s\ :sup:`-1`). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_flood_release`` | Flood release rate of reservoirs (kg s\ :sup:`-1`). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ +| ``res_storage`` | Water stored in reservoirs (kg). | np_rivers | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE and | | +| | :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++-------------------------------+-----------------------------------------------------------------------------------------------+------------+ | Output regridded to the JULES model grid | +-------------------------------+-----------------------------------------------------------------------------------------------+------------+ | ``rflow`` | Gridbox river flow rate (kg m\ :sup:`-2` s\ :sup:`-1`). | | diff --git a/rose-meta/jules-standalone/HEAD/rose-meta.conf b/rose-meta/jules-standalone/HEAD/rose-meta.conf index e30a27c5..5ac8bf83 100644 --- a/rose-meta/jules-standalone/HEAD/rose-meta.conf +++ b/rose-meta/jules-standalone/HEAD/rose-meta.conf @@ -2425,13 +2425,13 @@ sort-key=06 url=https://metoffice.github.io/jules/latest/namelists/initial_conditions.nml.html#JULES_INITIAL::var values='canht','canopy','cropcanht','cropdvi','cropharvc','croplai', ='cropreservec','croprootc','cs','frac','frac_agr_prev','frac_past_prev', - ='frac_biocrop_prev','gs','lai','n_inorg','nsnow','ns','rfm_bflowin_rp', - ='rfm_flowin_rp','rfm_substore_rp','rfm_surfstore_rp','rgrain','rgrainl', - ='rho_snow','rivers_sto_rp','rivers_outflow_rp','snow_depth','snow_ds', - ='snow_grnd','snow_ice','snow_liq','snow_tile','sthu_irr','sthuf', - ='sthzw','t_growth_gb','t_soil','tsnow','tstar_tile','tsoil_deep', - ='wood_prod_fast','wood_prod_med','wood_prod_slow','years_since_harvest', - ='zw','inland_outflow_rp' + ='frac_biocrop_prev','gs','lai','n_inorg','nsnow','ns','res_storage', + ='rfm_bflowin_rp','rfm_flowin_rp','rfm_substore_rp','rfm_surfstore_rp', + ='rgrain','rgrainl','rho_snow','rivers_sto_rp','rivers_outflow_rp', + ='snow_depth','snow_ds','snow_grnd','snow_ice','snow_liq','snow_tile', + ='sthu_irr','sthuf','sthzw','t_growth_gb','t_soil','tsnow','tstar_tile', + ='tsoil_deep','wood_prod_fast','wood_prod_med','wood_prod_slow', + ='years_since_harvest','zw','inland_outflow_rp' [namelist:jules_initial=var_name] compulsory=true @@ -4076,6 +4076,13 @@ value-titles=RFM,TRIP values=2,3 widget[rose-config-edit]=cylc8_compat.ComboBoxValueWidget +[namelist:jules_rivers=l_reservoirs] +compulsory=true +description=Switch to activate reservoirs +sort-key=a8 +type=logical +url=https://metoffice.github.io/jules/latest/namelists/jules_rivers.nml.html#JULES_RIVERS::l_reservoirs + [namelist:jules_rivers=l_riv_overbank] compulsory=true description=Enable overbank inundation @@ -4094,6 +4101,7 @@ trigger=namelist:jules_rivers=i_river_vn: .true.; = namelist:jules_rivers=l_riv_overbank: .true.; = namelist:jules_overbank: .true.; = namelist:jules_rivers_props: .true.; + = namelist:jules_rivers=l_reservoirs: .true.; type=logical url=https://metoffice.github.io/jules/latest/namelists/jules_rivers.nml.html#JULES_RIVERS::l_rivers diff --git a/rose-meta/jules-standalone/versions.py b/rose-meta/jules-standalone/versions.py index 1d4c41a9..2e442c0d 100644 --- a/rose-meta/jules-standalone/versions.py +++ b/rose-meta/jules-standalone/versions.py @@ -45,15 +45,16 @@ from .version81_82 import * -class vnYY_txxxx(MacroUpgrade): +class vn82_txxxx(MacroUpgrade): - """Upgrade macro from JULES by Author""" + """Upgrade macro from JULES by Helen Baron""" - BEFORE_TAG = "vnY.Y" - AFTER_TAG = "vnY.Y_txxxx" + BEFORE_TAG = "vn8.2" + AFTER_TAG = "vn8.2_txxxx" def upgrade(self, config, meta_config=None): """Upgrade a JULES runtime app configuration.""" # Add settings + self.add_setting(config, ["namelist:jules_rivers", "l_reservoirs"], ".false.") return config, self.reports diff --git a/src/control/shared/jules_rivers_mod.F90 b/src/control/shared/jules_rivers_mod.F90 index 38b019f9..c6c3836c 100644 --- a/src/control/shared/jules_rivers_mod.F90 +++ b/src/control/shared/jules_rivers_mod.F90 @@ -142,9 +142,13 @@ MODULE jules_rivers_mod ! - OASIS send field contains 'inland_outflow' ! - The diagnostic 'inland_outflow_rp' has been requested ! - UM-TRIP when l_inland = T - ,l_init_storage = .FALSE. + ,l_init_storage = .FALSE. & ! Set to true if an initial river storage ancillary file is ! used + ,l_reservoirs = .FALSE. + ! Switch for reservoirs. + ! .TRUE. = consider major reservoirs + ! .FALSE. = do not consider reservoirs INTEGER :: & nstep_rivers = imdi & @@ -238,7 +242,8 @@ MODULE jules_rivers_mod ! Single namelist definition for UM and standalone !------------------------------------------------------------------------------ NAMELIST /jules_rivers/ & - l_rivers, l_riv_overbank, l_adapt_timestep, l_sea_level, & + l_rivers, l_reservoirs, l_riv_overbank, & + l_adapt_timestep, l_sea_level, & l_vary_sea_level, i_river_vn, nstep_rivers, & trip_globe_shape, & cland, criver, cbland, cbriver, runoff_factor, retl, retr, & @@ -408,6 +413,27 @@ MODULE jules_rivers_mod river_manning_grid(:,:) ! Manning roughness coefficient for river channel (1). +!------------------------------------------------------------------------------ +! Ancillary arrays for reservoirs, defined on 2D rivers grid. +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), ALLOCATABLE :: & + res_capacity_grid(:,:), & + ! Storage capacity of reservoirs (kg). + res_catch_grid(:,:), & + ! Catchment area of reservoirs (m2). + res_year_grid(:,:), & + ! Initialisation year of reservoirs. + res_critical_grid(:,:), & + ! Critical storage of reservoir (kg). + res_flood_grid(:,:), & + ! Flood storage of reservoir (kg). + res_emergency_grid(:,:), & + ! Emergency storage of reservoir (kg). + res_normal_release_grid(:,:), & + ! Normal release rate of reservoir (kg s-1). + res_flood_release_grid(:,:) + ! Flood release rate of reservoir (kg s-1). + CHARACTER(LEN=*), PARAMETER, PRIVATE :: ModuleName='JULES_RIVERS_MOD' !------------------------------------------------------------------------------ @@ -596,6 +622,43 @@ MODULE jules_rivers_mod ! River outflow into the ocean (kg s-1) REAL(KIND=real_jlslsm), ALLOCATABLE :: inland_outflow_rp(:) ! Inland basin flow into soil moisture (kg m-2 s-1) + + !---------------------------------------------------------------------------- + ! Reservoir ancillary variables, defined on river points. + !---------------------------------------------------------------------------- + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_capacity(:) + ! Storage capacity of reservoirs (kg). + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_catch(:) + ! Catchment area of reservoirs (m). + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_year(:) + ! Initialisation year of reservoirs. + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_critical(:) + ! Critical storage of reservoir (kg). + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_flood(:) + ! Flood storage of reservoir (kg). + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_emergency(:) + ! Emergency storage of reservoir (kg). + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_normal_release(:) + ! Normal release rate of reservoir (kg s-1). + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_flood_release(:) + ! Flood release rate of reservoir (kg s-1). + + !---------------------------------------------------------------------------- + ! Reservoir prognostic variables, defined on river points. + !---------------------------------------------------------------------------- + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_storage(:) + ! Water stored in reservoirs (kg). + REAL(KIND=real_jlslsm), ALLOCATABLE :: res_cap_current(:) + ! Capacity of currently active reservoirs (kg). + + ! Variables for coupling with water resources. + !---------------------------------------------------------------------------- + REAL, ALLOCATABLE :: tot_abstracted_res_global(:) + ! Water abstracted from reservoirs over river + ! timestep, on global land points (kg). + REAL, ALLOCATABLE :: tot_net_abstracted_river_global(:) + ! Water abstracted from rivers over river timestep, + ! on global land points(kg m-2). END TYPE rivers_data_type TYPE :: rivers_type @@ -664,13 +727,26 @@ MODULE jules_rivers_mod REAL(KIND=real_jlslsm), POINTER :: rivers_boxareas_rp(:) REAL(KIND=real_jlslsm), POINTER :: rivers_outflow_rp(:) REAL(KIND=real_jlslsm), POINTER :: inland_outflow_rp(:) + REAL(KIND=real_jlslsm), POINTER :: res_capacity(:) + REAL(KIND=real_jlslsm), POINTER :: res_catch(:) + REAL(KIND=real_jlslsm), POINTER :: res_storage(:) + REAL(KIND=real_jlslsm), POINTER :: res_cap_current(:) + REAL(KIND=real_jlslsm), POINTER :: res_year(:) + REAL(KIND=real_jlslsm), POINTER :: res_critical(:) + REAL(KIND=real_jlslsm), POINTER :: res_flood(:) + REAL(KIND=real_jlslsm), POINTER :: res_emergency(:) + REAL(KIND=real_jlslsm), POINTER :: res_normal_release(:) + REAL(KIND=real_jlslsm), POINTER :: res_flood_release(:) + REAL(KIND=real_jlslsm), POINTER :: tot_abstracted_res_global(:) + REAL(KIND=real_jlslsm), POINTER :: tot_net_abstracted_river_global(:) END TYPE rivers_type CONTAINS !############################################################################## -SUBROUTINE jules_rivers_alloc(land_pts, t_i_length, t_j_length, rivers_data) +SUBROUTINE jules_rivers_alloc(land_pts, t_i_length, t_j_length, & + sw_river_source, l_water_resources, rivers_data) !No USE statements other than Dr Hook USE parkind1, ONLY: jprb, jpim @@ -680,6 +756,10 @@ SUBROUTINE jules_rivers_alloc(land_pts, t_i_length, t_j_length, rivers_data) !Arguments INTEGER, INTENT(IN) :: land_pts, t_i_length, t_j_length +INTEGER, INTENT(IN) :: sw_river_source + ! Index of river water in surface water source arrays. +LOGICAL, INTENT(IN) :: l_water_resources + ! Switch to select water resource management modelling. TYPE(rivers_data_type), INTENT(IN OUT) :: rivers_data !Local variables @@ -1174,6 +1254,8 @@ SUBROUTINE print_nlist_jules_rivers() WRITE(lineBuffer,*)' l_rivers = ',l_rivers CALL jules_print('jules_rivers',lineBuffer) +WRITE(lineBuffer,*)' l_reservoirs = ',l_reservoirs +CALL jules_print('jules_rivers',lineBuffer) WRITE(lineBuffer,*)' l_riv_overbank = ',l_riv_overbank CALL jules_print('jules_rivers',lineBuffer) WRITE(lineBuffer,*)' i_river_vn = ',i_river_vn @@ -1315,7 +1397,7 @@ SUBROUTINE read_nml_jules_rivers(unit_in) INTEGER, PARAMETER :: no_of_types = 3 INTEGER, PARAMETER :: n_int = 5 INTEGER, PARAMETER :: n_real = 12 -INTEGER, PARAMETER :: n_log = 5 +INTEGER, PARAMETER :: n_log = 6 TYPE :: my_namelist SEQUENCE @@ -1337,6 +1419,7 @@ SUBROUTINE read_nml_jules_rivers(unit_in) REAL(KIND=real_jlslsm) :: rivers_speed REAL(KIND=real_jlslsm) :: runoff_factor LOGICAL :: l_adapt_timestep + LOGICAL :: l_reservoirs LOGICAL :: l_riv_overbank LOGICAL :: l_rivers LOGICAL :: l_sea_level @@ -1375,6 +1458,7 @@ SUBROUTINE read_nml_jules_rivers(unit_in) my_nml % rivers_speed = rivers_speed my_nml % runoff_factor = runoff_factor my_nml % l_adapt_timestep = l_adapt_timestep + my_nml % l_reservoirs = l_reservoirs my_nml % l_riv_overbank = l_riv_overbank my_nml % l_rivers = l_rivers my_nml % l_sea_level = l_sea_level @@ -1402,6 +1486,7 @@ SUBROUTINE read_nml_jules_rivers(unit_in) rivers_speed = my_nml % rivers_speed runoff_factor = my_nml % runoff_factor l_adapt_timestep = my_nml % l_adapt_timestep + l_reservoirs = my_nml % l_reservoirs l_riv_overbank = my_nml % l_riv_overbank l_rivers = my_nml % l_rivers l_sea_level = my_nml % l_sea_level diff --git a/src/initialisation/rivers-standalone/allocate_river_arrays_mod.F90 b/src/initialisation/rivers-standalone/allocate_river_arrays_mod.F90 index 8c9bdf65..e67a3c87 100644 --- a/src/initialisation/rivers-standalone/allocate_river_arrays_mod.F90 +++ b/src/initialisation/rivers-standalone/allocate_river_arrays_mod.F90 @@ -95,6 +95,12 @@ SUBROUTINE allocate_river_arrays(psparms_data,ainfo_data, progs_data, & INTEGER, PARAMETER :: river_row_length_dum = 1 INTEGER, PARAMETER :: river_rows_dum = 1 +! Water resources is not available with rivers-standalone. Set values so that +! water resource variables are not allocated. +INTEGER, PARAMETER :: tmp_sw_river_source = 0 +LOGICAL, PARAMETER :: tmp_l_water_resources = .FALSE. +LOGICAL, PARAMETER :: tmp_l_reservoirs = .FALSE. + INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 REAL(KIND=jprb) :: zhook_handle @@ -152,7 +158,9 @@ SUBROUTINE allocate_river_arrays(psparms_data,ainfo_data, progs_data, & nice, nsoilt, ntype, & ainfo_data) -CALL jules_rivers_alloc(temp_size, t_i_length, t_j_length, rivers_data) +CALL jules_rivers_alloc(temp_size, t_i_length, t_j_length, & + tmp_sw_river_source, tmp_l_water_resources, & + rivers_data) CALL coastal_alloc(temp_size, t_i_length, t_j_length, & u_i_length, u_j_length, & diff --git a/src/initialisation/shared/allocate_jules_arrays.F90 b/src/initialisation/shared/allocate_jules_arrays.F90 index 533ff2f7..c6c716b1 100644 --- a/src/initialisation/shared/allocate_jules_arrays.F90 +++ b/src/initialisation/shared/allocate_jules_arrays.F90 @@ -31,8 +31,8 @@ SUBROUTINE allocate_jules_arrays(crop_vars_data,psparms_data,top_pdm_data, & lake_data, & forcing_data, & rivers_data, & - !veg3_parm_(data), & - !veg3_field_(data), & + !veg3_parm_(data), & + !veg3_field_(data), & chemvars_data, water_resources_data, & wtrac_jls_data & ) @@ -52,14 +52,20 @@ SUBROUTINE allocate_jules_arrays(crop_vars_data,psparms_data,top_pdm_data, & l_water_domestic, l_water_industry, & l_water_irrigation,l_water_livestock, & l_water_resources, l_water_transfers, & - n_sw_source, nwater_use + n_sw_source, nwater_use, sw_river_source USE jules_deposition_mod, ONLY: l_deposition, ndry_dep_species USE jules_water_tracers_mod, ONLY: l_wtrac_jls +USE jules_rivers_mod, ONLY: l_reservoirs !Variables- dimensions USE jules_surface_types_mod, ONLY: ncpft,nnpft USE jules_snow_mod, ONLY: nsmax, cansnowtile USE jules_surface_types_mod, ONLY: npft, nnvg, ntype +#if defined(UM_JULES) +USE atm_land_sea_mask, ONLY: global_land_pts => atmos_number_of_landpts +#else +USE model_grid_mod, ONLY: global_land_pts +#endif USE theta_field_sizes, ONLY: t_i_length, t_j_length, & u_i_length,u_j_length, & v_i_length,v_j_length @@ -142,6 +148,8 @@ SUBROUTINE allocate_jules_arrays(crop_vars_data,psparms_data,top_pdm_data, & USE Field_Types, ONLY: fld_type_r #endif +USE parallel_mod, ONLY: is_master_task + USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -265,7 +273,8 @@ SUBROUTINE allocate_jules_arrays(crop_vars_data,psparms_data,top_pdm_data, & nice,nsoilt,ntype, & ainfo_data) -CALL jules_rivers_alloc(land_pts, t_i_length, t_j_length, rivers_data) +CALL jules_rivers_alloc(land_pts, t_i_length, t_j_length, & + sw_river_source, l_water_resources, rivers_data) CALL forcing_alloc(t_i_length,t_j_length, u_i_length, u_j_length, & v_i_length, v_j_length, forcing_data) @@ -290,12 +299,13 @@ SUBROUTINE allocate_jules_arrays(crop_vars_data,psparms_data,top_pdm_data, & CALL deposition_species_alloc(ntype, ndry_dep_species, l_deposition) -CALL water_resources_alloc( land_pts, n_sw_source, nwater_use, & - l_have_groundwater, l_have_surface_water, & - l_water_domestic, l_water_industry, & - l_water_irrigation, l_water_livestock, & - l_water_resources, l_water_transfers, & - water_resources_data ) +CALL water_resources_alloc( global_land_pts, land_pts, n_sw_source, & + nwater_use, sw_river_source, l_have_groundwater, & + l_have_surface_water, is_master_task(), & + l_reservoirs, l_water_domestic, l_water_industry, & + l_water_irrigation, & + l_water_livestock, l_water_resources, & + l_water_transfers, water_resources_data ) ! Set up local river grid sizes (Note, water tracers only work in UM_JULES) #if defined(UM_JULES) diff --git a/src/initialisation/standalone/ancillaries/init_rivers_process_data_mod.F90 b/src/initialisation/standalone/ancillaries/init_rivers_process_data_mod.F90 index d2b50cfe..bfa55e2f 100644 --- a/src/initialisation/standalone/ancillaries/init_rivers_process_data_mod.F90 +++ b/src/initialisation/standalone/ancillaries/init_rivers_process_data_mod.F90 @@ -1029,9 +1029,11 @@ SUBROUTINE remap_ancil( nx_rivers, ny_rivers, rivers_dx, & USE jules_rivers_mod, ONLY: & channel_depth_grid, channel_width_grid, mean_sea_level_grid, & - river_distance_grid, river_elevation_grid, & - river_elevation_grid, river_length_grid, river_manning_grid, & - river_nextx_grid, river_nexty_grid, rivers_type + river_distance_grid, river_elevation_grid, river_length_grid, & + river_manning_grid, river_nextx_grid, river_nexty_grid, & + res_capacity_grid, res_catch_grid, res_year_grid, & + res_critical_grid, res_flood_grid, res_emergency_grid, & + res_normal_release_grid, res_flood_release_grid, rivers_type USE overbank_inundation_mod, ONLY: & hypsometric_quantiles_grid, logn_mean, logn_stdev, nquantile_hypso @@ -1085,7 +1087,7 @@ SUBROUTINE remap_ancil( nx_rivers, ny_rivers, rivers_dx, & SELECT CASE ( var ) !---------------------------------------------------------------------------- - ! Cases for river (not overbank) variables. + ! Cases for river variables. !---------------------------------------------------------------------------- CASE ( 'area' ) CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & @@ -1196,6 +1198,41 @@ SUBROUTINE remap_ancil( nx_rivers, ny_rivers, rivers_dx, & ! Extract the first column from the 2D field. rivers%rivers_ygrid(:) = tmp_2d_coord(1,:) + !---------------------------------------------------------------------------- + ! Cases for reservoir variables. + !---------------------------------------------------------------------------- +CASE ( 'res_capacity_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_capacity_grid ) + +CASE ( 'res_catch_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_catch_grid ) + +CASE ( 'res_year_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_year_grid ) + +CASE ( 'res_critical_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_critical_grid ) + +CASE ( 'res_flood_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_flood_grid ) + +CASE ( 'res_emergency_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_emergency_grid ) + +CASE ( 'res_normal_release_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_normal_release_grid ) + +CASE ( 'res_flood_release_grid' ) + CALL remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, l_reverse_y, & + rivers%rivers_xgrid, res_flood_release_grid ) + CASE DEFAULT CALL log_fatal( RoutineName, & "Do not recognise var: " // TRIM(var) ) @@ -1250,7 +1287,7 @@ SUBROUTINE remap_field( nx_rivers, ny_rivers, rivers_dx, l_shift_x, & !------------------------------------------------------------------------------ LOGICAL, OPTIONAL, INTENT(IN) :: & is_x_index_arg - ! Flag indicating if the field is an index of x location on grid.& + ! Flag indicating if the field is an index of x location on grid. & !------------------------------------------------------------------------------ ! Local scalar variables. @@ -2216,12 +2253,16 @@ SUBROUTINE set_river_point_values( rivers_x1_input, & USE jules_model_environment_mod, ONLY: l_oasis_rivers USE jules_rivers_mod, ONLY: & - a_thresh, channel_depth_grid, channel_width_grid, i_river_vn, l_sea_level, & - l_riv_overbank, mean_sea_level_grid, np_rivers, nseqmax, nx_rivers, & - ny_rivers, rfm_land, rfm_river, rivers_camaflood, river_distance_grid, & + a_thresh, channel_depth_grid, channel_width_grid, i_river_vn, & + l_sea_level, l_riv_overbank, mean_sea_level_grid, & + np_rivers, nseqmax, nx_rivers, ny_rivers, rfm_land, rfm_river, & + rivers_camaflood, river_distance_grid, & rivers_dx, rivers_dy, river_elevation_grid, river_length_grid, & river_manning_grid, rivers_rfm, rivers_trip, rivers_x1, l_outflow_per_river, & - l_init_storage, & + l_init_storage, l_reservoirs, & + res_capacity_grid, res_catch_grid, res_year_grid, & + res_critical_grid, res_flood_grid, res_emergency_grid, & + res_normal_release_grid, res_flood_release_grid, & ! types rivers_type @@ -2393,8 +2434,8 @@ SUBROUTINE set_river_point_values( rivers_x1_input, & rivers%channel_width(ip) = channel_width_grid(irx,iry) rivers%river_distance(ip) = river_distance_grid(irx,iry) rivers%river_elevation(ip) = river_elevation_grid(irx,iry) - rivers% river_length(ip) = river_length_grid(irx,iry) - rivers% river_manning(ip) = river_manning_grid(irx,iry) + rivers%river_length(ip) = river_length_grid(irx,iry) + rivers%river_manning(ip) = river_manning_grid(irx,iry) IF ( l_sea_level ) THEN rivers%mean_sea_level(ip) = mean_sea_level_grid(irx,iry) END IF @@ -2462,6 +2503,24 @@ SUBROUTINE set_river_point_values( rivers_x1_input, & END IF ! l_riv_overbank + !------------------------------------------------------------------------ + ! Set reservoir ancillary variables, if required. + !------------------------------------------------------------------------ + IF ( l_reservoirs ) THEN + + rivers%res_capacity(ip) = res_capacity_grid(irx,iry) + rivers%res_catch(ip) = res_catch_grid(irx,iry) + rivers%res_year(ip) = res_year_grid(irx,iry) + rivers%res_critical(ip) = res_critical_grid(irx,iry) + rivers%res_flood(ip) = res_flood_grid(irx,iry) + rivers%res_emergency(ip) = res_emergency_grid(irx,iry) + rivers%res_normal_release(ip) = & + res_normal_release_grid(irx,iry) + rivers%res_flood_release(ip) = & + res_flood_release_grid(irx,iry) + + END IF ! l_reservoirs + END IF ! point_mask END DO ! iy diff --git a/src/initialisation/standalone/ancillaries/init_rivers_props_mod.F90 b/src/initialisation/standalone/ancillaries/init_rivers_props_mod.F90 index 541870ba..4e3207ee 100644 --- a/src/initialisation/standalone/ancillaries/init_rivers_props_mod.F90 +++ b/src/initialisation/standalone/ancillaries/init_rivers_props_mod.F90 @@ -47,9 +47,11 @@ SUBROUTINE init_rivers_props(rivers, rivers_data) USE io_constants, ONLY: max_file_name_len, namelist_unit USE jules_rivers_mod, ONLY: & - i_river_vn, l_rivers, l_sea_level, rivers_camaflood, rivers_rfm, & + i_river_vn, l_rivers, l_sea_level, & + rivers_camaflood, rivers_rfm, & rivers_trip, y1_land_grid, l_riv_overbank, l_outflow_per_river, & l_init_storage, & + l_reservoirs, & ! types rivers_type, rivers_data_type @@ -307,6 +309,28 @@ SUBROUTINE init_rivers_props(rivers, rivers_data) required_vars(nvars_required) = 'rivers_outflow_number' END IF + !---------------------------------------------------------------------------- + ! Add variables for reservoirs. + !---------------------------------------------------------------------------- + IF ( l_reservoirs ) THEN + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_capacity_grid' + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_catch_grid' + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_year_grid' + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_critical_grid' + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_flood_grid' + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_emergency_grid' + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_normal_release_grid' + nvars_required = nvars_required + 1 + required_vars(nvars_required) = 'res_flood_release_grid' + END IF + nvars_optional = nvars_optional + 1 optional_vars(nvars_optional) = 'rivers_outflow_number' diff --git a/src/initialisation/standalone/ancillaries/jules_rivers_props_mod.F90 b/src/initialisation/standalone/ancillaries/jules_rivers_props_mod.F90 index d847d1a0..c2e01568 100644 --- a/src/initialisation/standalone/ancillaries/jules_rivers_props_mod.F90 +++ b/src/initialisation/standalone/ancillaries/jules_rivers_props_mod.F90 @@ -35,7 +35,7 @@ MODULE jules_rivers_props_mod !----------------------------------------------------------------------------- ! Scalars in namelist. -INTEGER, PARAMETER :: nvars_max = 9 +INTEGER, PARAMETER :: nvars_max = 17 ! The maximum possible number of ancillary variables that are read in. INTEGER :: nvars @@ -69,7 +69,7 @@ MODULE jules_rivers_props_mod l_ignore_ancil_rivers_check = .FALSE. ! Switch to check river routing & river number ancillary for compatibility. ! The current Rose stem test ancils are not compatible as the land sea - ! mask is different. Setting this to true maintains bit comparability & + ! mask is different. Setting this to true maintains bit comparability & ! protects the calculation of the coupling field needed for the LFRic ! coupled miniapp. @@ -165,11 +165,11 @@ MODULE jules_rivers_props_mod SUBROUTINE check_jules_rivers_props() USE jules_rivers_mod, ONLY: nx_rivers_in=>nx_rivers, ny_rivers_in=>ny_rivers, & - i_river_vn, land_dx, land_dy, l_riv_overbank, & - nx_land_grid, ny_land_grid, x1_land_grid, & - y1_land_grid, rivers_length, l_outflow_per_river, & - rivers_camaflood, rivers_regrid, rivers_rfm, & - l_init_storage + i_river_vn, land_dx, land_dy, l_init_storage, & + l_riv_overbank, nx_land_grid, ny_land_grid, & + x1_land_grid, y1_land_grid, rivers_length, & + l_outflow_per_river, rivers_camaflood, & + rivers_regrid, rivers_rfm, l_reservoirs USE model_grid_mod, ONLY: l_coord_latlon @@ -190,7 +190,7 @@ SUBROUTINE check_jules_rivers_props() INTEGER :: i ! Loop counter INTEGER, PARAMETER :: & - nvar_values = 9 ! Number of recognised identifiers + nvar_values = 17 ! Number of recognised identifiers CHARACTER(LEN=identifier_len) :: & var_values(nvar_values) ! Names of recognised identifiers @@ -206,15 +206,23 @@ SUBROUTINE check_jules_rivers_props() ny_rivers_in = ny_rivers ! Check that the requested identifiers are recognised. -var_values(:) = ['area ', & - 'direction ', & - 'sequence ', & - 'latitude_2d ', & - 'longitude_2d ', & - 'rivers_outflow_number', & - 'rivers_storage ', & - 'logn_mean ', & - 'logn_stdev '] +var_values(:) = ['area ', & + 'direction ', & + 'sequence ', & + 'latitude_2d ', & + 'longitude_2d ', & + 'rivers_outflow_number ', & + 'rivers_storage ', & + 'logn_mean ', & + 'logn_stdev ', & + 'res_capacity_grid ', & + 'res_catch_grid ', & + 'res_year_grid ', & + 'res_critical_grid ', & + 'res_flood_grid ', & + 'res_emergency_grid ', & + 'res_normal_release_grid', & + 'res_flood_release_grid '] DO i = 1, nvars IF ( .NOT. ANY( var_values(:) == var(i) ) ) THEN @@ -394,7 +402,12 @@ SUBROUTINE allocate_river_vars_grid( global_land_pts, nx_rivers, ny_rivers, & river_length_grid, river_manning_grid, & river_nextx_grid, river_nexty_grid, rivers_rfm, & rivers_trip, rivers_data_type, rivers_type, & - l_outflow_per_river, l_init_storage + l_outflow_per_river, l_init_storage, & + l_reservoirs, res_capacity_grid, res_catch_grid, & + res_year_grid, res_critical_grid, & + res_flood_grid, res_emergency_grid, & + res_normal_release_grid, & + res_flood_release_grid IMPLICIT NONE @@ -423,10 +436,10 @@ SUBROUTINE allocate_river_vars_grid( global_land_pts, nx_rivers, ny_rivers, & CHARACTER(LEN=*), PARAMETER :: RoutineName = 'ALLOCATE_RIVER_VARS_GRID' INTEGER :: & - ERROR, error_sum, & - ! Error values. - nx_size, ny_size - ! Sizes used for allocations. + ERROR, error_sum, & + ! Error values. + nx_size, ny_size + ! Sizes used for allocations. !end of header !------------------------------------------------------------------------------ @@ -548,7 +561,7 @@ SUBROUTINE allocate_river_vars_grid( global_land_pts, nx_rivers, ny_rivers, & ALLOCATE( rivers_data%rivers_outflow_number(nx_size,ny_size), STAT = ERROR ) error_sum = error_sum + ERROR - +! River storage. IF ( l_init_storage ) THEN nx_size = nx_rivers ny_size = ny_rivers @@ -559,6 +572,31 @@ SUBROUTINE allocate_river_vars_grid( global_land_pts, nx_rivers, ny_rivers, & ALLOCATE( rivers_data%rivers_storage(nx_size,ny_size), STAT = ERROR ) error_sum = error_sum + ERROR +! Reservoir 2D ancillary variables. +IF ( l_reservoirs ) THEN + nx_size = nx_rivers + ny_size = ny_rivers +ELSE + nx_size = 1 + ny_size = 1 +END IF +ALLOCATE( res_capacity_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( res_catch_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( res_year_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( res_critical_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( res_flood_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( res_emergency_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( res_normal_release_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( res_flood_release_grid(nx_size,ny_size), STAT = ERROR ) +error_sum = error_sum + ERROR + IF ( error_sum /= 0 ) THEN CALL log_fatal( RoutineName, "Error allocating for rivers arrays." ) END IF @@ -585,7 +623,15 @@ SUBROUTINE allocate_river_vars_grid( global_land_pts, nx_rivers, ny_rivers, & rivers_data%rivers_seq(:,:) = rmdi rivers_data%rivers_outflow_number(:,:) = rmdi rivers_data%rivers_storage(:,:) = rmdi -rivers_data%land_fraction_2d(:,:) = rmdi +rivers_data%land_fraction_2d(:,:) = rmdi +res_capacity_grid(:,:) = rmdi +res_catch_grid(:,:) = rmdi +res_year_grid(:,:) = rmdi +res_critical_grid(:,:) = rmdi +res_flood_grid(:,:) = rmdi +res_emergency_grid(:,:) = rmdi +res_normal_release_grid(:,:) = rmdi +res_flood_release_grid(:,:) = rmdi ! Associate pointers rivers%rivers_dir => rivers_data%rivers_dir @@ -608,16 +654,25 @@ SUBROUTINE allocate_rivers_vars_rp( np_rivers, rivers, rivers_data ) !------------------------------------------------------------------------------ ! Description: ! Allocate river point variables, initialise, and associate pointers. -! Also allocates at least one variable on land points. +! Also allocates related variables on land points. !------------------------------------------------------------------------------ USE ancil_info, ONLY: land_pts +#if defined(UM_JULES) +USE atm_land_sea_mask, ONLY: global_land_pts => atmos_number_of_landpts +#else +USE model_grid_mod, ONLY: global_land_pts +#endif + USE jules_model_environment_mod, ONLY: l_oasis_rivers USE jules_rivers_mod, ONLY: i_river_vn, l_sea_level, l_vary_sea_level, & rivers_camaflood, rivers_data_type, rivers_rfm, & - rivers_trip, rivers_type, l_outflow_per_river + rivers_trip, rivers_type, l_outflow_per_river, & + l_reservoirs + +USE jules_water_resources_mod, ONLY: l_water_resources USE missing_data_mod, ONLY: imdi, rmdi @@ -648,8 +703,10 @@ SUBROUTINE allocate_rivers_vars_rp( np_rivers, rivers, rivers_data ) ! Error flags. np_land_tmp, & ! Number of land points to allocate for. - np_rivers_tmp + np_rivers_tmp, & ! Number of river points to allocate for. + np_global_land_tmp + ! Number of global land points (across all tasks) to allocate for. !end of header !------------------------------------------------------------------------------ @@ -880,6 +937,62 @@ SUBROUTINE allocate_rivers_vars_rp( np_rivers, rivers, rivers_data ) ALLOCATE(rivers_data%rivers_outflow_number_rp(np_rivers_tmp), STAT = ERROR) error_sum = error_sum + ERROR +!------------------------------------------------------------------------------ +! Allocate reservoir variables. +!------------------------------------------------------------------------------ +IF ( l_reservoirs .AND. is_master_task() ) THEN + ! Full size. + np_rivers_tmp = np_rivers + np_global_land_tmp = global_land_pts +ELSE + ! Minimum size. + np_rivers_tmp = 1 + np_global_land_tmp = 1 +END IF + +! Reservoir ancillaries. +ALLOCATE( rivers_data%res_capacity(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_catch(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_year(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_critical(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_flood(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_emergency(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_normal_release(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_flood_release(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +! Reservoir prognostics. +ALLOCATE( rivers_data%res_storage(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +ALLOCATE( rivers_data%res_cap_current(np_rivers_tmp), STAT = ERROR ) +error_sum = error_sum + ERROR +! Reservoir coupling to rivers. +ALLOCATE( rivers_data%tot_abstracted_res_global(np_global_land_tmp), & + STAT = ERROR ) +error_sum = error_sum + ERROR + +!------------------------------------------------------------------------------ +! Allocate further variables for coupling to water resources. +!------------------------------------------------------------------------------ +IF ( l_water_resources .AND. is_master_task() ) THEN + ! Full size. + np_global_land_tmp = global_land_pts +ELSE + ! Minimum size. + np_global_land_tmp = 1 +END IF + +! Coupling abstractions from rivers. +ALLOCATE( rivers_data%tot_net_abstracted_river_global(np_global_land_tmp), & + STAT = ERROR ) +error_sum = error_sum + ERROR + IF ( error_sum /= 0 ) THEN CALL log_fatal( RoutineName, "Error allocating for routing point arrays." ) END IF @@ -956,6 +1069,23 @@ SUBROUTINE allocate_rivers_vars_rp( np_rivers, rivers, rivers_data ) rivers_data%sub_surf_roff_rp(:) = 0.0 rivers_data%surf_roff_rp(:) = 0.0 +! Initialise reservoir variables. Storage (the state variable) should be +! initialised to zero until better initialisation is possible. +rivers_data%res_capacity(:) = rmdi +rivers_data%res_catch(:) = rmdi +rivers_data%res_year(:) = rmdi +rivers_data%res_critical(:) = rmdi +rivers_data%res_flood(:) = rmdi +rivers_data%res_emergency(:) = rmdi +rivers_data%res_normal_release(:)= rmdi +rivers_data%res_flood_release(:) = rmdi +rivers_data%res_storage(:) = 0.0 +rivers_data%res_cap_current(:) = 0.0 + +! Initialise further variables to couple to water resources. +rivers_data%tot_abstracted_res_global(:) = rmdi +rivers_data%tot_net_abstracted_river_global(:) = rmdi + !------------------------------------------------------------------------------ ! Associate pointers !------------------------------------------------------------------------------ @@ -978,10 +1108,10 @@ SUBROUTINE allocate_rivers_vars_rp( np_rivers, rivers, rivers_data ) rivers%rrun_surf_rp => rivers_data%rrun_surf_rp rivers%inland_outflow_rp => rivers_data%inland_outflow_rp -! Associate pointers for regridding variable. +! Regridding variable. rivers%map_river_to_land_points => rivers_data%map_river_to_land_points -! Associate pointers for CaMa-Flood variables. +! CaMa-Flood variables. ! CaMa-Flood ancillary variables. rivers%channel_depth => rivers_data%channel_depth rivers%channel_width => rivers_data%channel_width @@ -1004,7 +1134,7 @@ SUBROUTINE allocate_rivers_vars_rp( np_rivers, rivers, rivers_data ) ! CaMa-Flood diagnostic variables. rivers%river_depth => rivers_data%river_depth -! Associate pointers for RFM variables. +! RFM variables. rivers%rfm_flowobs1_rp => rivers_data%rfm_flowobs1_rp rivers%rfm_iarea_rp => rivers_data%rfm_iarea_rp rivers%rfm_land_rp => rivers_data%rfm_land_rp @@ -1015,15 +1145,33 @@ SUBROUTINE allocate_rivers_vars_rp( np_rivers, rivers, rivers_data ) rivers%rfm_substore_rp => rivers_data%rfm_substore_rp rivers%rfm_surfstore_rp => rivers_data%rfm_surfstore_rp -! Associate pointers for TRIP variables. +! TRIP variables. rivers%rivers_seq_rp => rivers_data%rivers_seq_rp rivers%rivers_sto_rp => rivers_data%rivers_sto_rp -! Associate pointers for coupling variables. +! Coupling variables. rivers%rivers_outflow_number_rp => rivers_data%rivers_outflow_number_rp rivers%sub_surf_roff_rp => rivers_data%sub_surf_roff_rp rivers%surf_roff_rp => rivers_data%surf_roff_rp +! Further variables to couple to water resources. +rivers%tot_net_abstracted_river_global & + => rivers_data%tot_net_abstracted_river_global + +! Associate pointers for reservoir variables. +rivers%res_capacity => rivers_data%res_capacity +rivers%res_catch => rivers_data%res_catch +rivers%res_year => rivers_data%res_year +rivers%res_critical => rivers_data%res_critical +rivers%res_flood => rivers_data%res_flood +rivers%res_emergency => rivers_data%res_emergency +rivers%res_normal_release => rivers_data%res_normal_release +rivers%res_flood_release => rivers_data%res_flood_release +rivers%res_storage => rivers_data%res_storage +rivers%res_cap_current => rivers_data%res_cap_current +rivers%tot_abstracted_res_global & + => rivers_data%tot_abstracted_res_global + RETURN END SUBROUTINE allocate_rivers_vars_rp @@ -1037,15 +1185,18 @@ SUBROUTINE deallocate_river_props( ) !------------------------------------------------------------------------------ USE jules_rivers_mod, ONLY: & - channel_depth_grid, channel_width_grid, river_distance_grid, & - river_elevation_grid, river_length_grid, river_manning_grid, & - river_nextx_grid, river_nexty_grid + channel_depth_grid, channel_width_grid, & + river_distance_grid, river_elevation_grid, & + river_length_grid, river_manning_grid, river_nextx_grid, river_nexty_grid, & + res_capacity_grid, res_catch_grid, res_year_grid, & + res_critical_grid, res_flood_grid, res_emergency_grid, & + res_normal_release_grid, res_flood_release_grid IMPLICIT NONE !------------------------------------------------------------------------------ -! Variables that are only allocated on master task and in some configurations -! - hence first check if allocated. +! Variables that are only allocated in some configurations and/or on master +! task - hence first check if allocated. !------------------------------------------------------------------------------ IF ( ALLOCATED(direction_grid) ) THEN DEALLOCATE( direction_grid ) @@ -1075,7 +1226,7 @@ SUBROUTINE deallocate_river_props( ) DEALLOCATE( river_nexty_grid ) END IF -! CaMa-Flood variables that are only allocated on master task. +! CaMa-Flood variables. IF ( ALLOCATED( channel_depth_grid ) ) THEN DEALLOCATE( channel_depth_grid ) END IF @@ -1100,10 +1251,36 @@ SUBROUTINE deallocate_river_props( ) DEALLOCATE( river_manning_grid ) END IF +! Reservoir variables. +IF ( ALLOCATED( res_capacity_grid ) ) THEN + DEALLOCATE( res_capacity_grid ) +END IF +IF ( ALLOCATED( res_catch_grid ) ) THEN + DEALLOCATE( res_catch_grid ) +END IF +IF ( ALLOCATED( res_year_grid ) ) THEN + DEALLOCATE( res_year_grid ) +END IF +IF ( ALLOCATED( res_critical_grid ) ) THEN + DEALLOCATE( res_critical_grid ) +END IF +IF ( ALLOCATED( res_flood_grid ) ) THEN + DEALLOCATE( res_flood_grid ) +END IF +IF ( ALLOCATED( res_emergency_grid ) ) THEN + DEALLOCATE( res_emergency_grid ) +END IF +IF ( ALLOCATED( res_normal_release_grid ) ) THEN + DEALLOCATE( res_normal_release_grid ) +END IF +IF ( ALLOCATED( res_flood_release_grid ) ) THEN + DEALLOCATE( res_flood_release_grid ) +END IF + RETURN END SUBROUTINE deallocate_river_props !############################################################################## END MODULE jules_rivers_props_mod -#endif +#endif \ No newline at end of file diff --git a/src/initialisation/standalone/init_output_mod.F90 b/src/initialisation/standalone/init_output_mod.F90 index f35234c0..3d6d002c 100644 --- a/src/initialisation/standalone/init_output_mod.F90 +++ b/src/initialisation/standalone/init_output_mod.F90 @@ -344,10 +344,13 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & USE jules_water_resources_mod, ONLY: l_have_groundwater, l_have_renew_gwater, & l_have_surface_water, l_water_domestic, l_water_environment, & l_water_industry, l_water_irrigation, l_water_livestock, & - l_water_resources, l_water_transfers, no_model, nr_gwater_model + l_water_resources, l_water_transfers, no_model, nr_gwater_model, & + sw_river_source -USE jules_rivers_mod, ONLY: l_rivers, l_riv_overbank, l_outflow_per_river, & - i_river_vn, rivers_camaflood, rivers_rfm, rivers_trip, l_inland_outflow +USE jules_rivers_mod, ONLY: & + l_outflow_per_river, l_rivers, l_riv_overbank, i_river_vn, & + rivers_camaflood, rivers_rfm, rivers_trip, l_inland_outflow, & + l_reservoirs USE jules_deposition_mod, ONLY: l_deposition, l_deposition_flux @@ -664,10 +667,11 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & ELSE !------------------------------------------------------------------------- - ! Variables that are only allowed with river routing. + ! Rivers are not selected. + ! Check for variables that are only allowed with river routing. ! First we check for variables on the river output grid (which allow us to ! shortern the list by omitting those), then by name for any other river - ! variable.< + ! variable. !------------------------------------------------------------------------- IF ( is_river_var(var(j)) ) THEN remove_var = .TRUE. @@ -682,15 +686,19 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & END IF ! l_rivers - !--------------------------------------------------------------------------- + !---------------------------------------------------------------------------- ! River 2D ancillary variables (including for overbank inundation) - these ! can be read in but not output (because code for 2D river grid outputs does ! not exist). - !--------------------------------------------------------------------------- + !---------------------------------------------------------------------------- SELECT CASE ( var(j) ) CASE ( 'area', 'channel_depth_grid', 'channel_width_grid', & 'direction', 'hypso_quantiles_grid', 'latitude_2d', 'logn_mean', & 'logn_stdev', 'longitude_2d', 'mean_sea_level_grid', & + 'res_capacity_grid', 'res_catch_grid', & + 'res_year_grid', 'res_critical_grid', & + 'res_flood_grid', 'res_emergency_grid', & + 'res_normal_release_grid', 'res_flood_release_grid', & 'river_distance_grid', 'river_elevation_grid', 'river_length_grid', & 'river_manning_grid', 'river_nextx_grid', 'river_nexty_grid', & 'rivers_xgrid', 'rivers_ygrid', 'sequence' ) @@ -716,7 +724,8 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & ! Water resource variables that only require l_water_resources=T. IF ( .NOT. l_water_resources ) THEN SELECT CASE ( var(j) ) - CASE ( 'conv_loss_frac', 'water_demand', 'water_demand_unmet' ) + CASE ( 'conv_loss_frac', 'conveyance_loss', 'water_demand', & + 'water_demand_unmet', 'water_removed' ) remove_var = .TRUE. message = 'Water resources (l_water_resources) not selected.' END SELECT @@ -724,11 +733,11 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & ! Water resource variables that also require another switch. - ! These sector-specific switches are set to F if l_water_resources=F, - ! so we don't need to test l_water_resources here. + ! Sector-specific switches (such as l_water_domestic) are FALSE if + ! l_water_resources=F, so we don't need to test l_water_resources here. IF ( .NOT. l_water_domestic ) THEN SELECT CASE ( var(j) ) - CASE ( 'demand_domestic', 'demand_rate_domestic', 'domestic_unmet' ) + CASE ( 'demand_domestic', 'demand_rate_domestic', 'unmet_domestic' ) remove_var = .TRUE. message = 'l_water_domestic not selected.' END SELECT @@ -736,7 +745,7 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & IF ( .NOT. l_water_environment ) THEN SELECT CASE ( var(j) ) - CASE ( 'demand_environment', 'environment_unmet' ) + CASE ( 'demand_environment', 'unmet_environment' ) remove_var = .TRUE. message = 'l_water_environment not selected.' END SELECT @@ -744,7 +753,7 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & IF ( .NOT. l_water_industry ) THEN SELECT CASE ( var(j) ) - CASE ( 'demand_industry', 'demand_rate_industry', 'industry_unmet' ) + CASE ( 'demand_industry', 'demand_rate_industry', 'unmet_industry' ) remove_var = .TRUE. message = 'l_water_industry not selected.' END SELECT @@ -752,7 +761,7 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & IF ( .NOT. l_water_irrigation ) THEN SELECT CASE ( var(j) ) - CASE ( 'demand_irrigation', 'grid_area', 'irrigation_unmet' ) + CASE ( 'demand_irrigation', 'grid_area', 'unmet_irrigation' ) remove_var = .TRUE. message = 'l_water_irrigation not selected.' END SELECT @@ -760,7 +769,7 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & IF ( .NOT. l_water_livestock ) THEN SELECT CASE ( var(j) ) - CASE ( 'demand_livestock', 'demand_rate_livestock','livestock_unmet' ) + CASE ( 'demand_livestock', 'demand_rate_livestock','unmet_livestock' ) remove_var = .TRUE. message = 'l_water_livestock not selected.' END SELECT @@ -768,7 +777,7 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & IF ( .NOT. l_water_transfers ) THEN SELECT CASE ( var(j) ) - CASE ( 'demand_transfers', 'transfers_unmet' ) + CASE ( 'demand_transfers', 'unmet_transfers' ) remove_var = .TRUE. message = 'l_water_transfers not selected.' END SELECT @@ -782,7 +791,8 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & END SELECT END IF - IF ( nr_gwater_model == no_model ) THEN + IF ( .NOT. l_water_resources .OR. & + ( l_water_resources .AND. nr_gwater_model == no_model ) ) THEN SELECT CASE ( var(j) ) CASE ( 'gw_nr_abstracted' ) remove_var = .TRUE. @@ -806,6 +816,36 @@ SUBROUTINE check_output_vars( nvars, var, var_name, output_type, & END SELECT END IF + IF ( .NOT. l_water_resources .OR. .NOT. l_reservoirs ) THEN + SELECT CASE ( var(j) ) + CASE ( 'res_abstracted' ) + remove_var = .TRUE. + message = 'Water resources + reservoirs not used.' + END SELECT + END IF + + IF ( .NOT. l_water_resources .OR. sw_river_source == 0 ) THEN + ! Rivers aren't being used as a source of water for abstraction. + SELECT CASE ( var(j) ) + CASE ( 'river_abstracted' ) + remove_var = .TRUE. + message = 'Rivers not used for abstractions.' + END SELECT + END IF + + !---------------------------------------------------------------------------- + ! Variables that are only allowed with reservoirs. + !---------------------------------------------------------------------------- + IF ( .NOT. l_reservoirs ) THEN + SELECT CASE ( var(j) ) + CASE ( 'res_capacity', 'res_catch', 'res_storage', 'res_year', & + 'res_critical', 'res_flood', 'res_emergency', & + 'res_normal_release', 'res_flood_release' ) + remove_var = .TRUE. + message = 'Reservoirs not used.' + END SELECT + END IF ! l_reservoirs + !--------------------------------------------------------------------------- ! Variables that are only allowed with atmospheric deposition !--------------------------------------------------------------------------- diff --git a/src/initialisation/standalone/initial_conditions/init_ic_mod.F90 b/src/initialisation/standalone/initial_conditions/init_ic_mod.F90 index 15fc5381..39389d92 100644 --- a/src/initialisation/standalone/initial_conditions/init_ic_mod.F90 +++ b/src/initialisation/standalone/initial_conditions/init_ic_mod.F90 @@ -993,6 +993,10 @@ FUNCTION get_default_ic_values(total_snow) RESULT(defaults_dict) CASE ( 'lake_h_ice_gb' ) CALL dict_set(defaults_dict, 'lake_h_ice_gb' , lake_H_ice_0) + !Reservoirs + CASE ( 'res_storage' ) + CALL dict_set(defaults_dict, 'res_storage', 0.0) + !Snow variables ! CASE ( 'rho_snow' ) ! CASE ( 'snow_depth' ) diff --git a/src/io/dump/get_dim_info_mod.F90 b/src/io/dump/get_dim_info_mod.F90 index be189123..ecd42154 100644 --- a/src/io/dump/get_dim_info_mod.F90 +++ b/src/io/dump/get_dim_info_mod.F90 @@ -564,8 +564,7 @@ SUBROUTINE get_dim_info_rivers( l_reading, identifier, ndims, dim_sizes, & 'rivers_outflow_rp', 'outflow_per_river', 'inland_outflow_rp', & 'flood_flow', 'flood_flow_prev', 'flood_storage', 'flood_storage_prev', & 'river_channel_flow', 'river_channel_storage', 'river_depth_prev', & - 'river_flow_prev' & - ) + 'river_flow_prev','res_storage') ndims = 1 dim_names(1) = p_rivers_dim_name dim_sizes(1) = np_rivers diff --git a/src/io/dump/read_dump_mod.F90 b/src/io/dump/read_dump_mod.F90 index fcfd4493..7a216ee2 100644 --- a/src/io/dump/read_dump_mod.F90 +++ b/src/io/dump/read_dump_mod.F90 @@ -987,7 +987,8 @@ SUBROUTINE read_dump(file_name, identifiers) 'flood_storage_prev', 'rfm_bflowin_rp', 'rfm_flowin_rp', & 'rfm_substore_rp', 'rfm_surfstore_rp', 'river_channel_storage', & 'river_channel_flow', 'river_depth_prev', 'river_flow_prev', & - 'rivers_outflow_rp', 'rivers_sto_rp', 'inland_outflow_rp' ) + 'rivers_outflow_rp', 'rivers_sto_rp', 'inland_outflow_rp', & + 'res_storage' ) ! nothing to do !------------------------------------------------------------------------- diff --git a/src/io/dump/read_dump_var_rivers_mod.F90 b/src/io/dump/read_dump_var_rivers_mod.F90 index 2045b2b8..41fc4845 100644 --- a/src/io/dump/read_dump_var_rivers_mod.F90 +++ b/src/io/dump/read_dump_var_rivers_mod.F90 @@ -100,6 +100,9 @@ SUBROUTINE read_dump_var_rivers( identifier, FILE, var_id ) CASE ( 'inland_outflow_rp' ) CALL file_read_var(FILE, var_id, rivers%inland_outflow_rp) +CASE ( 'res_storage' ) + CALL file_read_var(FILE, var_id, rivers%res_storage) + CASE DEFAULT CALL log_fatal( RoutineName, & "Unexpected variable in dump - " // TRIM(identifier) ) diff --git a/src/io/dump/required_vars_for_rivers_mod.F90 b/src/io/dump/required_vars_for_rivers_mod.F90 index 8bbbe4a4..e88d9461 100644 --- a/src/io/dump/required_vars_for_rivers_mod.F90 +++ b/src/io/dump/required_vars_for_rivers_mod.F90 @@ -17,7 +17,7 @@ SUBROUTINE required_vars_for_rivers( nvars, identifiers, & USE jules_rivers_mod, ONLY: l_outflow_per_river, i_river_vn, rivers_camaflood, & rivers_rfm, rivers_trip, l_init_storage, & - l_inland_outflow + l_inland_outflow, l_reservoirs USE logging_mod, ONLY: log_warn @@ -139,6 +139,15 @@ SUBROUTINE required_vars_for_rivers( nvars, identifiers, & END SELECT +IF ( l_reservoirs ) THEN + IF ( read_or_write_dump_local ) THEN + CALL add_to_list( 'res_storage', nvars, identifiers ) + ELSE + CALL log_warn( RoutineName, & + "Storage in reservoirs will be initialised to zero.") + END IF +END IF + RETURN END SUBROUTINE required_vars_for_rivers diff --git a/src/io/dump/write_dump_mod.F90 b/src/io/dump/write_dump_mod.F90 index 55a56520..db69c742 100644 --- a/src/io/dump/write_dump_mod.F90 +++ b/src/io/dump/write_dump_mod.F90 @@ -676,14 +676,14 @@ SUBROUTINE write_dump() 'seed_rain', 'ch4_ppbv' ) ! Nothing to do - CASE ( 'rivers_lat_rp', 'rivers_lon_rp', 'rivers_sto_rp', & + CASE ( 'rivers_lat_rp', 'rivers_lon_rp', & 'rivers_x_coord_rp', 'rivers_y_coord_rp', & - 'rfm_surfstore_rp', 'rfm_substore_rp', & + 'rivers_sto_rp', 'rfm_surfstore_rp', 'rfm_substore_rp', & 'rfm_flowin_rp', 'rfm_bflowin_rp', 'rivers_outflow_rp', & 'inland_outflow_rp', & 'flood_flow', 'flood_flow_prev', 'flood_storage', & 'flood_storage_prev', 'river_channel_flow', 'river_channel_storage', & - 'river_depth_prev', 'river_flow_prev' ) + 'river_depth_prev', 'river_flow_prev', 'res_storage' ) ! Rivers are run on a single task, hence nothing to do here. !------------------------------------------------------------------------- diff --git a/src/io/dump/write_dump_var_rivers_mod.F90 b/src/io/dump/write_dump_var_rivers_mod.F90 index 20475240..ae2267c1 100644 --- a/src/io/dump/write_dump_var_rivers_mod.F90 +++ b/src/io/dump/write_dump_var_rivers_mod.F90 @@ -113,6 +113,9 @@ SUBROUTINE write_dump_var_rivers( identifier, FILE, var_id ) CASE ( 'inland_outflow_rp' ) CALL file_write_var(FILE, var_id, rivers%inland_outflow_rp) +CASE ( 'res_storage' ) + CALL file_write_var(FILE, var_id, rivers%res_storage) + CASE DEFAULT CALL log_fatal( RoutineName, & "Unrecognised variable for dump - " // TRIM(identifier) ) diff --git a/src/io/model_interface/extract_var.inc b/src/io/model_interface/extract_var.inc index 83077bc9..ef4ac36b 100644 --- a/src/io/model_interface/extract_var.inc +++ b/src/io/model_interface/extract_var.inc @@ -2730,6 +2730,33 @@ CASE ( 'frac_fplain_lp' ) cube_land = cube_from_array(frac_fplain_lp) cube = map_from_land(cube_land) +CASE ( 'res_capacity' ) + cube = cube_from_array(rivers%res_capacity) + +CASE ( 'res_catch' ) + cube = cube_from_array(rivers%res_catch) + +CASE ( 'res_year' ) + cube = cube_from_array(rivers%res_year) + +CASE ( 'res_critical' ) + cube = cube_from_array(rivers%res_critical) + +CASE ( 'res_flood' ) + cube = cube_from_array(rivers%res_flood) + +CASE ( 'res_emergency' ) + cube = cube_from_array(rivers%res_emergency) + +CASE ( 'res_normal_release' ) + cube = cube_from_array(rivers%res_normal_release) + +CASE ( 'res_flood_release' ) + cube = cube_from_array(rivers%res_flood_release) + +CASE ( 'res_storage' ) + cube = cube_from_array(rivers%res_storage) + CASE ( 'lake_depth' ) cube_land = cube_from_array(lake_vars%lake_depth_gb) cube = map_from_land(cube_land) @@ -3318,6 +3345,11 @@ CASE ( 'conv_loss_frac' ) cube_land = cube_from_array( water_resources%conv_loss_frac ) cube = map_from_land(cube_land) +CASE ( 'conveyance_loss' ) + ! Convert kg to kg s-1. + cube_land = cube_from_array(water_resources%conveyance_loss / timestep) + cube = map_from_land(cube_land) + CASE ( 'demand_rate_domestic' ) cube_land = cube_from_array(water_resources%demand_rate_domestic) cube = map_from_land(cube_land) @@ -3420,6 +3452,11 @@ CASE ( 'gw_nr_abstracted' ) cube_land = cube_from_array( water_resources%gw_nr_abstracted / timestep ) cube = map_from_land(cube_land) +CASE ( 'river_abstracted' ) + ! Convert kg to kg s-1. + cube_land = cube_from_array( water_resources%abstracted_river / timestep ) + cube = map_from_land(cube_land) + CASE ( 'sfc_water_frac' ) cube_land = cube_from_array( water_resources%sfc_water_frac ) cube = map_from_land(cube_land) @@ -3491,6 +3528,11 @@ CASE ( 'water_removed' ) cube_land = cube_from_array( water_resources%water_removed / timestep ) cube = map_from_land(cube_land) +CASE ( 'res_abstracted' ) + ! Convert kg to kg s-1. + cube_land = cube_from_array( water_resources%abstracted_res / timestep ) + cube = map_from_land(cube_land) + !--------------------------------------------------------------------------- ! Urban morphology on land points ! These are calculated when l_urban_empirical otherwise they will be the diff --git a/src/io/model_interface/model_interface_mod.F90 b/src/io/model_interface/model_interface_mod.F90 index b06b80a8..1d0c867f 100644 --- a/src/io/model_interface/model_interface_mod.F90 +++ b/src/io/model_interface/model_interface_mod.F90 @@ -176,7 +176,9 @@ MODULE model_interface_mod ! Array holding the metadata for all model variables that we can use for input ! or output. The CABLE land surface model adds 10 prognostics for tiled ! soil/snow prognostics. -INTEGER, PARAMETER :: n_vars = 723 + +INTEGER, PARAMETER :: n_vars = 743 + TYPE(var_metadata) :: metadata(n_vars) ! Include the metadata DATA statement diff --git a/src/io/model_interface/populate_var.inc b/src/io/model_interface/populate_var.inc index 5b1ca8af..1f715000 100644 --- a/src/io/model_interface/populate_var.inc +++ b/src/io/model_interface/populate_var.inc @@ -33,8 +33,12 @@ USE prognostics, ONLY: & USE jules_rivers_mod, ONLY: & ! 2D ancillary fields used during initialisation. channel_depth_grid, channel_width_grid, mean_sea_level_grid, & - river_distance_grid, river_elevation_grid, river_length_grid, & - river_manning_grid, river_nextx_grid, river_nexty_grid + river_distance_grid, & + river_elevation_grid, river_length_grid, river_manning_grid, & + river_nextx_grid, river_nexty_grid, & + res_capacity_grid, res_catch_grid, res_year_grid, & + res_critical_grid, res_flood_grid, res_emergency_grid, & + res_normal_release_grid, res_flood_release_grid USE jules_soil_mod, ONLY: & sm_levels, l_tile_soil, l_broadcast_ancils @@ -2487,6 +2491,14 @@ CASE ( 'sea_level_lp' ) rivers%sea_level_lp(:) = const_val END IF +CASE ( 'res_storage' ) + IF ( PRESENT(cube) ) THEN + cube_land = cube_create([ SIZE(cube%values) ]) + CALL cube_get_data(cube_land, rivers%res_storage) + ELSE + rivers%res_storage(:) = const_val + END IF + !----------------------------------------------------------------------------- ! Overbank inundation variables !----------------------------------------------------------------------------- @@ -2512,6 +2524,62 @@ CASE ( 'logn_stdev' ) logn_stdev(:,:) = const_val END IF +CASE ( 'res_capacity_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_capacity_grid) + ELSE + res_capacity_grid(:,:) = const_val + END IF + +CASE ( 'res_catch_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_catch_grid) + ELSE + res_catch_grid(:,:) = const_val + END IF + +CASE ( 'res_year_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_year_grid) + ELSE + res_year_grid(:,:) = const_val + END IF + +CASE ( 'res_critical_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_critical_grid) + ELSE + res_critical_grid(:,:) = const_val + END IF + +CASE ( 'res_flood_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_flood_grid) + ELSE + res_flood_grid(:,:) = const_val + END IF + +CASE ( 'res_emergency_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_emergency_grid) + ELSE + res_emergency_grid(:,:) = const_val + END IF + +CASE ( 'res_normal_release_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_normal_release_grid) + ELSE + res_normal_release_grid(:,:) = const_val + END IF + +CASE ( 'res_flood_release_grid' ) + IF ( PRESENT(cube) ) THEN + CALL cube_get_data(cube, res_flood_release_grid) + ELSE + res_flood_release_grid(:,:) = const_val + END IF + !----------------------------------------------------------------------------- ! Fire and metstats variables- land points only !----------------------------------------------------------------------------- diff --git a/src/io/model_interface/variable_metadata.inc b/src/io/model_interface/variable_metadata.inc index 8b7fa44c..956d2a00 100644 --- a/src/io/model_interface/variable_metadata.inc +++ b/src/io/model_interface/variable_metadata.inc @@ -7951,7 +7951,7 @@ DATA metadata(614) / var_metadata( & "kg m-2 s-1" & ) / !----------------------------------------------------------------------------- -! Metadata for grid_area +! Metadata for grid_area_rp !----------------------------------------------------------------------------- DATA metadata(615) / var_metadata( & ! String identifier @@ -9367,4 +9367,256 @@ DATA metadata(723) / var_metadata( & ! Units "kg m-2 s-1" & ) / +!----------------------------------------------------------------------------- +! Metadata for conveyance_loss +!----------------------------------------------------------------------------- +DATA metadata(724) / var_metadata( & +! String identifier + 'conveyance_loss', & +! Variable type + var_type_surface, & +! Long name + "Water that is lost during conveyance", & +! Units + "kg s-1" & + ) / +!----------------------------------------------------------------------------- +! Metadata for river_abstracted +!----------------------------------------------------------------------------- +DATA metadata(725) / var_metadata( & +! String identifier + 'river_abstracted', & +! Variable type + var_type_surface, & +! Long name + "Water abstracted from rivers", & +! Units + "kg s-1" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_capacity_grid +!----------------------------------------------------------------------------- +DATA metadata(726) / var_metadata( & +! String identifier + 'res_capacity_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_catch_grid +!----------------------------------------------------------------------------- +DATA metadata(727) / var_metadata( & +! String identifier + 'res_catch_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_year_grid +!----------------------------------------------------------------------------- +DATA metadata(728) / var_metadata( & +! String identifier + 'res_year_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_critical_grid +!----------------------------------------------------------------------------- +DATA metadata(729) / var_metadata( & +! String identifier + 'res_critical_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_flood_grid +!----------------------------------------------------------------------------- +DATA metadata(730) / var_metadata( & +! String identifier + 'res_flood_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_emergency_grid +!----------------------------------------------------------------------------- +DATA metadata(731) / var_metadata( & +! String identifier + 'res_emergency_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_normal_release_grid +!----------------------------------------------------------------------------- +DATA metadata(732) / var_metadata( & +! String identifier + 'res_normal_release_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_flood_release_grid +!----------------------------------------------------------------------------- +DATA metadata(733) / var_metadata( & +! String identifier + 'res_flood_release_grid', & +! Variable type + var_type_surface, & +! Variable is not available for output, so give dummy values for long name and +! units + "", "" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_capacity +!----------------------------------------------------------------------------- +DATA metadata(734) / var_metadata( & +! String identifier + 'res_capacity', & +! Variable type + var_type_rp, & +! Long name + "Storage capacity of reservoirs", & +! Units + "kg" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_catch +!----------------------------------------------------------------------------- +DATA metadata(735) / var_metadata( & +! String identifier + 'res_catch', & +! Variable type + var_type_rp, & +! Long name + "Catchment area of reservoirs", & +! Units + "m2" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_year +!----------------------------------------------------------------------------- +DATA metadata(736) / var_metadata( & +! String identifier + 'res_year', & +! Variable type + var_type_rp, & +! Long name + "Initialisation year of reservoirs", & +! Units + "1" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_critical +!----------------------------------------------------------------------------- +DATA metadata(737) / var_metadata( & +! String identifier + 'res_critical', & +! Variable type + var_type_rp, & +! Long name + "Critical storage of reservoirs", & +! Units + "kg" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_flood +!----------------------------------------------------------------------------- +DATA metadata(738) / var_metadata( & +! String identifier + 'res_flood', & +! Variable type + var_type_rp, & +! Long name + "Flood storage of reservoirs", & +! Units + "kg" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_emergency +!----------------------------------------------------------------------------- +DATA metadata(739) / var_metadata( & +! String identifier + 'res_emergency', & +! Variable type + var_type_rp, & +! Long name + "Emergency storage of reservoirs", & +! Units + "kg" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_normal_release +!----------------------------------------------------------------------------- +DATA metadata(740) / var_metadata( & +! String identifier + 'res_normal_release', & +! Variable type + var_type_rp, & +! Long name + "Normal release rate of reservoirs", & +! Units + "kg s-1" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_flood_release +!----------------------------------------------------------------------------- +DATA metadata(741) / var_metadata( & +! String identifier + 'res_flood_release', & +! Variable type + var_type_rp, & +! Long name + "Flood release rate of reservoirs", & +! Units + "kg s-1" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_storage +!----------------------------------------------------------------------------- +DATA metadata(742) / var_metadata( & +! String identifier + 'res_storage', & +! Variable type + var_type_rp, & +! Long name + "Water stored in reservoirs", & +! Units + "kg" & + ) / +!----------------------------------------------------------------------------- +! Metadata for res_abstracted +!----------------------------------------------------------------------------- +DATA metadata(743) / var_metadata( & +! String identifier + 'res_abstracted', & +! Variable type + var_type_surface, & +! Long name + "Water abstracted from reservoirs", & +! Units + "kg s-1" & + ) / #endif From 839e74c9ef883874a61a805e3b1209bf703a981a Mon Sep 17 00:00:00 2001 From: Helen Baron Date: Wed, 26 Aug 2026 14:22:42 +0100 Subject: [PATCH 2/4] Adjustments to water resource code --- doc/source/output-variables.rst | 41 +- .../shared/jules_water_resources_mod.F90 | 32 +- .../shared/water_resources_control_mod.F90 | 598 ++++++++++++++---- .../shared/water_resources_vars_mod.F90 | 213 +++++-- .../water_resources/calc_avail_water_mod.F90 | 47 +- .../water_resources/water_resources_drive.F90 | 176 +++--- 6 files changed, 772 insertions(+), 335 deletions(-) diff --git a/doc/source/output-variables.rst b/doc/source/output-variables.rst index 78bb6822..77910748 100644 --- a/doc/source/output-variables.rst +++ b/doc/source/output-variables.rst @@ -1337,10 +1337,7 @@ These variables are only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_r +---------------------------+--------------------------------------------------------------------------------+------------+ | Name | Description | Dimensions | +===========================+================================================================================+============+ -| ``conv_loss_frac`` | Fraction of abstracted water that is lost during conveyance from source to | | -| | user. | | -+---------------------------+--------------------------------------------------------------------------------+------------+ -| ``sfc_water_frac`` | Target for the fraction of demand to be met from surface water. | | +| Variables that quantify the demand for water: | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``demand_rate_domestic`` | Demand rate for water for domestic use (kg s\ :sup:`-1`). | | | | Only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_domestic` = TRUE. | | @@ -1354,9 +1351,6 @@ These variables are only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_r | ``demand_rate_transfers`` | Demand rate for water for transfers (kg s\ :sup:`-1`). | | | | Only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_transfers` = TRUE. | | +---------------------------+--------------------------------------------------------------------------------+------------+ -| ``water_demand`` | Demand for water across all water resource sectors (kg s\ :sup:`-1`), including| | -| | any allowance for conveyance loss. | | -+---------------------------+--------------------------------------------------------------------------------+------------+ | ``demand_domestic`` | Demand for water for domestic use (kg s\ :sup:`-1`), including any allowance | | | | for conveyance loss. | | | | Only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_domestic` = TRUE. | | @@ -1381,7 +1375,10 @@ These variables are only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_r | | conveyance loss. | | | | Only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_transfers` = TRUE. | | +---------------------------+--------------------------------------------------------------------------------+------------+ -| ``water_demand_unmet`` | The part of the total demand for water that is not satisfied (kg s\ :sup:`-1`).| | +| ``water_demand`` | Demand for water across all water resource sectors (kg s\ :sup:`-1`), including| | +| | any allowance for conveyance loss. | | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| Variables that quantify how much of the demand for water cannot be met: | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``unmet_domestic`` | The part of the demand for water for domestic use that is not satisfied | | | | (kg s\ :sup:`-1`). | | @@ -1407,18 +1404,42 @@ These variables are only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_r | | (kg s\ :sup:`-1`). | | | | Only available if :nml:mem:`JULES_WATER_RESOURCES::l_water_transfers` = TRUE. | | +---------------------------+--------------------------------------------------------------------------------+------------+ -| ``irrig_water`` | Water applied as irrigation (kg m\ :sup:`-2` s\ :sup:`-1`). | | -| | Only available if :nml:mem:`JULES_IRRIG::l_irrig_dmd` = TRUE. | | +| ``water_demand_unmet`` | The part of the total demand for water that is not satisfied (kg s\ :sup:`-1`).| | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| Variables that quantify how much water is available and abstracted: | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``gw_avail`` | Groundwater that is available for abstraction at start of timestep (kg) . | | +| | Only available if groundwater sources are modelled. | | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``sw_avail`` | Surface water that is available for abstraction at start of timestep (kg). | | +| | Only available if surface water sources are modelled. | | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``gw_abstracted`` | Water abstracted from renewable groundwater sources (kg s\ :sup:`-1`). | | +| | Only available if a suitable groundwater model is used. | | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``gw_nr_abstracted`` | Water abstracted from non-renewable groundwater sources (kg s\ :sup:`-1`). | | +| | Only available if a suitable groundwater model is used. | | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| ``res_abstracted`` | Water abstracted from reservoirs (kg s\ :sup:`-1`). | | +| | Only available if :nml:mem:`JULES_RIVERS::l_reservoirs` = TRUE. | | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| ``river_abstracted`` | Water abstracted from rivers (kg s\ :sup:`-1`). | | +| | Only available if :nml:mem:`JULES_RIVERS::l_rivers` = TRUE. | | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``sw_abstracted`` | Water abstracted from surface water sources (kg s\ :sup:`-1`). | | +| | Only available if surface water sources are modelled. | | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| ``irrig_water`` | Water applied as irrigation (kg m\ :sup:`-2` s\ :sup:`-1`). | | +| | Only available if :nml:mem:`JULES_IRRIG::l_irrig_dmd` = TRUE. | | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| Miscellaneous water resource variables: | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| ``conv_loss_frac`` | Fraction of abstracted water that is lost during conveyance from source to | | +| | user (an input to the model). | | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| ``conveyance_loss`` | Water that is lost during conveyance (kg s\ :sup:`-1`). | | ++---------------------------+--------------------------------------------------------------------------------+------------+ +| ``sfc_water_frac`` | Target for the fraction of demand to be met from surface water. | | +---------------------------+--------------------------------------------------------------------------------+------------+ | ``water_removed`` | Water that is removed from the system during use (kg s\ :sup:`-1`). | | +---------------------------+--------------------------------------------------------------------------------+------------+ diff --git a/src/control/shared/jules_water_resources_mod.F90 b/src/control/shared/jules_water_resources_mod.F90 index 2b44c3bb..3507e13e 100644 --- a/src/control/shared/jules_water_resources_mod.F90 +++ b/src/control/shared/jules_water_resources_mod.F90 @@ -200,17 +200,19 @@ MODULE jules_water_resources_mod ! Index of livestock use in multi-use arrays. use_transfers = 0, & ! Index of transfers in multi-use arrays. + sw_res_source = 0, & + ! Index of reservoirs in surface water source arrays. water_res_count = 0 ! Counter of timesteps done in current water resource timestep. LOGICAL :: & l_have_groundwater = .FALSE., & - ! Flag indicating if we have a model of groundwater (renewable or - ! non-renewable). + ! Flag indicating if groundwater sources (renewable or non-renewable) are + ! represented. l_have_renew_gwater = .FALSE., & ! Flag indicating if we have a model of renewable groundwater. l_have_surface_water = .FALSE. - ! Flag indicating if we have surface water represented (e.g. rivers). + ! Flag indicating if surface water sources (e.g. rivers) are represented. ! TRUE means n_sw_source > 0. CONTAINS @@ -224,6 +226,7 @@ SUBROUTINE check_jules_water_resources( l_top ) USE jules_print_mgr, ONLY: jules_message +! This dependency on jules_rivers_mod is not ideal. USE jules_rivers_mod, ONLY: l_rivers !------------------------------------------------------------------------------ @@ -548,16 +551,6 @@ SUBROUTINE check_jules_water_resources( l_top ) END IF ! l_have_surface_water .AND. l_have_groundwater - ! Check that any surface water sources include rivers - so that we know we - ! can add return flow to rivers. At present this is guaranteed because rivers - ! are the only possible surface water source that has been coded. - ! The dependency on jules_rivers_mod here is not ideal. - IF ( l_have_surface_water .AND. .NOT. l_rivers ) THEN - error_status = 101 ! a fatal error - CALL ereport ( RoutineName, error_status, & - "Rivers must be included in any surface water sources." ) - END IF - END IF ! l_water_resources END SUBROUTINE check_jules_water_resources @@ -572,7 +565,8 @@ SUBROUTINE set_jules_water_resources( l_top ) ! Sets values related to water resource code. !------------------------------------------------------------------------------ -USE jules_rivers_mod, ONLY: l_rivers +! This dependency on jules_rivers_mod is not ideal. +USE jules_rivers_mod, ONLY: l_reservoirs, l_rivers IMPLICIT NONE @@ -635,16 +629,18 @@ SUBROUTINE set_jules_water_resources( l_top ) !---------------------------------------------------------------------------- ! Set index for each available surface water source to show position in - ! surface water arrays. - ! At present there is only one possible source, but others will follow. + ! surface water arrays. If a source exists, it is used. !---------------------------------------------------------------------------- ! Initialise as no surface water sources. n_sw_source = 0 IF ( l_rivers ) THEN - ! If rivers are modelled, they are always included as a source of water. n_sw_source = n_sw_source + 1 sw_river_source = n_sw_source END IF + IF ( l_reservoirs ) THEN + n_sw_source = n_sw_source + 1 + sw_res_source = n_sw_source + END IF ! Set surface water flag if we have any surface water sources. IF ( n_sw_source > 0 ) THEN @@ -857,8 +853,8 @@ SUBROUTINE read_nml_jules_water_resources (unitnumber) nr_gwater_model = my_nml % nr_gwater_model nstep_water_res = my_nml % nstep_water_res partition_method = my_nml % partition_method - l_water_domestic = my_nml % l_water_domestic l_prioritise = my_nml % l_prioritise + l_water_domestic = my_nml % l_water_domestic l_water_environment = my_nml % l_water_environment l_water_industry = my_nml % l_water_industry l_water_irrigation = my_nml % l_water_irrigation diff --git a/src/control/shared/water_resources_control_mod.F90 b/src/control/shared/water_resources_control_mod.F90 index 335194d8..6627a21d 100644 --- a/src/control/shared/water_resources_control_mod.F90 +++ b/src/control/shared/water_resources_control_mod.F90 @@ -63,6 +63,8 @@ MODULE water_resources_control_mod ! Demands for water accumulated over the water resource timestep (kg). demand_unmet_global(:,:), & ! The part of the demand for water that is not satisfied (kg). + grid_area_global(:), & + ! Gridbox area (m2). gw_abstracted_global(:), & ! Water abstracted from renewable groundwater (kg). gw_avail_global(:), & @@ -70,10 +72,14 @@ MODULE water_resources_control_mod ! This does not include "non-renewable" groundwater. gw_nr_abstracted_global(:), & ! Water abstracted from non-renewable groundwater (kg). + res_storage_global(:), & + ! Water stored in reservoirs, on land points (kg). return_flow_gw_global(:), & ! Water that is returned to renewable groundwater after use (kg). - return_flow_sw_global(:), & + return_flow_river_global(:), & ! Water that is returned to rivers after use (kg). + river_storage_global(:), & + ! Water in rivers, on land points (kg). sfc_water_frac_global(:), & ! Fraction of demand to be met from surface water. supply_irrig_global(:), & @@ -89,30 +95,37 @@ MODULE water_resources_control_mod ! Water that is removed from the system during use, e.g. incorporated into ! manufactured goods (kg). +! Dr Hook variables +INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 +INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 +REAL(KIND=jprb) :: zhook_handle + CONTAINS !############################################################################## -SUBROUTINE water_resources_control( & +SUBROUTINE water_resources_control( global_land_pts, & global_land_index, land_index, map_river_to_land_points, & rivers_index_rp, con_rain_ij, con_snow_ij, & conv_loss_frac, demand_rate_domestic, demand_rate_industry, & demand_rate_livestock, demand_rate_transfers, dvi_cpft, & flandg, frac_irr_soilt, frac_soilt, & frac_surft, grid_area_ij, & - ls_rain_ij, ls_snow_ij, lw_down, smvccl_soilt, & + ls_rain_ij, ls_snow_ij, lw_down, & + res_storage, rfm_surfstore_rp, rivers_sto_rp, smvccl_soilt, & smvcst_soilt, smvcwt_soilt, sthf_soilt, & sw_surft, tl_1_ij, tstar_surft, & icntmax_gb, plant_n_gb, demand_accum, & prec_1_day_av_gb, prec_1_day_av_use_gb, & - rfm_surfstore_rp, rivers_sto_rp, rn_1_day_av_gb, & - rn_1_day_av_use_gb, sfc_water_frac, smcl_soilt, & + rn_1_day_av_gb,rn_1_day_av_use_gb, sfc_water_frac, smcl_soilt, & sthu_irr_soilt, sthu_soilt, & sthzw_soilt, sub_surf_roff, tl_1_day_av_gb, & - tl_1_day_av_use_gb, priority_order, demand_unmet, gw_abstracted, & - gw_avail_start, gw_nr_abstracted, & - irrig_water_gb, net_abstracted_river, sw_abstracted, & - sw_avail_total, water_removed ) + tl_1_day_av_use_gb, & + priority_order, abstracted_res_global, & + net_abstracted_river_global, gw_abstracted, irrig_water_gb, & + abstracted_res, abstracted_river, & + conveyance_loss, demand_unmet, gw_avail_start, gw_nr_abstracted, & + sw_abstracted, sw_avail_total, water_removed ) !------------------------------------------------------------------------------ ! Description: @@ -137,7 +150,7 @@ SUBROUTINE water_resources_control( & USE jules_irrig_mod, ONLY: irr_crop_doell, irr_crop -USE jules_rivers_mod, ONLY: l_rivers, np_rivers +USE jules_rivers_mod, ONLY: l_reservoirs, np_rivers USE jules_soil_mod, ONLY: sm_levels @@ -145,10 +158,8 @@ SUBROUTINE water_resources_control( & USE jules_water_resources_mod, ONLY: & l_have_groundwater, l_have_surface_water, l_water_irrigation, & - nstep_water_res, n_sw_source, nwater_use, sw_river_source, use_irrigation, & - water_res_count - -USE model_grid_mod, ONLY: global_land_pts + nstep_water_res, n_sw_source, nwater_use, & + sw_res_source, sw_river_source, use_irrigation, water_res_count USE parallel_mod, ONLY: is_master_task @@ -166,6 +177,13 @@ SUBROUTINE water_resources_control( & IMPLICIT NONE +!------------------------------------------------------------------------------ +! Scalar arguments with INTENT(IN) +!------------------------------------------------------------------------------ +INTEGER, INTENT(IN) :: & + global_land_pts + ! Number of land points (summed over all tasks). + !------------------------------------------------------------------------------ ! Array arguments with INTENT(IN) !------------------------------------------------------------------------------ @@ -212,6 +230,12 @@ SUBROUTINE water_resources_control( & ! Large-scale snowfall rate (kg m-2 s-1). lw_down(row_length,rows), & ! Surface downward longwave radiation (W m-2). + res_storage(np_rivers), & + ! Water stored in reservoirs (kg). + rfm_surfstore_rp(np_rivers), & + ! River surface storage (m3). + rivers_sto_rp(np_rivers), & + ! River water storage (kg). smvccl_soilt(land_pts,nsoilt,sm_levels), & ! Critical volumetric SMC (cubic m per cubic m of soil). smvcst_soilt(land_pts,nsoilt,sm_levels), & @@ -246,10 +270,6 @@ SUBROUTINE water_resources_control( & ! Average precipitation rate for the current day (kg m-2 s-1). prec_1_day_av_use_gb(land_pts,ndpy,nyav), & ! Daily average precipitation rate (kg m-2 s-1). - rfm_surfstore_rp(np_rivers), & - ! River surface storage (m3). - rivers_sto_rp(np_rivers), & - ! River water storage (kg). rn_1_day_av_gb(land_pts), & ! Average net radiation for the current day (W m-2). rn_1_day_av_use_gb(land_pts,ndpy,nyav), & @@ -282,20 +302,41 @@ SUBROUTINE water_resources_control( & ! Priorities of water demands at each gridpoint, in order of decreasing ! priority. Values are the index in multi-sector arrays. -! Diagnostics and coupling variables. +! Coupling variables and diagnostics. REAL(KIND=real_jlslsm), INTENT(OUT) :: & + !---------------------------------------------------------------------------- + ! Coupling variables + !---------------------------------------------------------------------------- + abstracted_res_global(global_land_pts), & + ! Water abstracted from reservoirs, on global land points (kg). + ! This is used to couple to reservoirs. + ! Note that this has reduced size if it is not required. + net_abstracted_river_global(global_land_pts), & + ! Net abstraction from river, on global land points (kg m-2). + ! Note that this has reduced size if it is not required. + !---------------------------------------------------------------------------- + ! Variables that are passed out for diagnostic purposes. + !---------------------------------------------------------------------------- + gw_abstracted(land_pts), & + ! Water abstracted from renewable groundwater (kg). This is for coupling + ! to a groundwater model - though that is currently done internally to + ! the water resources code. + irrig_water_gb(land_pts), & + ! Water added to soil via irrigation (kg m-2 s-1). This is for coupling + ! to a soil model - though that is currently done internally to the water + ! resources code. + abstracted_res(land_pts), & + ! Water abstracted from reservoirs (kg). Diagnostic only. + abstracted_river(land_pts), & + ! Water abstracted from rivers (kg). Diagnostic only. + conveyance_loss(land_pts), & + ! Water that is lost during conveyance (kg). demand_unmet(land_pts,nwater_use), & ! The part of the demand for water that is not satisfied (kg). - gw_abstracted(land_pts), & - ! Water abstracted from renewable groundwater (kg). gw_avail_start(land_pts), & ! Groundwater that is available for abstraction at start of timestep (kg). gw_nr_abstracted(land_pts), & ! Water abstracted from non-renewable groundwater (kg). - irrig_water_gb(land_pts), & - ! Water added to soil via irrigation (kg m-2 s-1). - net_abstracted_river(land_pts), & - ! Net abstraction from river (kg m-2). sw_abstracted(land_pts,n_sw_source), & ! Water that is abstracted from surface waters (kg). sw_avail_total(land_pts), & @@ -328,8 +369,6 @@ SUBROUTINE water_resources_control( & ! Local array variables !------------------------------------------------------------------------------ REAL(KIND=real_jlslsm) :: & - conveyance_loss(land_pts), & - ! Water that is lost during conveyance (kg). demand_irrig_layer(land_pts,nsoilt,sm_levels), & ! Demand for irrigation water for each soil layer (kg m-2). ! This has sm_levels layers (though we only need nlayer_irrig) for ease of @@ -346,18 +385,11 @@ SUBROUTINE water_resources_control( & ! Area of land in gridbox (m2). return_flow_gw(land_pts), & ! Water that is returned to renewable groundwater after use (kg). - return_flow_sw(land_pts), & + return_flow_river(land_pts), & ! Water that is returned to rivers after use (kg). - river_abstracted(np_rivers), & - ! Water abstracted from river store (kg). supply_irrig(land_pts) ! Water supplied for irrigation (kg). -! Dr Hook variables -INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 -INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 -REAL(KIND=jprb) :: zhook_handle - !------------------------------------------------------------------------------ !end of header IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) @@ -387,17 +419,35 @@ SUBROUTINE water_resources_control( & ! here because on timesteps on which the water resource model is not called ! they are potentially still used for coupling to other components and/or as ! diagnostics. As these fluxes are currently passed as increments (e.g. kg) -! rather than as rates (e.g. kg s-1) they need to be zero on these timesteps. -!------------------------------------------------------------------------------ -! Initialise coupling flux. -net_abstracted_river(:) = 0.0 +! rather than as rates (e.g. kg s-1) they need to be zero on these +! intermediate timesteps. +!------------------------------------------------------------------------------ +! Initialise coupling fluxes. +! These "global" arrays are only allocated at full size on the master task +! (see water_resources_alloc), so they must only be accessed there. +IF ( is_master_task() ) THEN + IF ( l_reservoirs ) THEN + abstracted_res_global(:) = 0.0 + END IF + IF ( sw_river_source > 0 ) THEN + net_abstracted_river_global(:) = 0.0 + END IF +END IF ! Initialise abstractions. gw_abstracted(:) = 0.0 gw_nr_abstracted(:) = 0.0 sw_abstracted(:,:) = 0.0 ! Initialise other fluxes. +conveyance_loss(:) = 0.0 demand_unmet(:,:) = 0.0 water_removed(:) = 0.0 +! Initialise diagnostics. +IF ( l_reservoirs ) THEN + abstracted_res(:) = 0.0 +END IF +IF ( sw_river_source > 0 ) THEN + abstracted_river(:) = 0.0 +END IF !------------------------------------------------------------------------------ ! Add to the accumulated demands (only for prescribed demands). @@ -471,33 +521,70 @@ SUBROUTINE water_resources_control( & ! Fields on land_pts need to be gathered into global equivalents. !---------------------------------------------------------------------------- CALL gather_global_water( priority_order, conv_loss_frac, demand_accum, & - gw_avail_start, sfc_water_frac ) + grid_area_lp, gw_avail_start, sfc_water_frac ) - !---------------------------------------------------------------------------- - ! Calculate the surface water available for abstraction. - ! As this can involve the river grid, we only do this on the master task. - ! Global arrays should have been allocated before this routine is called. - !---------------------------------------------------------------------------- - IF ( l_have_surface_water .AND. is_master_task() ) THEN - CALL calc_avail_surface_water( global_land_index, map_river_to_land_points,& - rivers_index_rp, rfm_surfstore_rp, & - rivers_sto_rp, sw_avail_global ) - ! Save the start of timestep total for diagnostic purposes. - sw_avail_total_start_global(:) = SUM(sw_avail_global,2) - END IF - - !---------------------------------------------------------------------------- - ! Call the top-level routine for the chosen model, on a single processor. - ! Initially only one model is available. - !---------------------------------------------------------------------------- IF ( is_master_task() ) THEN + + !---------------------------------------------------------------------------- + ! Calculate the surface water available for abstraction. + ! As this can involve the river grid, we only do this on the master task. + !---------------------------------------------------------------------------- + IF ( l_have_surface_water ) THEN + + !------------------------------------------------------------------------ + ! Regrid river variables onto land points. + ! Global arrays should have been allocated before this routine is called. + !------------------------------------------------------------------------ + CALL regrid_to_land( global_land_index, map_river_to_land_points, & + rivers_index_rp, res_storage, rfm_surfstore_rp, & + rivers_sto_rp, res_storage_global, & + river_storage_global ) + + ! Calculate water available for abstraction. + CALL calc_avail_surface_water( res_storage_global, & + river_storage_global, & + sw_avail_global ) + + ! Save the total at start of timestep for diagnostic purposes. + sw_avail_total_start_global(:) = SUM(sw_avail_global,2) + END IF + + !-------------------------------------------------------------------------- + ! Call the top-level routine for the chosen model. + ! Initially only one model is available. + !-------------------------------------------------------------------------- CALL water_resources_drive( global_land_pts, priority_order_global, & conv_loss_frac_global, demand_accum_global, & demand_unmet_global, gw_abstracted_global, gw_avail_global, & gw_nr_abstracted_global, sfc_water_frac_global, & sw_abstracted_global, sw_avail_global, water_removed_global, & conveyance_loss_global, return_flow_gw_global, & - return_flow_sw_global, supply_irrig_global ) + return_flow_river_global, supply_irrig_global ) + + ! If reservoirs are modelled, save abstraction in a new variable. + IF ( sw_res_source > 0 ) THEN + abstracted_res_global(:) & + = sw_abstracted_global(:,sw_res_source) + END IF + !-------------------------------------------------------------------------- + ! Decide where return flows should go. + ! This is essentially part of coupling to other parts of the model but the + ! current code structure requires this bit to be done now (before net + ! abstraction from rivers is calculated). + !-------------------------------------------------------------------------- + CALL redirect_return_flows( global_land_pts, return_flow_gw_global, & + return_flow_river_global ) + + !-------------------------------------------------------------------------- + ! Calculate the net abstraction from rivers. + !-------------------------------------------------------------------------- + IF ( sw_river_source > 0 ) THEN + CALL calc_river_flux( global_land_pts, grid_area_global, & + return_flow_river_global, & + sw_abstracted_global(:,sw_river_source), & + net_abstracted_river_global ) + END IF + END IF ! is_master_task !---------------------------------------------------------------------------- @@ -506,9 +593,10 @@ SUBROUTINE water_resources_control( & ! There are no prognostic variables to be scattered. !---------------------------------------------------------------------------- CALL scatter_global_water( conveyance_loss, demand_unmet, gw_abstracted, & - gw_nr_abstracted, return_flow_gw, return_flow_sw, & - sfc_water_frac, supply_irrig, sw_abstracted, & - sw_avail_total, water_removed ) + abstracted_res, abstracted_river, & + gw_nr_abstracted, return_flow_gw, & + return_flow_river, sfc_water_frac, supply_irrig, & + sw_abstracted, sw_avail_total, water_removed ) !---------------------------------------------------------------------------- ! Deallocate global arrays. @@ -516,15 +604,6 @@ SUBROUTINE water_resources_control( & l_allocate = .FALSE. CALL allocate_global_water( l_allocate ) - !---------------------------------------------------------------------------- - ! Calculate the net abstraction from rivers. - !---------------------------------------------------------------------------- - IF ( sw_river_source > 0 ) THEN - CALL calc_river_flux( land_pts, grid_area_lp, return_flow_sw, & - sw_abstracted(:,sw_river_source), & - net_abstracted_river ) - END IF - !---------------------------------------------------------------------------- ! Update soil and groundwater stores. !---------------------------------------------------------------------------- @@ -602,16 +681,10 @@ SUBROUTINE initialise_water_resources( priority_order ) i ! Loop counter. -! Dr Hook variables -INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 -INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 -REAL(KIND=jprb) :: zhook_handle - !------------------------------------------------------------------------------ !end of header IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) - IF ( l_prioritise ) THEN ! Set sector priorities at each location. ! At present these are the same at all locations and it is simply a case @@ -698,11 +771,6 @@ SUBROUTINE accumulate_demand( demand_rate_domestic, & !------------------------------------------------------------------------------ CHARACTER(LEN=*), PARAMETER :: RoutineName = 'ACCUMULATE_DEMAND' -! Dr Hook variables -INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 -INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 -REAL(KIND=jprb) :: zhook_handle - !------------------------------------------------------------------------------ !end of header IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) @@ -773,11 +841,6 @@ SUBROUTINE add_conveyance_loss( conv_loss_frac, demand_accum ) INTEGER :: i ! Loop counter. -! Dr Hook variables -INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 -INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 -REAL(KIND=jprb) :: zhook_handle - !------------------------------------------------------------------------------ !end of header IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) @@ -811,6 +874,8 @@ SUBROUTINE allocate_global_water( l_allocate ) USE ereport_mod, ONLY: ereport +USE jules_rivers_mod, ONLY: l_reservoirs, l_rivers + USE jules_water_resources_mod, ONLY: l_have_groundwater, l_have_surface_water, & l_prioritise, l_water_irrigation, n_sw_source, nwater_use @@ -846,21 +911,26 @@ SUBROUTINE allocate_global_water( l_allocate ) land_size, & ! Size for arrays. land_size_irrig, & - ! Size for arrays for irrigation variables. + ! Size for irrigation arrays. land_size_gw, & ! Size for groundwater arrays. - land_size_sw + land_size_sw, & ! Size for surface water arrays. + land_size_res, & + ! Size for reservoir arrays. + land_size_rivers + ! Size for river arrays. !------------------------------------------------------------------------------ !end of header +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) IF ( l_allocate ) THEN !---------------------------------------------------------------------------- ! Allocate arrays. ! These are allocated at full size only on the master task, and only if a - ! parameterisation is selected - otherwise alocate minimum size. + ! parameterisation is selected - otherwise allocate minimum size. !---------------------------------------------------------------------------- IF ( is_master_task() ) THEN land_size = global_land_pts @@ -874,6 +944,16 @@ SUBROUTINE allocate_global_water( l_allocate ) ELSE land_size_sw = 1 END IF + IF ( l_rivers ) THEN + land_size_rivers = global_land_pts + ELSE + land_size_rivers = 1 + END IF + IF ( l_reservoirs ) THEN + land_size_res = global_land_pts + ELSE + land_size_res = 1 + END IF IF ( l_water_irrigation ) THEN land_size_irrig = global_land_pts ELSE @@ -884,6 +964,8 @@ SUBROUTINE allocate_global_water( l_allocate ) land_size = 1 land_size_irrig = 1 land_size_gw = 1 + land_size_res = 1 + land_size_rivers = 1 land_size_sw = 1 END IF @@ -902,6 +984,8 @@ SUBROUTINE allocate_global_water( l_allocate ) error_sum = error_sum + ERROR ALLOCATE(gw_nr_abstracted_global(land_size_gw), STAT = ERROR) error_sum = error_sum + ERROR + ALLOCATE(res_storage_global(land_size_res), STAT = ERROR) + error_sum = error_sum + ERROR IF ( l_prioritise ) THEN ALLOCATE(priority_order_global(land_size,nwater_use), STAT = ERROR) @@ -910,11 +994,14 @@ SUBROUTINE allocate_global_water( l_allocate ) END IF error_sum = error_sum + ERROR - ! Both SW and GW return flows are needed at full size, even if either - ! source is not active. + ALLOCATE(river_storage_global(land_size_rivers), STAT = ERROR) + error_sum = error_sum + ERROR + + ! Return flows to groundater and rivers are both needed at full size, even + ! if either source is not active. ALLOCATE(return_flow_gw_global(land_size), STAT = ERROR) error_sum = error_sum + ERROR - ALLOCATE(return_flow_sw_global(land_size), STAT = ERROR) + ALLOCATE(return_flow_river_global(land_size), STAT = ERROR) error_sum = error_sum + ERROR IF ( l_have_groundwater .AND. l_have_surface_water ) THEN @@ -935,6 +1022,9 @@ SUBROUTINE allocate_global_water( l_allocate ) ALLOCATE(water_removed_global(land_size), STAT = ERROR) error_sum = error_sum + ERROR + ALLOCATE(grid_area_global(land_size_rivers), STAT = ERROR) + error_sum = error_sum + ERROR + IF ( error_sum == 0 ) THEN ! Initialise arrays. conveyance_loss_global(:) = 0.0 @@ -944,15 +1034,18 @@ SUBROUTINE allocate_global_water( l_allocate ) gw_abstracted_global(:) = 0.0 gw_avail_global(:) = 0.0 gw_nr_abstracted_global(:) = 0.0 + res_storage_global(:) = 0.0 priority_order_global(:,:) = 0 return_flow_gw_global(:) = 0.0 - return_flow_sw_global(:) = 0.0 + return_flow_river_global(:) = 0.0 + river_storage_global(:) = 0.0 sfc_water_frac_global(:) = 0.0 supply_irrig_global(:) = 0.0 sw_abstracted_global(:,:) = 0.0 sw_avail_global(:,:) = 0.0 sw_avail_total_start_global(:) = 0.0 water_removed_global(:) = 0.0 + grid_area_global(:) = 0.0 ELSE errorstatus = 10 CALL ereport( RoutineName, errorstatus, & @@ -965,6 +1058,7 @@ SUBROUTINE allocate_global_water( l_allocate ) ! l_allocate = .FALSE. ! Deallocate arrays, in opposite order to the allocation. !---------------------------------------------------------------------------- + IF ( ALLOCATED(grid_area_global) ) DEALLOCATE(grid_area_global) IF ( ALLOCATED(water_removed_global) ) DEALLOCATE(water_removed_global) IF ( ALLOCATED(sw_avail_total_start_global) ) THEN DEALLOCATE(sw_avail_total_start_global) @@ -973,9 +1067,15 @@ SUBROUTINE allocate_global_water( l_allocate ) IF ( ALLOCATED(sw_abstracted_global) ) DEALLOCATE(sw_abstracted_global) IF ( ALLOCATED(supply_irrig_global) ) DEALLOCATE(supply_irrig_global) IF ( ALLOCATED(sfc_water_frac_global) ) DEALLOCATE(sfc_water_frac_global) + IF ( ALLOCATED(river_storage_global) ) DEALLOCATE(river_storage_global) IF ( ALLOCATED(return_flow_gw_global) ) DEALLOCATE(return_flow_gw_global) - IF ( ALLOCATED(return_flow_sw_global) ) DEALLOCATE(return_flow_sw_global) + IF ( ALLOCATED(return_flow_river_global) ) THEN + DEALLOCATE(return_flow_river_global) + END IF IF ( ALLOCATED(priority_order_global) ) DEALLOCATE(priority_order_global) + IF ( ALLOCATED(res_storage_global) ) THEN + DEALLOCATE(res_storage_global) + END IF IF ( ALLOCATED(gw_nr_abstracted_global)) DEALLOCATE(gw_nr_abstracted_global) IF ( ALLOCATED(gw_avail_global) ) DEALLOCATE(gw_avail_global) IF ( ALLOCATED(gw_abstracted_global) ) DEALLOCATE(gw_abstracted_global) @@ -986,6 +1086,7 @@ SUBROUTINE allocate_global_water( l_allocate ) END IF ! l_allocate +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE allocate_global_water @@ -993,7 +1094,7 @@ END SUBROUTINE allocate_global_water !############################################################################## SUBROUTINE gather_global_water( priority_order, conv_loss_frac, demand_accum, & - gw_avail_start, sfc_water_frac ) + grid_area_lp, gw_avail_start, sfc_water_frac ) !------------------------------------------------------------------------------ ! Description: @@ -1003,7 +1104,7 @@ SUBROUTINE gather_global_water( priority_order, conv_loss_frac, demand_accum, & USE ancil_info, ONLY: land_pts ! for the current task USE jules_water_resources_mod, ONLY: l_have_groundwater, l_prioritise, & - nwater_use, partition_ancil, partition_method + nwater_use, partition_ancil, partition_method, sw_river_source USE model_grid_mod, ONLY: global_land_pts @@ -1024,11 +1125,18 @@ SUBROUTINE gather_global_water( priority_order, conv_loss_frac, demand_accum, & ! Fraction of water that is lost during conveyance from source to user. demand_accum(land_pts,nwater_use), & ! Demands for water accumulated over the water resource timestep (kg). + grid_area_lp(land_pts), & + ! Area of gridbox (m2). gw_avail_start(land_pts), & ! Groundwater that is available for abstraction at start of timestep (kg). sfc_water_frac(land_pts) ! Fraction of demand to be met from surface water. +!------------------------------------------------------------------------------ +! Local parameters +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'GATHER_GLOBAL_WATER' + !------------------------------------------------------------------------------ ! Local variables. !------------------------------------------------------------------------------ @@ -1043,6 +1151,7 @@ SUBROUTINE gather_global_water( priority_order, conv_loss_frac, demand_accum, & !------------------------------------------------------------------------------ !end of header +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) CALL gather_land_field( conv_loss_frac, conv_loss_frac_global ) @@ -1068,6 +1177,12 @@ SUBROUTINE gather_global_water( priority_order, conv_loss_frac, demand_accum, & CALL gather_land_field( sfc_water_frac, sfc_water_frac_global ) END IF +! We only need to gather grid_area_lp if there is abstraction from rivers. +IF ( sw_river_source > 0 ) THEN + CALL gather_land_field( grid_area_lp, grid_area_global ) +END IF + +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE gather_global_water @@ -1075,9 +1190,12 @@ END SUBROUTINE gather_global_water !############################################################################## SUBROUTINE scatter_global_water( conveyance_loss, demand_unmet, gw_abstracted, & + abstracted_res, & + abstracted_river, & gw_nr_abstracted, return_flow_gw, & - return_flow_sw, sfc_water_frac, supply_irrig, & - sw_abstracted, sw_avail_total, water_removed ) + return_flow_river, sfc_water_frac, & + supply_irrig, sw_abstracted, sw_avail_total, & + water_removed ) !------------------------------------------------------------------------------ ! Description: @@ -1086,9 +1204,12 @@ SUBROUTINE scatter_global_water( conveyance_loss, demand_unmet, gw_abstracted, & USE ancil_info, ONLY: land_pts ! for the current task +USE jules_rivers_mod, ONLY: l_reservoirs, l_rivers + USE jules_water_resources_mod, ONLY: l_have_groundwater, l_have_surface_water, & l_water_irrigation, n_sw_source, nwater_use, & - partition_calc_from_stores, partition_method + partition_calc_from_stores, partition_method, & + sw_res_source, sw_river_source USE model_grid_mod, ONLY: global_land_pts @@ -1108,9 +1229,13 @@ SUBROUTINE scatter_global_water( conveyance_loss, demand_unmet, gw_abstracted, & ! Water abstracted from renewable groundwater (kg). gw_nr_abstracted(land_pts), & ! Water abstracted from non-renewable groundwater (kg). + abstracted_res(land_pts), & + ! Water abstracted from reservoirs (kg). + abstracted_river(land_pts), & + ! Water abstracted from rivers (kg). return_flow_gw(land_pts), & ! Water that is returned to renewable groundwater after use (kg). - return_flow_sw(land_pts), & + return_flow_river(land_pts), & ! Water that is returned to rivers after use (kg). sfc_water_frac(land_pts), & ! Target fraction of demand to be met by surface water. @@ -1125,6 +1250,11 @@ SUBROUTINE scatter_global_water( conveyance_loss, demand_unmet, gw_abstracted, & ! Water that is removed from the system during use, e.g. incorporated into ! manufactured goods (kg). +!------------------------------------------------------------------------------ +! Local parameters +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'SCATTER_GLOBAL_WATER' + !------------------------------------------------------------------------------ ! Local scalar variables !------------------------------------------------------------------------------ @@ -1132,6 +1262,7 @@ SUBROUTINE scatter_global_water( conveyance_loss, demand_unmet, gw_abstracted, & !------------------------------------------------------------------------------ !end of header +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) ! Fields that are always required. CALL scatter_land_field( conveyance_loss_global, conveyance_loss ) @@ -1146,37 +1277,266 @@ SUBROUTINE scatter_global_water( conveyance_loss, demand_unmet, gw_abstracted, & CALL scatter_land_field( gw_nr_abstracted_global, gw_nr_abstracted ) END IF ! Note that return_flow_gw is always used, even if groundwater is not -! modelled. +! modelled (in that case it is added to runoff). CALL scatter_land_field( return_flow_gw_global, return_flow_gw ) +! Return flow to rivers is always calculated (but can be zero). However +! it is only required further if we have rivers. +IF ( l_rivers ) THEN + CALL scatter_land_field( return_flow_river_global, return_flow_river ) +END IF + ! Fields that are only required if we are modelling surface water sources. IF ( l_have_surface_water ) THEN + ! Surface water fraction only needs to be scattered if it was calculated by ! the master task. If it was read as an ancillary field, nothing to do here. IF ( partition_method == partition_calc_from_stores ) THEN CALL scatter_land_field( sfc_water_frac_global, sfc_water_frac ) END IF - ! Surface water always includes rivers, so we always calculate return flow. - CALL scatter_land_field( return_flow_sw_global, return_flow_sw ) + DO i = 1, n_sw_source CALL scatter_land_field( sw_abstracted_global(:,i), sw_abstracted(:,i) ) END DO CALL scatter_land_field( sw_avail_total_start_global, sw_avail_total ) -END IF + + IF ( l_reservoirs ) THEN + CALL scatter_land_field( sw_abstracted_global(:,sw_res_source), & + abstracted_res ) + END IF + + ! Save diagnostic of abstraction from rivers. + IF ( sw_river_source > 0 ) THEN + CALL scatter_land_field( sw_abstracted_global(:,sw_river_source), & + abstracted_river ) + END IF + +END IF ! l_have_surface_water ! Fields for irrigation. IF ( l_water_irrigation ) THEN CALL scatter_land_field( supply_irrig_global, supply_irrig ) END IF +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE scatter_global_water !############################################################################## !############################################################################## -SUBROUTINE calc_river_flux( land_pts, grid_area_lp, return_flow_sw, & - sw_abstracted_river, net_abstracted_river ) +SUBROUTINE regrid_to_land( global_land_index, map_river_to_land_points, & + rivers_index_rp, & + res_storage, rfm_surfstore_rp, rivers_sto_rp, & + res_storage_global, & + river_storage_global ) + +!------------------------------------------------------------------------------ +! Description: +! Regrid water resource-related variables from river grid to land grid. +!------------------------------------------------------------------------------ + +USE ereport_mod, ONLY: ereport + +USE jules_rivers_mod, ONLY: i_river_vn, l_reservoirs, & + l_rivers, np_rivers, rivers_rfm, rivers_trip + +USE model_grid_mod, ONLY: global_land_pts + +USE rivers_regrid_mod, ONLY: rivpts_to_landpts + +USE water_constants_mod, ONLY: rho_water + + +IMPLICIT NONE + +!------------------------------------------------------------------------------ +! Array arguments with INTENT(IN). +!------------------------------------------------------------------------------ +INTEGER, INTENT(IN) :: & + global_land_index(global_land_pts), & + ! List of indices for the land model grid. + map_river_to_land_points(np_rivers), & + ! List of coincident land point numbers, on river points. + rivers_index_rp(np_rivers) + ! Index of points where routing is calculated. + +REAL(KIND=real_jlslsm), INTENT(IN) :: & + res_storage(np_rivers), & + ! Water stored in reservoirs (kg). + rfm_surfstore_rp(np_rivers), & + ! River surface storage (m3). + rivers_sto_rp(np_rivers) + ! River water storage (kg). + +!------------------------------------------------------------------------------ +! Array arguments with INTENT(OUT). +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(OUT) :: & + res_storage_global(global_land_pts), & + ! Water stored in reservoirs, on land points (kg). + river_storage_global(global_land_pts) + ! Water in rivers, on land points (kg). + +!------------------------------------------------------------------------------ +! Local parameters. +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'REGRID_TO_LAND' + +!------------------------------------------------------------------------------ +! Local scalar variables. +!------------------------------------------------------------------------------ +INTEGER :: & + errorstatus + ! Error value. + +!------------------------------------------------------------------------------ +!end of header +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + +IF ( l_rivers ) THEN + + !---------------------------------------------------------------------------- + ! Convert river storage to a variable on land points. + ! Select code for the current river model as variables differ between models. + ! This is later used to define the river water available for abstraction, + ! which informs the approach taken below for rivers_rfm. + !---------------------------------------------------------------------------- + + ! Initialise. + river_storage_global(:) = 0.0 + + SELECT CASE ( i_river_vn ) + + CASE ( rivers_rfm ) + ! We will only include the surface store in the available water. + CALL rivpts_to_landpts( global_land_pts, np_rivers, & + map_river_to_land_points,global_land_index, & + rivers_index_rp, rfm_surfstore_rp, & + river_storage_global ) + ! Convert units from m3 to kg. + river_storage_global(:) = river_storage_global(:) * rho_water + + CASE ( rivers_trip ) + CALL rivpts_to_landpts( global_land_pts, np_rivers, & + map_river_to_land_points, global_land_index, & + rivers_index_rp, rivers_sto_rp, & + river_storage_global ) + + CASE DEFAULT + + errorstatus = 101 ! a fatal error + CALL ereport(RoutineName, errorstatus, 'Unknown value of i_river_vn.') + + END SELECT ! i_river_vn + +END IF ! l_rivers + +!------------------------------------------------------------------------------ +! Get reservoir storage onto land points. +!------------------------------------------------------------------------------ +IF ( l_reservoirs ) THEN + ! Initialise. + res_storage_global(:) = 0.0 + CALL rivpts_to_landpts( global_land_pts, np_rivers, & + map_river_to_land_points, global_land_index, & + rivers_index_rp, res_storage, & + res_storage_global ) +END IF + +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) +RETURN +END SUBROUTINE regrid_to_land + +!############################################################################## +!############################################################################## + +SUBROUTINE redirect_return_flows( global_land_pts, & + return_flow_gw_global, & + return_flow_river_global ) + +!------------------------------------------------------------------------------ +! Description: +! Change the destination of return flows depending on what sinks are +! available. +!------------------------------------------------------------------------------ + +USE jules_rivers_mod, ONLY: l_rivers +USE jules_water_resources_mod, ONLY: l_have_renew_gwater + +IMPLICIT NONE + +!------------------------------------------------------------------------------ +! Scalar arguments with INTENT(IN) +!------------------------------------------------------------------------------ +INTEGER, INTENT(IN) :: & + global_land_pts + ! The number of land points. + +!------------------------------------------------------------------------------ +! Array arguments with INTENT(IN OUT) +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(IN OUT) :: & + return_flow_gw_global(global_land_pts), & + ! Water that is returned to renewable groundwater after use (kg). + return_flow_river_global(global_land_pts) + ! Water that is returned to rivers after use (kg). + +!------------------------------------------------------------------------------ +! Local parameters +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'REDIRECT_RETURN_FLOWS' + +!------------------------------------------------------------------------------ +! Local variables. +!------------------------------------------------------------------------------ +INTEGER :: & + l ! Loop counter. + +!------------------------------------------------------------------------------ +!end of header +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + +!------------------------------------------------------------------------------ +! Return flows to groundwater and rivers have been calculated assuming those +! sinks are available. If the preferred sink is not modelled the other is used +! (if available). If neither is available we add to the groundwater return flow +! because later code will divert that to runoff. +! Note that if both rivers and renewable groundwater are modelled nothing is +! changed in this subroutine. +!------------------------------------------------------------------------------ +IF ( l_rivers .AND. .NOT. l_have_renew_gwater ) THEN + + ! River are modelled but renewable groundwater is not. + ! Direct all return flows to rivers. + DO l=1,global_land_pts + return_flow_river_global(l) = return_flow_river_global(l) & + + return_flow_gw_global(l) + return_flow_gw_global(l) = 0.0 + END DO + +ELSE IF ( .NOT. l_rivers ) THEN + + ! Rivers are not modelled. Groundwater might or might not be modelled. + ! Direct all return flows to groundwater (for now). + DO l=1,global_land_pts + return_flow_gw_global(l) = return_flow_gw_global(l) & + + return_flow_river_global(l) + return_flow_river_global(l) = 0.0 + END DO + +END IF + +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) +RETURN +END SUBROUTINE redirect_return_flows + +!############################################################################## +!############################################################################## +SUBROUTINE calc_river_flux( global_land_pts, grid_area_global, & + return_flow_river_global, & + sw_abstracted_river_global, & + net_abstracted_river_global ) !------------------------------------------------------------------------------ ! Description: @@ -1189,27 +1549,32 @@ SUBROUTINE calc_river_flux( land_pts, grid_area_lp, return_flow_sw, & ! Scalar arguments with INTENT(IN) !------------------------------------------------------------------------------ INTEGER, INTENT(IN) :: & - land_pts + global_land_pts ! The number of land points. !------------------------------------------------------------------------------ ! Array arguments with INTENT(IN) !------------------------------------------------------------------------------ REAL(KIND=real_jlslsm), INTENT(IN) :: & - grid_area_lp(land_pts), & + grid_area_global(global_land_pts), & ! Area of gridbox (m2). - return_flow_sw(land_pts), & - ! Water that is returned to surface waters after use (kg). - sw_abstracted_river(land_pts) + return_flow_river_global(global_land_pts), & + ! Water that is returned to rivers after use (kg). + sw_abstracted_river_global(global_land_pts) ! Water abstracted from river (kg). !------------------------------------------------------------------------------ ! Array arguments with INTENT(OUT) !------------------------------------------------------------------------------ REAL(KIND=real_jlslsm), INTENT(OUT) :: & - net_abstracted_river(land_pts) + net_abstracted_river_global(global_land_pts) ! Net abstraction from river (kg m-2). +!------------------------------------------------------------------------------ +! Local parameters +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'CALC_RIVER_FLUX' + !------------------------------------------------------------------------------ ! Local variables. !------------------------------------------------------------------------------ @@ -1218,13 +1583,16 @@ SUBROUTINE calc_river_flux( land_pts, grid_area_lp, return_flow_sw, & !------------------------------------------------------------------------------ !end of header +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) ! Calculate net abstraction and change units from kg to kg m-2. -DO l=1,land_pts - net_abstracted_river(l) = ( sw_abstracted_river(l) - return_flow_sw(l) ) & - / grid_area_lp(l) +DO l=1,global_land_pts + net_abstracted_river_global(l) = ( sw_abstracted_river_global(l) & + - return_flow_river_global(l) ) & + / grid_area_global(l) END DO +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE calc_river_flux diff --git a/src/control/shared/water_resources_vars_mod.F90 b/src/control/shared/water_resources_vars_mod.F90 index f6b9b8d5..d6cf77bf 100644 --- a/src/control/shared/water_resources_vars_mod.F90 +++ b/src/control/shared/water_resources_vars_mod.F90 @@ -9,7 +9,7 @@ ! [Met Office Ref SC0237] !******************************COPYRIGHT************************************** -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! Description: ! Contains variables and field types for water resource modelling. ! @@ -18,7 +18,7 @@ ! ! Code Description: ! Language: Fortran 90. -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ MODULE water_resources_vars_mod @@ -52,7 +52,7 @@ MODULE water_resources_vars_mod REAL(KIND=real_jlslsm), ALLOCATABLE :: & !-------------------------------------------------------------------------- - ! Ancillary fields. + ! Ancillary fields on land points. !-------------------------------------------------------------------------- conv_loss_frac(:), & ! Fraction of water that is lost during conveyance from source to user. @@ -72,13 +72,15 @@ MODULE water_resources_vars_mod demand_rate_transfers(:), & ! Demand for water for (explicit) transfers (kg s-1). !-------------------------------------------------------------------------- - ! Flux for coupling to rivers. + ! Fluxes for coupling. !-------------------------------------------------------------------------- - net_abstracted_river(:), & - ! Net abstraction from river (kg). + net_abstracted_river_global(:), & + ! Net abstraction from river, on global land points (kg m-2). + abstracted_res_global(:), & + ! Water abstracted from reservoirs, on global land points (kg). !-------------------------------------------------------------------------- - ! Diagnostics (also used internally). - ! All abstractions are gross (not net of returns) fluxes unless stated + ! Internal/work variables (some of which are also used as diagnostics). + ! All abstractions are gross fluxes (not net of returns) unless stated ! otherwise. !-------------------------------------------------------------------------- demand_accum(:,:), & @@ -89,7 +91,9 @@ MODULE water_resources_vars_mod demand_unmet(:,:), & ! The part of the demand for water that is not satisfied (kg). gw_abstracted(:), & - ! Water abstracted from renewable groundwater (kg). + ! Water abstracted from renewable groundwater (kg). This is for coupling + ! to a groundwater model - though that is currently done internally to + ! the water resources code. gw_avail(:), & ! Groundwater that is available for abstraction at the start of the ! timestep (kg). This does not include "non-renewable" groundwater. @@ -100,9 +104,18 @@ MODULE water_resources_vars_mod sw_avail_total(:), & ! Surface water that is available for abstraction at start of timestep, ! summed over sources (kg). - water_removed(:) + !-------------------------------------------------------------------------- + ! Diagnostics. + !-------------------------------------------------------------------------- + abstracted_river(:), & + ! Water abstracted from rivers (kg). + conveyance_loss(:), & + ! Water that is lost during conveyance (kg). + water_removed(:), & ! Water that is removed from the system during use, e.g. incorporated ! into manufactured goods (kg). + abstracted_res(:) + ! Water abstracted from reservoirs (kg). END TYPE @@ -117,7 +130,8 @@ MODULE water_resources_vars_mod REAL(KIND=real_jlslsm), POINTER :: demand_rate_industry(:) REAL(KIND=real_jlslsm), POINTER :: demand_rate_livestock(:) REAL(KIND=real_jlslsm), POINTER :: demand_rate_transfers(:) - REAL(KIND=real_jlslsm), POINTER :: net_abstracted_river(:) + REAL(KIND=real_jlslsm), POINTER :: net_abstracted_river_global(:) + REAL(KIND=real_jlslsm), POINTER :: abstracted_res_global(:) REAL(KIND=real_jlslsm), POINTER :: demand_accum(:,:) REAL(KIND=real_jlslsm), POINTER :: demand_unmet(:,:) REAL(KIND=real_jlslsm), POINTER :: gw_abstracted(:) @@ -125,7 +139,10 @@ MODULE water_resources_vars_mod REAL(KIND=real_jlslsm), POINTER :: gw_nr_abstracted(:) REAL(KIND=real_jlslsm), POINTER :: sw_abstracted(:,:) REAL(KIND=real_jlslsm), POINTER :: sw_avail_total(:) + REAL(KIND=real_jlslsm), POINTER :: abstracted_river(:) + REAL(KIND=real_jlslsm), POINTER :: conveyance_loss(:) REAL(KIND=real_jlslsm), POINTER :: water_removed(:) + REAL(KIND=real_jlslsm), POINTER :: abstracted_res(:) END TYPE @@ -137,11 +154,12 @@ MODULE water_resources_vars_mod !############################################################################## -SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & - l_have_groundwater, l_have_surface_water, & - l_water_domestic, l_water_industry, l_water_irrigation, & - l_water_livestock, l_water_resources, l_water_transfers, & - water_resources_data ) +SUBROUTINE water_resources_alloc( global_land_pts, land_pts, n_sw_source, & + nwater_use, sw_river_source, l_have_groundwater, & + l_have_surface_water, l_is_master_task, & + l_reservoirs, l_water_domestic, l_water_industry, & + l_water_irrigation, l_water_livestock, l_water_resources, & + l_water_transfers, water_resources_data ) !No USE statements other than Dr Hook USE parkind1, ONLY: jprb, jpim @@ -153,12 +171,16 @@ SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & ! Scalar arguments with INTENT(IN) !------------------------------------------------------------------------------ INTEGER, INTENT(IN) :: & + global_land_pts, & + ! Number of land points (total over all tasks). land_pts, & ! Number of land points (current task). n_sw_source, & ! Number of surface water sources. - nwater_use + nwater_use, & ! Number of water resource sectors that are considered. + sw_river_source + ! Indicates if rivers are modelled. LOGICAL, INTENT(IN) :: & l_have_groundwater, & @@ -166,7 +188,10 @@ SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & ! non-renewable). l_have_surface_water, & ! Flag indicating if we have surface water represented (e.g. rivers). - ! TRUE means n_sw_source > 0. + l_is_master_task, & + ! Flag indicating if this call is from the master task. + l_reservoirs, & + ! Switch for reservoirs. l_water_domestic, & ! Switch to consider demand for water for domestic use. l_water_industry, & @@ -189,7 +214,9 @@ SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & ! Local variables !------------------------------------------------------------------------------ INTEGER :: & - land_pts_dim, land_pts_gw, land_pts_sw, n_sw_source_dim, nwater_use_dim + global_land_pts_res, global_land_pts_rivers, & + land_pts_dim, land_pts_gw, land_pts_res, land_pts_sw, & + land_pts_rivers, n_sw_source_dim, nwater_use_dim ! Sizes used when allocating arrays. INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 @@ -203,46 +230,65 @@ SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! Arrays are always allocated, but with minimal size if the science is not ! selected. Decide on sizes. -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! Set default sizes that are used if water resources (or a particular part of ! the parameterisation) are not selected. -land_pts_dim = 1 -land_pts_gw = 1 -land_pts_sw = 1 -n_sw_source_dim = 1 -nwater_use_dim = 1 +global_land_pts_res = 1 +global_land_pts_rivers = 1 +land_pts_res = 1 +land_pts_dim = 1 +land_pts_gw = 1 +land_pts_rivers = 1 +land_pts_sw = 1 +n_sw_source_dim = 1 +nwater_use_dim = 1 IF ( l_water_resources ) THEN land_pts_dim = land_pts n_sw_source_dim = n_sw_source nwater_use_dim = nwater_use IF ( l_have_groundwater ) THEN - land_pts_gw = land_pts + land_pts_gw = land_pts END IF IF ( l_have_surface_water ) THEN - land_pts_sw = land_pts + land_pts_sw = land_pts END IF -END IF + IF ( sw_river_source > 0 ) THEN + land_pts_rivers = land_pts + ! The "global" variable only needs to be full size on master task. + IF ( l_is_master_task ) THEN + global_land_pts_rivers = global_land_pts + END IF + END IF + IF ( l_reservoirs ) THEN + land_pts_res = land_pts + ! The "global" variable only needs to be full size on master task. + IF ( l_is_master_task ) THEN + global_land_pts_res = global_land_pts + END IF + END IF +END IF ! l_water_resources !----------------------------------------------------------------------------- ! Priority order. !----------------------------------------------------------------------------- ALLOCATE( water_resources_data%priority_order(land_pts_dim,nwater_use_dim) ) -!----------------------------------------------------------------------------- -! Ancillary fields. -! Although these are not required for every configuration, they are commonly -! required and hence for convenience we always allocate them. -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ +! Ancillary fields on land points. +! Although these are not required for every configuration that includes water +! resources they are commonly required and hence for convenience we always +! allocate them. +!------------------------------------------------------------------------------ ALLOCATE( water_resources_data%conv_loss_frac(land_pts_dim) ) ALLOCATE( water_resources_data%sfc_water_frac(land_pts_dim) ) -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! Individual demands (which can be prescibed). ! We allocate a minimum size if a sector is not being used. -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ IF ( l_water_domestic ) THEN ALLOCATE( water_resources_data%demand_rate_domestic(land_pts_dim) ) ELSE @@ -267,17 +313,29 @@ SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & ALLOCATE( water_resources_data%demand_rate_transfers(1) ) END IF -!----------------------------------------------------------------------------- -! Coupling to rivers - which are always modelled if l_have_surface_water=T. -!----------------------------------------------------------------------------- -ALLOCATE( water_resources_data%net_abstracted_river(land_pts_sw) ) +!------------------------------------------------------------------------------ +! Coupling to reservoirs. +! Only allocated at full size if reservoirs are modelled. +! Allocated space for all (global) land points because both the water resource +! and river codes operate globally. +!------------------------------------------------------------------------------ +ALLOCATE( water_resources_data & + %abstracted_res_global(global_land_pts_res) ) -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ +! Coupling to rivers. +! Only allocated at full size if rivers are modelled. +! Allocated space for all (global) land points because both the water resource +! and river codes operate globally. +!------------------------------------------------------------------------------ +ALLOCATE( water_resources_data & + %net_abstracted_river_global(global_land_pts_rivers) ) + +!------------------------------------------------------------------------------ ! Other variables. -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ALLOCATE( water_resources_data%demand_accum(land_pts_dim,nwater_use_dim) ) ALLOCATE( water_resources_data%demand_unmet(land_pts_dim,nwater_use_dim) ) -ALLOCATE( water_resources_data%water_removed(land_pts_dim) ) ! Groundwater variables. ALLOCATE( water_resources_data%gw_abstracted(land_pts_gw) ) ALLOCATE( water_resources_data%gw_avail(land_pts_gw) ) @@ -285,10 +343,15 @@ SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & ! Surface water variables. ALLOCATE( water_resources_data%sw_abstracted(land_pts_sw,n_sw_source_dim) ) ALLOCATE( water_resources_data%sw_avail_total(land_pts_sw) ) +! Diagnostics. +ALLOCATE( water_resources_data%abstracted_res(land_pts_res) ) +ALLOCATE( water_resources_data%abstracted_river(land_pts_rivers) ) +ALLOCATE( water_resources_data%conveyance_loss(land_pts_dim) ) +ALLOCATE( water_resources_data%water_removed(land_pts_dim) ) -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! Initialise arrays. -!----------------------------------------------------------------------------- +!------------------------------------------------------------------------------ water_resources_data%priority_order(:,:) = 0 water_resources_data%conv_loss_frac(:) = 0.0 water_resources_data%sfc_water_frac(:) = 0.0 @@ -296,15 +359,19 @@ SUBROUTINE water_resources_alloc( land_pts, n_sw_source, nwater_use, & water_resources_data%demand_rate_industry(:) = 0.0 water_resources_data%demand_rate_livestock(:) = 0.0 water_resources_data%demand_rate_transfers(:) = 0.0 -water_resources_data%net_abstracted_river(:) = 0.0 +water_resources_data%abstracted_res_global(:) = 0.0 +water_resources_data%net_abstracted_river_global(:) = 0.0 water_resources_data%demand_accum(:,:) = 0.0 water_resources_data%demand_unmet(:,:) = 0.0 -water_resources_data%water_removed(:) = 0.0 water_resources_data%gw_abstracted(:) = 0.0 water_resources_data%gw_avail(:) = 0.0 water_resources_data%gw_nr_abstracted(:) = 0.0 water_resources_data%sw_abstracted(:,:) = 0.0 water_resources_data%sw_avail_total(:) = 0.0 +water_resources_data%abstracted_res(:) = 0.0 +water_resources_data%abstracted_river(:) = 0.0 +water_resources_data%conveyance_loss(:) = 0.0 +water_resources_data%water_removed(:) = 0.0 IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN @@ -340,22 +407,27 @@ SUBROUTINE water_resources_dealloc(water_resources_data) IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) -DEALLOCATE( water_resources_data%priority_order ) -DEALLOCATE( water_resources_data%conv_loss_frac ) -DEALLOCATE( water_resources_data%sfc_water_frac ) -DEALLOCATE( water_resources_data%demand_rate_domestic ) -DEALLOCATE( water_resources_data%demand_rate_industry ) -DEALLOCATE( water_resources_data%demand_rate_livestock ) -DEALLOCATE( water_resources_data%demand_rate_transfers ) -DEALLOCATE( water_resources_data%net_abstracted_river ) -DEALLOCATE( water_resources_data%demand_accum ) -DEALLOCATE( water_resources_data%demand_unmet ) +! Deallocate in reverse order to allocation. DEALLOCATE( water_resources_data%water_removed ) -DEALLOCATE( water_resources_data%gw_abstracted ) -DEALLOCATE( water_resources_data%gw_avail ) -DEALLOCATE( water_resources_data%gw_nr_abstracted ) -DEALLOCATE( water_resources_data%sw_abstracted ) +DEALLOCATE( water_resources_data%conveyance_loss ) +DEALLOCATE( water_resources_data%abstracted_river ) +DEALLOCATE( water_resources_data%abstracted_res ) DEALLOCATE( water_resources_data%sw_avail_total ) +DEALLOCATE( water_resources_data%sw_abstracted ) +DEALLOCATE( water_resources_data%gw_nr_abstracted ) +DEALLOCATE( water_resources_data%gw_avail ) +DEALLOCATE( water_resources_data%demand_unmet ) +DEALLOCATE( water_resources_data%gw_abstracted ) +DEALLOCATE( water_resources_data%demand_accum ) +DEALLOCATE( water_resources_data%net_abstracted_river_global ) +DEALLOCATE( water_resources_data%abstracted_res_global ) +DEALLOCATE( water_resources_data%demand_rate_transfers ) +DEALLOCATE( water_resources_data%demand_rate_livestock ) +DEALLOCATE( water_resources_data%demand_rate_industry ) +DEALLOCATE( water_resources_data%demand_rate_domestic ) +DEALLOCATE( water_resources_data%sfc_water_frac ) +DEALLOCATE( water_resources_data%conv_loss_frac ) +DEALLOCATE( water_resources_data%priority_order ) IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN @@ -412,16 +484,21 @@ SUBROUTINE water_resources_assoc(water_resources,water_resources_data) => water_resources_data%demand_rate_livestock water_resources%demand_rate_transfers & => water_resources_data%demand_rate_transfers -water_resources%net_abstracted_river & - => water_resources_data%net_abstracted_river +water_resources%abstracted_res_global & + => water_resources_data%abstracted_res_global +water_resources%net_abstracted_river_global & + => water_resources_data%net_abstracted_river_global water_resources%demand_accum => water_resources_data%demand_accum water_resources%demand_unmet => water_resources_data%demand_unmet -water_resources%water_removed => water_resources_data%water_removed water_resources%gw_abstracted => water_resources_data%gw_abstracted water_resources%gw_avail => water_resources_data%gw_avail water_resources%gw_nr_abstracted => water_resources_data%gw_nr_abstracted water_resources%sw_abstracted => water_resources_data%sw_abstracted water_resources%sw_avail_total => water_resources_data%sw_avail_total +water_resources%abstracted_res => water_resources_data%abstracted_res +water_resources%abstracted_river => water_resources_data%abstracted_river +water_resources%conveyance_loss => water_resources_data%conveyance_loss +water_resources%water_removed => water_resources_data%water_removed IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN @@ -461,15 +538,19 @@ SUBROUTINE water_resources_nullify(water_resources) NULLIFY( water_resources%demand_rate_industry ) NULLIFY( water_resources%demand_rate_livestock ) NULLIFY( water_resources%demand_rate_transfers ) -NULLIFY( water_resources%net_abstracted_river ) +NULLIFY( water_resources%abstracted_res_global ) +NULLIFY( water_resources%net_abstracted_river_global ) NULLIFY( water_resources%demand_accum ) NULLIFY( water_resources%demand_unmet ) -NULLIFY( water_resources%water_removed ) NULLIFY( water_resources%gw_abstracted ) NULLIFY( water_resources%gw_avail ) NULLIFY( water_resources%gw_nr_abstracted ) NULLIFY( water_resources%sw_abstracted ) NULLIFY( water_resources%sw_avail_total ) +NULLIFY( water_resources%abstracted_res ) +NULLIFY( water_resources%abstracted_river ) +NULLIFY( water_resources%conveyance_loss ) +NULLIFY( water_resources%water_removed ) IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN diff --git a/src/science/water_resources/calc_avail_water_mod.F90 b/src/science/water_resources/calc_avail_water_mod.F90 index 1c96bede..879737f9 100644 --- a/src/science/water_resources/calc_avail_water_mod.F90 +++ b/src/science/water_resources/calc_avail_water_mod.F90 @@ -162,10 +162,9 @@ END SUBROUTINE calc_avail_groundwater !############################################################################## !############################################################################## -SUBROUTINE calc_avail_surface_water( global_land_index, & - map_river_to_land_points, & - rivers_index_rp, rfm_surfstore_rp, & - rivers_sto_rp, sw_avail_global ) +SUBROUTINE calc_avail_surface_water( res_storage_global, & + river_storage_global, & + sw_avail_global ) !------------------------------------------------------------------------------ ! Description: @@ -175,28 +174,21 @@ SUBROUTINE calc_avail_surface_water( global_land_index, & USE model_grid_mod, ONLY: global_land_pts -USE jules_rivers_mod, ONLY: l_rivers, np_rivers +USE jules_rivers_mod, ONLY: l_reservoirs, l_rivers -USE jules_water_resources_mod, ONLY: n_sw_source, sw_river_source +USE jules_water_resources_mod, ONLY: n_sw_source, & + sw_res_source, sw_river_source IMPLICIT NONE !------------------------------------------------------------------------------ ! Array arguments with INTENT(IN). !------------------------------------------------------------------------------ -INTEGER, INTENT(IN) :: & - global_land_index(global_land_pts), & - ! List of indices for the land model grid. - map_river_to_land_points(np_rivers), & - ! List of coincident land point numbers, on river points. - rivers_index_rp(np_rivers) - ! Index of points where routing is calculated. - REAL(KIND=real_jlslsm), INTENT(IN) :: & - rfm_surfstore_rp(np_rivers), & - ! River surface storage (m3). - rivers_sto_rp(np_rivers) - ! River water storage (kg). + river_storage_global(global_land_pts), & + ! Water in rivers, on land points (kg). + res_storage_global(global_land_pts) + ! Water stored in reservoirs, on land points (kg). !------------------------------------------------------------------------------ ! Array arguments with INTENT(OUT). @@ -220,19 +212,20 @@ SUBROUTINE calc_avail_surface_water( global_land_index, & IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) !------------------------------------------------------------------------------ -! Initialise values. +! Initialise output values. !------------------------------------------------------------------------------ sw_avail_global(:,:) = 0.0 -!------------------------------------------------------------------------------ -! Calculate the available water in rivers. -! At present this is the only possible source of surface water. -!------------------------------------------------------------------------------ IF ( l_rivers ) THEN - CALL calc_avail_river_water( global_land_index, map_river_to_land_points, & - rivers_index_rp, rfm_surfstore_rp, & - rivers_sto_rp, & - sw_avail_global(:,sw_river_source) ) + + ! Include water available from rivers. + sw_avail_global(:,sw_river_source) = river_storage_global(:) + + IF ( l_reservoirs ) THEN + ! Include water stored in reservoirs. + sw_avail_global(:,sw_res_source) = res_storage_global(:) + END IF + END IF IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) diff --git a/src/science/water_resources/water_resources_drive.F90 b/src/science/water_resources/water_resources_drive.F90 index af7f0465..af6df544 100644 --- a/src/science/water_resources/water_resources_drive.F90 +++ b/src/science/water_resources/water_resources_drive.F90 @@ -34,7 +34,7 @@ MODULE water_resources_drive_mod PUBLIC water_resources_drive ! Module parameters. -CHARACTER(LEN=*), PARAMETER, PRIVATE :: & +CHARACTER(LEN=*), PARAMETER :: & ModuleName = 'WATER_RESOURCES_DRIVE_MOD' CONTAINS @@ -47,7 +47,7 @@ SUBROUTINE water_resources_drive( global_land_pts, priority_order, & gw_nr_abstracted, sfc_water_frac, & sw_abstracted, sw_avail, water_removed, & conveyance_loss, return_flow_gw, & - return_flow_sw, supply_irrig ) + return_flow_river, supply_irrig ) !------------------------------------------------------------------------------ ! Description: @@ -118,7 +118,7 @@ SUBROUTINE water_resources_drive( global_land_pts, priority_order, & ! Water that is lost during conveyance (kg). return_flow_gw(global_land_pts), & ! Water that is returned to renewable groundwater after use (kg). - return_flow_sw(global_land_pts), & + return_flow_river(global_land_pts), & ! Water that is returned to rivers after use (kg). supply_irrig(global_land_pts) ! Water supplied for irrigation (kg). @@ -139,7 +139,7 @@ SUBROUTINE water_resources_drive( global_land_pts, priority_order, & ! Local array variables !------------------------------------------------------------------------------ REAL(KIND=real_jlslsm) :: & - abstracted(global_land_pts,nwater_use), & + abstracted_per_use(global_land_pts,nwater_use), & ! Water abstracted for each use (kg). conveyance_loss_use(global_land_pts,nwater_use), & ! Water that is lost during conveyance, for each water use (kg). @@ -216,23 +216,24 @@ SUBROUTINE water_resources_drive( global_land_pts, priority_order, & !------------------------------------------------------------------------------ DO i = 1, nwater_use IF ( i == use_environment ) THEN - abstracted(:,i) = 0.0 + abstracted_per_use(:,i) = 0.0 ELSE - abstracted(:,i) = demand_accum(:,i) - demand_unmet(:,i) + abstracted_per_use(:,i) = demand_accum(:,i) - demand_unmet(:,i) END IF END DO !------------------------------------------------------------------------------ ! Calculate conveyance loss. !------------------------------------------------------------------------------ -CALL calc_conveyance_loss( global_land_pts, abstracted, conv_loss_frac, & +CALL calc_conveyance_loss( global_land_pts, abstracted_per_use, conv_loss_frac,& conveyance_loss, conveyance_loss_use ) -!--------------------------------------------------------------------------- +!------------------------------------------------------------------------------ ! Calculate return flows. !------------------------------------------------------------------------------ -CALL calc_return_flow( global_land_pts, abstracted, conveyance_loss_use, & - return_flow_gw, return_flow_sw, water_removed ) +CALL calc_return_flow( global_land_pts, abstracted_per_use, & + conveyance_loss_use, return_flow_gw, return_flow_river, & + water_removed ) !------------------------------------------------------------------------------ ! Calculate the water supplied for irrigation. @@ -240,7 +241,7 @@ SUBROUTINE water_resources_drive( global_land_pts, priority_order, & IF ( l_water_irrigation ) THEN DO l = 1, global_land_pts ! Water supplied is the water abstracted, minus conveyance loss. - supply_irrig(l) = abstracted(l,use_irrigation) & + supply_irrig(l) = abstracted_per_use(l,use_irrigation) & - conveyance_loss_use(l,use_irrigation) END DO END IF @@ -291,6 +292,16 @@ SUBROUTINE calc_target_split( global_land_pts, gw_avail, sw_avail, & sfc_water_frac(global_land_pts) ! Fraction of demand to be met from surface water. +!------------------------------------------------------------------------------ +! Local parameters. +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), PARAMETER :: & + water_min = 1.0e-10 + ! A minimum amount of water below which calculations are not performed + ! (kg). + +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'CALC_TARGET_SPLIT' + !------------------------------------------------------------------------------ ! Local scalar variables. !------------------------------------------------------------------------------ @@ -298,25 +309,12 @@ SUBROUTINE calc_target_split( global_land_pts, gw_avail, sw_avail, & l ! Loop counter. -!------------------------------------------------------------------------------ -! Local scalar variables. -!------------------------------------------------------------------------------ REAL(KIND=real_jlslsm) :: & denom, & ! Denominator (kg). surface_water ! Available surface water, multiplied by weighting factor (kg). -!------------------------------------------------------------------------------ -! Local parameters. -!------------------------------------------------------------------------------ -REAL(KIND=real_jlslsm), PARAMETER :: & - water_min = 1.0e-10 - ! A minimum amount of water below which calculations are not performed - ! (kg). - -CHARACTER(LEN=*), PARAMETER :: RoutineName = 'CALC_TARGET_SPLIT' - ! Dr Hook variables INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -386,6 +384,11 @@ SUBROUTINE split_demands( global_land_pts, demand_accum, sfc_water_frac, & demand_gw(global_land_pts,nwater_use) ! Demand for water from groundwater, for each water use (kg). +!------------------------------------------------------------------------------ +! Local parameters. +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'SPLIT_DEMANDS' + !------------------------------------------------------------------------------ ! Local scalar variables. !------------------------------------------------------------------------------ @@ -400,11 +403,6 @@ SUBROUTINE split_demands( global_land_pts, demand_accum, sfc_water_frac, & sw_frac(global_land_pts,nwater_use) ! Surface water fraction for each water use. -!------------------------------------------------------------------------------ -! Local parameters. -!------------------------------------------------------------------------------ -CHARACTER(LEN=*), PARAMETER :: RoutineName = 'SPLIT_DEMANDS' - ! Dr Hook variables INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -474,11 +472,12 @@ SUBROUTINE split_demands( global_land_pts, demand_accum, sfc_water_frac, & RETURN END SUBROUTINE split_demands -!############################################################################# -!############################################################################# +!############################################################################## +!############################################################################## -SUBROUTINE calc_conveyance_loss( global_land_pts, abstracted, conv_loss_frac, & - conveyance_loss, conveyance_loss_use ) +SUBROUTINE calc_conveyance_loss( global_land_pts, abstracted_per_use, & + conv_loss_frac, conveyance_loss, & + conveyance_loss_use ) !------------------------------------------------------------------------------ ! Description: @@ -504,7 +503,7 @@ SUBROUTINE calc_conveyance_loss( global_land_pts, abstracted, conv_loss_frac, & ! Array arguments with INTENT(IN) !------------------------------------------------------------------------------ REAL(KIND=real_jlslsm), INTENT(IN) :: & - abstracted(global_land_pts,nwater_use), & + abstracted_per_use(global_land_pts,nwater_use), & ! Water abstracted to meet each use (kg). conv_loss_frac(global_land_pts) ! Fraction of water that is lost during conveyance from source to user. @@ -518,6 +517,11 @@ SUBROUTINE calc_conveyance_loss( global_land_pts, abstracted, conv_loss_frac, & conveyance_loss_use(global_land_pts,nwater_use) ! Water that is lost during conveyance, for each water use (kg). +!------------------------------------------------------------------------------ +! Local parameters. +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'CALC_CONVEYANCE_LOSS' + !------------------------------------------------------------------------------ ! Local scalar variables. !------------------------------------------------------------------------------ @@ -532,11 +536,6 @@ SUBROUTINE calc_conveyance_loss( global_land_pts, abstracted, conv_loss_frac, & loss_frac(global_land_pts) ! Fraction of water that is lost during conveyance. -!------------------------------------------------------------------------------ -! Local parameters. -!------------------------------------------------------------------------------ -CHARACTER(LEN=*), PARAMETER :: RoutineName = 'CALC_CONVEYANCE_LOSS' - ! Dr Hook variables INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -568,7 +567,7 @@ SUBROUTINE calc_conveyance_loss( global_land_pts, abstracted, conv_loss_frac, & ! Calculate conveyance loss for this use and add to total. !--------------------------------------------------------------------------- DO l = 1, global_land_pts - conveyance_loss_use(l,i) = loss_frac(l) * abstracted(l,i) + conveyance_loss_use(l,i) = loss_frac(l) * abstracted_per_use(l,i) conveyance_loss(l) = conveyance_loss(l) + conveyance_loss_use(l,i) END DO @@ -578,22 +577,21 @@ SUBROUTINE calc_conveyance_loss( global_land_pts, abstracted, conv_loss_frac, & RETURN END SUBROUTINE calc_conveyance_loss -!############################################################################# -!############################################################################# +!############################################################################## +!############################################################################## -SUBROUTINE calc_return_flow( global_land_pts, abstracted, conveyance_loss_use, & - return_flow_gw, return_flow_sw, water_removed ) +SUBROUTINE calc_return_flow( global_land_pts, abstracted_per_use, & + conveyance_loss_use, return_flow_gw, & + return_flow_river, water_removed ) !------------------------------------------------------------------------------ ! Description: -! Calculate return flows (i.e. water that is returned after use) and water -! that is removed from the system during use. +! Partition the water delivered for each use into a return flow and water +! removed through that use. !------------------------------------------------------------------------------ -USE jules_rivers_mod, ONLY: l_rivers - USE jules_water_resources_mod, ONLY: & - l_have_groundwater, l_have_surface_water, nwater_use, rf_domestic, & + nwater_use, rf_domestic, & rf_livestock, rf_industry, use_environment, use_domestic, use_industry, & use_irrigation, use_livestock, use_transfers @@ -610,7 +608,7 @@ SUBROUTINE calc_return_flow( global_land_pts, abstracted, conveyance_loss_use, & ! Array arguments with INTENT(IN) !------------------------------------------------------------------------------ REAL(KIND=real_jlslsm), INTENT(IN) :: & - abstracted(global_land_pts,nwater_use), & + abstracted_per_use(global_land_pts,nwater_use), & ! Water abstracted to meet each use (kg). conveyance_loss_use(global_land_pts,nwater_use) ! Water that is lost during conveyance, for each water use (kg). @@ -621,14 +619,17 @@ SUBROUTINE calc_return_flow( global_land_pts, abstracted, conveyance_loss_use, & REAL(KIND=real_jlslsm), INTENT(OUT) :: & return_flow_gw(global_land_pts), & ! Water that is returned to renewable groundwater after use (kg). - ! If there is no renewable groundwater, this water is later added to - ! the runoff flux. - return_flow_sw(global_land_pts), & + return_flow_river(global_land_pts), & ! Water that is returned to rivers after use (kg). water_removed(global_land_pts) ! Water that is removed from the system during use, e.g. incorporated into ! manufactured goods (kg). +!------------------------------------------------------------------------------ +! Local parameters. +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'CALC_RETURN_FLOW' + !------------------------------------------------------------------------------ ! Local scalar variables. !------------------------------------------------------------------------------ @@ -645,14 +646,9 @@ SUBROUTINE calc_return_flow( global_land_pts, abstracted, conveyance_loss_use, & return_flow_frac_gw, & ! Fraction of the water delivered that is then returned to renewable ! groundwater. - return_flow_frac_sw + return_flow_frac_river ! Fraction of the water delivered that is then returned to rivers. -!------------------------------------------------------------------------------ -! Local parameters. -!------------------------------------------------------------------------------ -CHARACTER(LEN=*), PARAMETER :: RoutineName = 'CALC_RETURN_FLOW' - ! Dr Hook variables INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 @@ -663,50 +659,32 @@ SUBROUTINE calc_return_flow( global_land_pts, abstracted, conveyance_loss_use, & IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) ! Initialise return flows and water removed. -return_flow_gw(:) = 0.0 -return_flow_sw(:) = 0.0 -water_removed(:) = 0.0 +return_flow_gw(:) = 0.0 +return_flow_river(:) = 0.0 +water_removed(:) = 0.0 DO i = 1, nwater_use !--------------------------------------------------------------------------- - ! Set fraction of flow that is returned for this water use. Each use returns - ! water to either renewable groundwater or river water, but if that sink is - ! not modelled the other is used. We know that at least one of l_rivers and - ! l_have_groundwater is always TRUE when water resources are modelled. If - ! the only groundwater is "non-renewable", any groundwater return is instead - ! later added to runoff (not groundwater). + ! Set the fraction of flow that is returned for this water use. Each use + ! returns water to either groundwater or rivers. The complementary fraction + ! is removed from the system (e.g. by being incorporated into a product). !--------------------------------------------------------------------------- IF ( i == use_domestic ) THEN - ! Domestic water is returned to rivers, if those are modelled. - IF ( l_rivers ) THEN - return_flow_frac_gw = 0.0 - return_flow_frac_sw = rf_domestic - ELSE - return_flow_frac_gw = rf_domestic - return_flow_frac_sw = 0.0 - END IF - removed_frac = 1.0 - rf_domestic + ! Domestic water is returned to rivers. + return_flow_frac_gw = 0.0 + return_flow_frac_river = rf_domestic + removed_frac = 1.0 - rf_domestic ELSE IF ( i == use_industry ) THEN - ! Industrial water is returned to rivers, if those are modelled. - IF ( l_rivers ) THEN - return_flow_frac_gw = 0.0 - return_flow_frac_sw = rf_industry - ELSE - return_flow_frac_gw = rf_industry - return_flow_frac_sw = 0.0 - END IF - removed_frac = 1.0 - rf_industry + ! Industrial water is returned to rivers. + return_flow_frac_gw = 0.0 + return_flow_frac_river = rf_industry + removed_frac = 1.0 - rf_industry ELSE IF ( i == use_livestock ) THEN - ! Livestock water is returned to groundwater, if that is modelled. - IF ( l_have_groundwater ) THEN - return_flow_frac_gw = rf_livestock - return_flow_frac_sw = 0.0 - ELSE - return_flow_frac_gw = 0.0 - return_flow_frac_sw = rf_livestock - END IF - removed_frac = 1.0 - rf_livestock + ! Livestock water is returned to renewable groundwater. + return_flow_frac_gw = rf_livestock + return_flow_frac_river = 0.0 + removed_frac = 1.0 - rf_livestock ELSE ! For all other uses no water is returned and none is removed (in that the ! water remains in the system and is accounted for). There is nothing more @@ -721,12 +699,12 @@ SUBROUTINE calc_return_flow( global_land_pts, abstracted, conveyance_loss_use, & ! Calculate the water delivered to meet this demand in this gridbox. ! This is the water abstracted, minus conveyance loss. - delivered = abstracted(l,i) - conveyance_loss_use(l,i) + delivered = abstracted_per_use(l,i) - conveyance_loss_use(l,i) - ! Add to totals. Both GW and SW return flows are always calculated, though - ! either can be zero if those sources are not represented. + ! Add to totals return_flow_gw(l) = return_flow_gw(l) + return_flow_frac_gw * delivered - return_flow_sw(l) = return_flow_sw(l) + return_flow_frac_sw * delivered + return_flow_river(l) = return_flow_river(l) & + + return_flow_frac_river * delivered water_removed(l) = water_removed(l) + removed_frac * delivered END DO From 410fe2deac789088f1e5dcd7abeb0b0385788ec8 Mon Sep 17 00:00:00 2001 From: Helen Baron Date: Wed, 26 Aug 2026 14:22:42 +0100 Subject: [PATCH 3/4] Include reservoir routing --- CONTRIBUTORS.md | 1 + src/control/rivers-standalone/control_mod.F90 | 23 +- src/control/shared/surf_couple_extra_mod.F90 | 41 ++- src/control/shared/surf_couple_rivers_mod.F90 | 307 +++++++++++++----- .../river_routing/rivers_route_mod.F90 | 124 ++++++- .../river_routing/rivers_route_rfm_mod.F90 | 31 +- .../river_routing/rivers_route_trip_mod.F90 | 56 +++- .../river_routing/route_reservoirs_mod.F90 | 142 ++++++++ 8 files changed, 615 insertions(+), 110 deletions(-) create mode 100644 src/science/river_routing/route_reservoirs_mod.F90 diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index db1cc110..c15453c3 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -5,3 +5,4 @@ | yaswant | Yaswant Pradhan | Met Office | 2026-06-30 | | doucla | Douglas Clark | UKCEH | 2026-07-03 | | james-bruten-mo | James Bruten | Met Office | 2026-07-09 | +| HEBaron | Helen Baron | UKCEH | 2026-08-26 | diff --git a/src/control/rivers-standalone/control_mod.F90 b/src/control/rivers-standalone/control_mod.F90 index 14d32fa3..cc94f3c4 100644 --- a/src/control/rivers-standalone/control_mod.F90 +++ b/src/control/rivers-standalone/control_mod.F90 @@ -121,8 +121,14 @@ SUBROUTINE control ( & flowin(row_length, rows), & bflowin(row_length, rows), & twatstor(river_row_length, river_rows), & - net_abstracted_river(land_pts) - ! Net abstraction from rivers (kg m-2). + net_abstracted_river_global(land_pts), & + ! Net abstraction from rivers (kg m-2) + abstracted_res_global(land_pts), & + ! Water abstracted from reservoirs (kg). + tot_net_abstracted_river_global(land_pts), & + ! Water abstracted from rivers over river timestep (kg m-2). + tot_abstracted_res_global(land_pts) + ! Water abstracted from reservoirs over river timestep (kg). REAL(KIND=real_jlslsm) :: & !Passed between river routing and UM diagnostics_riv only @@ -137,8 +143,11 @@ SUBROUTINE control ( & !------------------------------------------------------------------------------ !End of header -! Abstraction of water is zero in rivers-only configuration. -net_abstracted_river(:) = 0.0 +! Abstractions of water are zero in rivers-only configuration. +abstracted_res_global(:) = 0.0 +net_abstracted_river_global(:) = 0.0 +tot_abstracted_res_global(:) = 0.0 +tot_net_abstracted_river_global(:) = 0.0 SELECT CASE ( lsm_id ) CASE ( rivers ) @@ -146,13 +155,15 @@ SUBROUTINE control ( & !INTEGER, INTENT(IN) land_pts, n_wtrac_jls, & !REAL, INTENT(IN) - net_abstracted_river, fluxes%sub_surf_roff_gb, fluxes%surf_roff_gb, & + abstracted_res_global, net_abstracted_river_global, & + fluxes%sub_surf_roff_gb, fluxes%surf_roff_gb, & wtrac_jls%sub_surf_roff_gb, wtrac_jls%surf_roff_gb, & !INTEGER, INTENT(INOUT) a_steps_since_riv, & !REAL, INTENT (INOUT) river%tot_surf_runoff_gb, river%tot_sub_runoff_gb, & - river%acc_lake_evap_gb, & + tot_abstracted_res_global, & + tot_net_abstracted_river_global, river%acc_lake_evap_gb, & wtrac_jls%tot_surf_runoff_gb, wtrac_jls%tot_sub_runoff_gb, & wtrac_jls%acc_lake_evap_gb, & !REAL, INTENT (OUT) diff --git a/src/control/shared/surf_couple_extra_mod.F90 b/src/control/shared/surf_couple_extra_mod.F90 index 4b52d3c1..3b9094ec 100644 --- a/src/control/shared/surf_couple_extra_mod.F90 +++ b/src/control/shared/surf_couple_extra_mod.F90 @@ -79,7 +79,8 @@ SUBROUTINE surf_couple_extra( & zw_soilt, sthzw_soilt, & ls_rainfrac_gb, & substore, surfstore, flowin, bflowin, & - tot_surf_runoff_gb, tot_sub_runoff_gb, acc_lake_evap_gb, twatstor, & + tot_surf_runoff_gb, tot_sub_runoff_gb, & + acc_lake_evap_gb, twatstor, & asteps_since_triffid, & inlandout_atm_gb, & !OUT @@ -122,7 +123,7 @@ SUBROUTINE surf_couple_extra( & USE water_resources_vars_mod, ONLY: water_resources_type USE jules_wtrac_type_mod, ONLY: jls_wtrac_type -! In general CABLE utilizes a required subset of tbe JULES types, however; +! In general CABLE utilizes a required subset of the JULES types, however; USE work_vars_mod_cbl, ONLY: work_vars_type ! and some kept thru timestep !Import interfaces to subroutines called @@ -248,7 +249,9 @@ SUBROUTINE surf_couple_extra( & #if defined(UM_JULES) && !defined(LFRIC) USE atm_fields_mod, ONLY: disturb_veg_prev -USE atm_step_local, ONLY: STASHwork19, STASHwork8, STASHwork26 +USE atm_land_sea_mask, ONLY: global_land_pts => atmos_number_of_landpts + +USE atm_step_local, ONLY: STASHwork19, STASHwork8, STASHwork26 USE model_domain_mod, ONLY: model_type, mt_single_column @@ -263,7 +266,7 @@ SUBROUTINE surf_couple_extra( & USE metstats_mod, ONLY: metstats_prog, l_metstats -USE model_grid_mod, ONLY: grid_area_ij +USE model_grid_mod, ONLY: global_land_pts, grid_area_ij USE model_time_mod, ONLY: current_time #endif @@ -841,7 +844,7 @@ SUBROUTINE surf_couple_extra( & #if !defined(UM_JULES) ! Water resources (standalone; not yet allowed in UM). IF ( l_water_resources ) THEN - CALL water_resources_control( & + CALL water_resources_control( global_land_pts, & rivers%global_land_index, ainfo%land_index, & rivers%map_river_to_land_points, rivers%rivers_index_rp, & forcing%con_rain_ij, forcing%con_snow_ij, & @@ -853,24 +856,29 @@ SUBROUTINE surf_couple_extra( & flandg, crop_vars%frac_irr_soilt, ainfo%frac_soilt, & ainfo%frac_surft, grid_area_ij, & forcing%ls_rain_ij, forcing%ls_snow_ij, forcing%lw_down_ij, & + rivers%res_storage, & + rivers%rfm_surfstore_rp, rivers%rivers_sto_rp, & psparms%smvccl_soilt, & psparms%smvcst_soilt, psparms%smvcwt_soilt, psparms%sthf_soilt, & fluxes%sw_surft, forcing%tl_1_ij, progs%tstar_surft, & crop_vars%icntmax_gb, crop_vars%plant_n_gb, & water_resources%demand_accum, & crop_vars%prec_1_day_av_gb, crop_vars%prec_1_day_av_use_gb, & - rivers%rfm_surfstore_rp, rivers%rivers_sto_rp, & crop_vars%rn_1_day_av_gb, & crop_vars%rn_1_day_av_use_gb, water_resources%sfc_water_frac, & progs%smcl_soilt, & crop_vars%sthu_irr_soilt, psparms%sthu_soilt, & sthzw_soilt, fluxes%sub_surf_roff_gb, crop_vars%tl_1_day_av_gb, & - crop_vars%tl_1_day_av_use_gb, water_resources%priority_order, & - water_resources%demand_unmet, water_resources%gw_abstracted, & - water_resources%gw_avail, water_resources%gw_nr_abstracted, & - crop_vars%irrig_water_gb, water_resources%net_abstracted_river, & - water_resources%sw_abstracted, water_resources%sw_avail_total, & - water_resources%water_removed ) + crop_vars%tl_1_day_av_use_gb, & + water_resources%priority_order, & + water_resources%abstracted_res_global, & + water_resources%net_abstracted_river_global, & + water_resources%gw_abstracted, crop_vars%irrig_water_gb, & + water_resources%abstracted_res, & + water_resources%abstracted_river, water_resources%conveyance_loss, & + water_resources%demand_unmet, water_resources%gw_avail, & + water_resources%gw_nr_abstracted, water_resources%sw_abstracted, & + water_resources%sw_avail_total, water_resources%water_removed ) END IF #endif @@ -879,13 +887,16 @@ SUBROUTINE surf_couple_extra( & !INTEGER, INTENT(IN) land_pts, n_wtrac_jls, & !REAL, INTENT(IN) - water_resources%net_abstracted_river, & + water_resources%abstracted_res_global, & + water_resources%net_abstracted_river_global, & fluxes%sub_surf_roff_gb, fluxes%surf_roff_gb, & wtrac_jls%sub_surf_roff_gb, wtrac_jls%surf_roff_gb, & !INTEGER, INTENT(INOUT) a_steps_since_riv, & !REAL, INTENT (INOUT) - tot_surf_runoff_gb, tot_sub_runoff_gb, acc_lake_evap_gb, & + tot_surf_runoff_gb, tot_sub_runoff_gb, & + rivers%tot_abstracted_res_global, & + rivers%tot_net_abstracted_river_global, acc_lake_evap_gb, & wtrac_jls%tot_surf_runoff_gb, wtrac_jls%tot_sub_runoff_gb, & wtrac_jls%acc_lake_evap_gb, & !REAL, INTENT (OUT) @@ -1420,4 +1431,4 @@ SUBROUTINE surf_couple_extra( & IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE surf_couple_extra -END MODULE surf_couple_extra_mod +END MODULE surf_couple_extra_mod \ No newline at end of file diff --git a/src/control/shared/surf_couple_rivers_mod.F90 b/src/control/shared/surf_couple_rivers_mod.F90 index 56f5740b..8893df68 100644 --- a/src/control/shared/surf_couple_rivers_mod.F90 +++ b/src/control/shared/surf_couple_rivers_mod.F90 @@ -18,12 +18,15 @@ SUBROUTINE surf_couple_rivers( & !INTEGER, INTENT(IN) land_pts, n_wtrac_jls, & !REAL, INTENT(IN) - net_abstracted_river, & + abstracted_res_global, & + net_abstracted_river_global, & sub_surf_roff, surf_roff, sub_surf_roff_wtrac, surf_roff_wtrac, & !INTEGER, INTENT(INOUT) a_steps_since_riv, & !REAL, INTENT (INOUT) - tot_surf_runoff_gb, tot_sub_runoff_gb, acc_lake_evap_gb, & + tot_surf_runoff_gb, tot_sub_runoff_gb, & + tot_abstracted_res_global, tot_net_abstracted_river_global, & + acc_lake_evap_gb, & tot_surf_runoff_gb_wtrac, tot_sub_runoff_gb_wtrac, & acc_lake_evap_gb_wtrac, & !REAL, INTENT (OUT) @@ -56,11 +59,13 @@ SUBROUTINE surf_couple_rivers( & USE jules_rivers_mod, ONLY: i_river_vn, nstep_rivers, & rivers_camaflood, rivers_rfm, & rivers_trip, rivers_call, np_rivers, & - rivers_type, l_riv_overbank, & + rivers_type, & + l_riv_overbank, & l_vary_sea_level, & ! UM only l_inland_outflow, rivers_um_trip, & - rivers_first + rivers_first, & + l_reservoirs USE jules_model_environment_mod, ONLY: l_oasis_rivers USE timestep_mod, ONLY: timestep @@ -80,6 +85,8 @@ SUBROUTINE surf_couple_rivers( & USE rivers_regrid_mod, ONLY: landpts_to_rivpts, rivpts_to_landpts +USE jules_water_resources_mod, ONLY: l_water_resources, sw_river_source + USE jules_water_tracers_mod, ONLY: l_wtrac_jls !Module imports - Variables required only in UM-mode @@ -145,7 +152,9 @@ SUBROUTINE surf_couple_rivers( & a_steps_since_riv REAL(KIND=real_jlslsm), INTENT(IN) :: & - net_abstracted_river(land_pts), & + abstracted_res_global(global_land_pts), & + ! Water abstracted from reservoirs (kg). + net_abstracted_river_global(global_land_pts), & ! Net abstraction from rivers (kg m-2). surf_roff(land_pts), & ! Surface runoff (kg m-2 s-1) @@ -185,6 +194,10 @@ SUBROUTINE surf_couple_rivers( & ! Average rate of surface runoff over river timestep (kg m-2 s-1). tot_sub_runoff_gb(land_pts), & ! Average rate of subsurface runoff over river timestep (kg m-2 s-1). + tot_abstracted_res_global(global_land_pts), & + ! Water abstracted from reservoirs over river timestep (kg). + tot_net_abstracted_river_global(global_land_pts), & + ! Water abstracted from rivers over river timestep (kg m-2). acc_lake_evap_gb(row_length,rows), & tot_surf_runoff_gb_wtrac(land_pts,n_wtrac_jls), & tot_sub_runoff_gb_wtrac(land_pts,n_wtrac_jls), & @@ -235,6 +248,11 @@ SUBROUTINE surf_couple_rivers( & REAL(KIND=real_jlslsm), ALLOCATABLE :: global_rrun(:) REAL(KIND=real_jlslsm), ALLOCATABLE :: global_rflow(:) REAL(KIND=real_jlslsm), ALLOCATABLE :: global_sea_level(:) +REAL(KIND=real_jlslsm), ALLOCATABLE :: abstracted_res_rp(:) + ! Water abstracted from reservoirs over river timestep, on river points + ! (kg). +REAL(KIND=real_jlslsm), ALLOCATABLE :: net_abstracted_river_rp(:) + ! Water abstracted from rivers over river timestep, on river points (kg m-2). #if defined(UM_JULES) REAL(KIND=real_jlslsm) :: & @@ -244,17 +262,17 @@ SUBROUTINE surf_couple_rivers( & REAL(KIND=real_jlslsm), ALLOCATABLE :: inlandout_atmos_wtrac(:,:,:) ! inlandout_atmos_gb_wtrac on (i,j) grid -REAL(KIND=real_jlslsm), ALLOCATABLE :: net_abstracted_river_wtrac(:,:) - ! Net abstraction of water tracers from rivers - ! (kg m-2). !Local variables INTEGER :: & ERROR, & - ! Error status from each call to ALLOCATE. + ! Error status from each call to ALLOCATE. error_sum, & - ! Accumulated error status. - l,i,j,ip,i_wt + ! Accumulated error status. + l,i,j,ip,i_wt, & + ! Counters + tmp_land_size, tmp_river_size + ! Sizes used in allocation. #if defined(UM_JULES) REAL(KIND=real_jlslsm) :: & @@ -283,10 +301,8 @@ SUBROUTINE surf_couple_rivers( & ! Allocate water tracer fields IF (l_wtrac_jls) THEN ALLOCATE(inlandout_atmos_wtrac(row_length,rows,n_wtrac_jls)) - ALLOCATE(net_abstracted_river_wtrac(land_pts,n_wtrac_jls)) ELSE ALLOCATE(inlandout_atmos_wtrac(1,1,1)) - ALLOCATE(net_abstracted_river_wtrac(1,1)) END IF #if defined(UM_JULES) @@ -323,9 +339,11 @@ SUBROUTINE surf_couple_rivers( & #endif -!Initialise the accumulated surface and subsurface runoff to zero -!at the beginning of river routing timestep +! Initialise accumulated fluxes to zero at the beginning of river routing +! timestep. IF ( a_steps_since_riv == 0 ) THEN + + ! Initialise the accumulated surface and subsurface runoff. tot_surf_runoff_gb(:) = 0.0 tot_sub_runoff_gb(:) = 0.0 acc_lake_evap_gb(:,:) = 0.0 @@ -348,7 +366,22 @@ SUBROUTINE surf_couple_rivers( & rivers%rrun_sub_surf_rp(:) = 0.0 rivers%rrun_surf_rp(:) = 0.0 END IF -END IF + + ! These "global" arrays are only allocated at full size on the master task, + ! so they must only be accessed there. + IF ( is_master_task() ) THEN + ! Initialise the accumulated abstraction from rivers. + IF ( sw_river_source > 0 ) THEN + tot_net_abstracted_river_global(:) = 0.0 + END IF + + ! Initialise the accumulated abstraction from reservoirs. + IF ( l_water_resources .AND. l_reservoirs ) THEN + tot_abstracted_res_global(:) = 0.0 + END IF + END IF + +END IF ! a_steps_since_riv == 0 ! Increment counters. a_steps_since_riv = a_steps_since_riv + 1 @@ -359,33 +392,48 @@ SUBROUTINE surf_couple_rivers( & rivers_call = .FALSE. END IF -!Accumulate the runoff as Kg/m2/s over the River Routing period (that is, -! between calls to river routing). +! Accumulate input fluxes over the river routing period (that is, between calls +! to river routing). IF (l_oasis_rivers) THEN - CALL accumulate_runoff(np_rivers, net_abstracted_river, & + ! Accumulate runoff (on river points). + CALL accumulate_runoff(np_rivers, & rivers%surf_roff_rp, rivers%sub_surf_roff_rp, & rivers%rrun_surf_rp, rivers%rrun_sub_surf_rp) ELSE - CALL accumulate_runoff(land_pts, net_abstracted_river, & + ! Accumulate runoff (on land points). + CALL accumulate_runoff(land_pts, & surf_roff, sub_surf_roff, & tot_surf_runoff_gb, tot_sub_runoff_gb) ! Repeat for water tracers IF (l_wtrac_jls) THEN - ! Water tracers are not yet represented in the water resources code, hence - ! set net abstraction of tracers to zero. - net_abstracted_river_wtrac(:,:) = 0.0 DO i_wt = 1, n_wtrac_jls - CALL accumulate_runoff(land_pts, net_abstracted_river_wtrac(:,i_wt), & + CALL accumulate_runoff(land_pts, & surf_roff_wtrac(:,i_wt), sub_surf_roff_wtrac(:,i_wt), & tot_surf_runoff_gb_wtrac(:,i_wt), & tot_sub_runoff_gb_wtrac(:,i_wt)) END DO END IF ! l_wtrac_jls + ! These "global" arrays are only allocated at full size on the master task, + ! so they must only be accessed there. + IF ( is_master_task() ) THEN + ! Accumulate net abstraction from rivers. + IF ( sw_river_source > 0 ) THEN + CALL accumulate_abstraction( global_land_pts, net_abstracted_river_global, & + tot_net_abstracted_river_global ) + END IF + + ! Accumulate abstraction from reservoirs. + IF ( l_water_resources .AND. l_reservoirs ) THEN + CALL accumulate_abstraction( global_land_pts, abstracted_res_global, & + tot_abstracted_res_global ) + END IF + END IF + END IF ! l_oasis_rivers) ! Could this UM ifdef be replaced with the i_rivers_vn switch? @@ -501,54 +549,86 @@ SUBROUTINE surf_couple_rivers( & END IF #endif !-------------------------------------------------------------------------- - ! Gather runoff information from all processors + ! Allocate coupling arrays. + ! These are global land point arrays or river point arrays. !-------------------------------------------------------------------------- - + ! Allocate runoff arrays. IF (.NOT. l_oasis_rivers .AND. is_master_task()) THEN - ALLOCATE(global_tot_sub_runoff(global_land_pts), STAT = ERROR) - error_sum = ERROR - ALLOCATE(global_tot_surf_runoff(global_land_pts), STAT = ERROR) - error_sum = error_sum + ERROR - ALLOCATE(global_rrun(global_land_pts), STAT = ERROR) - error_sum = error_sum + ERROR - ALLOCATE(global_rflow(global_land_pts), STAT = ERROR) - error_sum = error_sum + ERROR + tmp_land_size = global_land_pts ELSE - ALLOCATE(global_tot_sub_runoff(1), STAT = ERROR) - error_sum = ERROR - ALLOCATE(global_tot_surf_runoff(1), STAT = ERROR) - error_sum = error_sum + ERROR - ALLOCATE(global_rrun(1), STAT = ERROR) - error_sum = error_sum + ERROR - ALLOCATE(global_rflow(1), STAT = ERROR) - error_sum = error_sum + ERROR + ! Allocate at minimum size. + tmp_land_size = 1 END IF - + ALLOCATE(global_tot_sub_runoff(tmp_land_size), STAT = ERROR) + error_sum = ERROR + ALLOCATE(global_tot_surf_runoff(tmp_land_size), STAT = ERROR) + error_sum = error_sum + ERROR + ALLOCATE(global_rrun(tmp_land_size), STAT = ERROR) + error_sum = error_sum + ERROR + ALLOCATE(global_rflow(tmp_land_size), STAT = ERROR) + error_sum = error_sum + ERROR + + ! Allocate sea level array. IF ( .NOT. l_oasis_rivers .AND. l_vary_sea_level .AND. is_master_task() ) & THEN - ALLOCATE(global_sea_level(global_land_pts), STAT = ERROR) - error_sum = error_sum + ERROR + tmp_land_size = global_land_pts + ELSE + tmp_land_size = 1 + END IF + ALLOCATE(global_sea_level(tmp_land_size), STAT = ERROR) + error_sum = error_sum + ERROR + + ! Allocate abstraction from rivers. + IF ( .NOT. l_oasis_rivers .AND. sw_river_source > 0 & + .AND. is_master_task() ) THEN + ! Water resource code will calclate abstraction. + tmp_river_size = np_rivers ELSE - ALLOCATE(global_sea_level(1), STAT = ERROR) - error_sum = error_sum + ERROR + ! Allocate at minimum size. + tmp_river_size = 1 END IF + ALLOCATE(net_abstracted_river_rp(tmp_river_size), STAT = ERROR) + error_sum = error_sum + ERROR + + ! Allocate global abstraction from reservoirs. + IF ( .NOT. l_oasis_rivers .AND. l_reservoirs & + .AND. is_master_task() ) THEN + tmp_river_size = np_rivers + ELSE + ! Allocate at minimum size. + tmp_river_size = 1 + END IF + ALLOCATE(abstracted_res_rp(tmp_river_size), STAT = ERROR) + error_sum = error_sum + ERROR IF ( error_sum /= 0 ) THEN errorstatus = 10 CALL ereport( RoutineName, errorstatus, & - "Error related to allocation of runoff variables." ) + "Error related to allocation of abstraction variables." ) + ELSE + ! Initialise abstraction from reservoirs to zero - this value is + ! passed to subroutine rivers_route_rp and used if l_reservoirs=T but + ! l_water_resources=F. + abstracted_res_rp(:) = 0.0 + ! Initialise abstraction from rivers. + net_abstracted_river_rp(:) = 0.0 END IF - IF (.NOT. l_oasis_rivers) THEN + !-------------------------------------------------------------------------- + ! Gather fluxes from all processors + !-------------------------------------------------------------------------- + IF ( .NOT. l_oasis_rivers ) THEN + ! Gather runoff. CALL gather_land_field(tot_sub_runoff_gb, global_tot_sub_runoff, & rivers%global_land_index) CALL gather_land_field(tot_surf_runoff_gb, global_tot_surf_runoff, & rivers%global_land_index) + ! Gather sea level. IF ( l_vary_sea_level ) THEN CALL gather_land_field(rivers%sea_level_lp, global_sea_level, & rivers%global_land_index ) END IF - END IF + END IF ! .NOT. l_oasis_rivers ) THEN !------------------------------------------------------------------------- ! Call routing driver on single processor @@ -556,18 +636,24 @@ SUBROUTINE surf_couple_rivers( & IF ( is_master_task() ) THEN IF (l_oasis_rivers) THEN + ! No regridding is required between the land and river grids as this ! is done by the OASIS coupler. Only the science routine needs calling. - CALL rivers_route_rp( rivers ) + CALL rivers_route_rp( abstracted_res_rp, & + net_abstracted_river_rp, rivers ) + ELSE + ! Initialisation DO l = 1, global_land_pts global_rflow(l)= 0.0 global_rrun(l) = 0.0 END DO - ! Regrid surface and subsurface runoff from land points to rivers - ! points. + !---------------------------------------------------------------------- + ! Regrid coupling fields from land points to rivers. + !---------------------------------------------------------------------- + ! Regrid runoff. CALL landpts_to_rivpts( global_land_pts, np_rivers, & rivers%map_river_to_land_points, & rivers%global_land_index, & @@ -581,8 +667,8 @@ SUBROUTINE surf_couple_rivers( & rivers%rivers_index_rp, & global_tot_surf_runoff, rivers%rrun_surf_rp ) + ! Regrid sea level. IF ( l_vary_sea_level ) THEN - ! Regrid sea level from land points to river points. CALL landpts_to_rivpts( global_land_pts, np_rivers, & rivers%map_river_to_land_points, & rivers%global_land_index, & @@ -590,10 +676,31 @@ SUBROUTINE surf_couple_rivers( & global_sea_level, rivers%sea_level ) END IF + ! Regrid abstraction from rivers. + IF ( sw_river_source > 0 ) THEN + CALL landpts_to_rivpts( global_land_pts, np_rivers, & + rivers%map_river_to_land_points, & + rivers%global_land_index, & + rivers%rivers_index_rp, & + tot_net_abstracted_river_global, & + net_abstracted_river_rp ) + END IF + + ! Regrid abstraction from reservoirs. + IF ( l_reservoirs .AND. l_water_resources ) THEN + CALL landpts_to_rivpts( global_land_pts, np_rivers, & + rivers%map_river_to_land_points, & + rivers%global_land_index, & + rivers%rivers_index_rp, & + tot_abstracted_res_global, & + abstracted_res_rp ) + END IF + ! Call the routing science routine. - CALL rivers_route_rp( rivers ) + CALL rivers_route_rp( abstracted_res_rp, & + net_abstracted_river_rp, rivers ) - ! Regrid outputs from rivers to land grid + ! Regrid selected outputs from rivers to land grid. CALL rivpts_to_landpts( global_land_pts, np_rivers, & rivers%map_river_to_land_points, & rivers%global_land_index, & @@ -605,7 +712,8 @@ SUBROUTINE surf_couple_rivers( & rivers%global_land_index, & rivers%rivers_index_rp, & rivers%rrun_rp, global_rrun ) - END IF + + END IF ! l_oasis_rivers !----------------------------------------------------------------------- ! Compute overbank inundation @@ -613,6 +721,7 @@ SUBROUTINE surf_couple_rivers( & IF ( l_riv_overbank ) THEN CALL overbank_update(rivers%rfm_rivflow_rp) END IF + END IF ! end is_master !------------------------------------------------------------------------- @@ -655,6 +764,9 @@ SUBROUTINE surf_couple_rivers( & END IF ! .NOT. l_oasis_rivers #endif + ! Deallocate local variables (in reverse order to allocation). + DEALLOCATE(abstracted_res_rp) + DEALLOCATE(net_abstracted_river_rp) DEALLOCATE(global_sea_level) DEALLOCATE(global_rflow) DEALLOCATE(global_rrun) @@ -706,7 +818,6 @@ SUBROUTINE surf_couple_rivers( & END IF ! rivers_call ! Deallocate water tracer fields -DEALLOCATE(net_abstracted_river_wtrac) DEALLOCATE(inlandout_atmos_wtrac) IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) @@ -715,15 +826,13 @@ END SUBROUTINE surf_couple_rivers !--------------------------------------------------------------------------- -SUBROUTINE accumulate_runoff(npoints, net_abstracted_river, surf_runoff, & - sub_runoff, surf_runoff_accum, sub_runoff_accum) +SUBROUTINE accumulate_runoff(npoints, surf_runoff, sub_runoff, & + surf_runoff_accum, sub_runoff_accum) ! Accumulate the surface and subsurface runoff over the river routing period. -! This is a generic routine used for normal water and water tracers +! This is a generic routine used for normal water and water tracers. -USE jules_rivers_mod, ONLY: nstep_rivers -USE jules_water_resources_mod, ONLY: sw_river_source -USE timestep_mod, ONLY: timestep +USE jules_rivers_mod, ONLY: nstep_rivers USE parkind1, ONLY: jprb, jpim USE yomhook, ONLY: lhook, dr_hook @@ -733,8 +842,6 @@ SUBROUTINE accumulate_runoff(npoints, net_abstracted_river, surf_runoff, & INTEGER, INTENT(IN) :: npoints ! No. of points REAL(KIND=real_jlslsm), INTENT(IN) :: & - net_abstracted_river(npoints), & - ! Net abstraction from rivers (kg m-2). surf_runoff(npoints), & ! Surface runoff (kg m-2 s-1) sub_runoff(npoints) @@ -774,19 +881,69 @@ SUBROUTINE accumulate_runoff(npoints, net_abstracted_river, surf_runoff, & (sub_runoff(ip) / REAL(nstep_rivers)) END IF - ! Consider abstraction of river water by water resources. - IF ( sw_river_source > 0 ) THEN - ! Remove net abstraction from the accumulated surface runoff. - ! Convert units of abstraction from kg m-2 to kg m-2 s-1. - surf_runoff_accum(ip) = surf_runoff_accum(ip) - & - ( net_abstracted_river(ip) / & - (REAL(nstep_rivers) * timestep) ) - END IF - END DO IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) RETURN END SUBROUTINE accumulate_runoff +!############################################################################## + +SUBROUTINE accumulate_abstraction( npoints, abstracted, tot_abstracted ) + +! Accumulate an abstraction - essentially just adding to a sum. + +USE parkind1, ONLY: jprb, jpim +USE yomhook, ONLY: lhook, dr_hook + +IMPLICIT NONE + +!------------------------------------------------------------------------------ +! Scalar arguments with INTENT(IN) +!------------------------------------------------------------------------------ +INTEGER, INTENT(IN) :: npoints ! No. of points + +!------------------------------------------------------------------------------ +! Array arguments with INTENT(IN) +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(IN) :: & + abstracted(npoints) + ! Water abstracted (various units). + +!------------------------------------------------------------------------------ +! Array arguments with INTENT(IN OUT) +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(IN OUT) :: & + tot_abstracted(npoints) + ! Accumulated water abstracted (various units). + +!------------------------------------------------------------------------------ +! Local scalar variables +!------------------------------------------------------------------------------ +INTEGER :: ip ! Loop counter + +!------------------------------------------------------------------------------ +! Local scalar parameters +!------------------------------------------------------------------------------ +CHARACTER(LEN=*), PARAMETER :: RoutineName = 'ACCUMULATE_ABSTRACTION' + +!Dr Hook variables +INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0 +INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1 +REAL(KIND=jprb) :: zhook_handle +!end of header +!------------------------------------------------------------------------------ + +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_in,zhook_handle) + +DO ip = 1, npoints + + tot_abstracted(ip) = tot_abstracted(ip) + abstracted(ip) + +END DO + +IF (lhook) CALL dr_hook(ModuleName//':'//RoutineName,zhook_out,zhook_handle) +RETURN +END SUBROUTINE accumulate_abstraction + END MODULE surf_couple_rivers_mod diff --git a/src/science/river_routing/rivers_route_mod.F90 b/src/science/river_routing/rivers_route_mod.F90 index 9f94355b..c496a37a 100644 --- a/src/science/river_routing/rivers_route_mod.F90 +++ b/src/science/river_routing/rivers_route_mod.F90 @@ -28,7 +28,8 @@ MODULE rivers_route_mod !############################################################################## -SUBROUTINE rivers_route_rp(rivers) +SUBROUTINE rivers_route_rp( abstracted_res_rp, & + net_abstracted_river_rp, rivers ) !------------------------------------------------------------------------------ ! @@ -42,8 +43,11 @@ SUBROUTINE rivers_route_rp(rivers) USE jules_rivers_mod, ONLY: & ! imported scalars with intent(in) - i_river_vn, rivers_camaflood, rivers_rfm, rivers_trip, np_rivers, & - rivers_type + i_river_vn, l_reservoirs, np_rivers, & + nstep_rivers, rivers_camaflood, rivers_rfm, rivers_trip, rivers_type + +USE jules_water_resources_mod, ONLY: & + l_water_resources, sw_river_source USE rivers_route_camaflood_mod, ONLY: & ! imported procedures @@ -57,12 +61,28 @@ SUBROUTINE rivers_route_rp(rivers) ! imported procedures rivers_route_trip +USE timestep_mod, ONLY: & + timestep + USE jules_print_mgr, ONLY: & jules_message, & jules_print +USE time_info_mod, ONLY: current_model_time + IMPLICIT NONE +!------------------------------------------------------------------------------ +! Arguments with INTENT(IN) +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(IN) :: & + net_abstracted_river_rp(np_rivers), & + ! Water abstracted from rivers over river timestep, on river points + ! (kg m-2). + abstracted_res_rp(np_rivers) + ! Water abstracted from reservoirs over river timestep, on river + ! points (kg). + !------------------------------------------------------------------------------ ! Arguments with INTENT(IN OUT) !------------------------------------------------------------------------------ @@ -71,7 +91,14 @@ SUBROUTINE rivers_route_rp(rivers) !------------------------------------------------------------------------------ ! Local scalar variables !------------------------------------------------------------------------------ -INTEGER :: ip ! loop counter +INTEGER :: ip ! loop counter + +REAL(KIND=real_jlslsm) :: recip_timestep + ! Reciprocal of river timestep length (s-1). + +INTEGER :: year ! current model year +INTEGER, SAVE :: year_change = -1 +! year of last reservoir capacity update !------------------------------------------------------------------------------ ! Local array variables @@ -82,6 +109,21 @@ SUBROUTINE rivers_route_rp(rivers) !end of header +!------------------------------------------------------------------------------ +! Remove net abstraction from rivers from surface runoff - in effect the +! surface runoff variable becomes a more generic source/sink term for rivers. +! The resulting term can be negative. +!------------------------------------------------------------------------------ +IF ( l_water_resources .AND. sw_river_source > 0 ) THEN + ! Calculate reciprocal of timestep length. + recip_timestep = 1.0 / ( REAL(nstep_rivers) * timestep ) + DO ip = 1, np_rivers + rivers%rrun_surf_rp(ip) = rivers%rrun_surf_rp(ip) & + ! Convert abstraction from kg m-2 to kg m-2 s-1. + - net_abstracted_river_rp(ip) * recip_timestep + END DO +END IF + !------------------------------------------------------------------------------ ! Calculate total runoff diagnostic. !------------------------------------------------------------------------------ @@ -89,6 +131,19 @@ SUBROUTINE rivers_route_rp(rivers) rivers%rrun_rp(ip) = rivers%rrun_surf_rp(ip) + rivers%rrun_sub_surf_rp(ip) END DO +!------------------------------------------------------------------------------ +! Update reservoir capacities at the start and on each year change. +!------------------------------------------------------------------------------ +IF ( l_reservoirs ) THEN + CALL current_model_time(year=year) + IF ( year /= year_change ) THEN + CALL update_reservoirs( rivers%res_capacity, & + rivers%res_year, & + rivers%res_cap_current ) + year_change = year + END IF +END IF + !------------------------------------------------------------------------------ ! Call the routing science routine. !------------------------------------------------------------------------------ @@ -106,15 +161,25 @@ SUBROUTINE rivers_route_rp(rivers) CASE ( rivers_rfm ) CALL rivers_route_rfm( rivers%rrun_surf_rp, rivers%rrun_sub_surf_rp, & rivers%rflow_rp, baseflow, & + abstracted_res_rp, & ! imported river arrays rivers) CASE ( rivers_trip ) CALL rivers_route_trip( rivers%rrun_surf_rp, rivers%rrun_sub_surf_rp, & rivers%rflow_rp, baseflow, rivers%rivers_outflow_rp, & rivers%rivers_next_rp, rivers%rivers_seq_rp, & - rivers%rivers_sto_rp,rivers%rivers_boxareas_rp, & + rivers%rivers_sto_rp, rivers%rivers_boxareas_rp, & rivers%rivers_lat_rp, rivers%rivers_lon_rp, & - rivers%inland_outflow_rp, rivers%land_fraction_rp ) + rivers%inland_outflow_rp, rivers%land_fraction_rp, & + abstracted_res_rp, & + rivers%res_cap_current, & + rivers%res_catch, & + rivers%res_critical, rivers%res_flood, & + rivers%res_emergency, & + rivers%res_normal_release, & + rivers%res_flood_release, & + rivers%res_storage ) + CASE DEFAULT WRITE(jules_message,*) 'ERROR: rivers_drive: ' // & 'do not recognise i_river_vn=', i_river_vn @@ -330,4 +395,51 @@ END SUBROUTINE adjust_routestore !############################################################################## +SUBROUTINE update_reservoirs( res_capacity, res_year, & + res_cap_current ) + +USE jules_rivers_mod, ONLY: & +! imported scalars + np_rivers + +USE time_info_mod, ONLY: current_model_time + +IMPLICIT NONE + +!------------------------------------------------------------------------------ +! Array arguments with intent(in). +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(IN) :: & + res_capacity(np_rivers), & + ! Capacity of reservoirs (kg). + res_year(np_rivers) + ! Build year of reservoirs. + +!------------------------------------------------------------------------------ +! Array arguments with intent(inout). +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(INOUT) :: & + res_cap_current(np_rivers) + ! Capacity of currently active reservoirs (kg). + +!------------------------------------------------------------------------------ +! Local scalar variables. +!------------------------------------------------------------------------------ +INTEGER :: i ! Loop counter. +INTEGER :: year ! Current model year. + +!end of header +!------------------------------------------------------------------------------ +CALL current_model_time(year) + +DO i = 1, np_rivers + IF ( res_year(i) <= year ) THEN + res_cap_current(i) = res_capacity(i) + END IF +END DO + +END SUBROUTINE update_reservoirs + +!############################################################################## + END MODULE rivers_route_mod diff --git a/src/science/river_routing/rivers_route_rfm_mod.F90 b/src/science/river_routing/rivers_route_rfm_mod.F90 index 94b43b21..f333fbde 100644 --- a/src/science/river_routing/rivers_route_rfm_mod.F90 +++ b/src/science/river_routing/rivers_route_rfm_mod.F90 @@ -65,13 +65,13 @@ MODULE rivers_route_rfm_mod ! again to a flux density kg/m2/s SUBROUTINE rivers_route_rfm( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & - rivers ) + abstracted_res_rp, rivers ) USE jules_rivers_mod, ONLY: & ! imported scalars nstep_rivers, np_rivers, river_mouth & ,rivers_first, rivers_length, runoff_factor, cland, criver, cbland & - ,cbriver, retl, retr, rfm_land, rfm_river, rfm_sea, & + ,cbriver, retl, retr, rfm_land, rfm_river, rfm_sea, l_reservoirs, & ! imported type rivers_type @@ -87,6 +87,10 @@ SUBROUTINE rivers_route_rfm( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & !----------------------------------------------------------------------------- +USE route_reservoirs_mod, ONLY: & +! imported procedures + route_reservoirs + USE um_types, ONLY: real_jlslsm IMPLICIT NONE @@ -99,6 +103,8 @@ SUBROUTINE rivers_route_rfm( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & ! This includes any abstraction of water for water resources. REAL(KIND=real_jlslsm), INTENT(IN) :: sub_sfc_runoff(np_rivers) ! average rate of sub-surface runoff since last call (kg m-2 s-1) +REAL(KIND=real_jlslsm), INTENT(IN) :: abstracted_res_rp(np_rivers) + ! Water abstracted from reservoirs over river timestep (kg). REAL(KIND=real_jlslsm), INTENT(OUT) :: outflow(np_rivers) ! rate of channel surface flow leaving gridbox (kg m-2 s-1) @@ -125,8 +131,10 @@ SUBROUTINE rivers_route_rfm( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & ! initial river flow [m3/s] ,dt & ! river routing model timestep (s) - ,dx + ,dx & ! distance between midpoints of neighbouring cells (m) + ,reservoir_flow + ! flow in and out of reservoir (kg s-1) REAL(KIND=real_jlslsm) :: & substore_n(np_rivers) & @@ -237,6 +245,23 @@ SUBROUTINE rivers_route_rfm( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & rn = rivers%rivers_next_rp(ip) landtype = rivers%rfm_land_rp(ip) + !------------------------------------------------------------------------------- + ! If reservoirs are considered and capacity > 0, route through reservoirs. + !------------------------------------------------------------------------------- + IF (l_reservoirs .AND. rivers%res_cap_current(ip) > 0.0) THEN + reservoir_flow = rivers%rfm_flowin_rp(ip) * 1000.0 / dt + CALL route_reservoirs(dt, abstracted_res_rp(ip), & + rivers%res_cap_current(ip), & + rivers%res_catch(ip), & + rivers%res_critical(ip), & + rivers%res_flood(ip), & + rivers%res_emergency(ip), & + rivers%res_normal_release(ip), & + rivers%res_flood_release(ip), & + rivers%res_storage(ip), reservoir_flow) + rivers%rfm_flowin_rp(ip) = reservoir_flow * dt / 1000.0 + ENDIF + IF (landtype == rfm_land) THEN !Gridcell is land ! land surface (Lewis et al. 2018:eqnB4) diff --git a/src/science/river_routing/rivers_route_trip_mod.F90 b/src/science/river_routing/rivers_route_trip_mod.F90 index fb07095b..93e7048e 100644 --- a/src/science/river_routing/rivers_route_trip_mod.F90 +++ b/src/science/river_routing/rivers_route_trip_mod.F90 @@ -42,15 +42,18 @@ MODULE rivers_route_trip_mod SUBROUTINE rivers_route_trip( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & rivers_outflow_rp, rivers_next_rp, & - rivers_seq_rp, rivers_sto_rp , & - rivers_boxareas_rp, & - rivers_lat_rp, rivers_lon_rp, & - inland_outflow_rp, land_fraction_rp ) + rivers_seq_rp, rivers_sto_rp, & + rivers_boxareas_rp, rivers_lat_rp, rivers_lon_rp,& + inland_outflow_rp, land_fraction_rp, & + abstracted_res_rp,res_cap_current, res_catch, & + res_critical, res_flood,res_emergency, & + res_normal_release, res_flood_release, & + res_storage ) USE jules_rivers_mod, ONLY: & ! imported scalars with intent(in) np_rivers,nstep_rivers,nseqmax,river_mouth,rivers_meander,rivers_speed, & - inland_drainage, l_inland_outflow + inland_drainage, l_inland_outflow, l_reservoirs USE rivers_utils, ONLY: & ! imported procedures @@ -61,6 +64,10 @@ SUBROUTINE rivers_route_trip( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & USE missing_data_mod, ONLY: imdi +USE route_reservoirs_mod, ONLY: & +! imported procedures + route_reservoirs + USE um_types, ONLY: real_jlslsm IMPLICIT NONE @@ -98,6 +105,28 @@ SUBROUTINE rivers_route_trip( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & REAL, INTENT(IN OUT) :: rivers_sto_rp(:) REAL, INTENT(IN) :: land_fraction_rp(:) +REAL(KIND=real_jlslsm), INTENT(IN) :: & + abstracted_res_rp(np_rivers), & + ! Water abstracted from reservoirs over river timestep (kg). + res_cap_current(np_rivers), & + ! Storage capacity of reservoirs (kg). + res_catch(np_rivers), & + ! Upstream catchment area of reservoirs (m2). + res_critical(np_rivers), & + ! Critical storage threshold of reservoirs (kg). + res_flood(np_rivers), & + ! Flood storage threshold of reservoirs (kg). + res_emergency(np_rivers), & + ! Emergency storage threshold of reservoirs (kg). + res_normal_release(np_rivers), & + ! Normal release rate from reservoirs (kg s-1). + res_flood_release(np_rivers) + ! Flood release rate from reservoirs (kg s-1). + +REAL(KIND=real_jlslsm), INTENT(IN OUT) :: & + res_storage(np_rivers) + ! Water stored in reservoirs (kg). + INTEGER :: & ! local scalars (work/loop counters) ip, iseq & @@ -186,6 +215,21 @@ SUBROUTINE rivers_route_trip( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & ! Get index (location in rivers vector) of the point to consider. IF ( rivers_seq_rp(ip) == iseq ) THEN + !------------------------------------------------------------------------------- + ! If reservoirs are considered and capacity > 0, route through them. + !------------------------------------------------------------------------------- + IF (l_reservoirs .AND. res_cap_current(ip) > 0.0) THEN + CALL route_reservoirs(dt, abstracted_res_rp(ip), & + res_cap_current(ip), & + res_catch(ip), & + res_critical(ip), & + res_flood(ip), & + res_emergency(ip), & + res_normal_release(ip), & + res_flood_release(ip), & + res_storage(ip), inflow(ip)) + ENDIF + !----------------------------------------------------------------------- ! Calculate the coefficient "c" of the model. ! c=u/(d*r), where u is effective flow speed, @@ -211,6 +255,8 @@ SUBROUTINE rivers_route_trip( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & !----------------------------------------------------------------------- outflow(ip) = inflow(ip) + (store_old - rivers_sto_rp(ip)) / dt + + !----------------------------------------------------------------------- ! Add outflow to inflow of next downstream point. !----------------------------------------------------------------------- diff --git a/src/science/river_routing/route_reservoirs_mod.F90 b/src/science/river_routing/route_reservoirs_mod.F90 new file mode 100644 index 00000000..55bb3dce --- /dev/null +++ b/src/science/river_routing/route_reservoirs_mod.F90 @@ -0,0 +1,142 @@ +!******************************COPYRIGHT************************************** +! (c) UK Centre for Ecology & Hydrology. +! All rights reserved. +! +! This routine has been licensed to the other JULES partners for use and +! distribution under the JULES collaboration agreement, subject to the terms +! and conditions set out therein. +! +! [Met Office Ref SC0237] +!******************************COPYRIGHT************************************** + +MODULE route_reservoirs_mod + +USE um_types, ONLY: real_jlslsm + +IMPLICIT NONE + +PRIVATE ! Private scope by default +PUBLIC route_reservoirs + +CONTAINS + +!############################################################################## + +SUBROUTINE route_reservoirs( dt, abstracted_res_rp, & + res_cap_current, & + res_catch, & + res_critical, res_flood, & + res_emergency, & + res_normal_release, & + res_flood_release, & + res_storage, reservoir_flow ) + +! Route runoff through reservoirs. + +USE jules_rivers_mod, ONLY: dt_rivers, nstep_rivers + +USE timestep_mod, ONLY: timestep + +IMPLICIT NONE +!------------------------------------------------------------------------------ +! Scalar arguments with INTENT(IN) +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(IN) :: & + dt, & + ! Timestep of routing model (s) + abstracted_res_rp, & + ! Water abstracted from reservoirs (kg). + res_cap_current, & + ! Storage capacity of reservoirs (kg). + res_catch, & + ! Upstream catchment area of reservoirs (m2). + res_critical, & + ! Critical storage threshold of reservoirs (kg). + res_flood, & + ! Flood storage threshold of reservoirs (kg). + res_emergency, & + ! Emergency storage threshold of reservoirs (kg). + res_normal_release, & + ! Normal release rate from reservoirs (kg s-1). + res_flood_release + ! Flood release rate from reservoirs (kg s-1). + +!------------------------------------------------------------------------------ +! Scalar arguments with INTENT(IN OUT) +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm), INTENT(IN OUT) :: & + res_storage, & + ! Water stored in reservoirs (kg). + reservoir_flow + ! Flow in and out of reservoir (kg s-1) + +!------------------------------------------------------------------------------ +! Local scalar variables +!------------------------------------------------------------------------------ +REAL(KIND=real_jlslsm) :: & + k + ! Release coefficient. + +!end of header +!------------------------------------------------------------------------------ + +! Update reservoir storage with inflow and abstractions +res_storage = res_storage + reservoir_flow * dt & + - abstracted_res_rp + +k = MAX( 0., 1. - (res_cap_current - res_storage) / 1000. & + / ( res_catch * 0.2 ) ) + +! Calculate reservoir release based on storage and inflow + +! Release calculation if inflow is not extremely high +IF (reservoir_flow < res_flood_release) THEN + + ! Calculate outflow depending on storage level + IF (res_storage <= res_critical) THEN + reservoir_flow = res_normal_release * & + res_storage / res_flood + ELSE IF (res_storage <= res_emergency) THEN + reservoir_flow = res_normal_release / 2. + & + ( (res_storage - res_critical) / & + (res_emergency - res_critical) ) ** 2. * & + ( res_flood_release - res_normal_release ) + ELSE + reservoir_flow = res_flood_release + END IF + +! Release calculation if inflow is extremely high +ELSE + + ! Calculate outflow depending on storage level + IF (res_storage <= res_critical) THEN + reservoir_flow = res_normal_release * & + res_storage / res_flood + ELSE IF (res_storage <= res_flood) THEN + reservoir_flow = res_normal_release / 2. + & + (res_storage - res_critical) / & + (res_flood - res_critical) * & + ( res_flood_release - res_normal_release ) + ELSE IF (res_storage <= res_emergency) THEN + reservoir_flow = res_flood_release + k * & + (res_storage - res_flood) / & + (res_emergency - res_flood) * & + ( reservoir_flow - res_flood_release ) + ELSE + reservoir_flow = reservoir_flow + END IF + +END IF + +! Ensure that release never exceeds storage +reservoir_flow = MIN( reservoir_flow, res_storage/dt ) + +! Update reservoir storage with release +res_storage = res_storage - reservoir_flow * dt + +RETURN +END SUBROUTINE route_reservoirs + +!############################################################################## + +END MODULE route_reservoirs_mod From 23dc657bb0c92356218a9c3e2d85bb45a5737c20 Mon Sep 17 00:00:00 2001 From: Helen Baron Date: Fri, 28 Aug 2026 11:04:57 +0100 Subject: [PATCH 4/4] fix to prevent indexing errors, and update to update- macro --- rose-meta/jules-standalone/versions.py | 4 +-- .../river_routing/rivers_route_rfm_mod.F90 | 26 ++++++++++--------- .../river_routing/rivers_route_trip_mod.F90 | 22 +++++++++------- 3 files changed, 28 insertions(+), 24 deletions(-) diff --git a/rose-meta/jules-standalone/versions.py b/rose-meta/jules-standalone/versions.py index 2e442c0d..279de4ea 100644 --- a/rose-meta/jules-standalone/versions.py +++ b/rose-meta/jules-standalone/versions.py @@ -45,12 +45,12 @@ from .version81_82 import * -class vn82_txxxx(MacroUpgrade): +class vn82_t148(MacroUpgrade): """Upgrade macro from JULES by Helen Baron""" BEFORE_TAG = "vn8.2" - AFTER_TAG = "vn8.2_txxxx" + AFTER_TAG = "vn8.2_t148" def upgrade(self, config, meta_config=None): """Upgrade a JULES runtime app configuration.""" diff --git a/src/science/river_routing/rivers_route_rfm_mod.F90 b/src/science/river_routing/rivers_route_rfm_mod.F90 index f333fbde..bcb9377c 100644 --- a/src/science/river_routing/rivers_route_rfm_mod.F90 +++ b/src/science/river_routing/rivers_route_rfm_mod.F90 @@ -248,18 +248,20 @@ SUBROUTINE rivers_route_rfm( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & !------------------------------------------------------------------------------- ! If reservoirs are considered and capacity > 0, route through reservoirs. !------------------------------------------------------------------------------- - IF (l_reservoirs .AND. rivers%res_cap_current(ip) > 0.0) THEN - reservoir_flow = rivers%rfm_flowin_rp(ip) * 1000.0 / dt - CALL route_reservoirs(dt, abstracted_res_rp(ip), & - rivers%res_cap_current(ip), & - rivers%res_catch(ip), & - rivers%res_critical(ip), & - rivers%res_flood(ip), & - rivers%res_emergency(ip), & - rivers%res_normal_release(ip), & - rivers%res_flood_release(ip), & - rivers%res_storage(ip), reservoir_flow) - rivers%rfm_flowin_rp(ip) = reservoir_flow * dt / 1000.0 + IF (l_reservoirs) THEN + IF (rivers%res_cap_current(ip) > 0.0) THEN + reservoir_flow = rivers%rfm_flowin_rp(ip) * 1000.0 / dt + CALL route_reservoirs(dt, abstracted_res_rp(ip), & + rivers%res_cap_current(ip), & + rivers%res_catch(ip), & + rivers%res_critical(ip), & + rivers%res_flood(ip), & + rivers%res_emergency(ip), & + rivers%res_normal_release(ip), & + rivers%res_flood_release(ip), & + rivers%res_storage(ip), reservoir_flow) + rivers%rfm_flowin_rp(ip) = reservoir_flow * dt / 1000.0 + ENDIF ENDIF IF (landtype == rfm_land) THEN !Gridcell is land diff --git a/src/science/river_routing/rivers_route_trip_mod.F90 b/src/science/river_routing/rivers_route_trip_mod.F90 index 93e7048e..f90b0ec5 100644 --- a/src/science/river_routing/rivers_route_trip_mod.F90 +++ b/src/science/river_routing/rivers_route_trip_mod.F90 @@ -218,16 +218,18 @@ SUBROUTINE rivers_route_trip( sfc_runoff, sub_sfc_runoff, outflow, baseflow, & !------------------------------------------------------------------------------- ! If reservoirs are considered and capacity > 0, route through them. !------------------------------------------------------------------------------- - IF (l_reservoirs .AND. res_cap_current(ip) > 0.0) THEN - CALL route_reservoirs(dt, abstracted_res_rp(ip), & - res_cap_current(ip), & - res_catch(ip), & - res_critical(ip), & - res_flood(ip), & - res_emergency(ip), & - res_normal_release(ip), & - res_flood_release(ip), & - res_storage(ip), inflow(ip)) + IF (l_reservoirs) THEN + IF (res_cap_current(ip) > 0.0) THEN + CALL route_reservoirs(dt, abstracted_res_rp(ip), & + res_cap_current(ip), & + res_catch(ip), & + res_critical(ip), & + res_flood(ip), & + res_emergency(ip), & + res_normal_release(ip), & + res_flood_release(ip), & + res_storage(ip), inflow(ip)) + ENDIF ENDIF !-----------------------------------------------------------------------