adding new geothermal bedrock

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 |
+| eleanorgb       | Eleanor Burke    | Met Office  | 2026-08-24 |

doc/source/namelists/jules_soil.nml.rst

@@ -113,7 +113,7 @@ This file sets the soil options and parameters. It contains one namelist called
 
    Switch for using a thermal bedrock column beneath the soil column. The bedrock has no hydrological processes - diffusion of heat is the only process represented.
 
-   Properties of the bedrock can be set using :nml:mem:`ns_deep`, :nml:mem:`hcapdeep`, :nml:mem:`hcondeep` and :nml:mem:`dzdeep`.
+   Properties of the bedrock can be set using :nml:mem:`ns_deep`, :nml:mem:`hcapdeep`, :nml:mem:`hcondeep`, :nml:mem:`dzdeep` and :nml:mem:`hflux_geo`.
 
    TRUE
        An additional bedrock column is used below the soil column.
@@ -147,9 +147,9 @@ This file sets the soil options and parameters. It contains one namelist called
    .. nml:member:: hcondeep
 
       :type: real
-      :default: 8.6
+      :default: 3.0
 
-      The heat conductivity of the bedrock (W m\ :sup:`-2` K\ :sup:`-1` ).
+      The heat conductivity of the bedrock (W m\ :sup:`-1` K\ :sup:`-1` ).
 
 
    .. nml:member:: dzdeep
@@ -159,6 +159,13 @@ This file sets the soil options and parameters. It contains one namelist called
 
       The thickness of the bedrock layers (m).
 
+   .. nml:member:: hflux_geo
+
+      :type: real
+      :default: 0.067
+
+      The geothermal heat flux (W m\ :sup:`-2` ).
+
 
 .. nml:member:: cs_min
 

rose-meta/jules-standalone/HEAD/rose-meta.conf

@@ -4700,21 +4700,30 @@ url=https://metoffice.github.io/jules/latest/namelists/jules_soil.nml.html#JULES
 
 [namelist:jules_soil=hcondeep]
 compulsory=true
-description=Thermal conductivity of bedrock (W m-2 K-1)
+description=Thermal conductivity of bedrock (W m-1 K-1)
 range=0.4:12.0
 sort-key=Panel-E06c
 type=real
 url=https://metoffice.github.io/jules/latest/namelists/jules_soil.nml.html#JULES_SOIL::hcondeep
 
+[namelist:jules_soil=hflux_geo]
+compulsory=true
+description=Geothermal heat flux (W/m2)
+range=-0.052:15.6
+sort-key=Panel-E06e
+type=real
+url=https://metoffice.github.io/jules/latest/namelists/jules_soil.nml.html#JULES_SOIL::hflux_geo
+
 [namelist:jules_soil=l_bedrock]
 compulsory=true
-description=Include bedrock below soil column with thermal diffusion; does not include hydrology
+description=Include bedrock below soil column with thermal diffusion - does not include hydrology
            =NOT AVAILABLE TO UM
 sort-key=Panel-E06
 trigger=namelist:jules_soil=ns_deep: .true.;
        = namelist:jules_soil=hcapdeep: .true.;
        = namelist:jules_soil=hcondeep: .true.;
        = namelist:jules_soil=dzdeep: .true.;
+       = namelist:jules_soil=hflux_geo: .true.;
 type=logical
 url=https://metoffice.github.io/jules/latest/namelists/jules_soil.nml.html#JULES_SOIL::l_bedrock
 

rose-meta/jules-standalone/versions.py

@@ -45,15 +45,17 @@
 from .version81_82 import *
 
 
-class vnYY_txxxx(MacroUpgrade):
+class vn82_t76(MacroUpgrade):
 
     """Upgrade macro from JULES by Author"""
 
-    BEFORE_TAG = "vnY.Y"
-    AFTER_TAG = "vnY.Y_txxxx"
+    BEFORE_TAG = "vn8.2"
+    AFTER_TAG = "vn8.2_t76"
 
     def upgrade(self, config, meta_config=None):
         """Upgrade a JULES runtime app configuration."""
 
         # Add settings
+        self.add_setting(config, ["namelist:jules_soil", "hflux_geo"], "0.067")
+        self.change_setting_value(config, ["namelist:jules_soil", "hcondeep"], "3.0")
         return config, self.reports

src/control/shared/jules_soil_mod.F90

@@ -7,6 +7,7 @@
 MODULE jules_soil_mod
 
 USE max_dimensions, ONLY: sm_levels_max
