diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 00c6e6edb..383f9863f 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -34,4 +34,5 @@ | thomasmelvin | Thomas Melvin | Met Office | 2026-01-15 | | tinyendian | Wolfgang Hayek | Earth Sciences New Zealand | 2026-02-02 | | DanStoneMO | Daniel Stone | Met Office | 2026-02-26 | -| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | \ No newline at end of file +| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | +| mo-cjsmith | Chris Smith | Met Office | 2026-05-22 | diff --git a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 index 0907509b0..0fb3c8b46 100644 --- a/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 +++ b/interfaces/physics_schemes_interface/source/algorithm/pres_lev_diags_alg_mod.x90 @@ -12,6 +12,7 @@ module pres_lev_diags_alg_mod use io_config_mod, only: write_diag, use_xios_io use timing_mod, only: start_timing, stop_timing, tik, LPROF use constants_mod, only: r_def, i_def, l_def, str_def + use dycore_constants_mod, only: get_geopotential use initialise_diagnostics_mod, only: init_diag => init_diagnostic_field, & diag_samp => diagnostic_to_be_sampled use lfric_xios_diag_mod, only: get_axis_dimension, get_axis_values @@ -21,6 +22,8 @@ module pres_lev_diags_alg_mod use mr_indices_mod, only: nummr, imr_v use moist_dyn_mod, only: num_moist_factors, total_mass use planet_config_mod, only: p_zero, kappa, gravity, cp + use extrusion_config_mod, only: planet_radius + use formulation_config_mod, only: shallow use planet_constants_mod, only: ex_power implicit none @@ -65,8 +68,8 @@ contains type( field_type ), pointer :: v_in_w3 type( field_type ), pointer :: w_in_wth type( field_type ), pointer :: wetrho_in_wth - type( field_type ), pointer :: height_w3 type( field_type ), pointer :: height_wth + type( field_type ), pointer :: geopotential ! Define fields type( field_type ) :: plev_heaviside, plev_temp, plev_u, plev_v, & @@ -261,12 +264,13 @@ contains plev_geopot_flag = init_diag(plev_geopot, 'plev__geopot', activate=plev_geopot_clim_flag) if ((plev_geopot_flag .or. plev_geopot_clim_flag) .and. use_xios_io) then - height_w3 => get_height_fv(W3, exner_w3%get_mesh_id()) - call invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, & - height_wth, exner_wth, & - nplev, plevs, plev_geopot, & - p_zero, kappa, cp, gravity, & - ex_power) + geopotential => get_geopotential( exner_w3%get_mesh_id() ) + call invoke_geo_on_pres_kernel_type(geopotential, exner_w3, theta_wth, & + height_wth, exner_wth, & + nplev, plevs, plev_geopot, & + p_zero, kappa, cp, gravity, & + planet_radius, ex_power, & + shallow) if (plev_geopot_flag) call plev_geopot%write_field() if (plev_geopot_clim_flag) then diff --git a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 index da1bf2dcf..8aa2ae93a 100644 --- a/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 +++ b/interfaces/physics_schemes_interface/source/kernel/geo_on_pres_kernel_mod.F90 @@ -7,16 +7,16 @@ module geo_on_pres_kernel_mod - use argument_mod, only: arg_type, & - GH_FIELD, GH_SCALAR, & - GH_READ, GH_WRITE, & - GH_INTEGER, & - GH_REAL, CELL_COLUMN, & - ANY_DISCONTINUOUS_SPACE_1, & - ANY_DISCONTINUOUS_SPACE_2 - use fs_continuity_mod, only: WTHETA - use constants_mod, only: r_def, i_def - use kernel_mod, only: kernel_type + use argument_mod, only: arg_type, & + GH_FIELD, GH_SCALAR, & + GH_READ, GH_WRITE, & + GH_INTEGER, GH_LOGICAL, & + GH_REAL, CELL_COLUMN, & + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2 + use fs_continuity_mod, only: WTHETA + use constants_mod, only: r_def, i_def, l_def + use kernel_mod, only: kernel_type implicit none @@ -25,20 +25,21 @@ module geo_on_pres_kernel_mod !> Kernel metadata for PSyclone type, public, extends(kernel_type) :: geo_on_pres_kernel_type private - type(arg_type) :: meta_args(12) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & - arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & - arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & - arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & - arg_type(GH_SCALAR,GH_INTEGER, GH_READ), & + type(arg_type) :: meta_args(13) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & + arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), & + arg_type(GH_SCALAR,GH_INTEGER, GH_READ), & ! arg_type(GH_SCALAR_ARRAY,GH_REAL, GH_READ, 1), see PSyclone issue #1312 - arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & - arg_type(GH_SCALAR,GH_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ), & - arg_type(GH_SCALAR,GH_REAL, GH_READ) & + arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_REAL, GH_READ), & + arg_type(GH_SCALAR,GH_LOGICAL, GH_READ) & /) integer :: operates_on = CELL_COLUMN contains @@ -65,7 +66,9 @@ module geo_on_pres_kernel_mod !> @param[in] kappa Rd / cp !> @param[in] cp Specific heat at constant pressure !> @param[in] gravity Acceleration due to gravity + !> @param[in] planet_radius Radius of the planet !> @param[in] ex_power (cp * lapse_rate ) / g + !> @param[in] shallow If true, gravity is constant !> @param[in] ndf_in Number of degrees of freedom per cell for in fields !> @param[in] undf_in Number of total degrees of freedom for in fields !> @param[in] map_in Dofmap for the cell at the base of the column for in fields @@ -85,7 +88,9 @@ subroutine geo_on_pres_code(nlayers, & kappa, & cp, & gravity, & + planet_radius, & ex_power, & + shallow, & ndf_in, undf_in, map_in, & ndf_wth, undf_wth, map_wth, & ndf_out, undf_out, map_out) @@ -110,12 +115,15 @@ subroutine geo_on_pres_code(nlayers, & real(kind=r_def), intent(inout), dimension(undf_out) :: data_out ! Constants passed explicitly from algorithm - real(kind=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power + real(kind=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power, & + planet_radius real(kind=r_def), intent(in), dimension(nplev) :: plevs + logical(kind=l_def), intent(in) :: shallow ! Internal variables integer(kind=i_def) :: k, level_above, top_df, kp, level_extrap real(kind=r_def) :: desired_ex, h_ref_lev + real(kind=r_def) :: geo_ht_extrap, geo_ht_lowest real(kind=r_def), parameter:: extrap_height = 2000.0_r_def do kp = 1, nplev @@ -138,9 +146,9 @@ subroutine geo_on_pres_code(nlayers, & if (level_above == -1_i_def) then ! Desired level is above model top, extrapolate up - data_out(map_out(1)+kp-1) = data_in(map_in(1)+top_df) - & - (desired_ex - ex_at_data(map_in(1)+top_df)) & - * cp * theta(map_wth(1)+nlayers) / gravity + data_out(map_out(1)+kp-1) = ( data_in(map_in(1)+top_df) - & + (desired_ex - ex_at_data(map_in(1)+top_df)) & + * cp * theta(map_wth(1)+nlayers) ) / gravity else if (level_above == 0_i_def) then ! Desired level is below surface, extrapolate down do k = 1, nlayers @@ -150,10 +158,18 @@ subroutine geo_on_pres_code(nlayers, & exit end if end do - h_ref_lev = (height_wth(map_wth(1)+level_extrap) - data_in(map_in(1))) & + if ( shallow ) then + geo_ht_extrap = height_wth(map_wth(1)+level_extrap) + else + geo_ht_extrap = height_wth(map_wth(1)+level_extrap) / & + ( 1.0_r_def + height_wth(map_wth(1)+level_extrap) / & + planet_radius ) + end if + geo_ht_lowest = data_in(map_in(1)) / gravity + h_ref_lev = ( geo_ht_extrap - geo_ht_lowest ) & / (1.0_r_def - (exner_wth(map_wth(1)+level_extrap) / & ex_at_data(map_in(1)))**ex_power ) - data_out(map_out(1)+kp-1) = data_in(map_in(1)+level_above) & + data_out(map_out(1)+kp-1) = geo_ht_lowest & + h_ref_lev * (1.0_r_def - & (desired_ex/ex_at_data(map_in(1)))**ex_power) else @@ -165,7 +181,7 @@ subroutine geo_on_pres_code(nlayers, & ex_at_data(map_in(1)+level_above)) * & data_in(map_in(1)+level_above-1) ) / & (ex_at_data(map_in(1)+level_above) - & - ex_at_data(map_in(1)+level_above-1)) + ex_at_data(map_in(1)+level_above-1)) / gravity end if end do diff --git a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 index a064523c8..b66bae444 100644 --- a/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 +++ b/interfaces/physics_schemes_interface/source/psy/psykal_lite_phys_mod.F90 @@ -10,7 +10,7 @@ module psykal_lite_phys_mod - use constants_mod, only : i_def, r_def + use constants_mod, only : i_def, r_def, l_def use field_mod, only : field_type, field_proxy_type use integer_field_mod, only : integer_field_type, integer_field_proxy_type use mesh_mod, only : mesh_type @@ -991,14 +991,18 @@ END SUBROUTINE invoke_pres_interp_kernel_type !> see https://github.com/stfc/PSyclone/issues/1312 !> Hence this module could be removed once the PSyclone ticket is !> completed - SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height_wth, exner_wth, nplev, plevs, plev_geopot, & -&p_zero, kappa, cp, gravity, ex_power) + SUBROUTINE invoke_geo_on_pres_kernel_type(geopot_w3, exner_w3, & + theta_wth, height_wth, exner_wth, & + nplev, plevs, plev_geopot, & + p_zero, kappa, cp, gravity, planet_radius, ex_power, shallow) USE geo_on_pres_kernel_mod, ONLY: geo_on_pres_code USE mesh_mod, ONLY: mesh_type - REAL(KIND=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power + REAL(KIND=r_def), intent(in) :: p_zero, kappa, cp, gravity, ex_power, & + planet_radius INTEGER(KIND=i_def), intent(in) :: nplev + LOGICAL(KIND=l_def), intent(in) :: shallow REAL(KIND=r_def), intent(in) :: plevs(nplev) - TYPE(field_type), intent(in) :: height_w3, exner_w3, theta_wth, height_wth, exner_wth, plev_geopot + TYPE(field_type), intent(in) :: geopot_w3, exner_w3, theta_wth, height_wth, exner_wth, plev_geopot INTEGER(KIND=i_def) cell INTEGER(KIND=i_def) loop0_start, loop0_stop INTEGER(KIND=i_def) nlayers @@ -1007,19 +1011,19 @@ SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height REAL(KIND=r_def), pointer, dimension(:) :: height_wth_data => null() REAL(KIND=r_def), pointer, dimension(:) :: theta_wth_data => null() REAL(KIND=r_def), pointer, dimension(:) :: exner_w3_data => null() - REAL(KIND=r_def), pointer, dimension(:) :: height_w3_data => null() - TYPE(field_proxy_type) height_w3_proxy, exner_w3_proxy, theta_wth_proxy, height_wth_proxy, exner_wth_proxy, plev_geopot_proxy - INTEGER(KIND=i_def), pointer :: map_adspc1_height_w3(:,:) => null(), map_adspc2_plev_geopot(:,:) => null(), & + REAL(KIND=r_def), pointer, dimension(:) :: geopot_w3_data => null() + TYPE(field_proxy_type) geopot_w3_proxy, exner_w3_proxy, theta_wth_proxy, height_wth_proxy, exner_wth_proxy, plev_geopot_proxy + INTEGER(KIND=i_def), pointer :: map_adspc1_geopot_w3(:,:) => null(), map_adspc2_plev_geopot(:,:) => null(), & &map_wtheta(:,:) => null() - INTEGER(KIND=i_def) ndf_adspc1_height_w3, undf_adspc1_height_w3, ndf_wtheta, undf_wtheta, ndf_adspc2_plev_geopot, & + INTEGER(KIND=i_def) ndf_adspc1_geopot_w3, undf_adspc1_geopot_w3, ndf_wtheta, undf_wtheta, ndf_adspc2_plev_geopot, & &undf_adspc2_plev_geopot INTEGER(KIND=i_def) max_halo_depth_mesh TYPE(mesh_type), pointer :: mesh => null() ! ! Initialise field and/or operator proxies ! - height_w3_proxy = height_w3%get_proxy() - height_w3_data => height_w3_proxy%data + geopot_w3_proxy = geopot_w3%get_proxy() + geopot_w3_data => geopot_w3_proxy%data exner_w3_proxy = exner_w3%get_proxy() exner_w3_data => exner_w3_proxy%data theta_wth_proxy = theta_wth%get_proxy() @@ -1033,23 +1037,23 @@ SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height ! ! Initialise number of layers ! - nlayers = height_w3_proxy%vspace%get_nlayers() + nlayers = geopot_w3_proxy%vspace%get_nlayers() ! ! Create a mesh object ! - mesh => height_w3_proxy%vspace%get_mesh() + mesh => geopot_w3_proxy%vspace%get_mesh() max_halo_depth_mesh = mesh%get_halo_depth() ! ! Look-up dofmaps for each function space ! - map_adspc1_height_w3 => height_w3_proxy%vspace%get_whole_dofmap() + map_adspc1_geopot_w3 => geopot_w3_proxy%vspace%get_whole_dofmap() map_wtheta => theta_wth_proxy%vspace%get_whole_dofmap() map_adspc2_plev_geopot => plev_geopot_proxy%vspace%get_whole_dofmap() ! - ! Initialise number of DoFs for adspc1_height_w3 + ! Initialise number of DoFs for adspc1_geopot_w3 ! - ndf_adspc1_height_w3 = height_w3_proxy%vspace%get_ndf() - undf_adspc1_height_w3 = height_w3_proxy%vspace%get_undf() + ndf_adspc1_geopot_w3 = geopot_w3_proxy%vspace%get_ndf() + undf_adspc1_geopot_w3 = geopot_w3_proxy%vspace%get_undf() ! ! Initialise number of DoFs for wtheta ! @@ -1070,10 +1074,14 @@ SUBROUTINE invoke_geo_on_pres_kernel_type(height_w3, exner_w3, theta_wth, height ! DO cell=loop0_start,loop0_stop ! - CALL geo_on_pres_code(nlayers, height_w3_data, exner_w3_data, theta_wth_data, height_wth_data, exner_wth_data, nplev, & -&plevs, plev_geopot_data, p_zero, kappa, cp, gravity, ex_power, ndf_adspc1_height_w3, undf_adspc1_height_w3, & -&map_adspc1_height_w3(:,cell), ndf_wtheta, undf_wtheta, map_wtheta(:,cell), ndf_adspc2_plev_geopot, undf_adspc2_plev_geopot, & -&map_adspc2_plev_geopot(:,cell)) + CALL geo_on_pres_code(nlayers, geopot_w3_data, exner_w3_data, & + theta_wth_data, height_wth_data, exner_wth_data, & + nplev, plevs, plev_geopot_data, & + p_zero, kappa, cp, gravity, planet_radius, ex_power, shallow, & + ndf_adspc1_geopot_w3, undf_adspc1_geopot_w3, map_adspc1_geopot_w3(:,cell), & + ndf_wtheta, undf_wtheta, map_wtheta(:,cell), & + ndf_adspc2_plev_geopot, undf_adspc2_plev_geopot, & + map_adspc2_plev_geopot(:,cell)) END DO ! ! Set halos dirty/clean for fields modified in the above loop diff --git a/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf b/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf index 69d9d0086..18a184a10 100644 --- a/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf +++ b/interfaces/physics_schemes_interface/unit-test/kernel/geo_on_pres_kernel_mod_test.pf @@ -8,7 +8,7 @@ module geo_on_pres_kernel_mod_test use funit - use constants_mod, only : r_def, i_def + use constants_mod, only : r_def, i_def, l_def implicit none @@ -83,12 +83,14 @@ contains class(geo_on_pres_test_type), intent(inout) :: this - real(r_def), parameter :: tol = 1.0e-14_r_def - real(r_def), parameter :: p_zero = 100000.0_r_def - real(r_def), parameter :: kappa = 1.0_r_def - real(r_def), parameter :: ex_power = 1.0_r_def - real(r_def), parameter :: cp = 1000.0_r_def - real(r_def), parameter :: gravity = 10.0_r_def + real(r_def), parameter :: tol = 1.0e-14_r_def + real(r_def), parameter :: p_zero = 100000.0_r_def + real(r_def), parameter :: kappa = 1.0_r_def + real(r_def), parameter :: ex_power = 1.0_r_def + real(r_def), parameter :: cp = 1000.0_r_def + real(r_def), parameter :: gravity = 10.0_r_def + real(r_def), parameter :: planet_radius = 6.4e6_r_def + logical(l_def), parameter :: shallow = .true. ! Input arrays this%map_in = 1_i_def @@ -105,6 +107,8 @@ contains this%data_in(1) = 0.0_r_def this%data_in(2) = 100.0_r_def this%data_in(3) = 1000.0_r_def + ! Conver heights to geopotential + this%data_in(:) = gravity * this%data_in(:) this%height(1) = 0.0_r_def this%height(2) = 100.0_r_def this%height(3) = 3000.0_r_def @@ -129,7 +133,9 @@ contains kappa, & cp, & gravity, & + planet_radius, & ex_power, & + shallow, & ndf_in, & undf_in, & this%map_in, & diff --git a/rose-stem/app/lfric_atm/opt/rose-app-regrav.conf b/rose-stem/app/lfric_atm/opt/rose-app-regrav.conf new file mode 100644 index 000000000..fa2c9058e --- /dev/null +++ b/rose-stem/app/lfric_atm/opt/rose-app-regrav.conf @@ -0,0 +1,5 @@ +[namelist:formulation] +shallow=.false. + +[namelist:initialization] +regravitate='geopot' diff --git a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc index a49293d58..02837ff0b 100644 --- a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc +++ b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc @@ -50,16 +50,16 @@ "panel_decomp": "auto_nonuniform", }) %} -{% elif task_ns.conf_name == "nwp_gal9_deep-C48_MG" %} +{% elif task_ns.conf_name == "nwp_gal9_regrav-C48_MG" %} {% do task_dict.update({ - "opt_confs": ["um_dump","deep","test_diags","physics_segmentation"], + "opt_confs": ["um_dump","regrav","test_diags","physics_segmentation"], "resolution": "C48_MG", "DT": 1800, "tsteps": 36, "mpi_parts": 24, "kgo_checks": ["checksum"], - "plot_str": "plot_map.py $NODAL_DATA_DIR/lfric_diag.nc $PLOT_DIR", + "plot_str": "plot_map.py $NODAL_DATA_DIR/lfric_diagnostics.nc $PLOT_DIR", }) %} {% elif task_ns.conf_name == "nwp_gal9_nomg-C48" %} diff --git a/rose-stem/site/meto/groups/groups_lfric_atm.cylc b/rose-stem/site/meto/groups/groups_lfric_atm.cylc index ea319d080..7e73477b1 100644 --- a/rose-stem/site/meto/groups/groups_lfric_atm.cylc +++ b/rose-stem/site/meto/groups/groups_lfric_atm.cylc @@ -19,6 +19,7 @@ "lfric_atm_nwp_azspice_extra": [ "lfric_atm_nwp_gal9-C48_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9_debug-C48_MG_azspice_gnu_full-debug-32bit", + "lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit", "lfric_atm_nwp_gal9_coarse_aero-C48_MG_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9_da-C12_azspice_gnu_fast-debug-32bit", "lfric_atm_nwp_gal9_eda-C12_azspice_gnu_fast-debug-32bit", diff --git a/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt new file mode 100644 index 000000000..652a56c0d --- /dev/null +++ b/rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_gal9_regrav-C48_MG_azspice_gnu_full-debug-64bit.txt @@ -0,0 +1,9 @@ +Inner product checksum rho = 411AFEB6C87F3787 +Inner product checksum theta = 42718B765093CD28 +Inner product checksum u = 4552C9D6D84CAF9C +Inner product checksum mr1 = 40399BB2B2F333E2 +Inner product checksum mr2 = 3F39BBA239439B78 +Inner product checksum mr3 = 3EF535A3F74F3252 +Inner product checksum mr4 = 3F2BE0BCF15A0966 +Inner product checksum mr5 = 0 +Inner product checksum mr6 = 0 diff --git a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf index 4380951b7..208b7fdfa 100644 --- a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf +++ b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf @@ -1971,6 +1971,7 @@ help=When true, the gravitational field is constant, as would be the under the s =from planet centre. !kind=default sort-key=Panel-A04 +trigger=namelist:initialization=regravitate: this == ".false."; type=logical [namelist:formulation=si_momentum_equation] @@ -3220,6 +3221,27 @@ help=Horizontal winds are read in on a W2h space from a start dump. If sort-key=02b type=logical +[namelist:initialization=regravitate] +compulsory=true +description=Method for switching gravity +!enumeration=true +help=Method for switching gravity from constant to height-varying. + =Used when the shallow switch is false and reading a start dump + =that has been created by a constant gravity model, such as the UM. + =All methods calculate an Exner field in hydrostatic balance with + =height-varying gravity. + =Isothermodynamic: Preserve absolute temperature and pressure on + = model levels by allowing the heights of levels + = to change. Absolute temperature is then interpolated + = back onto the original model levels. + =Geopotential remap: + = Physical height and geopotential height are the same for data in the + = start dump, due to gravity being constant. +!kind=default +sort-key=Panel-A04b +value-titles=Geopotential remap, Isothermal vn1, Isothermal vn2, Isothermal vn3 +values='geopot', 'isotherm1', 'isotherm2', 'isotherm3' + [namelist:initialization=sea_ice_source] compulsory=true description=Where to read the sea ice fields from diff --git a/science/gungho/rose-meta/lfric-gungho/versions.py b/science/gungho/rose-meta/lfric-gungho/versions.py index 01798ad2b..6667b0aea 100644 --- a/science/gungho/rose-meta/lfric-gungho/versions.py +++ b/science/gungho/rose-meta/lfric-gungho/versions.py @@ -18,16 +18,13 @@ def __repr__(self): __str__ = __repr__ -""" -Copy this template and complete to add your macro +class vn31_t218(MacroUpgrade): + """Upgrade macro for issue #218 by Chris Smith.""" -class vnXX_txxx(MacroUpgrade): - # Upgrade macro for by - - BEFORE_TAG = "vnX.X" - AFTER_TAG = "vnX.X_txxx" + BEFORE_TAG = "vn3.1" + AFTER_TAG = "vn3.1_t218" def upgrade(self, config, meta_config=None): - # Add settings + self.add_setting(config, ["namelist:initialization", "regravitate"], "'isotherm1'") return config, self.reports -""" + diff --git a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 index 8f81cc672..af6cdd9c0 100644 --- a/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 +++ b/science/gungho/source/algorithm/physics/map_fd_to_prognostics_alg_mod.x90 @@ -51,7 +51,9 @@ module map_fd_to_prognostics_alg_mod use random_perturb_field_alg_mod, only: random_perturb_field use mr_indices_mod, only: imr_v, imr_cl, imr_ci, imr_r, & imr_s, imr_g - use moist_dyn_mod, only: num_moist_factors, gas_law + use moist_dyn_mod, only: num_moist_factors, & + gas_law, & + total_mass use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg ! Configuration modules @@ -59,8 +61,16 @@ module map_fd_to_prognostics_alg_mod geometry_spherical, & topology, & topology_fully_periodic + use extrusion_config_mod, only: planet_radius use formulation_config_mod, only: rotating, shallow use finite_element_config_mod, only: coord_system, coord_system_native + use initialization_config_mod, only: model_eos_height, & + read_w2h_wind, & + regravitate, & + regravitate_geopot, & + regravitate_isotherm1, & + regravitate_isotherm2, & + regravitate_isotherm3 use physics_config_mod, only: sample_physics_winds, & sample_physics_winds_correction use planet_config_mod, only: gravity, p_zero, kappa, rd, cp @@ -114,15 +124,6 @@ contains type(integer_field_type), pointer :: face_selector_ew type(integer_field_type), pointer :: face_selector_ns - integer(i_def) :: model_eos_height - logical(l_def) :: read_w2h_wind - type(namelist_type), pointer :: initialization_nml - - initialization_nml => configuration%get_namelist('initialization') - call initialization_nml%get_value( 'model_eos_height', model_eos_height ) - call initialization_nml%get_value( 'read_w2h_wind', read_w2h_wind ) - nullify( initialization_nml ) - ! FD (source) fields if ( read_w2h_wind ) then call fd_fields%get_field('h_wind', h_wind_in_w2h) @@ -347,8 +348,8 @@ contains use sci_enforce_bc_kernel_mod, only : enforce_bc_kernel_type use hydrostatic_coriolis_kernel_mod, & only : hydrostatic_coriolis_kernel_type - use hydro_shallow_to_deep_kernel_mod, & - only : hydro_shallow_to_deep_kernel_type + use regrav_isotherm_kernel_mod, only : regrav_isotherm_kernel_type + use regrav_geopot_kernel_mod, only : regrav_geopot_kernel_type use moist_dyn_mod, only : num_moist_factors use matrix_vector_kernel_mod, only : matrix_vector_kernel_type use sci_nodal_xyz_coordinates_kernel_mod, & @@ -374,6 +375,7 @@ contains type(field_type ) :: c_vert_term, r_c type(field_type ) :: exner_in_wth type(field_type ) :: temperature + type(field_type ) :: geopotential_g_const type(field_type), pointer :: geopotential type(field_type), pointer :: chi(:) type(field_type), pointer :: panel_id @@ -383,7 +385,7 @@ contains type(mesh_type), pointer :: mesh type(field_type), pointer :: mask type(field_type), target :: dummy_mask - integer(i_def) :: eos_level + integer(i_def) :: eos_level, nlayers nullify(geopotential, chi, panel_id, height_w3, & height_wth, vert_coriolis, mesh, mask ) @@ -427,7 +429,8 @@ contains ! Get the index of the level closest to the height, ignoring the top eta ! level as eos_level is to be used for a W3 space rather than Wtheta eos_level = get_nearest_level( mesh, eos_height ) - if ( eos_level == (mesh%get_nlayers() + 1) ) eos_level = eos_level - 1 + nlayers = mesh%get_nlayers() + if ( eos_level == nlayers + 1 ) eos_level = nlayers if (shallow) then ! Initialize the vertical profile by starting with the equation of state at @@ -449,7 +452,12 @@ contains ! by rederiving potential temperature (theta) and density (rho). call theta%copy_field_properties(exner_in_wth) call theta%copy_field_properties(temperature) - call invoke( & + call geopotential%copy_field_properties(geopotential_g_const) + call invoke( a_times_X( geopotential_g_const, gravity, height_w3 ) ) + + select case( regravitate ) + case( regravitate_isotherm2 ) + call invoke( & ! Compute exner from EoS sample_eos_pressure_kernel_type( exner, rho, & theta, moist_dyn(gas_law) ), & @@ -460,17 +468,86 @@ contains ! Compute absolute temperature T = theta * exner X_times_Y( temperature, exner_in_wth, theta), & + ! Find exner at level 1 by integrating hydrostatic equation + ! downwards from eos_level, preserving theta. Gravity + ! assumed constant. + hydrostatic_coriolis_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, geopotential_g_const, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + ! Find exner at all levels by integrating hydrostatic + ! equation upwards from level 1, preserving absolute + ! temperature. Gravity is height-dependent. + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & + mask, cp ) ) + case( regravitate_isotherm3 ) + call invoke( & + ! Find exner below eos_level by integrating hydrostatic + ! equation downwards, assuming constant gravity. + hydrostatic_coriolis_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, geopotential_g_const, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + + ! Find exner in Wtheta space + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & - ! Rebalance exner in vertical, using hydrostatic balance. - ! This integrates down using g (shallow) from eos_level and - ! then integrates up using phi (deep). - hydro_shallow_to_deep_kernel_type( exner, temperature, & - c_vert_term, & - moist_dyn, geopotential, & - height_w3, height_wth, & - mask, gravity, & - cp, eos_level ), & - + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + ! Find exner at all levels by integrating hydrostatic + ! equation upwards from level 1, preserving absolute + ! temperature. Gravity is height-dependent. + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & + mask, cp ) ) + case( regravitate_geopot ) + call invoke( & + ! Find exner below eos_level by integrating hydrostatic + ! equation downwards, assuming constant gravity. + hydrostatic_coriolis_kernel_type( exner, rho, theta, & + c_vert_term, & + moist_dyn, geopotential_g_const, & + height_w3, mask, & + p_zero, kappa, & + rd, cp, & + eos_level ), & + ! Find exner in Wtheta space + sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & + height_wth, height_w3 ), & + ! Compute absolute temperature T = theta * exner + X_times_Y( temperature, exner_in_wth, theta), & + ! Interpolate temperature from heights corresponding to + ! geopotential height onto heights of the model grid. + regrav_geopot_kernel_type( temperature, height_wth, & + planet_radius ), & + ! Find exner at all levels by integrating hydrostatic + ! equation upwards from level 1, preserving absolute + ! temperature. Gravity is height-dependent. + regrav_isotherm_kernel_type( exner, temperature, & + c_vert_term, & + moist_dyn(gas_law), & + moist_dyn(total_mass), & + geopotential, & + height_w3, height_wth, & + mask, cp ) ) + end select + call invoke( & ! Recompute exner in Wtheta (required to calculate theta) sample_w3_to_wtheta_kernel_type( exner_in_wth, exner, & height_wth, height_w3 ), & diff --git a/science/gungho/source/kernel/initialisation/hydro_shallow_to_deep_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hydro_shallow_to_deep_kernel_mod.F90 deleted file mode 100644 index a8b733462..000000000 --- a/science/gungho/source/kernel/initialisation/hydro_shallow_to_deep_kernel_mod.F90 +++ /dev/null @@ -1,211 +0,0 @@ -!----------------------------------------------------------------------------- -! (C) Crown copyright 2025 Met Office. All rights reserved. -! The file LICENCE, distributed with this code, contains details of the terms -! under which the code may be used. -!----------------------------------------------------------------------------- - -!> @brief Defines exner pressure from vertical balance with a deep gravity, -!! but using input data based on a shallow gravity. -!> @details Using the input exner pressure at level eos_index, integrate down -!! to the surface using hydrostatic balance, plus Coriolis terms, -!! using a discretisation using absolute temperature and constant -!! gravity g (shallow definition). Then integrate from the surface to -!! the top but using the model geopotential (may be defined as shallow -!! or deep depending on configuration). -module hydro_shallow_to_deep_kernel_mod - -use argument_mod, only : arg_type, func_type, & - GH_FIELD, GH_REAL, & - GH_SCALAR, GH_INTEGER, & - GH_READ, GH_READWRITE, & - ANY_SPACE_9, ANY_SPACE_1, & - GH_BASIS, CELL_COLUMN, GH_EVALUATOR -use constants_mod, only : r_def, i_def -use fs_continuity_mod, only : Wtheta, W3 -use kernel_mod, only : kernel_type -use reference_element_mod, only : B - -implicit none - -private - -!------------------------------------------------------------------------------- -! Public types -!------------------------------------------------------------------------------- -type, public, extends(kernel_type) :: hydro_shallow_to_deep_kernel_type - private - type(arg_type) :: meta_args(11) = (/ & - arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD*3, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & - arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_REAL, GH_READ), & - arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & - /) - type(func_type) :: meta_funcs(2) = (/ & - func_type(W3, GH_BASIS), & - func_type(Wtheta, GH_BASIS) & - /) - integer :: operates_on = CELL_COLUMN - integer :: gh_shape = GH_EVALUATOR -contains - procedure, nopass :: hydro_shallow_to_deep_code -end type - -!------------------------------------------------------------------------------- -! Contained functions/subroutines -!------------------------------------------------------------------------------- -public :: hydro_shallow_to_deep_code - -contains - -!! @param[in] nlayers Number of layers -!! @param[in,out] exner Exner pressure field -!! @param[in] temperature Absolute temperature field -!! @param[in] coriolis_term Vertical component of the coriolis term -!! @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon -!! @param[in] moist_dyn_tot Total mass factor 1 + sum m_x -!! @param[in] moist_dyn_fac Water factor -!! @param[in] phi Geopotential field -!! @param[in] height_w3 Height coordinate in w3 -!! @param[in] height_wth Height coordinate in wth -!! @param[in] w3_mask LBC mask or Dummy mask for w3 space -!! @param[in] gravity The planet gravity -!! @param[in] cp Specific heat of dry air at constant pressure -!! @param[in] eos_index Vertical level at which the equation of state -!! is satisfied -!! @param[in] ndf_w3 Number of degrees of freedom per cell for w3 -!! @param[in] undf_w3 Total number of degrees of freedom for w3 -!! @param[in] map_w3 Dofmap for the cell at column base for w3 -!! @param[in] basis_w3 Basis functions evaluated at w3 nodes -!! @param[in] ndf_wt Number of degrees of freedom per cell for wtheta -!! @param[in] undf_wt Total number of degrees of freedom for wtheta -!! @param[in] map_wt Dofmap for the cell at column base for wt -!! @param[in] basis_wt Basis functions evaluated at wt nodes -subroutine hydro_shallow_to_deep_code( & - nlayers, & - exner, & - temperature, & - coriolis_term, & - moist_dyn_gas, & - moist_dyn_tot, & - moist_dyn_fac, & - phi, & - height_w3, & - height_wth, & - w3_mask, & - gravity, & - cp, & - eos_index, & - ndf_w3, & - undf_w3, & - map_w3, & - basis_w3, & - ndf_wt, & - undf_wt, & - map_wt, & - basis_wt ) - - implicit none - - ! Arguments - integer(kind=i_def), intent(in) :: nlayers, & - ndf_w3, & - undf_w3, & - ndf_wt, & - undf_wt - integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 - integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt - integer(kind=i_def), intent(in) :: eos_index - - real(kind=r_def), dimension(undf_w3), intent(inout) :: exner - real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & - w3_mask, & - phi - real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & - moist_dyn_tot, & - moist_dyn_fac - real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & - height_wth - real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term - real(kind=r_def), dimension(1,ndf_w3,ndf_w3), intent(in) :: basis_w3 - real(kind=r_def), dimension(1,ndf_wt,ndf_w3), intent(in) :: basis_wt - real(kind=r_def), intent(in) :: gravity - real(kind=r_def), intent(in) :: cp - - ! Internal variables - integer(kind=i_def) :: k, eos_index_m1 - real(kind=r_def) :: temp_moist, dz, & - weight1, weight2 - - ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) - ! setting exner to 1 - if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then - do k = 0, nlayers-1 - exner(map_w3(1) + k ) = 1.0_r_def - enddo - return - end if - - ! Levels in the kernel are numbered 0 to nlayers-1, rather than 1 to nlayers - eos_index_m1 = eos_index - 1 - - ! Integrate from eos_level to surface using shallow gravity - do k = eos_index_m1, 1, -1 - - ! Absolute temperature at which dry air would have same pressure and density - ! T_moist = T (1 + mr/eps)/(1 + sum mr) - temp_moist = moist_dyn_gas( map_wt(1) + k ) * & - temperature( map_wt(1) + k ) / & - moist_dyn_tot( map_wt(1) + k ) - - ! Discretisation weights e.g. weight 1 = w * dz where w = (z_k - z_k-1/2 ) / dz - ! and weight 2 = (1 - w) * dz, dz = z_k - z_k-1 - ! - ! z_k W3 k - ! z_k-1/2 Wtheta k - ! z_k-1 W3 k-1 - weight1 = height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) - weight2 = height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) - - ! Vertical momentum equation in hydrostatic balance (plus coriolis) with T = theta Pi - ! cp T d Pi /dz = Pi ( F - d Phi/ dz ) - ! Discretised, using d Phi /dz = g as - ! cp T_{k-1/2} ( Pi_k - Pi_k-1 ) / dz = ( w Pi_k + (1-w) Pi_k-1 ) ( F_k-1/2 - g) - ! Rearranged as - ! Pi_k-1 = Pi_k * ( cp T_{k-1/2} + (g - F_k-1/2) * (1-w) * dz ) / - ! ( cp T_{k-1/2} - (g - F_k-1/2) * w * dz ) - exner( map_w3 (1) + k-1 ) = exner( map_w3(1) + k ) * & - ( cp * temp_moist + ( gravity - coriolis_term( map_wt(1) + k ) ) * weight2 ) / & - ( cp * temp_moist - ( gravity - coriolis_term( map_wt(1) + k ) ) * weight1 ) - - end do - - ! Integrate from surface to top using model geopotential phi - do k = 1, nlayers-1 - - temp_moist = moist_dyn_gas( map_wt(1) + k ) * temperature( map_wt(1) + k ) / & - moist_dyn_tot( map_wt(1) + k ) - - dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) - weight1 = ( height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) ) / dz - weight2 = ( height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) ) /dz - - ! Pi_k = Pi_k-1 * ( cp T_k-1/2 - ( (phi_k - phi_k-1) - F_k-1/2 * dz) * (1-w) ) / - ! ( cp T_k-1/2 + ( (phi_k - phi_k-1) - F_k-1/2 * dz) * w ) - exner( map_w3(1) + k ) = exner( map_w3 (1) + k-1 ) * & - ( cp * temp_moist - ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & - - coriolis_term( map_wt(1) + k ) * dz ) * weight1 ) / & - ( cp * temp_moist + ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & - - coriolis_term( map_wt(1) + k ) * dz ) * weight2 ) - - end do - -end subroutine hydro_shallow_to_deep_code - -end module hydro_shallow_to_deep_kernel_mod diff --git a/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 new file mode 100644 index 000000000..498b323ab --- /dev/null +++ b/science/gungho/source/kernel/initialisation/regrav_geopot_kernel_mod.F90 @@ -0,0 +1,130 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!> @brief Convert temperature on geopotential height to temperature on height +!> @details The input temperature is assumed to have been generated by a model +!> with constant gravity. For such a model physical height and +!> geopotential height are identical. In this kernel gravity is assumed +!> to vary correctly with height, so that physical and geopotential +!> heights are no longer identical. The physical height, height_phys, +!> corresponding to height_wth, interpreted as geopotential height, is +!> calculated. Temperature is then interpolated from height_phys onto +!> model levels height_wth. +module regrav_geopot_kernel_mod + +use argument_mod, only: arg_type, & + GH_FIELD, GH_REAL, & + GH_READ, GH_READWRITE, & + GH_SCALAR, & + CELL_COLUMN +use constants_mod, only: r_def, i_def +use fs_continuity_mod, only: Wtheta +use kernel_mod, only: kernel_type + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +type, public, extends(kernel_type) :: regrav_geopot_kernel_type + private + type(arg_type) :: meta_args(3) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_SCALAR, GH_REAL, GH_READ) & + /) + integer :: operates_on = CELL_COLUMN +contains + procedure, nopass :: regrav_geopot_code +end type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- +public :: regrav_geopot_code + +contains + +!> @param[in] nlayers Number of layers +!> @param[in,out] temperature Absolute temperature field +!> @param[in] height_wth Height coordinate in wth +!> @param[in] planet_radius Radius of the planet +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +subroutine regrav_geopot_code( nlayers, & + temperature, & + height_wth, & + planet_radius, & + ndf_wt, & + undf_wt, & + map_wt ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + + real(kind=r_def), dimension(undf_wt), intent(inout) :: temperature + real(kind=r_def), dimension(undf_wt), intent(in) :: height_wth + real(kind=r_def), intent(in) :: planet_radius + + ! Internal variables + integer(kind=i_def) :: k + real(kind=r_def) :: height_phys(0:nlayers) + real(kind=r_def) :: temp(0:nlayers) + + ! Height corresponding to height_wth when latter represents geopotential height + do k = 0, nlayers + height_phys(k) = planet_radius * height_wth(map_wt(1)+k) / & + ( planet_radius - height_wth(map_wt(1)+k) ) + end do + +! Temporary storage for T + do k = 0, nlayers + temp(k) = temperature( map_wt(1) + k ) + end do + + ! Interpolate temperature from newly computed grid back onto current model grid + call interp( nlayers+1, height_phys, temp, nlayers+1, height_wth(map_wt(1)), & + temperature(map_wt(1)) ) + +end subroutine regrav_geopot_code + +subroutine interp( nin, zin, fin, nout, zout, fout ) + +implicit none + +integer(kind=i_def), intent(in) :: nin, nout +real(kind=r_def), intent(in) :: zin(nin), fin(nin), zout(nout) +real(kind=r_def), intent(inout) :: fout(nout) + +integer(kind=i_def) :: kout, kin +real(kind=r_def) :: xi + +do kout = 1, nout + if ( zout(kout) <= zin(1) ) then + fout(kout) = fin(1) + else if ( zout(kout) >= zin(nin) ) then + fout(kout) = fin(nin) + else + do kin = 1, nin - 1 + if ( zin(kin+1) > zout(kout) ) then + xi = ( zout(kout) - zin(kin) ) / ( zin(kin+1) - zin(kin) ) + fout(kout) = ( 1.0_r_def - xi ) * fin(kin) + xi * fin(kin+1) + end if + end do + end if +end do + +end subroutine interp + +end module regrav_geopot_kernel_mod diff --git a/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 new file mode 100644 index 000000000..556750dc3 --- /dev/null +++ b/science/gungho/source/kernel/initialisation/regrav_isotherm_kernel_mod.F90 @@ -0,0 +1,152 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!> @brief Defines exner pressure from vertical balance with a deep gravity, +!> but using input data based on a shallow gravity. +!> @details Using the input exner pressure at level eos_index, integrate down +!> to the surface using hydrostatic balance, plus Coriolis terms, +!> using a discretisation using absolute temperature and constant +!> gravity g (shallow definition). Then integrate from the surface to +!> the top but using the model geopotential (may be defined as shallow +!> or deep depending on configuration). +module regrav_isotherm_kernel_mod + +use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_SCALAR, GH_INTEGER, & + GH_READ, GH_READWRITE, & + CELL_COLUMN +use constants_mod, only : r_def, i_def +use fs_continuity_mod, only : Wtheta, W3 +use kernel_mod, only : kernel_type + +implicit none + +private + +!------------------------------------------------------------------------------- +! Public types +!------------------------------------------------------------------------------- +type, public, extends(kernel_type) :: regrav_isotherm_kernel_type + private + type(arg_type) :: meta_args(10) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_SCALAR, GH_REAL, GH_READ) + /) + integer :: operates_on = CELL_COLUMN +contains + procedure, nopass :: regrav_isotherm_code +end type + +!------------------------------------------------------------------------------- +! Contained functions/subroutines +!------------------------------------------------------------------------------- +public :: regrav_isotherm_code + +contains + +!> @param[in] nlayers Number of layers +!> @param[in,out] exner Exner pressure field +!> @param[in] temperature Absolute temperature field +!> @param[in] coriolis_term Vertical component of the coriolis term +!> @param[in] moist_dyn_gas Gas factor 1+ m_v/epsilon +!> @param[in] moist_dyn_tot Total mass factor 1 + sum m_x +!> @param[in] phi Geopotential field +!> @param[in] height_w3 Height coordinate in w3 +!> @param[in] height_wth Height coordinate in wth +!> @param[in] w3_mask LBC mask or Dummy mask for w3 space +!> @param[in] ndf_w3 Number of degrees of freedom per cell for w3 +!> @param[in] undf_w3 Total number of degrees of freedom for w3 +!> @param[in] map_w3 Dofmap for the cell at column base for w3 +!> @param[in] ndf_wt Number of degrees of freedom per cell for wtheta +!> @param[in] undf_wt Total number of degrees of freedom for wtheta +!> @param[in] map_wt Dofmap for the cell at column base for wt +subroutine regrav_isotherm_code( nlayers, & + exner, & + temperature, & + coriolis_term, & + moist_dyn_gas, & + moist_dyn_tot, & + phi, & + height_w3, & + height_wth, & + w3_mask, & + cp, & + ndf_w3, & + undf_w3, & + map_w3, & + ndf_wt, & + undf_wt, & + map_wt ) + + implicit none + + ! Arguments + + integer(kind=i_def), intent(in) :: nlayers, & + ndf_w3, & + undf_w3, & + ndf_wt, & + undf_wt + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt + + real(kind=r_def), dimension(undf_w3), intent(inout) :: exner + real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3, & + w3_mask, & + phi + real(kind=r_def), dimension(undf_wt), intent(in) :: moist_dyn_gas, & + moist_dyn_tot + real(kind=r_def), dimension(undf_wt), intent(in) :: temperature, & + height_wth + real(kind=r_def), dimension(undf_wt), intent(in) :: coriolis_term + real(kind=r_def), intent(in) :: cp + + ! Internal variables + integer(kind=i_def) :: k + real(kind=r_def) :: temp_virtual + real(kind=r_def) :: dz, weight1, weight2 + + ! Return if the mask is 0 (with tolerance of 0.5 as mask is real, 0 or 1) + ! setting exner to 1 + if ( w3_mask( map_w3(1) ) < 0.5_r_def ) then + do k = 0, nlayers-1 + exner(map_w3(1) + k ) = 1.0_r_def + enddo + return + end if + + ! Integrate from surface to top using model geopotential phi + do k = 1, nlayers-1 + + temp_virtual = moist_dyn_gas( map_wt(1) + k ) * temperature( map_wt(1) + k ) / & + moist_dyn_tot( map_wt(1) + k ) + + dz = height_w3( map_w3(1) + k ) - height_w3( map_w3(1) + k - 1 ) + weight1 = ( height_w3( map_w3(1) + k ) - height_wth( map_wt(1) + k ) ) / dz + weight2 = ( height_wth( map_wt(1) + k ) - height_w3( map_w3(1) + k - 1 ) ) /dz + + ! Pi_k = Pi_k-1 * ( cp T_k-1/2 - ( (phi_k - phi_k-1) - F_k-1/2 * dz) * (1-w) ) / + ! ( cp T_k-1/2 + ( (phi_k - phi_k-1) - F_k-1/2 * dz) * w ) + exner( map_w3(1) + k ) = exner( map_w3 (1) + k-1 ) * & + ( cp * temp_virtual - ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & + - coriolis_term( map_wt(1) + k ) * dz ) * weight1 ) / & + ( cp * temp_virtual + ( phi( map_w3(1) + k ) - phi( map_w3(1) + k - 1 ) & + - coriolis_term( map_wt(1) + k ) * dz ) * weight2 ) + + end do + +end subroutine regrav_isotherm_code + +end module regrav_isotherm_kernel_mod diff --git a/science/gungho/unit-test/kernel/initialisation/hydro_shallow_to_deep_kernel_mod_test.pf b/science/gungho/unit-test/kernel/initialisation/hydro_shallow_to_deep_kernel_mod_test.pf deleted file mode 100644 index 4610f6287..000000000 --- a/science/gungho/unit-test/kernel/initialisation/hydro_shallow_to_deep_kernel_mod_test.pf +++ /dev/null @@ -1,238 +0,0 @@ -!----------------------------------------------------------------------------- -! (C) Crown copyright 2025 Met Office. All rights reserved. -! The file LICENCE, distributed with this code, contains details of the terms -! under which the code may be used. -!----------------------------------------------------------------------------- - -!> Test the pressure initialization computation. -module hydro_shallow_to_deep_kernel_mod_test - - use constants_mod, only : i_def, r_def - use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w0_m3x3_q3x3x3_size, & - get_w3_m3x3_q3x3x3_size, & - get_wtheta_m3x3_q3x3x3_size - use get_unit_test_w3nodal_basis_mod, only : get_wtheta_w3nodal_basis - use get_unit_test_m3x3_dofmap_mod, only : get_w0_m3x3_dofmap, & - get_w3_m3x3_dofmap, & - get_wtheta_m3x3_dofmap - use funit - - implicit none - - private - public :: hydro_shallow_to_deep_test_type, test_all - - @TestCase - type, extends(TestCase) :: hydro_shallow_to_deep_test_type - private - integer(i_def) :: ndf_w3, undf_w3, ndf_wt, undf_wt - integer(i_def), allocatable :: map_w3(:,:), map_wt(:,:) - real(r_def), allocatable :: basis_w3(:,:,:), basis_wt(:,:,:) - real(r_def), allocatable :: exner_data(:) - real(r_def), allocatable :: temperature_data(:) - real(r_def), allocatable :: coriolis_data(:) - real(r_def), allocatable :: phi_data(:) - real(r_def), allocatable :: moist_dyn_data(:,:) - real(r_def), allocatable :: height_w3_data(:) - real(r_def), allocatable :: height_wth_data(:) - real(r_def), allocatable :: mask_w3_data(:) - contains - procedure setUp - procedure tearDown - procedure test_all - end type hydro_shallow_to_deep_test_type - - real(r_def), parameter :: gravity = 10.0_r_def - real(r_def), parameter :: radius = 6000000.0_r_def - real(r_def), parameter :: omega = 8.0E-5_r_def - real(r_def), parameter :: cp = 1000.0_r_def - real(r_def), parameter :: scaling = 1.0_r_def - integer(i_def), parameter :: eos_level = 2 - -contains - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine setUp( this ) - - use base_mesh_config_mod, only : geometry_planar, & - topology_fully_periodic - use finite_element_config_mod, only : cellshape_quadrilateral - use feign_config_mod, only : feign_base_mesh_config, & - feign_planet_config - - implicit none - - class(hydro_shallow_to_deep_test_type), intent(inout) :: this - - real(r_def), parameter :: dx = 6000.0_r_def, & - dy = 1000.0_r_def, & - dz = 2000.0_r_def - - integer(i_def) :: ncells, icell - integer(i_def) :: dim_space, dim_space_diff - integer(i_def) :: nqp_h, nqp_v - integer(i_def) :: nlayers, i, j, k - - call feign_base_mesh_config & - ( file_prefix='foo', & - prime_mesh_name='unit_test', & - geometry=geometry_planar, & - prepartitioned=.false., & - topology=topology_fully_periodic, & - fplane=.false., & - f_lat_deg=0.0_r_def ) - - nlayers=3 - - call get_wtheta_w3nodal_basis(this%basis_wt) - allocate(this%basis_w3(1,1,1)) - this%basis_w3(:,:,:) = 1.0_r_def - - call get_w3_m3x3_dofmap( this%map_w3, nlayers ) - call get_w3_m3x3_q3x3x3_size( this%ndf_w3, this%undf_w3, ncells, & - dim_space, dim_space_diff, & - nqp_h, nqp_v, & - nlayers ) - - call get_wtheta_m3x3_dofmap( this%map_wt, nlayers ) - call get_wtheta_m3x3_q3x3x3_size( this%ndf_wt, this%undf_wt, ncells, & - dim_space, dim_space_diff, & - nqp_h, nqp_v, & - nlayers ) - - allocate(this%height_w3_data(this%undf_w3)) - icell = 1 - do j = 1,3 - do i = 1,3 - do k = 0,nlayers-1 - this%height_w3_data(this%map_w3(1,icell)+k) = real(k+0.5_r_def)*dz - end do - icell = icell + 1 - end do - end do - - allocate(this%height_wth_data(this%undf_wt)) - icell = 1 - do j = 1,3 - do i = 1,3 - do k = 0,nlayers - this%height_wth_data(this%map_wt(1,icell)+k) = real(k)*dz - end do - icell = icell + 1 - end do - end do - - ! Create the data - allocate(this%exner_data(this%undf_w3)) - allocate(this%phi_data(this%undf_w3)) - allocate(this%temperature_data(this%undf_wt)) - allocate(this%coriolis_data(this%undf_wt)) - allocate(this%moist_dyn_data(this%undf_wt, 3)) - allocate(this%mask_w3_data(this%undf_w3)) - - end subroutine setUp - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine tearDown( this ) - - use configuration_mod, only : final_configuration - - implicit none - - class(hydro_shallow_to_deep_test_type), intent(inout) :: this - - deallocate(this%basis_w3) - deallocate(this%basis_wt) - deallocate(this%map_w3) - deallocate(this%map_wt) - deallocate(this%exner_data) - deallocate(this%phi_data) - deallocate(this%temperature_data) - deallocate(this%coriolis_data) - deallocate(this%moist_dyn_data) - deallocate(this%height_w3_data) - deallocate(this%mask_w3_data) - - call final_configuration() - - end subroutine tearDown - - !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - @test - subroutine test_all( this ) - - use hydro_shallow_to_deep_kernel_mod, only : hydro_shallow_to_deep_code - - implicit none - - class(hydro_shallow_to_deep_test_type), intent(inout) :: this - - real(r_def), parameter :: tol = 1.0e-6_r_def - real(r_def), parameter :: mv = 0.01_r_def ! set 10 g/kg of vapour mixing ratio - real(r_def) :: recip_epsilon - real(r_def) :: answer - integer(i_def) :: nlayers, cell, k - - nlayers=3 - cell=1 - - recip_epsilon = 461.51_r_def/300.0_r_def - ! With u=50m/s, theta=800K, rho=10^-5 - ! coriolis term = 2 \Omega u cos(lat) = 2 x 8e-5 x 50 x cos 45 = 5e-3 - this%exner_data(:) = 1.0_r_def - this%temperature_data(:) = 300.0_r_def - this%coriolis_data(:) = 5.0e-3_r_def - this%mask_w3_data(:) = 1.0_r_def - this%moist_dyn_data(:,1) = 1.0_r_def + recip_epsilon*mv - this%moist_dyn_data(:,2) = 1.0_r_def + mv - this%moist_dyn_data(:,3) = 1.0_r_def - this%phi_data(:) = 1.5_r_def * gravity * this%height_w3_data(:) - call hydro_shallow_to_deep_code( nlayers, & - this%exner_data, & - this%temperature_data, & - this%coriolis_data, & - this%moist_dyn_data(:,1), & - this%moist_dyn_data(:,2), & - this%moist_dyn_data(:,3), & - this%phi_data, & - this%height_w3_data, & - this%height_wth_data, & - this%mask_w3_data, & - gravity, & - cp, & - eos_level, & - this%ndf_w3, this%undf_w3, & - this%map_w3(:,cell), this%basis_w3, & - this%ndf_wt, this%undf_wt, & - this%map_wt(:,cell), this%basis_wt & - ) - - ! Checks hydrostatic plus coriolis balance using g (bottom level) - ! T_moist = T ( 1+ mv * recip_epsilon )/(1 + mv) - ! = 300 * (1 + (461.51/300) *0.01 )/ (1 + 0.01) - ! = 301.59910891089106 - ! exner_{k} = exner_{k+1} * ( cp T - (g - coriolis) * (1-w) ) / - ! ( cp T + (g - coriolis) * w ) - ! = 1 * ( (1000 * 301.59910891089106 ) + (10 - 0.005) * 1000 )/ - ! ( (1000 * 301.59910891089106 ) - (10 - 0.005) * 1000 ) - ! = 1.0685518461130759 - - k=0 - answer = 1.068551846113075_r_def - @assertEqual(answer, this%exner_data(this%map_w3(1, cell)+k), tol) - - ! Checks hydrostatic plus coriolis balance using phi (next level up) - ! exner_{k+1} = exner_{k} * ( cp T + (1.5g - coriolis) * (w) ) / - ! ( cp T - (1.5g - coriolis) * (1-w) ) - ! = 1.0685518461130759 * - ! ( (1000 * 301.59910891089106 ) - (15 - 0.005) * 1000 )/ - ! ( (1000 * 301.59910891089106 ) + (15 - 0.005) * 1000 ) - ! = 0.9673311696602781 - - k=1 - answer = 0.9673311696602781_r_def - @assertEqual(answer, this%exner_data(this%map_w3(1, cell)+k), tol) - - end subroutine test_all - -end module hydro_shallow_to_deep_kernel_mod_test diff --git a/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf b/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf new file mode 100644 index 000000000..e950751f4 --- /dev/null +++ b/science/gungho/unit-test/kernel/initialisation/regrav_geopot_kernel_mod_test.pf @@ -0,0 +1,104 @@ +!----------------------------------------------------------------------------- +! (c) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- +!> Test of interpolation from constant gravity geopotential height to true +!> geopotentail height. +!> +module regrav_geopot_kernel_mod_test + +use constants_mod, only: i_def, r_def +use, intrinsic :: iso_fortran_env, only: real64 +use funit + +implicit none + +private +public :: regrav_geopot_test_type, test_all + +@TestCase +type, extends(TestCase) :: regrav_geopot_test_type + private +contains + procedure setUp + procedure tearDown + procedure test_all +end type regrav_geopot_test_type + +integer(i_def), parameter :: nlayers = 3_i_def +integer(i_def), parameter :: ndf_wt = 2_i_def +integer(i_def), parameter :: undf_wt = nlayers + 1_i_def +real(r_def), parameter :: height_domain = 40.0e3_r_def +real(r_def), parameter :: T_surf = 300.0_r_def +real(r_def), parameter :: dTdz = -2.0e-3_r_def +real(r_def), parameter :: planet_radius = 6.4e6_r_def + +integer(i_def), allocatable :: map_wt(:) +real(r_def), allocatable :: temperature(:) +real(r_def), allocatable :: height_wt(:) +real(r_def), allocatable :: height_phys(:) + +contains + + subroutine setUp( this ) + + implicit none + + class(regrav_geopot_test_type), intent(inout) :: this + + integer(i_def) :: k + + allocate( map_wt(ndf_wt) ) + allocate( temperature(undf_wt) ) + allocate( height_wt(undf_wt), height_phys(undf_wt) ) + + ! Create the model grid + map_wt(:) = (/ 1_i_def, 2_i_def /) + + do k = 0, nlayers + height_wt(map_wt(1) + k) = k * height_domain / nlayers + height_phys(map_wt(1) + k) = height_wt(map_wt(1) + k) / & + ( 1.0_r_def - height_wt(map_wt(1) + k) / & + planet_radius) + temperature(map_wt(1) + k) = T_surf + dTdz * height_phys(map_wt(1) + k) + end do + + end subroutine setUp + + subroutine tearDown( this ) + + implicit none + + class(regrav_geopot_test_type), intent(inout) :: this + + deallocate( map_wt ) + deallocate( temperature ) + deallocate( height_wt, height_phys ) + + end subroutine tearDown + + @Test + subroutine test_all( this ) + + use regrav_geopot_kernel_mod, only: regrav_geopot_code + + implicit none + + class(regrav_geopot_test_type), intent(inout) :: this + + integer(i_def) :: k + real(r_def) :: target_value, tolerance + + call regrav_geopot_code( nlayers, temperature, height_wt, planet_radius, & + ndf_wt, undf_wt, map_wt ) + + do k = 0, nlayers + target_value = T_surf + height_wt(map_wt(1)+k) * dTdz + tolerance = 10.0_r_def * spacing(target_value) + @assertEqual(target_value, temperature(map_wt(1)+k), tolerance) + end do + + end subroutine test_all + +end module regrav_geopot_kernel_mod_test