+USE ancil_info,     ONLY: nsoilt
 
 USE missing_data_mod, ONLY: rmdi, imdi
 
@@ -129,12 +130,13 @@ MODULE jules_soil_mod
 
 REAL(KIND=real_jlslsm) ::                                                      &
   hcapdeep = rmdi,                                                             &
-      ! Heat capacity of bedrock
+      ! Heat capacity of bedrock (J/K/m3)
   hcondeep = rmdi,                                                             &
-      ! Thermal conductivity of bedrock
-  dzdeep = rmdi
-      ! Thickness of bedrock layers
-
+      ! Thermal conductivity of bedrock (W/m/K)
+  dzdeep = rmdi,                                                               &
+      ! Thickness of bedrock layers (m)
+  hflux_geo = rmdi
+      ! Geothermal heat flux (W/m2)
 
 !-----------------------------------------------------------------------------
 ! Variable length arrays that can be set using the namelist
@@ -172,7 +174,7 @@ MODULE jules_soil_mod
     l_holdwater, l_tile_soil,                                                  &
 ! Parameters
     cs_min, zsmc, zst, confrac, ns_deep, hcapdeep, hcondeep,                   &
-    dzdeep, dzsoil_io, dzsoil_elev
+    dzdeep, dzsoil_io, dzsoil_elev, hflux_geo
 
 
 
@@ -211,7 +213,6 @@ SUBROUTINE check_jules_soil()
 
 ! Check that sm_levels has been set
 IF ( sm_levels < 1 ) THEN
-  errorstatus = 101
   CALL ereport("check_jules_soil", errorstatus,                                &
                "sm_levels must be between 1 and sm_levels_max")
 END IF
@@ -280,33 +281,44 @@ SUBROUTINE check_jules_soil()
 END IF
 
 IF ( l_bedrock ) THEN
-  ! Check ns_deep, dzdeep, hcondeep. hcapdeep are set and valid
-  IF ( ABS( ns_deep - rmdi ) < EPSILON(1.0) ) THEN
-    CALL ereport(RoutineName, errorstatus, 'ns_deep not found')
+  ! Check ns_deep, dzdeep, hcondeep. hcapdeep, hflux_geo are set and valid
+  IF ( ns_deep == imdi ) THEN
+    CALL ereport(RoutineName, errorstatus, 'ns_deep not specified')
   ELSE IF ( ns_deep < 1 ) THEN
     CALL ereport(RoutineName, errorstatus,                                     &
                  'Bedrock must have at least one layer')
+  ELSE IF ( nsoilt > 1 ) THEN
+    CALL ereport(RoutineName, errorstatus,                                     &
+              'Bedrock not yet compatible with soil tiling')
   END IF
 
   IF ( ABS( dzdeep - rmdi ) < EPSILON(1.0) ) THEN
     CALL ereport(RoutineName, errorstatus, 'dzdeep not found')
   ELSE IF ( dzdeep < 0.01 ) THEN
     CALL ereport(RoutineName, errorstatus,                                     &
-                 'layer thickness for bedrock (dzdeep) must be > 0.01')
+                 'layer thickness for bedrock (dzdeep) must be > 0.01 m')
   END IF
 
   IF ( ABS( hcondeep - rmdi ) < EPSILON(1.0) ) THEN
     CALL ereport(RoutineName, errorstatus, 'hcondeep not found')
   ELSE IF ( hcondeep < 0.4 .OR. hcondeep > 12.0 ) THEN
     CALL ereport(RoutineName, errorstatus,                                     &
-                   'hcondeep must lie in the range 0.4 to 12.0')
+                   'hcondeep must lie in the range 0.4 to 12.0 W/m/K')
   END IF
 
   IF ( ABS( hcapdeep - rmdi ) < EPSILON(1.0) ) THEN
     CALL ereport(RoutineName, errorstatus, 'hcapdeep not found')
   ELSE IF (hcapdeep < 100000.0 .OR. hcapdeep > 8000000.0 ) THEN
     CALL ereport(RoutineName, errorstatus,                                     &
-                   'hcapdeep must lie in the range 100000 to 8000000')
+                   'hcapdeep must lie in the range 100000 to 8000000 J/K/m3')
+  END IF
+
+  ! For continents, areal average 0.067 (F.LUCAZEAU, 2019)
+  IF ( ABS( hflux_geo - rmdi ) < EPSILON(1.0) ) THEN
+    CALL ereport(RoutineName, errorstatus, 'hflux_geo not found')
+  ELSE IF (hflux_geo < -0.052 .OR. hflux_geo > 15.6 ) THEN
+    CALL ereport(RoutineName, errorstatus,                                     &
+                   'hflux_geo must lie in the range -0.052 to 15.6 W/m2')
   END IF
 END IF  ! end if l_bedrock
 
@@ -372,6 +384,9 @@ SUBROUTINE print_nlist_jules_soil()
 WRITE(lineBuffer, *) '  dzdeep = ', dzdeep
 CALL jules_print('jules_soil', lineBuffer)
 
+WRITE(lineBuffer, *) '  hflux_geo = ', hflux_geo
+CALL jules_print('jules_soil', lineBuffer)
+
 WRITE(lineBuffer, *) '  dzsoil_io = ', dzsoil_io(1:sm_levels)
 CALL jules_print('jules_soil', lineBuffer)
 
@@ -418,7 +433,7 @@ SUBROUTINE read_nml_jules_soil (unitnumber)
 ! set number of each type of variable in my_namelist type
 INTEGER, PARAMETER :: no_of_types = 3
 INTEGER, PARAMETER :: n_int = 3
-INTEGER, PARAMETER :: n_real = 8 + sm_levels_max
+INTEGER, PARAMETER :: n_real = 9 + sm_levels_max
 INTEGER, PARAMETER :: n_log = 6
 
 TYPE :: my_namelist
@@ -433,6 +448,7 @@ SUBROUTINE read_nml_jules_soil (unitnumber)
   REAL(KIND=real_jlslsm) :: hcapdeep
   REAL(KIND=real_jlslsm) :: hcondeep
   REAL(KIND=real_jlslsm) :: dzdeep
+  REAL(KIND=real_jlslsm) :: hflux_geo
   REAL(KIND=real_jlslsm) ::dzsoil_io(sm_levels_max)
   REAL(KIND=real_jlslsm) ::dzsoil_elev
   LOGICAL :: l_vg_soil
@@ -467,6 +483,7 @@ SUBROUTINE read_nml_jules_soil (unitnumber)
   my_nml % confrac         = confrac
   my_nml % hcapdeep        = hcapdeep
   my_nml % hcondeep        = hcondeep
+  my_nml % hflux_geo       = hflux_geo
   my_nml % dzdeep          = dzdeep
   my_nml % dzsoil_io       = dzsoil_io
   my_nml % dzsoil_elev     = dzsoil_elev
@@ -493,6 +510,7 @@ SUBROUTINE read_nml_jules_soil (unitnumber)
   hcapdeep        = my_nml % hcapdeep
   hcondeep        = my_nml % hcondeep
   dzdeep          = my_nml % dzdeep
+  hflux_geo       = my_nml % hflux_geo
   dzsoil_io       = my_nml % dzsoil_io
   dzsoil_elev     = my_nml % dzsoil_elev
   l_vg_soil       = my_nml % l_vg_soil

src/control/shared/prognostics.F90

@@ -48,6 +48,8 @@ MODULE prognostics
 
   REAL(KIND=real_jlslsm), ALLOCATABLE :: tsoil_deep_gb(:,:)
     ! Deep soil temperatures (K)
+  REAL(KIND=real_jlslsm), ALLOCATABLE :: tsoil_deep_corr_acc(:,:)
+     ! Accumulated correction in deep soil (bedrock) temperature (K).
   REAL(KIND=real_jlslsm), ALLOCATABLE :: sice_surft(:,:,:)
     ! Snow layer ice mass on tiles (kg/m2)
   REAL(KIND=real_jlslsm), ALLOCATABLE :: sliq_surft(:,:,:)
@@ -201,6 +203,7 @@ MODULE prognostics
   INTEGER, POINTER :: years_since_harvest(:,:)
 
   REAL(KIND=real_jlslsm), POINTER :: tsoil_deep_gb(:,:)
+  REAL(KIND=real_jlslsm), POINTER :: tsoil_deep_corr_acc(:,:)
   REAL(KIND=real_jlslsm), POINTER :: sice_surft(:,:,:)
   REAL(KIND=real_jlslsm), POINTER ::   sliq_surft(:,:,:)
   REAL(KIND=real_jlslsm), POINTER :: snowdepth_surft(:,:)
@@ -409,10 +412,13 @@ SUBROUTINE prognostics_alloc(land_pts, t_i_length, t_j_length,                 &
 ! Only allocate the bedrock tsoil_deep_gb if bedrock is being used
 IF ( l_bedrock ) THEN
   ALLOCATE(progs_data%tsoil_deep_gb(land_pts,ns_deep))
-  progs_data%tsoil_deep_gb(:,:) = 0.0
+  ALLOCATE(progs_data%tsoil_deep_corr_acc(land_pts,ns_deep))
 ELSE
   ALLOCATE(progs_data%tsoil_deep_gb(1,1))
+  ALLOCATE(progs_data%tsoil_deep_corr_acc(1,1))
 END IF
+progs_data%tsoil_deep_corr_acc(:,:) = 0.0
+progs_data%tsoil_deep_gb(:,:) = 0.0
 
 ! Prognostics for microbial methane scheme
 ALLOCATE(progs_data%substr_ch4(land_pts,dim_ch4layer))
@@ -566,6 +572,7 @@ SUBROUTINE prognostics_dealloc(progs_data)
 END IF
 
 DEALLOCATE(progs_data%tsoil_deep_gb)
+DEALLOCATE(progs_data%tsoil_deep_corr_acc)
 
 IF ( ALLOCATED(progs_data%t_home_gb) ) THEN
   DEALLOCATE(progs_data%t_home_gb)
@@ -662,6 +669,7 @@ SUBROUTINE prognostics_assoc(progs,progs_data)
 progs%frac_biocrop_prev_gb => progs_data%frac_biocrop_prev_gb
 progs%triffid_co2_gb => progs_data%triffid_co2_gb
 progs%tsoil_deep_gb => progs_data%tsoil_deep_gb
+progs%tsoil_deep_corr_acc => progs_data%tsoil_deep_corr_acc
 progs%lai_pft => progs_data%lai_pft
 progs%canht_pft => progs_data%canht_pft
 progs%smcl_soilt => progs_data%smcl_soilt
@@ -752,6 +760,7 @@ SUBROUTINE prognostics_nullify(progs)
 NULLIFY(progs%frac_biocrop_prev_gb)
 NULLIFY(progs%triffid_co2_gb)
 NULLIFY(progs%tsoil_deep_gb)
+NULLIFY(progs%tsoil_deep_corr_acc)
 NULLIFY(progs%lai_pft)
 NULLIFY(progs%canht_pft)
 NULLIFY(progs%smcl_soilt)

src/control/shared/surf_couple_extra_mod.F90

@@ -778,6 +778,7 @@ SUBROUTINE surf_couple_extra(                                                  &
       progs%canopy_surft, wtrac_jls%canopy_surft,                              &
       progs%smcl_soilt, psparms%sthf_soilt,                                    &
       psparms%sthu_soilt, crop_vars%sthu_irr_soilt, progs%tsoil_deep_gb,       &
+      progs%tsoil_deep_corr_acc,                                               &
       progs%t_soil_soilt, progs%t_soil_soilt_acc, progs%tsurf_elev_surft,      &
       wtrac_jls%smcl_soilt, wtrac_jls%sthf_soilt, wtrac_jls%sthu_soilt,        &
       fsat_soilt, fwetl_soilt, sthzw_soilt, zw_soilt, wtrac_jls%sthzw_soilt,   &

src/initialisation/standalone/init_soil.F90

@@ -29,11 +29,11 @@ SUBROUTINE init_soil(nml_dir)
 USE jules_soil_mod, ONLY: jules_soil, sm_levels, dzsoil, dzsoil_elev,          &
                           l_vg_soil, l_soil_sat_down, soilhc_method,           &
                           l_bedrock, l_tile_soil, l_broadcast_ancils,          &
-                          check_jules_soil
+                          check_jules_soil, ns_deep
 
 USE jules_surface_mod, ONLY: l_elev_land_ice
 
-USE logging_mod, ONLY: log_info, log_fatal
+USE logging_mod, ONLY: log_info, log_fatal, log_warn
 
 USE errormessagelength_mod, ONLY: errormessagelength
 
@@ -111,9 +111,15 @@ SUBROUTINE init_soil(nml_dir)
                 "soilHc_method = 3 - Chadburn et al. (2015)")
 END IF
 
-IF ( l_bedrock )                                                               &
+IF ( l_bedrock ) THEN
   CALL log_info("init_soil", "Bedrock will be included at base of soil")
 
+  IF ( ns_deep ==  1) THEN
+    CALL log_warn("init_soil", "Bedrock will be modelled as a single layer")
+    CALL log_warn("init_soil", "There is no geothermal heat flux applied")
+  END IF
+END IF
+
 IF ( l_tile_soil ) THEN
   CALL log_info("init_soil",                                                   &
     "l_tile_soil = T. Soil tiling is switched on: nsoilt = nsurft")

src/io/dump/get_dim_info_mod.F90

@@ -229,7 +229,7 @@ SUBROUTINE get_dim_info( l_reading, identifier, ndims, dim_sizes, dim_names,   &
   dim_names(1:ndims) = [ land_dim_name, ch4layer_dim_name ]
   dim_sizes(1:ndims) = [ global_land_pts, dim_ch4layer ]
 
-CASE ( 'tsoil_deep' )
+CASE ( 'tsoil_deep', 'tsoil_deep_corr_acc' )
   ndims = 2
   dim_names(1:ndims) = [ land_dim_name, bedrock_dim_name ]
   dim_sizes(1:ndims) = [ global_land_pts, ns_deep ]

src/io/dump/read_dump_mod.F90

@@ -237,7 +237,7 @@ SUBROUTINE read_dump(file_name, identifiers)
       CALL file_read_var(FILE, var_ids(i),                                     &
                          global_data_2d(:,1:dim_ch4layer))
 
-    CASE ( 'tsoil_deep' )
+    CASE ( 'tsoil_deep', 'tsoil_deep_corr_acc' )
       CALL file_read_var(FILE, var_ids(i), global_data_2d(:,1:ns_deep))
 
     CASE ( 'canopy', 'nsnow', 'rgrain', 'rho_snow', 'snow_tile',               &
@@ -706,6 +706,12 @@ SUBROUTINE read_dump(file_name, identifiers)
       CALL scatter_land_field(global_data_2d(:,n), progs%tsoil_deep_gb(:,n))
     END DO
 
+  CASE ( 'tsoil_deep_corr_acc' )
+    DO n = 1,ns_deep
+      CALL scatter_land_field(global_data_2d(:,n), progs%tsoil_deep_corr_acc(:,n))
+    END DO
+
+
   CASE ( 'sthu_irr' )
     IF ( l_tile_soil .AND. l_broadcast_soilt ) THEN
       DO m = 1,nsoilt

src/io/dump/required_vars_for_configuration_mod.F90

@@ -296,6 +296,7 @@ SUBROUTINE required_vars_for_configuration(nvars, identifiers,                 &
 ! Additional deep soil temperature if bedrock is on
 IF ( l_bedrock ) THEN
   CALL add_to_list( 'tsoil_deep', nvars, identifiers )
+  CALL add_to_list( 'tsoil_deep_corr_acc', nvars, identifiers )
 END IF
 
 ! FLake prognostics

src/io/dump/write_dump_mod.F90

@@ -455,6 +455,12 @@ SUBROUTINE write_dump()
       CALL gather_land_field(progs%tsoil_deep_gb(:,n), global_data_2d(:,n))
     END DO
 
+  CASE ( 'tsoil_deep_corr_acc' )
+    DO n = 1,ns_deep
+      CALL gather_land_field(progs%tsoil_deep_corr_acc(:,n), global_data_2d(:,n))
+    END DO
+
+
     !Case if nsoilt == 1, so it is OK to hardwire the 2nd dimension to 1
   CASE ( 'sthu_irr' )
     DO n = 1,sm_levels
@@ -1013,7 +1019,7 @@ SUBROUTINE write_dump()
       CALL file_write_var(FILE, var_ids(i),                                    &
                           global_data_2d(:,1:dim_ch4layer))
 
-    CASE ( 'tsoil_deep' )
+    CASE ( 'tsoil_deep', 'tsoil_deep_corr_acc' )
       CALL file_write_var(FILE, var_ids(i), global_data_2d(:,1:ns_deep))
 
     CASE ( 'canopy', 'nsnow', 'rgrain', 'rho_snow', 'snow_tile',               &

src/io/model_interface/extract_var.inc

@@ -2218,6 +2218,10 @@ CASE ( 'tsoil_deep' )
   cube_land = cube_from_array(progs%tsoil_deep_gb(:,1:ns_deep))
   cube = map_from_land(cube_land)
 
+CASE ( 'tsoil_deep_corr_acc' )
+  cube_land = cube_from_array(progs%tsoil_deep_corr_acc(:,1:ns_deep))
+  cube = map_from_land(cube_land)
+
 CASE ( 'alb_tile_1' )
   cube_land = cube_from_array(fluxes%alb_surft(:,:,1))
   cube = map_from_land(cube_land)