diff --git a/applications/gungho_model/example/configuration.nml b/applications/gungho_model/example/configuration.nml
index ea60b8aa1..39b5d77f4 100644
--- a/applications/gungho_model/example/configuration.nml
+++ b/applications/gungho_model/example/configuration.nml
@@ -71,6 +71,7 @@ rotating=.true.,
shallow=.false.,
si_momentum_equation=.false.,
theta_moist_source=.false.
+traditional_coriolis=.true.
use_multires_coupling=.false.,
use_physics=.false.,
use_wavedynamics=.true.,
diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml
index 007ee86df..78813faa4 100644
--- a/applications/lfric_atm/metadata/field_def_diags.xml
+++ b/applications/lfric_atm/metadata/field_def_diags.xml
@@ -20,6 +20,7 @@
+
@@ -43,6 +44,7 @@
+
zoomed__orography
diff --git a/applications/lfric_atm/metadata/field_def_initial_diags.xml b/applications/lfric_atm/metadata/field_def_initial_diags.xml
index 707675af2..1c0b905f9 100644
--- a/applications/lfric_atm/metadata/field_def_initial_diags.xml
+++ b/applications/lfric_atm/metadata/field_def_initial_diags.xml
@@ -9,6 +9,7 @@
+
@@ -19,8 +20,10 @@
+
+
diff --git a/applications/lfric_atm/metadata/lfric_dictionary.xml b/applications/lfric_atm/metadata/lfric_dictionary.xml
index aa5a1efc5..5de27f625 100644
--- a/applications/lfric_atm/metadata/lfric_dictionary.xml
+++ b/applications/lfric_atm/metadata/lfric_dictionary.xml
@@ -12,6 +12,7 @@
+
diff --git a/dependencies.yaml b/dependencies.yaml
index caab36487..3e284764a 100644
--- a/dependencies.yaml
+++ b/dependencies.yaml
@@ -30,8 +30,8 @@ lfric_apps:
ref:
lfric_core:
- source: git@github.com:MetOffice/lfric_core.git
- ref: 2025.12.1
+ source: git@github.com:tommbendall/lfric_core.git
+ ref: TBendall/vn3.0_dcmip
moci:
source: git@github.com:MetOffice/moci.git
diff --git a/rose-stem/app/gungho_model/file/file_def_diags_nwp.xml b/rose-stem/app/gungho_model/file/file_def_diags_nwp.xml
index a9cfe35e5..80903075e 100644
--- a/rose-stem/app/gungho_model/file/file_def_diags_nwp.xml
+++ b/rose-stem/app/gungho_model/file/file_def_diags_nwp.xml
@@ -5,12 +5,16 @@
+
+
+
+
diff --git a/rose-stem/app/gungho_model/file/file_def_initial.xml b/rose-stem/app/gungho_model/file/file_def_initial.xml
index 2b1050253..291907dc8 100644
--- a/rose-stem/app/gungho_model/file/file_def_initial.xml
+++ b/rose-stem/app/gungho_model/file/file_def_initial.xml
@@ -11,8 +11,10 @@
+
+
diff --git a/rose-stem/app/gungho_model/opt/rose-app-breaking_gw.conf b/rose-stem/app/gungho_model/opt/rose-app-breaking_gw.conf
new file mode 100644
index 000000000..51838223b
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-breaking_gw.conf
@@ -0,0 +1,46 @@
+[namelist:external_forcing]
+hs_random=.false.
+wind_forcing='rayleigh_friction'
+
+[namelist:extrusion]
+domain_height=40000.0
+method='uniform'
+number_of_layers=88
+
+[namelist:formulation]
+shallow=.false.
+use_physics=.true.
+
+[namelist:idealised]
+test='breaking_gw'
+
+[namelist:initial_wind]
+profile='breaking_gw'
+
+[namelist:initialization]
+zero_w2v_wind=.true.
+
+[namelist:io]
+diagnostic_frequency=48
+
+[namelist:orography]
+orog_init_option='analytic'
+profile='bell'
+
+[namelist:orography_bell_spherical]
+lambda_centre1=72.0
+lambda_centre2=140.0
+mountain_height=3000.0
+phi_centre1=45.0
+phi_centre2=45.0
+radius_lat=20.0
+radius_lon=3.5
+
+[namelist:physics]
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:section_choice]
+external_forcing=.true.
diff --git a/rose-stem/app/gungho_model/opt/rose-app-mountain_gap.conf b/rose-stem/app/gungho_model/opt/rose-app-mountain_gap.conf
new file mode 100644
index 000000000..e4ec0a846
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-mountain_gap.conf
@@ -0,0 +1,55 @@
+[namelist:damping_layer]
+dl_base=10000.0
+dl_str=0.15
+dl_type='standard'
+
+[namelist:external_forcing]
+hs_random=.false.
+wind_forcing='rayleigh_friction'
+
+[namelist:extrusion]
+domain_height=20000.0
+method='dcmip_mountain'
+number_of_layers=57
+
+[namelist:formulation]
+dlayer_on=.false.
+shallow=.true.
+traditional_coriolis=.true.
+use_physics=.true.
+
+[namelist:idealised]
+test='horizontal_mountain'
+
+[namelist:initial_pressure]
+method='balanced'
+
+[namelist:initial_wind]
+profile='horizontal_mountain'
+
+[namelist:initialization]
+zero_w2v_wind=.true.
+
+[namelist:io]
+diagnostic_frequency=120
+
+[namelist:orography]
+orog_init_option='analytic'
+profile='dcmip_gap'
+
+[namelist:physics]
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:planet]
+cp=1004.64
+gravity=9.80616
+omega=7.29211E-5
+p_zero=100000.0
+rd=287.04
+scaling_factor=20.0
+
+[namelist:section_choice]
+external_forcing=.true.
diff --git a/rose-stem/app/gungho_model/opt/rose-app-mountain_gap_fail.conf b/rose-stem/app/gungho_model/opt/rose-app-mountain_gap_fail.conf
new file mode 100644
index 000000000..7db7980c6
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-mountain_gap_fail.conf
@@ -0,0 +1,50 @@
+[namelist:external_forcing]
+hs_random=.false.
+wind_forcing='rayleigh_friction'
+
+[namelist:extrusion]
+domain_height=20000.0
+method='dcmip_mountain'
+number_of_layers=57
+
+[namelist:formulation]
+dlayer_on=.false.
+shallow=.true.
+traditional_coriolis=.true.
+use_physics=.true.
+
+[namelist:idealised]
+test='horizontal_mountain'
+
+[namelist:initial_pressure]
+method='balanced'
+
+[namelist:initial_wind]
+profile='horizontal_mountain'
+
+[namelist:initialization]
+zero_w2v_wind=.true.
+
+[namelist:io]
+diagnostic_frequency=160
+
+[namelist:orography]
+orog_init_option='analytic'
+profile='dcmip_gap'
+
+[namelist:physics]
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:planet]
+cp=1004.64
+gravity=9.80616
+omega=7.29211E-5
+p_zero=100000.0
+rd=287.04
+scaling_factor=20.0
+
+[namelist:section_choice]
+external_forcing=.true.
diff --git a/rose-stem/app/gungho_model/opt/rose-app-mountain_gap_no_rotation.conf b/rose-stem/app/gungho_model/opt/rose-app-mountain_gap_no_rotation.conf
new file mode 100644
index 000000000..471b307d7
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-mountain_gap_no_rotation.conf
@@ -0,0 +1,50 @@
+[namelist:external_forcing]
+hs_random=.false.
+wind_forcing='rayleigh_friction'
+
+[namelist:extrusion]
+domain_height=20000.0
+method='dcmip_mountain'
+number_of_layers=57
+
+[namelist:formulation]
+dlayer_on=.false.
+rotating=.false.
+shallow=.true.
+use_physics=.true.
+
+[namelist:idealised]
+test='horizontal_mountain'
+
+[namelist:initial_pressure]
+method='balanced'
+
+[namelist:initial_wind]
+profile='horizontal_mountain'
+
+[namelist:initialization]
+zero_w2v_wind=.true.
+
+[namelist:io]
+diagnostic_frequency=160
+
+[namelist:orography]
+orog_init_option='analytic'
+profile='dcmip_gap'
+
+[namelist:physics]
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:planet]
+cp=1004.64
+gravity=9.80616
+omega=0.0
+p_zero=100000.0
+rd=287.04
+scaling_factor=20.0
+
+[namelist:section_choice]
+external_forcing=.true.
diff --git a/rose-stem/app/gungho_model/opt/rose-app-mountain_vortex.conf b/rose-stem/app/gungho_model/opt/rose-app-mountain_vortex.conf
new file mode 100644
index 000000000..01f7338eb
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-mountain_vortex.conf
@@ -0,0 +1,50 @@
+[namelist:external_forcing]
+hs_random=.false.
+wind_forcing='rayleigh_friction'
+
+[namelist:extrusion]
+domain_height=20000.0
+method='dcmip_mountain'
+number_of_layers=57
+
+[namelist:formulation]
+dlayer_on=.false.
+shallow=.true.
+traditional_coriolis=.true.
+use_physics=.true.
+
+[namelist:idealised]
+test='horizontal_mountain'
+
+[namelist:initial_pressure]
+method='balanced'
+
+[namelist:initial_wind]
+profile='horizontal_mountain'
+
+[namelist:initialization]
+zero_w2v_wind=.true.
+
+[namelist:io]
+diagnostic_frequency=160
+
+[namelist:orography]
+orog_init_option='analytic'
+profile='dcmip_vortex'
+
+[namelist:physics]
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:planet]
+cp=1004.64
+gravity=9.80616
+omega=7.29211E-5
+p_zero=100000.0
+rd=287.04
+scaling_factor=20.0
+
+[namelist:section_choice]
+external_forcing=.true.
diff --git a/rose-stem/app/gungho_model/opt/rose-app-mountain_vortex_no_rotation.conf b/rose-stem/app/gungho_model/opt/rose-app-mountain_vortex_no_rotation.conf
new file mode 100644
index 000000000..c79b84106
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-mountain_vortex_no_rotation.conf
@@ -0,0 +1,50 @@
+[namelist:external_forcing]
+hs_random=.false.
+wind_forcing='rayleigh_friction'
+
+[namelist:extrusion]
+domain_height=20000.0
+method='dcmip_mountain'
+number_of_layers=57
+
+[namelist:formulation]
+dlayer_on=.false.
+rotating=.false.
+shallow=.true.
+use_physics=.true.
+
+[namelist:idealised]
+test='horizontal_mountain'
+
+[namelist:initial_pressure]
+method='balanced'
+
+[namelist:initial_wind]
+profile='horizontal_mountain'
+
+[namelist:initialization]
+zero_w2v_wind=.true.
+
+[namelist:io]
+diagnostic_frequency=160
+
+[namelist:orography]
+orog_init_option='analytic'
+profile='dcmip_vortex'
+
+[namelist:physics]
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:planet]
+cp=1004.64
+gravity=9.80616
+omega=0.0
+p_zero=100000.0
+rd=287.04
+scaling_factor=20.0
+
+[namelist:section_choice]
+external_forcing=.true.
diff --git a/rose-stem/app/gungho_model/opt/rose-app-squall_line.conf b/rose-stem/app/gungho_model/opt/rose-app-squall_line.conf
new file mode 100644
index 000000000..32716a324
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-squall_line.conf
@@ -0,0 +1,55 @@
+[namelist:damping_layer]
+dl_base=18000.0
+dl_str=0.15
+dl_type='standard'
+
+[namelist:extrusion]
+domain_height=20000.0
+method='uniform'
+number_of_layers=40
+
+[namelist:formulation]
+dlayer_on=.true.
+moisture_formulation='traditional'
+rotating=.false.
+shallow=.true.
+theta_moist_source=.true.
+traditional_coriolis=.true.
+use_physics=.true.
+
+[namelist:idealised]
+test='squall_line'
+
+[namelist:initial_wind]
+profile='squall_line'
+
+[namelist:io]
+diagnostic_frequency=30
+
+[namelist:physics]
+evap_condense_placement='fast'
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:planet]
+cp=1004.64
+gravity=9.80616
+omega=7.29211E-5
+p_zero=100000.0
+rd=287.04
+scaling_factor=60.0
+
+[namelist:section_choice]
+aerosol='none'
+boundary_layer='none'
+chemistry='none'
+cloud='kessler'
+methane_oxidation=.false.
+microphysics='none'
+orographic_drag='none'
+radiation='none'
+spectral_gwd='none'
+stochastic_physics='none'
+surface='none'
diff --git a/rose-stem/app/gungho_model/opt/rose-app-supercell.conf b/rose-stem/app/gungho_model/opt/rose-app-supercell.conf
new file mode 100644
index 000000000..39cfa14f0
--- /dev/null
+++ b/rose-stem/app/gungho_model/opt/rose-app-supercell.conf
@@ -0,0 +1,55 @@
+[namelist:damping_layer]
+dl_base=18000.0
+dl_str=0.15
+dl_type='standard'
+
+[namelist:extrusion]
+domain_height=20000.0
+method='uniform'
+number_of_layers=40
+
+[namelist:formulation]
+dlayer_on=.true.
+moisture_formulation='traditional'
+rotating=.false.
+shallow=.true.
+theta_moist_source=.true.
+traditional_coriolis=.true.
+use_physics=.true.
+
+[namelist:idealised]
+test='supercell'
+
+[namelist:initial_wind]
+profile='supercell'
+
+[namelist:io]
+diagnostic_frequency=60
+
+[namelist:physics]
+evap_condense_placement='fast'
+limit_drag_incs=.false.
+sample_physics_scalars=.true.
+sample_physics_winds=.true.
+sample_physics_winds_correction=.false.
+
+[namelist:planet]
+cp=1004.64
+gravity=9.80616
+omega=7.29211E-5
+p_zero=100000.0
+rd=287.04
+scaling_factor=120.0
+
+[namelist:section_choice]
+aerosol='none'
+boundary_layer='none'
+chemistry='none'
+cloud='kessler'
+methane_oxidation=.false.
+microphysics='none'
+orographic_drag='none'
+radiation='none'
+spectral_gwd='none'
+stochastic_physics='none'
+surface='none'
diff --git a/rose-stem/app/gungho_model/rose-app.conf b/rose-stem/app/gungho_model/rose-app.conf
index d4bc76fa1..cfd48f276 100644
--- a/rose-stem/app/gungho_model/rose-app.conf
+++ b/rose-stem/app/gungho_model/rose-app.conf
@@ -413,6 +413,7 @@ rotating=.true.
shallow=.true.
si_momentum_equation=.false.
!!theta_moist_source=.true.
+traditional_coriolis=.false.
use_multires_coupling=.false.
use_physics=.false.
use_wavedynamics=.true.
diff --git a/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc b/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc
index 17c6843c5..7d825bd37 100644
--- a/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc
+++ b/rose-stem/site/common/gungho_model/tasks_gungho_model.cylc
@@ -112,6 +112,20 @@
"kgo_checks": [],
}) %}
+{% elif task_ns.conf_name == "breaking_gw-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["breaking_gw"],
+ "resolution": "C96_MG",
+ "DT": 450,
+ "mpi_parts": 216,
+ "tsteps": 1920,
+ "wallclock": 60,
+ "use_xios": true,
+ "log_level": "info",
+ "kgo_checks": [],
+ }) %}
+
{% elif task_ns.conf_name == "bryan_fritsch-dry-BiP200x10-100x100" %}
{% do task_dict.update({
@@ -207,6 +221,131 @@
"plot_str": "sbr.py $NODAL_DATA_DIR/lfric_diag.nc $PLOT_DIR um70 30 71 lfric"
}) %}
+{% elif task_ns.conf_name == "mountain_gap-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_gap"],
+ "resolution": "C96_MG",
+ "DT": 120,
+ "mpi_parts": 216,
+ "tsteps": 960,
+ "wallclock": 60,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_gap-C192_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_gap"],
+ "resolution": "C192_MG",
+ "DT": 60,
+ "mpi_parts": 864,
+ "tsteps": 720,
+ "wallclock": 60,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_gap_fail-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_gap_fail"],
+ "resolution": "C96_MG",
+ "DT": 120,
+ "mpi_parts": 216,
+ "tsteps": 960,
+ "wallclock": 60,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_gap_no_rotation-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_gap_no_rotation"],
+ "resolution": "C96_MG",
+ "DT": 45,
+ "mpi_parts": 216,
+ "tsteps": 960,
+ "wallclock": 60,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_gap_no_rotation-C192_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_gap_no_rotation"],
+ "resolution": "C192_MG",
+ "DT": 45,
+ "mpi_parts": 864,
+ "tsteps": 960,
+ "wallclock": 60,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_vortex-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_vortex"],
+ "resolution": "C96_MG",
+ "DT": 45,
+ "mpi_parts": 216,
+ "tsteps": 1920,
+ "wallclock": 60,
+ "crun": 1,
+ "crun_compare": false,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_vortex-C192_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_vortex"],
+ "resolution": "C192_MG",
+ "DT": 45,
+ "mpi_parts": 864,
+ "tsteps": 1920,
+ "wallclock": 60,
+ "crun": 1,
+ "crun_compare": false,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_vortex_no_rotation-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_vortex_no_rotation"],
+ "resolution": "C96_MG",
+ "DT": 45,
+ "mpi_parts": 216,
+ "tsteps": 1920,
+ "wallclock": 60,
+ "crun": 1,
+ "crun_compare": false,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "mountain_vortex_no_rotation-C192_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["mountain_vortex_no_rotation"],
+ "resolution": "C192_MG",
+ "DT": 45,
+ "mpi_parts": 864,
+ "tsteps": 1920,
+ "wallclock": 60,
+ "crun": 1,
+ "crun_compare": false,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
{% elif task_ns.conf_name == "no-timestep-method-C12" %}
{% do task_dict.update({
@@ -451,6 +590,32 @@
"plot_str": "gw_high_order_cart_plot.py diagGungho $NODAL_DATA_DIR 300 8 theta T000000:T000100:T000200:T000300 $PLOT_DIR"
}) %}
+{% elif task_ns.conf_name == "squall_line-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["squall_line"],
+ "resolution": "C96_MG",
+ "DT": 20,
+ "mpi_parts": 216,
+ "tsteps": 540,
+ "wallclock": 60,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
+{% elif task_ns.conf_name == "supercell-C96_MG" %}
+
+ {% do task_dict.update({
+ "opt_confs": ["supercell"],
+ "resolution": "C96_MG",
+ "DT": 10,
+ "mpi_parts": 216,
+ "tsteps": 720,
+ "wallclock": 60,
+ "use_xios": true,
+ "kgo_checks": [],
+ }) %}
+
{% elif task_ns.conf_name == "straka_200m-BiP256x8-200x200" %}
{% do task_dict.update({
diff --git a/rose-stem/site/meto/groups/groups_gungho_model.cylc b/rose-stem/site/meto/groups/groups_gungho_model.cylc
index 0605786b9..4a6351f13 100644
--- a/rose-stem/site/meto/groups/groups_gungho_model.cylc
+++ b/rose-stem/site/meto/groups/groups_gungho_model.cylc
@@ -111,6 +111,20 @@
"gungho_model_tidally-locked-earth-C24_MG_ex1a_gnu_fast-debug-64bit",
"gungho_model_ex1a_canned",
],
+ "gungho_model_ex1a_dcmip":[
+ "gungho_model_breaking_gw-C96_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_mountain_gap-C192_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_mountain_gap_fail-C96_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_mountain_gap_no_rotation-C192_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_mountain_vortex-C192_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_mountain_vortex_no_rotation-C192_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_squall_line-C96_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_supercell-C96_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ ],
+ "gungho_model_ex1a_squall":[
+ "gungho_model_squall_line-C96_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ "gungho_model_supercell-C96_MG_ex1a_gnu_fast-debug-64bit-rtran32",
+ ],
"gungho_model_ex1a_C256_performance_perftools": [
"gungho_model_gh_profile_omp1_3nodes-C256_MG_ex1a_perftools-gnu_production-64bit",
"gungho_model_gh_profile_omp1_6nodes-C256_MG_ex1a_perftools-gnu_production-64bit",
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 15f3d9684..80d6f63f1 100644
--- a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf
+++ b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf
@@ -1126,8 +1126,9 @@ trigger=namelist:external_forcing=filtering_order: this == "'filtering'";
=namelist:external_forcing=diffusion_order: this == "'diffusion'";
=namelist:external_forcing=diffusion_coefficient: this == "'diffusion'";
=namelist:wind_forcing: this == "'profile'";
-value-titles=None, Held-Suarez, Filtering, Diffusion, Profile
-values='none','held_suarez', 'filtering', 'diffusion', 'profile'
+value-titles=None, Held-Suarez, Filtering, Diffusion, Profile, Rayleigh Friction
+ =Breaking Gravity Wave, Horizontal DCMIP Mountains
+values='none','held_suarez', 'filtering', 'diffusion', 'profile', 'rayleigh_friction'
[namelist:files]
compulsory=true
@@ -1997,6 +1998,16 @@ help=Adds a pointwise Crank-Nicolson discretisation of the moist source term
sort-key=Panel-A02
type=logical
+[namelist:formulation=traditional_coriolis]
+compulsory=true
+description=Whether the traditional approximation is used for the Coriolis terms
+help=If true, the Coriolis force is neglected in the vertical velocity equation,
+ =and contributions from the vertical velocity to the Coriolis force are
+ =neglected in the horizontal velocity equations.
+!kind=default
+sort-key=Panel-A11
+type=logical
+
[namelist:formulation=use_multires_coupling]
compulsory=true
description=Enable use of multiple resolution meshes for coupling
@@ -2421,6 +2432,10 @@ fail-if=this == "'yz_cosine_hill'" and namelist:base_mesh=geometry == "
=this == "'bryan_fritsch'" and namelist:base_mesh=geometry == "'spherical'";
=this == "'grabowski_clark'" and namelist:base_mesh=geometry == "'spherical'";
=this == "'grabowski_clark'" and namelist:formulation=moisture_formulation == "'dry'";
+ =this == "'breaking_gw'" and namelist:base_mesh=geometry == "'planar'";
+ =this == "'horizontal_mountain'" and namelist:base_mesh=geometry == "'planar'";
+ =this == "'squall_line'" and namelist:base_mesh=geometry == "'planar'";
+ =this == "'supercell'" and namelist:base_mesh=geometry == "'planar'";
help=gravity_wave : Gravity wave test (planar|spherical)
=isot_atm : Isothermal atmosphere
=isot_cold_atm : Isothermal atmosphere, cold ground temperature
@@ -2454,6 +2469,10 @@ help=gravity_wave : Gravity wave test (planar|spherical)
=specified_profiles : Use profile data provided by initial namelists
=bryan_fritsch : Rising bubble test of Bryan and Fritsch (2002)
=grabowski_clark : Rising bubble test of Grabowski and Clark (1991)
+ =breaking_gw : Breaking gravity wave test from DCMIP 2025
+ =horizontal_mountain : Horizontal mountain flow from DCMIP 2025
+ =squall_line : Squall line test from DCMIP 2025
+ =supercell : Supercell test from DCMIP 2016
sort-key=Panel-A01
trigger=namelist:initial_temperature=profile_data: 'specified_profiles';
=namelist:initial_temperature=profile_heights: 'specified_profiles';
@@ -2472,6 +2491,7 @@ value-titles=None, Gravity Wave, Isothermal atmosphere, Cold ground temperature,
=Uniform translation, Vertical solid body rotation, Cylinder for vertical slice,
=Specified profiles, Bryan and Fritsch (2002) rising bubble,
=Grabowski and Clark (1991) moist rising bubble in an unsaturated atmosphere
+ =Breaking gravity wave, Horizontal mountain flow, Squall line, Supercell
# initial temp namelists
# XS
# trigger=namelist:initial_density: 'none' ;
@@ -2484,7 +2504,8 @@ values='none', 'gravity_wave', 'isot_atm', 'isot_cold_atm', 'isot_dry_atm', 'col
='solid_body_rotation_alt', 'deep_baroclinic_wave', 'isentropic',
='cos_phi','cosine_bubble',
='eternal_fountain', 'curl_free_reversible', 'div_free_reversible', 'translational',
- ='rotational','vertical_cylinder', 'specified_profiles', 'bryan_fritsch', 'grabowski_clark'
+ ='rotational','vertical_cylinder', 'specified_profiles', 'bryan_fritsch', 'grabowski_clark',
+ ='breaking_gw', 'horizontal_mountain', 'squall_line', 'supercell'
#==============================================================================
# INITIAL DENSITY
@@ -2844,6 +2865,7 @@ value-titles=None, Solid body rotation, Solid body rotation-alt, Constant uv, Co
=Solid body rotation with vertical wind, DCMIP 2012 Test 1.1, Vertical Deformational Flow,
=Solid body rotation around sphere in four parts,
=Piecewise linear vertical profile
+ =Breaking gravity wave, Horizontal mountain, Squall line, Supercell
# PLEASE add coments to indicate which of these values should allowed given
# the value of namelist:idealised=test
# the value of the mesh geometry
@@ -2852,7 +2874,8 @@ values='none', 'solid_body_rotation', 'solid_body_rotation_alt', 'constant_uv',
='NL_case_1', 'NL_case_2', 'NL_case_3', 'NL_case_4', 'xy_NL_case_1', 'yz_NL_case_1',
='hadley_like_dcmip', 'vortex', 'sbr_streamfunction', 'dcmip301_streamfunction', 'eternal_fountain',
='curl_free_reversible', 'div_free_reversible', 'rotational',
- ='sbr_with_vertical','dcmip_101', 'vertical_deformation', 'four_part_sbr', 'pw_linear'
+ ='sbr_with_vertical','dcmip_101', 'vertical_deformation', 'four_part_sbr', 'pw_linear',
+ ='breaking_gw', 'horizontal_mountain', 'squall_line', 'supercell'
[namelist:initial_wind=profile_data_u]
!bounds=namelist:initial_wind=profile_size_uv
@@ -4086,11 +4109,14 @@ description=orography profile
!enumeration=true
fail-if=this == "'dcmip200'" and namelist:base_mesh=geometry != "'spherical'" ;
=this == "'none'" and namelist:orography=orog_init_option == "'analytic'" ;
-help=none : No orography (flat surface)
- =schar : Schar mountain profile
- =agnesi : Witch-of-Agnesi mountain profile
- =dcmip200 : DCMIP 2.0.0 mountain orography
- =bell : Bell-shaped mountain profile
+help=none : No orography (flat surface)
+ =schar : Schar mountain profile
+ =agnesi : Witch-of-Agnesi mountain profile
+ =dcmip200 : DCMIP 2.0.0 mountain orography
+ =bell : Bell-shaped mountain profile
+ =dcmip_breaking_gw : DCMIP Breaking gravity wave mountain
+ =dcmip_gap : DCMIP Mountain gap
+ =dcmip_vortex : DCMIP vortex shedding mountain
sort-key=Panel-A02
trigger=namelist:orography_agnesi_cartesian: 'agnesi' ;
=namelist:orography_agnesi_spherical: 'agnesi' ;
@@ -4099,8 +4125,9 @@ trigger=namelist:orography_agnesi_cartesian: 'agnesi' ;
=namelist:orography_dcmip200_spherical: 'dcmip200' ;
=namelist:orography_schar_cartesian: 'schar' ;
=namelist:orography_schar_spherical: 'schar' ;
-value-titles=None, Agnesi, Bell-shaped, Schar, DCMIP 2.0.0
-values='none', 'agnesi', 'bell', 'schar', 'dcmip200'
+value-titles=None, Agnesi, Bell-shaped, Schar, DCMIP 2.0.0,
+ =DCMIP Breaking GW, DCMIP Gap, DCMIP Vortex
+values='none', 'agnesi', 'bell', 'schar', 'dcmip200', 'dcmip_breaking_gw', 'dcmip_gap', 'dcmip_vortex'
#==============================================================================
# OROGRAPHY (AGNESI CARTESIAN)
@@ -5109,8 +5136,8 @@ sort-key=Panel-A03
trigger=namelist:cloud: this == "'um'" ;
=namelist:physics=evap_condense_placement: this == "'evap_condense'" ;
=namelist:section_choice=microphysics: this != "'none'" ;
-value-titles=None, Unified Model, Evaporation-Condensation
-values='none', 'um', 'evap_condense'
+value-titles=None, Unified Model, Evaporation-Condensation, Kessler
+values='none', 'um', 'evap_condense', 'kessler'
[namelist:section_choice=convection]
compulsory=true
@@ -5222,8 +5249,8 @@ sort-key=Panel-A06
trigger=namelist:physics=microphysics_placement: this != "'none'" ;
=namelist:section_choice=electric: this!= "'none'";
=namelist:microphysics: this == "'um'" ;
-value-titles=None, Unified Model
-values='none', 'um'
+value-titles=None, Unified Model, Kessler
+values='none', 'um', 'kessler'
[namelist:section_choice=orographic_drag]
compulsory=true
diff --git a/science/gungho/source/algorithm/diagnostics/diagnostic_alg_mod.x90 b/science/gungho/source/algorithm/diagnostics/diagnostic_alg_mod.x90
index 402ac004b..2e3345af8 100644
--- a/science/gungho/source/algorithm/diagnostics/diagnostic_alg_mod.x90
+++ b/science/gungho/source/algorithm/diagnostics/diagnostic_alg_mod.x90
@@ -36,6 +36,8 @@ module diagnostic_alg_mod
LOG_LEVEL_INFO, &
LOG_LEVEL_ERROR
use initialise_diagnostics_mod, only: init_diag => init_diagnostic_field
+ use mr_indices_mod, only: nummr, imr_v, imr_cl, imr_r, &
+ mr_names
implicit none
@@ -51,6 +53,7 @@ module diagnostic_alg_mod
public :: column_total_diagnostics_alg
public :: calc_wbig_diagnostic_alg
public :: pressure_diag_alg
+ public :: dbz_diagnostic_alg
contains
@@ -96,6 +99,37 @@ contains
end subroutine divergence_diagnostic_alg
+ !> @details Calculates the l2 error norm for wind divergence
+ !> @param[in,out] divergence Divergence field
+ !> @param[in,out] l2 l2 norm
+ !> @param[in] u 3D wind field
+ !> @param[in] mesh Mesh
+ subroutine dbz_diagnostic_alg(dbz, exner_in_wth, theta, rho, mr, mesh)
+ use dbz_kernel_mod, only: dbz_kernel_type
+ use planet_config_mod, only: recip_epsilon, kappa, p_zero
+ use physics_mappings_alg_mod, only: map_physics_scalars
+
+ implicit none
+ type(field_type), intent(inout) :: dbz
+ type(field_type), intent(in) :: exner_in_wth
+ type(field_type), intent(in) :: theta
+ type(field_type), intent(in) :: rho
+ type(field_type), intent(in) :: mr(nummr)
+ type(mesh_type), intent(in), pointer :: mesh
+
+ type(field_type) :: rho_in_wth
+
+ call rho_in_wth%initialise( theta%get_function_space() )
+ call dbz%initialise( theta%get_function_space() )
+
+ call map_physics_scalars(rho_in_wth, rho)
+
+ call invoke( dbz_kernel_type( dbz, exner_in_Wth, theta, rho_in_wth, &
+ mr(imr_v), mr(imr_cl), mr(imr_r), &
+ recip_epsilon, kappa, p_zero ) )
+
+ end subroutine dbz_diagnostic_alg
+
!=============================================================================!
!> @details An algorithm for calculating the l2 error for the density field.
!> It is capable of calculating the error at each timestep by using
@@ -592,7 +626,6 @@ contains
use extrusion_config_mod, only: planet_radius
use sci_geometric_constants_mod, only: get_height_fe, get_height_fv
use io_config_mod, only: write_diag, use_xios_io
- use mr_indices_mod, only: nummr, mr_names
use function_space_mod, only: function_space_type
use physics_mappings_alg_mod, only: map_physics_scalars
use planet_config_mod, only: cv, gravity
diff --git a/science/gungho/source/algorithm/initialisation/init_gungho_prognostics_alg_mod.x90 b/science/gungho/source/algorithm/initialisation/init_gungho_prognostics_alg_mod.x90
index 6244df745..399bce38a 100644
--- a/science/gungho/source/algorithm/initialisation/init_gungho_prognostics_alg_mod.x90
+++ b/science/gungho/source/algorithm/initialisation/init_gungho_prognostics_alg_mod.x90
@@ -12,7 +12,7 @@ module init_gungho_prognostics_alg_mod
LOG_LEVEL_INFO
use constants_mod, only: r_def, i_def, radians_to_degrees
use domain_mod, only: domain_type
- use extrusion_config_mod, only: planet_radius
+ use planet_config_mod, only: scaled_radius
use sci_mass_matrix_solver_alg_mod, only: mass_matrix_solver_alg
use sample_init_u_alg_mod, only: sample_init_u_alg
use sci_geometric_constants_mod, only: get_coordinates, &
@@ -30,9 +30,12 @@ module init_gungho_prognostics_alg_mod
eos_method_projected, &
moisture_formulation, &
moisture_formulation_dry
- use idealised_config_mod, only: test, &
- test_bryan_fritsch, &
- test_grabowski_clark
+ use idealised_config_mod, only: test, &
+ test_bryan_fritsch, &
+ test_grabowski_clark, &
+ test_squall_line, &
+ test_supercell, &
+ test_horizontal_mountain
use initial_pressure_config_mod, only: method, &
method_sampled, &
method_projected, &
@@ -51,6 +54,9 @@ module init_gungho_prognostics_alg_mod
use planet_config_mod, only: gravity, cp, rd, p_zero, kappa, &
recip_epsilon
+ use log_mod, only: LOG_LEVEL_INFO
+ use sci_field_minmax_alg_mod, only: log_field_minmax
+
! PsyKAl-lite kernels
use sci_field_bundle_builtins_mod, only: set_bundle_scalar
@@ -66,6 +72,7 @@ module init_gungho_prognostics_alg_mod
use project_eos_rho_kernel_mod, only: project_eos_rho_kernel_type
use sample_eos_rho_kernel_mod, only: sample_eos_rho_kernel_type
use hydrostatic_exner_kernel_mod, only: hydrostatic_exner_kernel_type
+ use hydrostatic_exner_squall_kernel_mod, only: hydrostatic_exner_squall_kernel_type
use initial_mr_kernel_mod, only: initial_mr_kernel_type
use initial_streamfunc_kernel_mod, only: initial_streamfunc_kernel_type
use strong_curl_kernel_mod, only: strong_curl_kernel_type
@@ -90,6 +97,7 @@ module init_gungho_prognostics_alg_mod
public :: init_u_field
public :: init_theta_field
public :: init_exner_field
+ public :: init_exner_squall
public :: init_rho_field
public :: init_mr_fields
@@ -231,12 +239,14 @@ contains
integer(kind=i_def) :: mesh_id
type(field_type), pointer :: chi(:) => null()
type(field_type), pointer :: panel_id => null()
+ type(field_type), pointer :: height_wt => null()
mesh_id = theta%get_mesh_id()
chi => get_coordinates(mesh_id)
panel_id => get_panel_id(mesh_id)
+ height_wt => get_height_fe(Wtheta, mesh_id)
- call invoke( initial_theta_kernel_type( theta, chi, panel_id ) )
+ call invoke( initial_theta_kernel_type( theta, height_wt, chi, panel_id ) )
nullify( chi, panel_id )
@@ -250,6 +260,8 @@ contains
!> @param[in] time time, used for setting analytic profiles
subroutine init_exner_field(exner, theta, moist_dyn, time)
+ use physics_mappings_alg_mod, only: map_physics_scalars
+
implicit none
! Arguments
@@ -267,6 +279,7 @@ contains
type(field_type), pointer :: height_wth => null()
type(field_type), pointer :: height_w3 => null()
type(quadrature_xyoz_type), pointer :: qr => null()
+ type(field_type) :: exner_surf, exner_surf_wt
! Get element orders
element_order_h = exner%get_element_order_h()
@@ -289,17 +302,30 @@ contains
! Initialise Exner pressure according to formulation
select case(method)
case(method_balanced)
- ! Initialise exner assuming hydrostatic balance
- call invoke( hydrostatic_exner_kernel_type( exner, theta, moist_dyn, &
- height_wth, height_w3, &
- chi, panel_id, gravity, &
- p_zero, rd, cp, planet_radius ) )
+ if (test == test_horizontal_mountain) then
+ ! Calculate surface pressure
+ call exner_surf%initialise( exner%get_function_space() )
+ call exner_surf_wt%initialise( theta%get_function_space() )
+ call invoke( initial_exner_sample_kernel_type( exner_surf, height_wth, &
+ chi, panel_id, time ) )
+ call map_physics_scalars(exner_surf_wt, exner_surf)
+ call invoke( hydrostatic_exner_squall_kernel_type( exner, theta, moist_dyn, &
+ exner_surf_wt, height_wth, height_w3, &
+ chi, panel_id, gravity, &
+ rd, cp, scaled_radius ) )
+ else
+ ! Initialise exner assuming hydrostatic balance
+ call invoke( hydrostatic_exner_kernel_type( exner, theta, moist_dyn, &
+ height_wth, height_w3, &
+ chi, panel_id, gravity, &
+ p_zero, rd, cp, scaled_radius ) )
+ end if
case(method_sampled)
! Initialise exner through sampling analytic solution
- call invoke (initial_exner_sample_kernel_type( exner, chi, panel_id, time ) )
+ call invoke (initial_exner_sample_kernel_type( exner, height_wth, chi, panel_id, time ) )
case(method_projected)
! Initialise exner by projecting analytic solution
- call invoke( set_exner_kernel_type ( exner, chi, panel_id, time, qr ) )
+ call invoke( set_exner_kernel_type ( exner, height_wth, chi, panel_id, time, qr ) )
case default
call log_event( "Gungho: Unrecognised method in initial_pressure", LOG_LEVEL_ERROR )
end select
@@ -308,6 +334,39 @@ contains
end subroutine init_exner_field
+ subroutine init_exner_squall(exner, exner_surf, theta, moist_dyn, time)
+
+ implicit none
+
+ ! Arguments
+ type(field_type), intent(inout) :: exner
+ type(field_type), intent(in) :: exner_surf
+ type(field_type), intent(in) :: theta
+ type(field_type), intent(in) :: moist_dyn(num_moist_factors)
+ real(kind=r_def), intent(in) :: time
+
+ ! Internal variables
+ integer(kind=i_def) :: mesh_id
+ type(field_type), pointer :: chi(:) => null()
+ type(field_type), pointer :: panel_id => null()
+ type(field_type), pointer :: height_wth => null()
+ type(field_type), pointer :: height_w3 => null()
+
+ ! Point to objects for this mesh
+ mesh_id = exner%get_mesh_id()
+ chi => get_coordinates(mesh_id)
+ panel_id => get_panel_id(mesh_id)
+ height_w3 => get_height_fe(W3, mesh_id)
+ height_wth => get_height_fe(Wtheta, mesh_id)
+
+ call invoke( hydrostatic_exner_squall_kernel_type( exner, theta, moist_dyn, &
+ exner_surf, height_wth, height_w3, &
+ chi, panel_id, gravity, &
+ rd, cp, scaled_radius ) )
+ nullify( chi, panel_id, height_wth, height_w3 )
+
+ end subroutine init_exner_squall
+
!> @details This routine initialises the Exner pressure field. It is callable
!> so that it can be used by other modules.
@@ -424,7 +483,9 @@ contains
call moist_dyn_factors_alg( moist_dyn, mr )
if ( (test /= test_grabowski_clark) .and. &
- (test /= test_bryan_fritsch) ) then
+ (test /= test_bryan_fritsch) .and. &
+ (test /= test_squall_line) .and. &
+ (test /= test_supercell) ) then
! Recompute hydrostatic balance following the addition of moisture
! TODO: This should not be done here -- will be fixed by #2877
! This should be done by separate routines from the driver level
@@ -439,7 +500,7 @@ contains
call invoke( hydrostatic_exner_kernel_type( exner, theta, moist_dyn, &
height_wth, height_w3, &
chi, panel_id, gravity, &
- p_zero, rd, cp, planet_radius ) )
+ p_zero, rd, cp, scaled_radius ) )
call init_rho_field( rho, theta, exner, moist_dyn, initial_time )
end if
diff --git a/science/gungho/source/algorithm/initialisation/init_saturated_profile_alg_mod.x90 b/science/gungho/source/algorithm/initialisation/init_saturated_profile_alg_mod.x90
index f52599b02..f49887fc5 100644
--- a/science/gungho/source/algorithm/initialisation/init_saturated_profile_alg_mod.x90
+++ b/science/gungho/source/algorithm/initialisation/init_saturated_profile_alg_mod.x90
@@ -17,7 +17,8 @@ use formulation_config_mod, only: moisture_formulation, &
moisture_formulation_dry, &
theta_moist_source
use function_space_mod, only: function_space_type
-use sci_geometric_constants_mod, only: get_coordinates, get_panel_id
+use sci_geometric_constants_mod, only: get_coordinates, get_panel_id, &
+ get_height_fv
use idealised_config_mod, only: test, test_bryan_fritsch
use init_gungho_prognostics_alg_mod, only: init_exner_field, &
init_rho_field, &
@@ -37,6 +38,7 @@ use planet_config_mod, only: rd, p_zero, kappa, cp, recip_epsilon
use theta_e_kernel_mod, only: theta_e_kernel_type
use transport_config_mod, only: theta_variable, &
theta_variable_moist
+use fs_continuity_mod, only: Wtheta
implicit none
@@ -90,6 +92,7 @@ subroutine init_saturated_profile_alg( theta, mr, exner, rho, moist_dyn )
type(field_type), allocatable :: latest_mr(:) ! latest guess for mr fields
type(field_type), pointer :: chi(:) => null()
+ type(field_type), pointer :: height_wth => null()
type(field_type), pointer :: panel_id => null()
type(function_space_type), pointer :: wt_fs => null()
@@ -123,6 +126,7 @@ subroutine init_saturated_profile_alg( theta, mr, exner, rho, moist_dyn )
chi => get_coordinates( mesh%get_id() )
panel_id => get_panel_id( mesh%get_id() )
wt_fs => theta%get_function_space()
+ height_wth => get_height_fv(Wtheta, mesh%get_id())
call specified_theta_e%initialise( vector_space = wt_fs )
call latest_theta_e%initialise( vector_space = wt_fs )
@@ -274,7 +278,7 @@ subroutine init_saturated_profile_alg( theta, mr, exner, rho, moist_dyn )
! Apply perturbation to theta_e
!----------------------------------------------------------------------------!
- call invoke( initial_theta_kernel_type( theta_pert, chi, panel_id ), &
+ call invoke( initial_theta_kernel_type( theta_pert, height_wth, chi, panel_id ), &
inc_X_times_Y( theta, theta_pert ) )
! Return to saturation
diff --git a/science/gungho/source/algorithm/initialisation/init_squall_line_profile_alg_mod.x90 b/science/gungho/source/algorithm/initialisation/init_squall_line_profile_alg_mod.x90
new file mode 100644
index 000000000..36c7d3b6f
--- /dev/null
+++ b/science/gungho/source/algorithm/initialisation/init_squall_line_profile_alg_mod.x90
@@ -0,0 +1,201 @@
+!-----------------------------------------------------------------------------
+! (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 Compute theta and moisture mixing ratios for the squall line.
+
+module init_squall_line_profile_alg_mod
+
+use constants_mod, only: r_def, i_def
+use sci_field_bundle_builtins_mod, only: set_bundle_scalar, clone_bundle
+use field_mod, only: field_type
+use formulation_config_mod, only: moisture_formulation, &
+ moisture_formulation_dry
+use function_space_mod, only: function_space_type
+use sci_geometric_constants_mod, only: get_coordinates, get_panel_id
+use idealised_config_mod, only: test, test_squall_line, test_supercell
+use init_gungho_prognostics_alg_mod, only: init_exner_squall, &
+ init_rho_field, &
+ init_mr_fields
+use initial_theta_kernel_mod, only: initial_theta_kernel_type
+use log_mod, only: LOG_LEVEL_ERROR, &
+ LOG_LEVEL_INFO, &
+ log_event, &
+ log_scratch_space
+use mesh_mod, only: mesh_type
+use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg
+use moist_dyn_mod, only: num_moist_factors
+use mr_indices_mod, only: nummr, imr_v, imr_cl
+use norm_alg_mod, only: rel_l2_error_alg
+use physics_mappings_alg_mod, only: map_physics_scalars
+use planet_config_mod, only: rd, p_zero, kappa, cp, recip_epsilon
+use theta_v_kernel_mod, only: theta_v_kernel_type
+use relative_humidity_kernel_mod, only: relative_humidity_kernel_type
+use init_unsaturated_profile_alg_mod, only: init_unsaturated_profile_alg
+use sci_geometric_constants_mod, only: get_coordinates, get_panel_id, get_height_fv
+use initial_theta_pert_kernel_mod, only: initial_theta_pert_kernel_type
+use initial_rel_humidity_kernel_mod, only: initial_rel_humidity_kernel_type
+use initial_theta_v_kernel_mod, only: initial_theta_v_kernel_type
+use fs_continuity_mod, only: Wtheta
+use sci_field_minmax_alg_mod, only: log_field_minmax
+use moist_dyn_mod, only: gas_law, total_mass
+
+implicit none
+
+private
+
+!-------------------------------------------------------------------------------
+! Contained functions/subroutines
+!-------------------------------------------------------------------------------
+public :: init_squall_line_profile_alg
+
+contains
+
+!> @brief Compute theta and moisture mixing ratios for a saturated atmosphere.
+!> @details The initial theta and moisture mixing ratio fields are computed for
+!> a saturated atmosphere, whose thermodynamics are specified in terms
+!> of a wet equivalent potential temperature, and a total moisture
+!> field. The atmosphere is saturated, so the water vapour field is set
+!> to the saturation mixing ratio, and the remainder of the moisture is
+!> cloud liquid.
+!>
+!> The algebraic relationship between these variables is complex. To
+!> compute our prognostic fields we use three nested loops, relying on
+!> Picard iteration to convergence to the solution. The Exner and rho
+!> fields are computed by taking hydrostatic balance.
+!>
+!> @param[in,out] theta potential temperature field
+!> @param[in,out] mr array of moisture mixing ratio fields
+!> @param[in,out] exner Exner pressure field
+!> @param[in,out] rho dry density field
+!> @param[in,out] moist_dyn array of fields with moist dynamics factors
+subroutine init_squall_line_profile_alg( theta, mr, exner, rho, moist_dyn )
+
+ implicit none
+
+ ! Arguments
+ type(field_type), intent(inout) :: theta
+ type(field_type), intent(inout) :: mr(nummr)
+ type(field_type), intent(inout) :: exner
+ type(field_type), intent(inout) :: rho
+ type(field_type), intent(inout) :: moist_dyn(num_moist_factors)
+
+ ! Internal variables
+ type(field_type) :: theta_v
+ type(field_type) :: theta_pert ! potential temperature perturbation
+ type(field_type) :: theta_eq ! equatorial theta field
+ type(field_type) :: theta_v_eq ! equatorial virtual theta field
+ type(field_type) :: exner_eq ! equatorial Exner field
+ type(field_type) :: exner_eq_wt ! Exner field at Wtheta points
+ type(field_type) :: exner_wt ! Exner field at Wtheta points
+ type(field_type) :: rho_eq ! equatorial dry density field
+ type(field_type) :: mr_eq(nummr) ! equatorial moisture mixing ratios
+
+ type(field_type), pointer :: chi(:) => null()
+ type(field_type), pointer :: panel_id => null()
+ type(field_type), pointer :: height_wth => null()
+ type(function_space_type), pointer :: wt_fs => null()
+ type(function_space_type), pointer :: w3_fs => null()
+
+ type(mesh_type), pointer :: mesh => null()
+
+ real(kind=r_def) :: initial_time
+ integer(kind=i_def) :: i
+
+ !----------------------------------------------------------------------------!
+ ! Initialise fields
+ !----------------------------------------------------------------------------!
+
+ initial_time = 0.0_r_def
+
+ mesh => theta%get_mesh()
+ chi => get_coordinates( mesh%get_id() )
+ panel_id => get_panel_id( mesh%get_id() )
+ height_wth => get_height_fv(Wtheta, mesh%get_id())
+ wt_fs => theta%get_function_space()
+ w3_fs => rho%get_function_space()
+
+ call exner_wt%initialise( vector_space = wt_fs )
+ call exner_eq_wt%initialise( vector_space = wt_fs )
+ call theta_pert%initialise( vector_space = wt_fs )
+ call theta_eq%initialise( vector_space = wt_fs )
+ call theta_v%initialise( vector_space = wt_fs )
+ call theta_v_eq%initialise( vector_space = wt_fs )
+ call exner_eq%initialise( vector_space = w3_fs )
+ call rho_eq%initialise( vector_space = w3_fs )
+ do i = 1, nummr
+ call mr_eq(i)%initialise( vector_space = wt_fs )
+ end do
+
+ !----------------------------------------------------------------------------!
+ ! FIRST: compute equatorial theta / m_v using unsaturated initialisation
+ !----------------------------------------------------------------------------!
+
+ ! Set initial moisture fields to be zero
+ call set_bundle_scalar(0.0_r_def, mr, nummr)
+
+ ! Set equatorial profiles
+ call invoke( setval_X(theta_eq, theta), &
+ setval_X(exner_eq, exner), &
+ setval_X(rho_eq, rho) )
+
+ ! Iterative procedure to get theta_v, mr_v and Exner at equator, given the
+ ! specified theta_d and relative humidity
+ call init_unsaturated_profile_alg( theta_eq, mr_eq, exner_eq, rho_eq, moist_dyn )
+
+ !----------------------------------------------------------------------------!
+ ! SECOND: set moisture mixing ratio everywhere to be equatorial value
+ !----------------------------------------------------------------------------!
+
+ call invoke( setval_X(mr(imr_v), mr_eq(imr_v)) )
+ call moist_dyn_factors_alg( moist_dyn, mr )
+
+ !----------------------------------------------------------------------------!
+ ! THIRD: Iterative procedure to get theta_v and Exner everywhere
+ !----------------------------------------------------------------------------!
+
+ call invoke( X_times_Y(theta_v_eq, theta_eq, moist_dyn(gas_law)), &
+ inc_X_divideby_Y(theta_v_eq, moist_dyn(total_mass)) )
+
+ call map_physics_scalars( exner_eq_wt, exner_eq )
+
+ call invoke( initial_theta_v_kernel_type(theta_v, exner_wt, &
+ theta_v_eq, exner_eq_wt, &
+ height_wth, chi, panel_id, test) )
+
+ !----------------------------------------------------------------------------!
+ ! FOURTH: Obtain initial theta from theta_v and moisture
+ !----------------------------------------------------------------------------!
+
+ call invoke( X_divideby_Y(theta, theta_v, moist_dyn(gas_law)), &
+ inc_X_times_Y(theta, moist_dyn(total_mass)) )
+
+ !----------------------------------------------------------------------------!
+ ! FIFTH: Determine Exner from hydrostatic balance
+ !----------------------------------------------------------------------------!
+
+ ! call init_exner_squall( exner, exner_surf, theta, moist_dyn, initial_time )
+
+ call map_physics_scalars( exner, exner_wt )
+
+ !----------------------------------------------------------------------------!
+ ! SIXTH: Apply perturbation to theta_v
+ !----------------------------------------------------------------------------!
+
+ if ( test == test_squall_line .or. test == test_supercell ) then
+ call theta_pert%initialise( vector_space = wt_fs )
+ call invoke( initial_theta_pert_kernel_type(theta_pert, chi, panel_id), &
+ inc_X_plus_Y(theta, theta_pert))
+ end if
+
+ !----------------------------------------------------------------------------!
+ ! SEVENTH: find the dry density, presumably from equation of state
+ !----------------------------------------------------------------------------!
+
+ call init_rho_field( rho, theta, exner, moist_dyn, initial_time )
+
+end subroutine init_squall_line_profile_alg
+
+end module init_squall_line_profile_alg_mod
diff --git a/science/gungho/source/algorithm/initialisation/init_unsaturated_profile_alg_mod.x90 b/science/gungho/source/algorithm/initialisation/init_unsaturated_profile_alg_mod.x90
index b23d043aa..00c2cb8f7 100644
--- a/science/gungho/source/algorithm/initialisation/init_unsaturated_profile_alg_mod.x90
+++ b/science/gungho/source/algorithm/initialisation/init_unsaturated_profile_alg_mod.x90
@@ -12,7 +12,8 @@ use constants_mod, only: r_def, i_def
use sci_field_bundle_builtins_mod, only: set_bundle_scalar, clone_bundle
use field_mod, only: field_type
use function_space_mod, only: function_space_type
-use idealised_config_mod, only: test, test_grabowski_clark
+use idealised_config_mod, only: test, test_grabowski_clark, &
+ test_squall_line, test_supercell
use init_gungho_prognostics_alg_mod, only: init_exner_field, &
init_mr_fields, &
init_rho_field
@@ -29,6 +30,11 @@ use planet_config_mod, only: rd, p_zero, kappa, &
recip_epsilon
use relative_humidity_kernel_mod, only: relative_humidity_kernel_type
+use sci_geometric_constants_mod, only: get_coordinates, get_panel_id
+use initial_rel_humidity_kernel_mod, only: initial_rel_humidity_kernel_type
+
+use sci_field_minmax_alg_mod, only: log_field_minmax
+
implicit none
private
@@ -50,7 +56,7 @@ contains
!> pressure is used to calculate the vapour saturation field. We use
!> an iterative procedure to arrive at appropriate values.
!>
-!> @param[in] theta potential temperature field
+!> @param[in,out] theta potential temperature field
!> @param[in,out] mr array of moisture mixing ratio fields
!> @param[in,out] exner Exner pressure field
!> @param[in,out] rho dry density field
@@ -73,6 +79,9 @@ subroutine init_unsaturated_profile_alg( theta, mr, exner, rho, moist_dyn )
type(field_type) :: latest_rel_hum ! latest relative humidity field
type(field_type) :: latest_mr_v ! latest water vapour mixing ratio field
+ type(field_type), pointer :: chi(:)
+ type(field_type), pointer :: panel_id
+
type(function_space_type), pointer :: wt_fs => null()
integer(kind=i_def) :: outer_count, inner_count
@@ -103,29 +112,35 @@ subroutine init_unsaturated_profile_alg( theta, mr, exner, rho, moist_dyn )
call latest_exner%initialise( vector_space = exner%get_function_space() )
call exner_at_wt%initialise( vector_space = wt_fs )
+ chi => get_coordinates(theta%get_mesh_id())
+ panel_id => get_panel_id(theta%get_mesh_id())
+
!----------------------------------------------------------------------------!
! Specify background relative humidity field
!----------------------------------------------------------------------------!
+ ! Set initial moisture fields to be zero
+ ! This will not be changed if moisture_formulation = 'dry'
+ call set_bundle_scalar(0.0_r_def, mr, nummr)
+
select case( test )
case( test_grabowski_clark )
- call invoke( setval_c(specified_rel_hum, 0.2_r_def ) )
+ call invoke( setval_c(specified_rel_hum, 0.2_r_def ), &
+ setval_c( mr(imr_v), 0.02_r_def) )
+
+ case ( test_squall_line, test_supercell )
+ call invoke( initial_rel_humidity_kernel_type(specified_rel_hum, chi, panel_id), &
+ setval_c( mr(imr_v), 1.e-5_r_def) )
- case default
- call log_event( "Gungho: Your idealised test is not compatible " // &
- "with the unsaturated profile routine", LOG_LEVEL_ERROR )
+ case default
+ call log_event( "Gungho: Your idealised test is not compatible " // &
+ "with the unsaturated profile routine", LOG_LEVEL_ERROR )
end select
!----------------------------------------------------------------------------!
! Do nested loops to compute background prognostic fields
!----------------------------------------------------------------------------!
- ! Set initial moisture fields to be zero
- ! This will not be changed if moisture_formulation = 'dry'
- call set_bundle_scalar(0.0_r_def, mr, nummr)
-
- ! First guess for mr_v is 0.02
- call invoke( setval_c( mr(imr_v), 0.02_r_def) )
! Give first guess for Exner
call moist_dyn_factors_alg( moist_dyn, mr )
@@ -196,14 +211,20 @@ subroutine init_unsaturated_profile_alg( theta, mr, exner, rho, moist_dyn )
call init_exner_field( exner, theta, moist_dyn, initial_time )
!----------------------------------------------------------------------------!
- ! Apply perturbation to water vapour
+ ! Apply perturbation to water vapour / theta
!----------------------------------------------------------------------------!
- ! Now apply perturbation to water vapour field
- call init_mr_fields( mr, theta, exner, rho, moist_dyn )
+ if ( test == test_grabowski_clark ) then
+ ! Now apply perturbation to water vapour field
+ call init_mr_fields( mr, theta, exner, rho, moist_dyn )
+
+ ! Adjust the dry density field so that the pressure is kept the same
+ call init_rho_field( rho, theta, exner, moist_dyn, initial_time )
+ end if
- ! Adjust the dry density field so that the pressure is kept the same
- call init_rho_field( rho, theta, exner, moist_dyn, initial_time )
+ call log_field_minmax(LOG_LEVEL_INFO, 'init_unsaturated_end: mr_v', mr(imr_v))
+ call log_field_minmax(LOG_LEVEL_INFO, 'init_unsaturated_end: exner', exner)
+ call log_field_minmax(LOG_LEVEL_INFO, 'init_unsaturated_end: rho', rho)
!----------------------------------------------------------------------------!
! Finalise internal fields and pointers
diff --git a/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90 b/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90
index 1a0eb304d..042f8892e 100644
--- a/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90
+++ b/science/gungho/source/algorithm/physics/external_forcing_alg_mod.X90
@@ -24,6 +24,7 @@ module external_forcing_alg_mod
wind_forcing_filtering, &
wind_forcing_diffusion, &
wind_forcing_profile, &
+ wind_forcing_rayleigh_friction, &
vertadvect_forcing, &
geostrophic_forcing, &
filtering_order, &
@@ -34,6 +35,7 @@ module external_forcing_alg_mod
LOG_LEVEL_INFO, LOG_LEVEL_ERROR
use lfric_xios_write_mod, only: write_field_generic
use field_parent_mod, only: write_interface
+ use sci_field_minmax_alg_mod, only: log_field_minmax
use held_suarez_fv_kernel_mod, only: held_suarez_fv_kernel_type
use held_suarez_fv_wind_kernel_mod, only: held_suarez_fv_wind_kernel_type
@@ -42,6 +44,7 @@ module external_forcing_alg_mod
use shallow_hot_jupiter_kernel_mod, only: shallow_hot_jupiter_kernel_type
use deep_hot_jupiter_kernel_mod, only: deep_hot_jupiter_kernel_type
use momentum_smagorinsky_kernel_mod, only: momentum_smagorinsky_kernel_type
+ use rayleigh_friction_kernel_mod, only: rayleigh_friction_kernel_type
#ifdef UM_PHYSICS
use pc2_force_response_alg_mod, only: pc2_force_response_alg
#endif
@@ -61,10 +64,11 @@ module external_forcing_alg_mod
use mr_indices_mod, only: nummr
use io_config_mod, only: subroutine_timers, write_diag, &
use_xios_io
- use planet_config_mod, only: kappa
+ use planet_config_mod, only: kappa, p_zero
use section_choice_config_mod, only: cloud, cloud_um
use cloud_config_mod, only: scheme, scheme_pc2
-
+ use idealised_config_mod, only: test, test_horizontal_mountain, &
+ test_breaking_gw
use reference_element_mod, only: T
use sci_set_any_dof_kernel_mod, only: set_any_dof_kernel_type
@@ -137,6 +141,7 @@ contains
! Temporary field to unpack from field collection
type( field_type ), pointer :: u_physics => null()
+ type( field_type ), pointer :: u_ref => null()
type( field_type ), pointer :: w2_rmultiplicity => null()
type( field_type ), pointer :: dA => null()
type( field_type ), pointer :: chi(:) => null()
@@ -150,6 +155,7 @@ contains
integer( kind=i_def ) :: order, nlayers
real( kind=r_def ) :: diffusion_coefficient_scaled, dx_sum, dx_size
+ real( kind=r_def ) :: p_cutoff, tau
type( mesh_type ), pointer :: mesh => null()
procedure(write_interface), pointer :: write_behaviour => null()
@@ -351,6 +357,36 @@ contains
call wind_forcing_profile_alg( du_forcing, exner, time_now )
+ case ( wind_forcing_rayleigh_friction )
+
+ w2_rmultiplicity => get_rmultiplicity_fv( W2, mesh%get_id() ) ! 1/multiplicity of w2
+ call derived_fields%get_field('u_ref', u_ref)
+ call du_forcing%initialise( vector_space = u%get_function_space() )
+
+ if (test == test_horizontal_mountain) then
+ ! For the horizontal mountain test we use a different p_zero and p_cutoff
+ tau = 0.1_r_def*24.0_r_def*60.0_r_def*60.0_r_def
+ p_cutoff = 30500.0_r_def
+ else if (test == test_breaking_gw) then
+ ! For the horizontal mountain test we use a different p_zero and p_cutoff
+ tau = 5.0_r_def*24.0_r_def*60.0_r_def*60.0_r_def
+ p_cutoff = 16500.0_r_def
+ else
+ call log_event('Rayleigh friction not implemented for this test', &
+ LOG_LEVEL_ERROR)
+ end if
+
+ call invoke( setval_c(du_forcing, 0.0_r_def), &
+ rayleigh_friction_kernel_type(du_forcing, u, u_ref, &
+ w2_rmultiplicity, &
+ exner, exner_in_wth, &
+ kappa, p_zero, p_cutoff, &
+ tau, dt))
+
+ call log_field_minmax(LOG_LEVEL_INFO, 'Rayleigh friction u', u)
+ call log_field_minmax(LOG_LEVEL_INFO, 'Rayleigh friction u_ref', u_ref)
+ call log_field_minmax(LOG_LEVEL_INFO, 'Rayleigh friction du', du_forcing)
+
case default
call log_event( 'slow_physics: Incorrect wind_forcing option', &
diff --git a/science/gungho/source/algorithm/physics/fast_physics_alg_mod.X90 b/science/gungho/source/algorithm/physics/fast_physics_alg_mod.X90
index 34bc852f4..5c426f2d7 100644
--- a/science/gungho/source/algorithm/physics/fast_physics_alg_mod.X90
+++ b/science/gungho/source/algorithm/physics/fast_physics_alg_mod.X90
@@ -11,6 +11,7 @@ module fast_physics_alg_mod
surface, surface_jules, &
cloud, cloud_um, &
cloud_evap_condense, &
+ cloud_kessler, &
stochastic_physics, &
stochastic_physics_um
@@ -27,6 +28,8 @@ module fast_physics_alg_mod
use io_config_mod, only: subroutine_timers, &
write_conservation_diag
use timer_mod, only: timer
+ use extrusion_mod, only : TWOD
+ use fs_continuity_mod, only : W3, Wtheta
use physics_config_mod, only: blayer_placement, blayer_placement_fast, &
convection_placement, &
@@ -45,6 +48,8 @@ module fast_physics_alg_mod
use convection_config_mod,only: cv_scheme, cv_scheme_lambert_lewis, &
cv_scheme_gregory_rowntree, &
cv_scheme_comorph
+ use sci_geometric_constants_mod, only: get_height_fv
+ use physics_mappings_alg_mod, only: map_physics_scalars
use cloud_config_mod, only: scheme, scheme_pc2, cloud_call_b4_conv
use stochastic_physics_config_mod, &
only: blpert_type, blpert_type_off
@@ -67,6 +72,11 @@ module fast_physics_alg_mod
only : outer_iterations
use evap_condense_kernel_mod, &
only: evap_condense_kernel_type
+ use kessler_kernel_mod, only: kessler_kernel_type
+ use mesh_collection_mod, only: mesh_collection
+ use function_space_collection_mod, only: function_space_collection
+ use function_space_mod, only: function_space_type
+ use mesh_mod, only: mesh_type
#ifdef UM_PHYSICS
use cld_alg_mod, only: cld_alg
! UM BL scheme
@@ -162,7 +172,7 @@ contains
! ...for boundary layer
type( field_type ) :: dtheta_bl, du_bl
! ...for evaporation condensation
- type( field_type ) :: dtheta_cld
+ type( field_type ) :: dtheta_cld, rho_in_wth, rain_2d
type( field_type ) :: dmr_cld(nummr)
#ifdef UM_PHYSICS
@@ -191,6 +201,9 @@ contains
type( field_type), pointer :: theta_star => null()
type( field_type), pointer :: u_star => null()
type( field_type), pointer :: exner_in_wth => null()
+ type( field_type), pointer :: height_in_wth
+ type( function_space_type ), pointer :: w3_2d
+ type(mesh_type), pointer :: twod_mesh
! Flags to indicate whether scheme has run and increments need adding
logical(kind=l_def) :: boundary_layer_done, stph_done, blpert_done, &
@@ -359,6 +372,23 @@ contains
exner_in_wth, &
Rd, Rv, cp, cpv, cl, p_zero) )
evap_condense_done = .true.
+ else if (evap_condense_placement == evap_condense_placement_fast .and. &
+ cloud == cloud_kessler) then
+ call derived_fields%get_field('exner_in_wth', exner_in_wth)
+ call clone_bundle(mr, dmr_cld, nummr)
+ twod_mesh => mesh_collection%get_mesh(theta%get_mesh(), TWOD)
+ w3_2d => function_space_collection%get_fs( twod_mesh, 0, 0, W3 )
+ call rain_2d%initialise(w3_2d)
+ call dtheta%copy_field_properties(dtheta_cld)
+ call dtheta%copy_field_properties(rho_in_wth)
+ call map_physics_scalars(rho_in_wth, rho)
+ height_in_wth => get_height_fv(Wtheta, theta%get_mesh_id())
+ call invoke( kessler_kernel_type(dtheta_cld, theta, &
+ dmr_cld, &
+ mr, &
+ exner_in_wth, rho_in_wth, &
+ height_in_wth, rain_2d, dt) )
+ evap_condense_done = .true.
end if
!=====================================================================
diff --git a/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90 b/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90
index f2bbe547f..5ddaac13c 100644
--- a/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90
+++ b/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90
@@ -10,6 +10,7 @@ module slow_physics_alg_mod
only: heat_capacity_h2o_vapour, &
heat_capacity_h2o, &
gas_constant_h2o
+ use sci_geometric_constants_mod, only: get_height_fv
use energy_correction_config_mod, &
only: encorr_usage, &
encorr_usage_apply
@@ -65,6 +66,7 @@ module slow_physics_alg_mod
surface, surface_jules, &
convection, convection_um, &
cloud, cloud_um, cloud_evap_condense, &
+ cloud_kessler, &
external_forcing
! Derived Types
@@ -94,6 +96,7 @@ module slow_physics_alg_mod
use timer_mod, only: timer
use external_forcing_alg_mod, only: external_forcing_alg
use evap_condense_kernel_mod, only: evap_condense_kernel_type
+ use kessler_kernel_mod, only: kessler_kernel_type
use physics_mappings_alg_mod, only: map_physics_scalars
use planet_config_mod, only: p_zero, kappa, cp, Rd
! Physics schemes called
@@ -324,6 +327,7 @@ contains
type( field_type ), pointer :: exner_in_wth => null()
type( field_type ), pointer :: rho_in_wth => null()
type( field_type ), pointer :: exner_wth_n => null()
+ type( field_type ), pointer :: height_in_wth => null()
! Fields required by PC2
type( field_type ), pointer :: liquid_fraction => null()
type( field_type ), pointer :: frozen_fraction => null()
@@ -331,7 +335,8 @@ contains
type( field_type ), pointer :: cf_liq_n
type( field_type ), pointer :: cf_fro_n
type( field_type ), pointer :: cf_bulk_n
-
+ type( field_type ) :: rain_2d
+ type(function_space_type), pointer :: w3_2d
integer( kind=i_def ) :: i_mr
! Flags to indicate whether scheme has run and increments need adding
@@ -1003,6 +1008,19 @@ contains
exner_in_wth, &
Rd, Rv, cp, cpv, cl, p_zero) )
evap_condense_done = .true.
+ else if (evap_condense_placement == evap_condense_placement_slow .and. &
+ cloud == cloud_kessler) then
+
+ call clone_bundle(mr, dmr_cld, nummr)
+ w3_2d => function_space_collection%get_fs( twod_mesh, 0, 0, W3 )
+ call rain_2d%initialise(w3_2d)
+ height_in_wth => get_height_fv(Wtheta, mesh%get_id())
+ call invoke( kessler_kernel_type(dtheta_cld, theta, &
+ dmr_cld, &
+ mr, &
+ exner_in_wth, rho_in_wth, &
+ height_in_wth, rain_2d, dt) )
+ evap_condense_done = .true.
end if
!=====================================================================
diff --git a/science/gungho/source/algorithm/runtime_constants/dycore_constants_mod.x90 b/science/gungho/source/algorithm/runtime_constants/dycore_constants_mod.x90
index 9421aba5b..713ce0d40 100644
--- a/science/gungho/source/algorithm/runtime_constants/dycore_constants_mod.x90
+++ b/science/gungho/source/algorithm/runtime_constants/dycore_constants_mod.x90
@@ -90,6 +90,9 @@ contains
use base_mesh_config_mod, only: f_lat
use compute_coriolis_matrix_kernel_mod, &
only: compute_coriolis_matrix_kernel_type
+ use compute_trad_coriolis_matrix_kernel_mod, &
+ only: compute_trad_coriolis_matrix_kernel_type
+ use formulation_config_mod, only: traditional_coriolis
use planet_config_mod, only: scaled_omega
implicit none
@@ -123,10 +126,17 @@ contains
call coriolis_inventory%add_operator(coriolis, w2_fs, w2_fs, mesh)
- call invoke( compute_coriolis_matrix_kernel_type(coriolis, &
- chi, panel_id, &
- scaled_omega, &
- f_lat, qr) )
+ if (traditional_coriolis) then
+ call invoke( compute_trad_coriolis_matrix_kernel_type(coriolis, &
+ chi, panel_id, &
+ scaled_omega, &
+ f_lat, qr) )
+ else
+ call invoke( compute_coriolis_matrix_kernel_type(coriolis, &
+ chi, panel_id, &
+ scaled_omega, &
+ f_lat, qr) )
+ end if
if ( subroutine_timers ) call timer('runtime_constants.dycore')
end if
@@ -143,7 +153,10 @@ contains
use base_mesh_config_mod, only: f_lat
use compute_vert_coriolis_matrix_kernel_mod, &
only: compute_vert_coriolis_matrix_kernel_type
+ use formulation_config_mod, only: traditional_coriolis
use planet_config_mod, only: scaled_omega
+ use sci_operator_algebra_kernel_mod, &
+ only: operator_setval_c_kernel_type
implicit none
@@ -178,10 +191,14 @@ contains
call vert_coriolis_inventory%add_operator(vert_coriolis, wt_fs, w2_fs, mesh)
- call invoke( compute_vert_coriolis_matrix_kernel_type(vert_coriolis, &
- chi, panel_id, &
- scaled_omega, &
- f_lat, qr) )
+ if (traditional_coriolis) then
+ call invoke( operator_setval_c_kernel_type(vert_coriolis, 0.0_r_def) )
+ else
+ call invoke( compute_vert_coriolis_matrix_kernel_type(vert_coriolis, &
+ chi, panel_id, &
+ scaled_omega, &
+ f_lat, qr) )
+ end if
if ( subroutine_timers ) call timer('runtime_constants.dycore')
end if
diff --git a/science/gungho/source/configuration/orography_control_mod.F90 b/science/gungho/source/configuration/orography_control_mod.F90
index 4c19841e2..2ba65ac4c 100644
--- a/science/gungho/source/configuration/orography_control_mod.F90
+++ b/science/gungho/source/configuration/orography_control_mod.F90
@@ -15,17 +15,21 @@
module orography_control_mod
use constants_mod, only : i_def, str_short
- use base_mesh_config_mod, only : geometry, &
- geometry_spherical, &
+ use base_mesh_config_mod, only : geometry, &
+ geometry_spherical, &
geometry_planar
- use orography_config_mod, only : profile, &
- profile_schar, &
- profile_agnesi, &
- profile_bell, &
- profile_dcmip200
- use log_mod, only : log_event, &
- log_scratch_space, &
- LOG_LEVEL_INFO
+ use orography_config_mod, only : profile, &
+ profile_schar, &
+ profile_agnesi, &
+ profile_bell, &
+ profile_dcmip200, &
+ profile_dcmip_breaking_gw, &
+ profile_dcmip_gap, &
+ profile_dcmip_vortex
+ use log_mod, only : log_event, &
+ log_scratch_space, &
+ LOG_LEVEL_INFO, &
+ LOG_LEVEL_ERROR
use analytic_orography_mod, only : orography_profile
implicit none
@@ -86,6 +90,30 @@ subroutine set_orography_option()
! coordinates and initialise the corresponding type
call set_orography_dcmip200_spherical()
end if
+ ! DCMIP 2025 orography
+ case( profile_dcmip_breaking_gw )
+ call set_orography_dcmip_breaking_gw()
+
+ ! DCMIP 2025 orography
+ case( profile_dcmip_gap )
+ if ( geometry == geometry_spherical ) then
+ ! Read parameters for DCMIP 2025 mountain in spherical
+ ! coordinates and initialise the corresponding type
+ call set_orography_dcmip_gap()
+ else
+ call log_event('DCMIP 2025 orography only appropriate for sphere', LOG_LEVEL_ERROR)
+ end if
+
+ ! DCMIP 2025 orography
+ case( profile_dcmip_vortex )
+ if ( geometry == geometry_spherical ) then
+ ! Read parameters for DCMIP 2025 mountain in spherical
+ ! coordinates and initialise the corresponding type
+ call set_orography_dcmip_vortex()
+ else
+ call log_event('DCMIP 2025 orography only appropriate for sphere', LOG_LEVEL_ERROR)
+ end if
+
! Bell-shaped orography
case( profile_bell )
if ( geometry == geometry_planar ) then
@@ -263,6 +291,50 @@ subroutine set_orography_dcmip200_spherical()
return
end subroutine set_orography_dcmip200_spherical
+ subroutine set_orography_dcmip_breaking_gw()
+
+ use dcmip_2025_orography_mod, only : dcmip_breaking_gw_type
+
+ implicit none
+
+ allocate( orography_profile, source=dcmip_breaking_gw_type() )
+
+ write(log_scratch_space,'(A,A)') "set_orography_dcmip_breaking_gw: "// &
+ "Set analytic orography type to dcmip breaking GW mountain."
+ call log_event(log_scratch_space, LOG_LEVEL_INFO)
+
+ return
+ end subroutine set_orography_dcmip_breaking_gw
+
+ subroutine set_orography_dcmip_gap()
+
+ use dcmip_2025_orography_mod, only : dcmip_gap_type
+
+ implicit none
+
+ allocate( orography_profile, source = dcmip_gap_type() )
+
+ write(log_scratch_space,'(A,A)') "set_orography_dcmip_gap: "// &
+ "Set analytic orography type to dcmip gap mountain."
+ call log_event(log_scratch_space, LOG_LEVEL_INFO)
+
+ return
+ end subroutine set_orography_dcmip_gap
+
+ subroutine set_orography_dcmip_vortex()
+
+ use dcmip_2025_orography_mod, only : dcmip_vortex_type
+
+ implicit none
+
+ allocate( orography_profile, source=dcmip_vortex_type() )
+
+ write(log_scratch_space,'(A,A)') "set_orography_dcmip_vortex: "// &
+ "Set analytic orography type to dcmip vortex mountain."
+ call log_event(log_scratch_space, LOG_LEVEL_INFO)
+
+ return
+ end subroutine set_orography_dcmip_vortex
!=============================================================================
!> @brief Initialises analytic orography type for bell-shaped mountain in
diff --git a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90 b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90
index f1787a8bf..7457230d0 100755
--- a/science/gungho/source/diagnostics/diagnostics_calc_mod.F90
+++ b/science/gungho/source/diagnostics/diagnostics_calc_mod.F90
@@ -16,7 +16,8 @@ module diagnostics_calc_mod
use diagnostic_alg_mod, only: divergence_diagnostic_alg, &
hydbal_diagnostic_alg, &
vorticity_diagnostic_alg, &
- potential_vorticity_diagnostic_alg
+ potential_vorticity_diagnostic_alg, &
+ dbz_diagnostic_alg
use io_config_mod, only: use_xios_io, &
nodal_output_on_w3
use files_config_mod, only: diag_stem_name
@@ -43,13 +44,16 @@ module diagnostics_calc_mod
use mesh_mod, only: mesh_type
use sci_geometric_constants_mod, only: get_coordinates, &
get_panel_id
+ use mr_indices_mod, only: nummr
implicit none
private
public :: write_divergence_diagnostic, &
write_pv_diagnostic, &
write_hydbal_diagnostic, &
- write_vorticity_diagnostic
+ write_vorticity_diagnostic, &
+ write_dbz_diagnostic, &
+ write_exner_wth_diagnostic
contains
@@ -101,6 +105,68 @@ subroutine write_divergence_diagnostic(u_field, clock, mesh)
end subroutine write_divergence_diagnostic
+subroutine write_dbz_diagnostic(exner_in_wth, theta, rho, mr, clock, mesh)
+ implicit none
+
+ type(field_type), intent(in) :: exner_in_wth
+ type(field_type), intent(in) :: theta
+ type(field_type), intent(in) :: rho
+ type(field_type), intent(in) :: mr(nummr)
+ class(model_clock_type), intent(in) :: clock
+ type(mesh_type), intent(in), pointer :: mesh
+
+ type(field_type) :: dbz_field
+
+ procedure(write_interface), pointer :: tmp_write_ptr
+
+ if (diagnostic_to_be_sampled('dbz')) then
+
+ ! Create the divergence diagnostic
+ call dbz_diagnostic_alg( dbz_field, exner_in_Wth, theta, rho, mr, mesh )
+
+ if (use_xios_io) then
+ !If using XIOS, we need to set a field I/O method appropriately
+ tmp_write_ptr => write_field_generic
+ call dbz_field%set_write_behaviour(tmp_write_ptr)
+ end if
+
+ call write_scalar_diagnostic( 'dbz', dbz_field, &
+ clock, mesh, .false. )
+
+ nullify(tmp_write_ptr)
+
+ end if
+
+end subroutine write_dbz_diagnostic
+
+
+subroutine write_exner_wth_diagnostic(exner_wth, clock, mesh)
+ implicit none
+
+ type(field_type), intent(inout) :: exner_wth
+ class(model_clock_type), intent(in) :: clock
+ type(mesh_type), intent(in), pointer :: mesh
+
+ procedure(write_interface), pointer :: tmp_write_ptr
+
+ if (diagnostic_to_be_sampled('exner_in_wth') .or. &
+ diagnostic_to_be_sampled('init_exner_in_wth') ) then
+
+ if (use_xios_io) then
+ !If using XIOS, we need to set a field I/O method appropriately
+ tmp_write_ptr => write_field_generic
+ call exner_wth%set_write_behaviour(tmp_write_ptr)
+ end if
+
+ call write_scalar_diagnostic( 'exner_in_wth', exner_wth, &
+ clock, mesh, .false. )
+
+ nullify(tmp_write_ptr)
+
+ end if
+
+end subroutine write_exner_wth_diagnostic
+
!-------------------------------------------------------------------------------
!> @brief Handles hydrostatic balance diagnostic processing
!!
diff --git a/science/gungho/source/driver/create_physics_prognostics_mod.F90 b/science/gungho/source/driver/create_physics_prognostics_mod.F90
index 649d501bd..beccb2c31 100644
--- a/science/gungho/source/driver/create_physics_prognostics_mod.F90
+++ b/science/gungho/source/driver/create_physics_prognostics_mod.F90
@@ -62,10 +62,13 @@ module create_physics_prognostics_mod
orographic_drag_um, &
convection, convection_um, &
stochastic_physics, &
- stochastic_physics_um
+ stochastic_physics_um, &
+ external_forcing
use cloud_config_mod, only : scheme, &
scheme_pc2
use convection_config_mod, only : cv_scheme, cv_scheme_comorph
+ use external_forcing_config_mod, only : wind_forcing, &
+ wind_forcing_rayleigh_friction
use microphysics_config_mod, only : microphysics_casim
use multires_coupling_config_mod, only : coarse_rad_aerosol
@@ -212,6 +215,10 @@ subroutine process_physics_prognostics(processor)
call processor%apply(make_spec('u_star', main%derived, W2))
call processor%apply(make_spec('wetrho_in_w2', main%derived, W2))
+ if ( external_forcing .and. wind_forcing == wind_forcing_rayleigh_friction ) then
+ call processor%apply(make_spec('u_ref', main%derived, W2))
+ end if
+
! W2H fields
call processor%apply(make_spec('u_in_w2h', main%derived, W2H))
call processor%apply(make_spec('v_in_w2h', main%derived, W2H))
diff --git a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90 b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90
index 12cc745a4..6f6d01499 100644
--- a/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90
+++ b/science/gungho/source/driver/gungho_diagnostics_driver_mod.F90
@@ -22,7 +22,9 @@ module gungho_diagnostics_driver_mod
use diagnostics_calc_mod, only : write_divergence_diagnostic, &
write_hydbal_diagnostic, &
write_vorticity_diagnostic, &
- write_pv_diagnostic
+ write_pv_diagnostic, &
+ write_dbz_diagnostic, &
+ write_exner_wth_diagnostic
use initialise_diagnostics_mod, only : diagnostic_to_be_sampled
use field_array_mod, only : field_array_type
use field_collection_iterator_mod, &
@@ -35,7 +37,7 @@ module gungho_diagnostics_driver_mod
use formulation_config_mod, only : use_physics, &
moisture_formulation, &
moisture_formulation_dry
- use fs_continuity_mod, only : W3, Wtheta
+ use fs_continuity_mod, only : W3, Wtheta, W2H
use integer_field_mod, only : integer_field_type
use initialization_config_mod, only : ls_option, &
ls_option_analytic, &
@@ -52,6 +54,7 @@ module gungho_diagnostics_driver_mod
use timer_mod, only : timer
use transport_config_mod, only : transport_ageofair
use driver_modeldb_mod, only : modeldb_type
+ use physics_mappings_alg_mod, only : map_physics_scalars
#ifdef UM_PHYSICS
use pres_lev_diags_alg_mod, only : pres_lev_diags_alg
@@ -103,6 +106,7 @@ subroutine gungho_diagnostics_driver( modeldb, &
type( field_type), pointer :: panel_id => null()
type( field_type), pointer :: height_w3 => null()
type( field_type), pointer :: height_wth => null()
+ type( field_type), pointer :: height_w2h => null()
type( field_type), pointer :: lbc_u => null()
type( field_type), pointer :: lbc_theta => null()
type( field_type), pointer :: lbc_rho => null()
@@ -115,6 +119,7 @@ subroutine gungho_diagnostics_driver( modeldb, &
type( field_type), pointer :: ageofair => null()
type( field_type), pointer :: exner_in_wth => null()
type( field_type), pointer :: dA => null()
+ type( field_type ) :: exner_in_wth_init
type(field_array_type), pointer :: mr_array => null()
type(field_array_type), pointer :: moist_dyn_array => null()
@@ -169,10 +174,12 @@ subroutine gungho_diagnostics_driver( modeldb, &
! Get the finite element height
height_w3 => get_height_fe(W3, mesh%get_id())
height_wth => get_height_fe(Wtheta, mesh%get_id())
+ height_w2h => get_height_fe(W2H, mesh%get_id())
else
! Get the finite volume height
height_w3 => get_height_fv(W3, mesh%get_id())
height_wth => get_height_fv(Wtheta, mesh%get_id())
+ height_w2h => get_height_fv(W2H, mesh%get_id())
end if
! Scalar fields
@@ -186,6 +193,8 @@ subroutine gungho_diagnostics_driver( modeldb, &
modeldb%clock, mesh, nodal_output_on_w3)
call write_scalar_diagnostic('height_wth', height_wth, &
modeldb%clock, mesh, nodal_output_on_w3)
+ call write_scalar_diagnostic('height_w2h', height_w2h, &
+ modeldb%clock, mesh, nodal_output_on_w3)
if (transport_ageofair) then
call con_tracer_last_outer%get_field('ageofair',ageofair)
@@ -332,6 +341,23 @@ subroutine gungho_diagnostics_driver( modeldb, &
! Don't output for the tangent linear model
call write_divergence_diagnostic( u, modeldb%clock, mesh )
call write_hydbal_diagnostic( theta, moist_dyn, exner, mesh )
+
+ if ( .not. modeldb%clock%is_initialisation() ) then
+ call derived_fields%get_field('exner_in_wth', exner_in_wth)
+ call write_exner_wth_diagnostic( exner_in_wth, modeldb%clock, mesh )
+ else if (use_xios_io .and. modeldb%clock%is_initialisation() &
+ .and. diagnostic_to_be_sampled("init_exner_in_wth") ) then
+ call exner_in_wth_init%initialise(theta%get_function_space())
+ call map_physics_scalars(exner_in_wth_init, exner)
+ tmp_write_ptr => write_field_generic
+ call exner_in_wth_init%set_write_behaviour(tmp_write_ptr)
+ call exner_in_wth_init%write_field("init_exner_in_wth")
+ end if
+
+ if ( moisture_formulation /= moisture_formulation_dry .and. .not. modeldb%clock%is_initialisation() ) then
+ call derived_fields%get_field('exner_in_wth', exner_in_wth)
+ call write_dbz_diagnostic( exner_in_wth, theta, rho, mr, modeldb%clock, mesh )
+ end if
end if
if ( subroutine_timers ) call timer('gungho_diagnostics_driver')
diff --git a/science/gungho/source/driver/gungho_init_fields_mod.X90 b/science/gungho/source/driver/gungho_init_fields_mod.X90
index db9e0d3db..6c59f7c1a 100644
--- a/science/gungho/source/driver/gungho_init_fields_mod.X90
+++ b/science/gungho/source/driver/gungho_init_fields_mod.X90
@@ -10,6 +10,7 @@
!>
module gungho_init_fields_mod
+ use copy_field_alg_mod, only : copy_field
use mr_indices_mod, only : nummr, mr_names
use moist_dyn_mod, only : num_moist_factors, moist_dyn_names
use field_array_mod, only : field_array_type
@@ -88,6 +89,9 @@ module gungho_init_fields_mod
use derived_config_mod, only : l_esm_couple
use ageofair_alg_mod, only : init_ageofair
use transport_config_mod, only : transport_ageofair
+ use section_choice_config_mod, only : external_forcing
+ use external_forcing_config_mod, only : wind_forcing, &
+ wind_forcing_rayleigh_friction
#ifdef COUPLED
use coupler_mod, only : zero_coupling_send_fields
#endif
@@ -553,6 +557,7 @@ subroutine create_model_data( modeldb, &
type( field_type ), pointer :: exner => null()
type( field_type ), pointer :: accumulated_fluxes => null()
type( field_type ), pointer :: ageofair => null()
+ type( field_type ), pointer :: u_ref => null()
type( io_value_type ), pointer :: temp_corr_io_value
type( io_value_type ), pointer :: random_seed_io_value
@@ -815,6 +820,11 @@ subroutine create_model_data( modeldb, &
call prognostic_fields%get_field('exner',exner)
call prognostic_fields%get_field('rho',rho)
+ if (external_forcing .and. wind_forcing == wind_forcing_rayleigh_friction) then
+ call derived_fields%get_field('u_ref', u_ref)
+ call copy_field(u, u_ref)
+ end if
+
if ( encorr_usage /= encorr_usage_none ) then
if ( geometry /= geometry_spherical ) then
write(log_scratch_space, '(a)') &
diff --git a/science/gungho/source/driver/gungho_init_prognostics_driver_mod.f90 b/science/gungho/source/driver/gungho_init_prognostics_driver_mod.f90
index c70dfee16..2a6d4b0b9 100644
--- a/science/gungho/source/driver/gungho_init_prognostics_driver_mod.f90
+++ b/science/gungho/source/driver/gungho_init_prognostics_driver_mod.f90
@@ -18,11 +18,14 @@ module gungho_init_prognostics_driver_mod
init_mr_fields
use init_saturated_profile_alg_mod, only: init_saturated_profile_alg
use init_unsaturated_profile_alg_mod, only: init_unsaturated_profile_alg
+ use init_squall_line_profile_alg_mod, only: init_squall_line_profile_alg
use init_thermo_profile_alg_mod, only: init_thermo_profile_alg
use idealised_config_mod, only: test, &
test_specified_profiles, &
test_bryan_fritsch, &
test_grabowski_clark, &
+ test_squall_line, &
+ test_supercell, &
perturb_init, &
perturb_magnitude
use random_perturb_field_alg_mod, only: random_perturb_field
@@ -93,7 +96,9 @@ subroutine init_gungho_prognostics(prognostic_fields, mr, moist_dyn)
if ( test == test_grabowski_clark ) then
call init_unsaturated_profile_alg( theta, mr, exner, rho, moist_dyn )
- else if (test /= test_bryan_fritsch .and. test /= test_specified_profiles) then
+ else if ( test == test_squall_line .or. test == test_supercell ) then
+ call init_squall_line_profile_alg( theta, mr, exner, rho, moist_dyn )
+ else if ( test /= test_bryan_fritsch .and. test /= test_specified_profiles ) then
call init_mr_fields( mr, theta, exner, rho, moist_dyn )
end if
diff --git a/science/gungho/source/initialisation/analytic_density_profiles_mod.F90 b/science/gungho/source/initialisation/analytic_density_profiles_mod.F90
index 1d4c7c7f5..bc25b8e9f 100644
--- a/science/gungho/source/initialisation/analytic_density_profiles_mod.F90
+++ b/science/gungho/source/initialisation/analytic_density_profiles_mod.F90
@@ -31,6 +31,7 @@ module analytic_density_profiles_mod
test_solid_body_rotation, &
test_solid_body_rotation_alt, &
test_deep_baroclinic_wave, &
+ test_breaking_gw, &
test_isentropic, &
test_isot_atm, &
test_isot_cold_atm, &
@@ -53,6 +54,7 @@ module analytic_density_profiles_mod
use analytic_temperature_profiles_mod, &
only : analytic_temperature
use deep_baroclinic_wave_mod, only : deep_baroclinic_wave
+use breaking_gravity_wave_mod, only : breaking_gravity_wave
use extrusion_config_mod, only : domain_height
implicit none
@@ -243,7 +245,7 @@ function analytic_density(chi, choice, time) result(density)
case( test_solid_body_rotation, &
test_solid_body_rotation_alt )
t0 = 280.0_r_def
- temperature = analytic_temperature(chi, choice)
+ temperature = analytic_temperature(chi, choice, 0.0_r_def)
pressure = t0 / temperature
density = p_zero * pressure**( (1.0_r_def - kappa )/ kappa ) &
/ (Rd * temperature)
@@ -251,6 +253,11 @@ function analytic_density(chi, choice, time) result(density)
call deep_baroclinic_wave(long, lat, radius-scaled_radius, &
pressure, temperature, density, &
u, v, w, mr_v)
+ case( test_breaking_gw )
+ call breaking_gravity_wave(long, lat, radius-scaled_radius, &
+ pressure, temperature, density, &
+ u, v, w, mr_v)
+
case( test_cosine_bubble )
l1 = sqrt( ((chi(1) - x1)/r1)**2 + ((chi(3) - y1)/r2)**2 )
if ( l1 < 1.0_r_def ) then
diff --git a/science/gungho/source/initialisation/analytic_moisture_profiles_mod.F90 b/science/gungho/source/initialisation/analytic_moisture_profiles_mod.F90
index b3329459b..884818832 100644
--- a/science/gungho/source/initialisation/analytic_moisture_profiles_mod.F90
+++ b/science/gungho/source/initialisation/analytic_moisture_profiles_mod.F90
@@ -15,7 +15,9 @@ module analytic_moisture_profiles_mod
LOG_LEVEL_ERROR
use idealised_config_mod, only : test_grabowski_clark, &
test_deep_baroclinic_wave, &
- test_isot_dry_atm
+ test_isot_dry_atm, &
+ test_squall_line, &
+ test_supercell
use physics_common_mod, only : qsaturation
use planet_config_mod, only : recip_epsilon, scaled_radius
use coord_transform_mod, only : xyz2llr, central_angle
@@ -29,6 +31,7 @@ module analytic_moisture_profiles_mod
private
public :: analytic_moisture
+public :: analytic_rel_humidity
contains
@@ -96,6 +99,13 @@ function analytic_moisture(chi, temperature, pressure, choice) result(moisture)
moisture = rel_hum * mr_sat / &
( 1.0_r_def + (1.0_r_def - rel_hum)*mr_sat*recip_epsilon)
+ case ( test_squall_line )
+
+ rel_hum = analytic_rel_humidity(chi, choice)
+ moisture = rel_hum * mr_sat / &
+ ( 1.0_r_def + (1.0_r_def - rel_hum)*mr_sat*recip_epsilon)
+
+
! Isothermal **dry** atmosphere
case( test_isot_dry_atm )
moisture = 0.0_r_def
@@ -108,4 +118,38 @@ function analytic_moisture(chi, temperature, pressure, choice) result(moisture)
end function analytic_moisture
+
+function analytic_rel_humidity(chi, choice) result(rel_hum)
+
+ implicit none
+ real(kind=r_def), intent(in) :: chi(3)
+ integer(kind=i_def), intent(in) :: choice
+ real(kind=r_def) :: rel_hum
+ real(kind=r_def) :: z, long, lat, radius
+ real(kind=r_def), parameter :: z_trop = 12000.0_r_def
+
+ if ( geometry == geometry_spherical ) then
+ call xyz2llr(chi(1), chi(2), chi(3), long, lat, radius)
+ z = MAX(radius - scaled_radius, 0.0_r_def)
+ else
+ z = chi(3)
+ end if
+
+ select case( choice )
+
+ case ( test_squall_line, test_supercell )
+ if (z < z_trop) then
+ rel_hum = 0.999_r_def - 0.749_r_def * (z / z_trop)**1.25_r_def
+ else
+ rel_hum = 0.25_r_def
+ end if
+
+ case default
+ ! In other cases, relative humidity is set to zero
+ rel_hum = 0.0_r_def
+
+ end select
+
+end function analytic_rel_humidity
+
end module analytic_moisture_profiles_mod
diff --git a/science/gungho/source/initialisation/analytic_pressure_profiles_mod.F90 b/science/gungho/source/initialisation/analytic_pressure_profiles_mod.F90
index b88c07755..6dc819a3c 100644
--- a/science/gungho/source/initialisation/analytic_pressure_profiles_mod.F90
+++ b/science/gungho/source/initialisation/analytic_pressure_profiles_mod.F90
@@ -34,15 +34,21 @@ module analytic_pressure_profiles_mod
test_cosine_bubble, &
test_specified_profiles, &
test_bryan_fritsch, &
- test_grabowski_clark
+ test_grabowski_clark, &
+ test_horizontal_mountain, &
+ test_breaking_gw, &
+ test_squall_line, &
+ test_supercell
use initial_density_config_mod, only : r1, x1, y1, z1, r2, x2, y2, z2, &
density_max, density_background
use base_mesh_config_mod, only : geometry, &
geometry_spherical
-use planet_config_mod, only : p_zero, Rd, kappa, scaled_radius
+use planet_config_mod, only : p_zero, Rd, kappa, scaled_radius, &
+ scaled_omega, gravity, cp
use reference_profile_mod, only : reference_profile
use analytic_temperature_profiles_mod, only: analytic_temperature
use deep_baroclinic_wave_mod, only : deep_baroclinic_wave
+use breaking_gravity_wave_mod, only : breaking_gravity_wave
implicit none
@@ -109,12 +115,13 @@ end function vortex_field
!> @param[in] chi Position in physical coordinates
!> @param[in] choice Integer defining which specified formula to use
!> @result pressure The result pressure field
- function analytic_pressure(chi, choice, time) result(pressure)
+ function analytic_pressure(chi, choice, time, surface_height) result(pressure)
implicit none
real(kind=r_def), intent(in) :: chi(3)
integer, intent(in) :: choice
real(kind=r_def), intent(in) :: time
+ real(kind=r_def), intent(in) :: surface_height
real(kind=r_def) :: pressure
real(kind=r_def) :: l, dt
@@ -127,6 +134,8 @@ function analytic_pressure(chi, choice, time) result(pressure)
real(kind=r_def) :: density, temperature, mr_v
real(kind=r_def) :: t0
real(kind=r_def) :: u, v, w
+ real(kind=r_def) :: p_surf, Phi_s, Nsq
+ real(kind=r_def) :: psp, u0
integer :: id
@@ -163,10 +172,27 @@ function analytic_pressure(chi, choice, time) result(pressure)
!> No perturbation needed for warm bubble tests so just use background
!> (isentropic) value
- case (test_warm_bubble, test_warm_bubble_3d, &
- test_bryan_fritsch, test_grabowski_clark )
+ case (test_warm_bubble, test_warm_bubble_3d, test_supercell, &
+ test_bryan_fritsch, test_grabowski_clark, test_squall_line)
call reference_profile(pressure, density, temperature, chi, choice)
+ case ( test_horizontal_mountain )
+ T0 = 288.0_r_def
+ Nsq = gravity**2/(cp*T0)
+ u0 = 10.0_r_def
+ psp = 100000.0_r_def
+ Phi_s = gravity*surface_height
+
+ ! Specify surface pressure
+ p_surf = psp * EXP( &
+ -scaled_radius*Nsq*u0/(2.0_r_def*gravity**2*kappa) &
+ * (u0 / scaled_radius + 2.0_r_def*scaled_omega) * (SIN(lat))**2 &
+ - Nsq*Phi_s/(gravity**2*kappa) &
+ )
+
+ ! Obtain vertical pressure given isothermal atmosphere
+ pressure = (p_surf * EXP( -gravity*(radius-scaled_radius) / (Rd*T0) ) / p_zero) ** kappa
+
case (test_constant_field)
pressure = density_background
@@ -184,7 +210,7 @@ function analytic_pressure(chi, choice, time) result(pressure)
case (test_solid_body_rotation, &
test_solid_body_rotation_alt)
t0 = 280.0_r_def
- temperature = analytic_temperature(chi, choice)
+ temperature = analytic_temperature(chi, choice, 0.0_r_def)
pressure = t0/temperature
density = p_zero/(Rd*temperature) * pressure**( (1.0_r_def - kappa )/ kappa )
case (test_deep_baroclinic_wave)
@@ -192,6 +218,11 @@ function analytic_pressure(chi, choice, time) result(pressure)
pressure, temperature, density, &
u, v, w, mr_v)
+ case (test_breaking_gw)
+ call breaking_gravity_wave(long, lat, radius-scaled_radius, &
+ pressure, temperature, density, &
+ u, v, w, mr_v)
+
case( test_cos_phi )
pressure = density_max*cos(lat)**4
diff --git a/science/gungho/source/initialisation/analytic_temperature_profiles_mod.F90 b/science/gungho/source/initialisation/analytic_temperature_profiles_mod.F90
index 84cf5b808..fbf6bae7c 100644
--- a/science/gungho/source/initialisation/analytic_temperature_profiles_mod.F90
+++ b/science/gungho/source/initialisation/analytic_temperature_profiles_mod.F90
@@ -9,7 +9,7 @@
!! point based upon a specified analytic formula
module analytic_temperature_profiles_mod
-use constants_mod, only : r_def, pi
+use constants_mod, only : r_def, pi, i_def, l_def
use log_mod, only : log_event, &
log_scratch_space, &
LOG_LEVEL_ERROR
@@ -25,7 +25,11 @@ module analytic_temperature_profiles_mod
test_cos_phi, &
test_cosine_bubble, &
test_bryan_fritsch, &
- test_grabowski_clark
+ test_grabowski_clark, &
+ test_squall_line, &
+ test_supercell, &
+ test_horizontal_mountain, &
+ test_breaking_gw
use initial_density_config_mod, only : r1, x1, y1, r2, x2, y2, &
density_max, density_background
use initial_pressure_config_mod, only : surface_pressure
@@ -37,6 +41,7 @@ module analytic_temperature_profiles_mod
use generate_global_gw_fields_mod, only : generate_global_gw_pert
use initial_wind_config_mod, only : u0, sbr_angle_lat
use deep_baroclinic_wave_mod, only : deep_baroclinic_wave
+use breaking_gravity_wave_mod, only : breaking_gravity_wave
use formulation_config_mod, only : shallow
use extrusion_config_mod, only : domain_height
@@ -45,6 +50,9 @@ module analytic_temperature_profiles_mod
private
public :: analytic_temperature
+public :: analytic_theta_pert
+public :: integrate_theta_v
+public :: integrate_exner_surf
contains
@@ -52,15 +60,16 @@ module analytic_temperature_profiles_mod
!> @param[in] chi Position in physical coordinates
!> @param[in] choice Integer defining which specified formula to use
!> @result temperature The result temperature field
-function analytic_temperature(chi, choice) result(temperature)
+function analytic_temperature(chi, choice, surface_height) result(temperature)
implicit none
real(kind=r_def), intent(in) :: chi(3)
integer, intent(in) :: choice
+ real(kind=r_def), intent(in) :: surface_height
real(kind=r_def) :: temperature
real(kind=r_def) :: l, dt
- real(kind=r_def), parameter :: THETA0 = 0.01_r_def
+ real(kind=r_def) :: theta0
real(kind=r_def), parameter :: XC = 0.0_r_def
real(kind=r_def), parameter :: YC = 0.0_r_def
real(kind=r_def), parameter :: A = 5000.0_r_def
@@ -72,10 +81,14 @@ function analytic_temperature(chi, choice) result(temperature)
ZR = 2000.0_r_def
real(kind=r_def) :: long, lat, radius, z
real(kind=r_def) :: l1, l2
- real(kind=r_def) :: pressure, density, mr_v
+ real(kind=r_def) :: pressure, density, mr_v, exner
real(kind=r_def) :: s, u00, f_sb, t0
real(kind=r_def) :: r_on_a
real(kind=r_def) :: u, v, w
+ real(kind=r_def) :: theta_trop, T_trop, theta_eq
+ real(kind=r_def) :: z_trop
+ real(kind=r_def) :: p_surf, Phi_s, Nsq
+ real(kind=r_def) :: psp, u0
if ( geometry == geometry_spherical ) then
call xyz2llr(chi(1),chi(2),chi(3),long,lat,radius)
@@ -100,6 +113,7 @@ function analytic_temperature(chi, choice) result(temperature)
temperature = temperature &
+ generate_global_gw_pert(long,lat,radius-scaled_radius)
else
+ theta0 = 0.01_r_def
temperature = temperature + THETA0 * sin ( PI * chi(3) / H ) &
/ ( 1.0_r_def + ( chi(1) - XC )**2/A**2 )
end if
@@ -152,6 +166,44 @@ function analytic_temperature(chi, choice) result(temperature)
s = 1.3e-5_r_def ! stability, in m^{-1}
temperature = t0*exp(s*chi(3))
+ ! Test from DCMIP 2025
+ case ( test_squall_line, test_supercell )
+
+ ! Specify potential temperature at equator
+ z = MAX(radius - scaled_radius, 0.0_r_def)
+ theta_trop = 343.0_r_def
+ theta0 = 300.0_r_def
+ T_trop = 213.0_r_def
+ z_trop = 12000.0_r_def
+ if (z < z_trop) then
+ theta_eq = theta0 + (theta_trop - theta0)*(z/z_trop)**1.25_r_def
+ else
+ theta_eq = theta_trop*EXP(gravity*(z-z_trop)/(cp*T_trop))
+ end if
+
+ temperature = theta_eq
+
+ case ( test_horizontal_mountain )
+ T0 = 288.0_r_def
+ Nsq = gravity**2/(cp*T0)
+ u0 = 10.0_r_def
+ psp = 100000.0_r_def
+ Phi_s = gravity*surface_height
+
+ ! Specify surface pressure
+ p_surf = psp * EXP( &
+ -scaled_radius*Nsq*u0/(2.0_r_def*gravity**2*kappa) &
+ * (u0 / scaled_radius + 2.0_r_def*scaled_omega) * (SIN(lat))**2 &
+ - Nsq*Phi_s/(gravity**2*kappa) &
+ )
+
+ ! Obtain vertical pressure given isothermal atmosphere
+ pressure = p_surf * EXP( -gravity*(radius-scaled_radius) / (Rd*T0) )
+
+ ! Obtain potential temperature
+ exner = (pressure / p_zero) ** kappa
+ temperature = T0 / exner
+
!> Test from Kelly & Giraldo
case( test_warm_bubble_3d )
! Set default value
@@ -204,6 +256,11 @@ function analytic_temperature(chi, choice) result(temperature)
pressure, temperature, density, &
u, v, w, mr_v)
+ case( test_breaking_gw )
+ call breaking_gravity_wave(long, lat, radius-scaled_radius, &
+ pressure, temperature, density, &
+ u, v, w, mr_v)
+
case( test_cos_phi )
temperature = density_max*cos(lat)**4
@@ -223,4 +280,237 @@ function analytic_temperature(chi, choice) result(temperature)
end function analytic_temperature
+function analytic_theta_pert(chi, choice) result(theta_pert)
+
+ use extrusion_config_mod, only: planet_radius
+
+ implicit none
+ real(kind=r_def), intent(in) :: chi(3)
+ integer(kind=i_def), intent(in) :: choice
+ real(kind=r_def) :: theta_pert
+ real(kind=r_def) :: long, lat, radius, l, z
+ real(kind=r_def) :: lon_p, lat_p, Rtheta, ds, r_p
+ real(kind=r_def), parameter :: lon_c = 2.0_r_def*PI/3.0_r_def
+ real(kind=r_def), parameter :: lat_c = 0.0_r_def
+ real(kind=r_def), parameter :: z_c = 1500.0_r_def
+ real(kind=r_def), parameter :: z_p = 1500.0_r_def
+ real(kind=r_def), parameter :: dtheta = 3.0_r_def
+ integer(kind=i_def), parameter :: num_bubbles = 9
+
+ integer(kind=i_def) :: i
+
+ if ( geometry == geometry_spherical ) then
+ call xyz2llr(chi(1), chi(2), chi(3), long, lat, radius)
+ else
+ long = chi(1)
+ lat = chi(2)
+ end if
+
+ select case( choice )
+
+ case ( test_squall_line )
+
+ z = radius - scaled_radius
+ theta_pert = 0.0_r_def
+ ds = 10000.0_r_def
+ r_p = 5000.0_r_def
+
+ do i = 1, num_bubbles
+ lon_p = lon_c
+ lat_p = lat_c + (real(i, r_def) - 0.5_r_def * real(num_bubbles + 1, r_def)) * ds / scaled_radius
+ call central_angle(long, lat, lon_p, lat_p, l)
+ Rtheta = SQRT((l * scaled_radius / r_p)**2 + ((z - z_c) / z_p)**2)
+
+ if (Rtheta <= 1.0_r_def) then
+ theta_pert = theta_pert + dtheta * (COS(PI/2.0_r_def*Rtheta))**2
+ end if
+ end do
+
+ case ( test_supercell )
+
+ r_p = 10000.0_r_def
+
+ z = radius - scaled_radius
+ theta_pert = 0.0_r_def
+
+ call central_angle(long, lat, lon_c, lat_c, l)
+ Rtheta = SQRT((l * scaled_radius / r_p)**2 + ((z - z_c) / z_p)**2)
+
+ if (Rtheta <= 1.0_r_def) then
+ theta_pert = theta_pert + dtheta * (COS(PI/2.0_r_def*Rtheta))**2
+ end if
+
+ case default
+ ! In other cases, relative humidity is set to zero
+ theta_pert = 0.0_r_def
+
+ end select
+
+end function analytic_theta_pert
+
+
+subroutine integrate_theta_v(theta_v, theta_v_prev, dtheta_v_dz, theta_eq, &
+ dusq_eq_dz, height_wth, lat_points, nl, num_quad_points, test)
+
+ use planet_config_mod, only: cp, gravity
+
+ implicit none
+
+ integer(kind=i_def), intent(in) :: nl, test
+ integer(kind=i_def), intent(in) :: num_quad_points
+ real(kind=r_def), intent(in) :: theta_v_prev(nl, num_quad_points)
+ real(kind=r_def), intent(in) :: theta_eq(nl)
+ real(kind=r_def), intent(in) :: dtheta_v_dz(nl, num_quad_points)
+ real(kind=r_def), intent(in) :: height_wth(nl)
+ real(kind=r_def), intent(in) :: lat_points(num_quad_points+1)
+ real(kind=r_def), intent(inout) :: theta_v(nl, num_quad_points)
+ real(kind=r_def), intent(in) :: dusq_eq_dz(nl)
+
+
+ real(kind=r_def) :: u_eq(nl)
+ real(kind=r_def) :: dusq_eq_dz_2(nl)
+ real(kind=r_def) :: z, zS, dzU, Us, Uc
+ real(kind=r_def) :: A, B, C
+ real(kind=r_def) :: integrand(nl, num_quad_points)
+ real(kind=r_def) :: integral(nl, num_quad_points)
+ integer(kind=i_def) :: k, j
+
+ ! Squall line test case from DCMIP 2025
+ if ( test == test_squall_line ) then
+ zS = 2500.0_r_def
+ dzU = 1000.0_r_def
+ Us = 12.0_r_def
+ Uc = 5.0_r_def
+
+ ! Supercell test case from DCMIP 2016
+ else if ( test == test_supercell ) then
+ zS = 5000.0_r_def
+ dzU = 1000.0_r_def
+ Us = 30.0_r_def
+ Uc = 15.0_r_def
+ end if
+
+ ! Zonal solid body rotation but with wind shear
+ A = 0.25_r_def * (-zs**2 + 2.0_r_def*zS*dzU - dzU**2) / (zS * dzU)
+ B = 0.5_r_def * (zS + dzU) / dzU
+ C = - 0.25_r_def * (zS / dzU)
+
+ do j = 1, num_quad_points
+ do k = 1, nl
+ ! Evaluate analytic wind
+ z = height_wth(k)
+ if (z < zS - dzU + 1.0E-6_r_def) then
+ u_eq(k) = Us*z/zS - Uc
+ dusq_eq_dz_2(k) = 2.0_r_def*u_eq(k)*Us/zs
+ else if (ABS(z - zS) < dzU + 1.0E-6_r_def) then
+ u_eq(k) = (A + B*z/zS + C*(z/zS)**2)*Us - Uc
+ dusq_eq_dz_2(k) = 2.0_r_def*u_eq(k)*Us*(B/zs + 2.0_r_def*C*z/(zs**2))
+ else
+ u_eq(k) = Us - Uc
+ dusq_eq_dz_2(k) = 0.0_r_def
+ end if
+ ! Evaluate integrand
+ integrand(k,j) = ( &
+ 0.5_r_def * SIN(2.0_r_def*lat_points(j)) / gravity &
+ * (u_eq(k)**2 * dtheta_v_dz(k,j) - theta_v_prev(k,j) * dusq_eq_dz_2(k)) &
+ )
+ end do
+ end do
+
+ ! Integrate the integrand over the latitude points
+ integral(:,:) = 0.0_r_def
+ do j = 2, num_quad_points
+ do k = 1, nl
+ ! Use trapezoidal rule for integration
+ integral(k,j) = integral(k,j-1) &
+ + 0.5_r_def * (integrand(k,j) + integrand(k,j-1)) &
+ * (lat_points(j) - lat_points(j-1))
+ end do
+ end do
+
+ do j = 1, num_quad_points
+ theta_v(:,j) = theta_eq(:) + integral(:,j)
+ end do
+
+end subroutine integrate_theta_v
+
+subroutine integrate_exner_surf(exner_wt, theta_v, exner_eq, &
+ height_wth, lat_points, nl, num_quad_points, test)
+
+ use planet_config_mod, only: cp, gravity
+
+ implicit none
+
+ integer(kind=i_def), intent(in) :: nl, test
+ integer(kind=i_def), intent(in) :: num_quad_points
+ real(kind=r_def), intent(in) :: exner_eq(nl)
+ real(kind=r_def), intent(in) :: lat_points(num_quad_points+1)
+ real(kind=r_def), intent(in) :: theta_v(nl, num_quad_points)
+ real(kind=r_def), intent(in) :: height_wth(nl)
+ real(kind=r_def), intent(inout) :: exner_wt(nl, num_quad_points)
+
+ real(kind=r_def) :: u_eq(nl)
+ real(kind=r_def) :: z, zS, dzU, Us, Uc
+ real(kind=r_def) :: A, B, C
+ real(kind=r_def) :: integrand(nl, num_quad_points)
+ real(kind=r_def) :: integral(nl, num_quad_points)
+ integer(kind=i_def) :: j, k
+
+ ! Squall line test case from DCMIP 2025
+ if ( test == test_squall_line ) then
+ zS = 2500.0_r_def
+ dzU = 1000.0_r_def
+ Us = 12.0_r_def
+ Uc = 5.0_r_def
+
+ ! Supercell test case from DCMIP 2016
+ else if ( test == test_supercell ) then
+ zS = 5000.0_r_def
+ dzU = 1000.0_r_def
+ Us = 30.0_r_def
+ Uc = 15.0_r_def
+ end if
+
+ ! Zonal solid body rotation but with wind shear
+ A = 0.25_r_def * (-zs**2 + 2.0_r_def*zS*dzU - dzU**2) / (zS * dzU)
+ B = 0.5_r_def * (zS + dzU) / dzU
+ C = - 0.25_r_def * (zS / dzU)
+
+ do j = 1, num_quad_points
+ do k = 1, nl
+ ! Evaluate analytic wind
+ z = height_wth(k)
+ if (z < zS - dzU) then
+ u_eq(k) = Us*z/zS - Uc
+ else if (ABS(z - zS) < dzU) then
+ u_eq(k) = (A + B*z/zS + C*(z/zS)**2)*Us - Uc
+ else
+ u_eq(k) = Us - Uc
+ end if
+
+ ! Evaluate integrand
+ integrand(k,j) = ( &
+ u_eq(k)**2 * SIN(lat_points(j)) * COS(lat_points(j)) &
+ / (cp * theta_v(k,j)) &
+ )
+ end do
+ end do
+
+ ! Integrate the integrand over the latitude points
+ integral(:,:) = 0.0_r_def
+ do j = 2, num_quad_points
+ do k = 1, nl
+ ! Use trapezoidal rule for integration
+ integral(k,j) = integral(k,j-1) &
+ + 0.5_r_def * (integrand(k,j) + integrand(k,j-1)) &
+ * (lat_points(j) - lat_points(j-1))
+ end do
+ end do
+
+ do j = 1, num_quad_points
+ exner_wt(:,j) = exner_eq(:) - integral(:,j)
+ end do
+
+end subroutine integrate_exner_surf
+
end module analytic_temperature_profiles_mod
diff --git a/science/gungho/source/initialisation/analytic_wind_profiles_mod.F90 b/science/gungho/source/initialisation/analytic_wind_profiles_mod.F90
index f4392c006..0fda840d3 100644
--- a/science/gungho/source/initialisation/analytic_wind_profiles_mod.F90
+++ b/science/gungho/source/initialisation/analytic_wind_profiles_mod.F90
@@ -32,12 +32,17 @@ module analytic_wind_profiles_mod
profile_dcmip_101, &
profile_vertical_deformation, &
profile_four_part_sbr, &
+ profile_breaking_gw, &
+ profile_horizontal_mountain, &
+ profile_squall_line, &
+ profile_supercell, &
wind_time_period
use planet_config_mod, only: scaled_radius
use log_mod, only: log_event, &
log_scratch_space, &
LOG_LEVEL_ERROR
use deep_baroclinic_wave_mod, only: deep_baroclinic_wave
+use breaking_gravity_wave_mod, only: breaking_gravity_wave
use formulation_config_mod, only: shallow
implicit none
@@ -454,7 +459,9 @@ function analytic_wind( chi, time, choice, num_options, &
real(kind=r_def), parameter :: p0 = 100000.0_r_def ! surface pressure for DCMIP 1.1
real(kind=r_def), parameter :: Rd = 287.0_r_def ! gas constant for DCMIP 1.1
real(kind=r_def), parameter :: T0 = 300.0_r_def ! ref temperature for DCMIP 1.1
- real(kind=r_def), parameter :: g = 9.80616_r_def ! gravity for DCMIP 1.1
+ real(kind=r_def), parameter :: g = 9.80616_r_def ! gravity for DCMIP 1.1
+ real(kind=r_def) :: A, C ! variables for DCMIP tests
+ real(kind=r_def) :: zS, dzU, Us, Uc ! variables for squall line and supercell
integer(kind=i_def) :: num_periods
if ( .not. present(option_arg) ) then
@@ -484,7 +491,7 @@ function analytic_wind( chi, time, choice, num_options, &
- cos(lat_pole) * cos(chi(1)-lon_pole) * sin(chi(2)) )
u(2) = s * option(1) * cos(lat_pole) * sin(chi(1)-lon_pole)
u(3) = 0.0_r_def
- case ( profile_solid_body_rotation_alt)
+ case ( profile_solid_body_rotation_alt )
! Rotated pole version of equation (74) of Staniforth & White (2007)
! with m=1, A=0, n therefore arbitrary.
! In shallow geometry the r/a factor is replaced by 1 (see section 6
@@ -520,6 +527,47 @@ function analytic_wind( chi, time, choice, num_options, &
- cos(lat_pole) * cos(chi(1)-lon_pole) * sin(chi(2)) )
u(2) = option(1) * cos(lat_pole) * sin(chi(1)-lon_pole)
u(3) = 0.0_r_def
+
+ case ( profile_horizontal_mountain )
+ ! Zonal solid body rotation, but constant with height
+ u0 = 10.0_r_def
+ u(1) = u0*COS(chi(2))
+ u(2) = 0.0_r_def
+ u(3) = 0.0_r_def
+
+ case ( profile_squall_line, profile_supercell )
+ ! Squall line test case from DCMIP 2025
+ if ( choice == profile_squall_line ) then
+ zS = 2500.0_r_def
+ dzU = 1000.0_r_def
+ Us = 12.0_r_def
+ Uc = 5.0_r_def
+
+ ! Supercell test case from DCMIP 2016
+ else if ( choice == profile_supercell ) then
+ zS = 5000.0_r_def
+ dzU = 1000.0_r_def
+ Us = 30.0_r_def
+ Uc = 15.0_r_def
+ end if
+
+ ! Zonal solid body rotation but with wind shear
+ z = chi(3) - scaled_radius
+ A = 0.25_r_def * (-zs**2 + 2.0_r_def*zS*dzU - dzU**2) / (zS * dzU)
+ B = 0.5_r_def * (zS + dzU) / dzU
+ C = - 0.25_r_def * (zS / dzU)
+
+ if (z < zS - dzU) then
+ u(1) = (Us*z/zS - Uc)*COS(chi(2))
+ else if (ABS(z - zS) < dzU) then
+ u(1) = ((A + B*z/zS + C*(z/zS)**2)*Us - Uc)*COS(chi(2))
+ else
+ u(1) = (Us - Uc)*COS(chi(2))
+ end if
+
+ u(2) = 0.0_r_def
+ u(3) = 0.0_r_def
+
case ( profile_constant_uv )
u(1) = option(1)
u(2) = option(2)
@@ -538,6 +586,10 @@ function analytic_wind( chi, time, choice, num_options, &
call deep_baroclinic_wave(chi(1), chi(2), chi(3)-scaled_radius, &
pressure, temperature, density, &
u(1), u(2), u(3), mr_v)
+ case ( profile_breaking_gw )
+ call breaking_gravity_wave(chi(1), chi(2), chi(3)-scaled_radius, &
+ pressure, temperature, density, &
+ u(1), u(2), u(3), mr_v)
case ( profile_vortex )
u = vortex_wind(chi(2),chi(1),chi(3))
case ( profile_NL_case_1 )
diff --git a/science/gungho/source/initialisation/breaking_gravity_wave_mod.F90 b/science/gungho/source/initialisation/breaking_gravity_wave_mod.F90
new file mode 100644
index 000000000..282d0f662
--- /dev/null
+++ b/science/gungho/source/initialisation/breaking_gravity_wave_mod.F90
@@ -0,0 +1,180 @@
+!-----------------------------------------------------------------------------
+! Copyright (c) 2017, Met Office, on behalf of HMSO and Queen's Printer
+! For further details please refer to the file LICENCE.original which you
+! should have received as part of this distribution.
+!-----------------------------------------------------------------------------
+
+!> @brief Functions to compute the analytic profiles of the deep atmosphere
+!> baroclinic wave of Ullrich et.al. QJRMS2013
+module breaking_gravity_wave_mod
+
+ use constants_mod, only : r_def, i_def, pi
+ use planet_config_mod, only : scaled_radius, gravity, Rd, cp, &
+ scaled_omega, kappa, p_zero
+ use formulation_config_mod, only : shallow
+ use initial_wind_config_mod, only : profile, &
+ profile_breaking_gw
+ implicit none
+
+ private
+
+ public :: breaking_gravity_wave, T0E, T0P, B, KK, lapse, pertu0, &
+ pertr, pertlon, pertlat, pertz, dxepsilon
+
+!=======================================================================
+! Deep baroclinic wave parameters
+!=======================================================================
+ real(kind=r_def), parameter :: &
+ T0E = 310.0_r_def, & ! temperature at equatorial surface (K)
+ T0P = 240.0_r_def, & ! temperature at polar surface (K)
+ B = 2.0_r_def , & ! jet half-width parameter
+ KK = 3.0_r_def , & ! jet width parameter
+ lapse = 0.005_r_def, & ! lapse rate parameter
+ z_strat = 20000.0_r_def, & ! Tropopause height (m)
+ lapse_strat = -0.005_r_def ! Stratosphere lapse rate (K m^-1)
+
+ real(kind=r_def), parameter :: &
+ pertu0 = 0.5_r_def, & ! Perturbation wind velocity (m/s)
+ pertr = 1.0_r_def/6.0_r_def, & ! Perturbation radius (Earth radii)
+ pertlon = pi/9.0_r_def, & ! Perturbation longitude
+ pertlat = 2.0_r_def*pi/9.0_r_def,& ! Perturbation latitude
+ pertz = 15000.0_r_def , & ! Perturbation height cap
+ dxepsilon = 1.0e-5_r_def ! Small value for numerical derivatives
+
+
+contains
+
+ subroutine breaking_gravity_wave(lon, lat, z, &
+ exner, theta, rho, &
+ u, v, w, mr_v)
+
+ implicit none
+
+!-----------------------------------------------------------------------
+! input/output params parameters at given location
+!-----------------------------------------------------------------------
+ real(kind=r_def), intent(in) :: &
+ lon, & ! Longitude (radians)
+ lat, & ! Latitude (radians)
+ z ! Altitude (m)
+
+
+ real(kind=r_def), intent(out) :: &
+ exner, & ! Exner pressure
+ theta, & ! Potential Temperature (K)
+ rho, & ! density (kg m^-3)
+ u, & ! Zonal wind (m s^-1)
+ v, & ! Meridional wind (m s^-1)
+ w, & ! Vertical Velocity (m s^-1)
+ mr_v ! Water vapour mixing ratio (kg kg^-1)
+
+
+ real(kind=r_def) :: p, t, q, eta, zeta
+ real(kind=r_def) :: T0, constA, constB, constC, constH, scaledZ
+ real(kind=r_def) :: tau1, tau2, inttau1, inttau2, intzeta
+ real(kind=r_def) :: rratio, inttermT, inttermU, bigU, rcoslat, omegarcoslat
+
+!-----------------------------------------------------------------------
+! constants / parameters
+!-----------------------------------------------------------------------
+
+ T0 = 0.5_r_def * (T0E + T0P)
+ constA = 1.0_r_def / lapse
+ constB = (T0 - T0P) / (T0 * T0P)
+ constC = 0.5_r_def * (KK + 2.0_r_def) * (T0E - T0P) / (T0E * T0P)
+ constH = Rd * T0 / gravity
+ scaledZ = z / (B * constH)
+
+!-----------------------------------------------------------------------
+! tau values
+!-----------------------------------------------------------------------
+ if (z < z_strat) then
+ zeta = constA * lapse / T0 * exp(lapse * z / T0)
+ intzeta = constA * (exp(lapse * z / T0) - 1.0_r_def)
+ else
+ zeta = constA * lapse / T0 * exp(lapse * z_strat / T0) &
+ + constA * lapse / T0 * (exp(lapse_strat * (z - z_strat) / T0) - 1.0_r_def)
+ intzeta = constA * (exp(lapse * z_strat / T0) - 1.0_r_def) &
+ + 1.0_r_def / T0 * (z-z_strat) * (exp(lapse * z_strat / T0) - 1.0_r_def) &
+ + 1.0_r_def / lapse_strat * (exp(lapse_strat * (z - z_strat) / T0) - 1.0_r_def)
+ end if
+ tau1 = zeta + constB * (1.0_r_def - 2.0_r_def * scaledZ**2) * exp(- scaledZ**2)
+ tau2 = constC * (1.0_r_def - 2.0_r_def * scaledZ**2) * exp(- scaledZ**2)
+
+ inttau1 = intzeta + constB * z * exp(- scaledZ**2)
+ inttau2 = constC * z * exp(- scaledZ**2)
+
+!-----------------------------------------------------------------------
+! radius ratio
+!-----------------------------------------------------------------------
+ if ( shallow ) then
+ rratio = 1.0_r_def
+ else
+ rratio = (z + scaled_radius) / scaled_radius;
+ end if
+
+!-----------------------------------------------------------------------
+! interior term on temperature expression
+!-----------------------------------------------------------------------
+ inttermT = (rratio * cos(lat))**int(KK,i_def) &
+ - KK / (KK + 2.0_r_def) &
+ * (rratio * cos(lat))**int(KK + 2.0_r_def,i_def)
+
+!-----------------------------------------------------------------------
+! temperature
+!-----------------------------------------------------------------------
+ t = 1.0_r_def / (rratio**2 * (tau1 - tau2 * inttermT))
+
+!-----------------------------------------------------------------------
+! hydrostatic pressure
+!-----------------------------------------------------------------------
+ p = p_zero * exp(- gravity / Rd * (inttau1 - inttau2 * inttermT))
+
+!-----------------------------------------------------------------------
+! density (via ideal gas law)
+!-----------------------------------------------------------------------
+ rho = p / (Rd * t)
+!-----------------------------------------------------------------------
+! Compute exner pressure
+!-----------------------------------------------------------------------
+ exner = (p/p_zero)**kappa
+
+!-----------------------------------------------------------------------
+! Compute potential temperaure
+!-----------------------------------------------------------------------
+ theta = t/exner
+
+!-----------------------------------------------------------------------
+! velocity field
+!-----------------------------------------------------------------------
+ inttermU = (rratio * cos(lat))**int(KK - 1.0_r_def,i_def) &
+ - (rratio * cos(lat))**int(KK + 1.0_r_def,i_def)
+ bigU = gravity / scaled_radius * KK * inttau2 * inttermU * t
+ if ( shallow ) then
+ rcoslat = scaled_radius * cos(lat)
+ else
+ rcoslat = (z + scaled_radius) * cos(lat)
+ end if
+
+ omegarcoslat = scaled_omega * rcoslat
+
+ u = - omegarcoslat + sqrt(omegarcoslat**2 + rcoslat * bigU)
+ v = 0.0_r_def
+ w = 0.0_r_def
+
+
+ eta = p/p_zero
+
+ if ( eta > 1.0e4_r_def/p_zero) then
+ q = 1.8e-2_r_def*exp(-(lat/(2.0_r_def*pi/9.0_r_def))**4.0_r_def) &
+ *exp(-((eta-1.0_r_def)*p_zero/340.0e2_r_def)**2.0_r_def )
+ else
+ q = 1.0e-12_r_def
+ end if
+
+ mr_v = q/(1.0_r_def - q)
+
+
+ end subroutine breaking_gravity_wave
+
+end module breaking_gravity_wave_mod
diff --git a/science/gungho/source/kernel/core_dynamics/compute_trad_coriolis_matrix_kernel_mod.F90 b/science/gungho/source/kernel/core_dynamics/compute_trad_coriolis_matrix_kernel_mod.F90
new file mode 100644
index 000000000..07dc9981a
--- /dev/null
+++ b/science/gungho/source/kernel/core_dynamics/compute_trad_coriolis_matrix_kernel_mod.F90
@@ -0,0 +1,218 @@
+!-----------------------------------------------------------------------------
+! (C) Crown copyright 2018 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 Compute the Traditional Coriolis operator
+!> @details The form of the Tradiational Coriolis operator is:
+!! \f[ < v, H[2\Omega \times H[u]] > \f] where v is the test function
+!! for the velocity space W2, u corresponds to the wind field in W2 to
+!! which this operator will be applied, and Omega is the rotation
+!! vector of the domain.
+!! Multiplying this operator by a wind field will assemble the weak
+!! form of the Traditional Coriolis term.
+
+module compute_trad_coriolis_matrix_kernel_mod
+
+use constants_mod, only: i_def, r_def
+use kernel_mod, only: kernel_type
+use argument_mod, only: arg_type, func_type, &
+ GH_OPERATOR, GH_FIELD, &
+ GH_READ, GH_WRITE, &
+ GH_REAL, GH_SCALAR, &
+ ANY_SPACE_9, &
+ ANY_DISCONTINUOUS_SPACE_3, &
+ GH_BASIS, GH_DIFF_BASIS, &
+ CELL_COLUMN, GH_QUADRATURE_XYoZ
+use fs_continuity_mod, only: W2
+
+use sci_coordinate_jacobian_mod, only: coordinate_jacobian
+use base_mesh_config_mod, only: geometry, &
+ geometry_spherical
+use rotation_vector_mod, only: rotation_vector_fplane, &
+ rotation_vector_sphere, &
+ vert_vector_sphere
+use cross_product_mod, only: cross_product
+
+implicit none
+private
+
+!-------------------------------------------------------------------------------
+! Public types
+!-------------------------------------------------------------------------------
+
+type, public, extends(kernel_type) :: compute_trad_coriolis_matrix_kernel_type
+ private
+ type(arg_type) :: meta_args(5) = (/ &
+ arg_type(GH_OPERATOR, GH_REAL, GH_WRITE, W2, W2), &
+ arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &
+ arg_type(GH_SCALAR, GH_REAL, GH_READ), &
+ arg_type(GH_SCALAR, GH_REAL, GH_READ) &
+ /)
+ type(func_type) :: meta_funcs(2) = (/ &
+ func_type(W2, GH_BASIS), &
+ func_type(ANY_SPACE_9, GH_BASIS, GH_DIFF_BASIS) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+ integer :: gh_shape = GH_QUADRATURE_XYoZ
+contains
+ procedure, nopass :: compute_trad_coriolis_matrix_code
+end type
+
+!-------------------------------------------------------------------------------
+! Contained functions/subroutines
+!-------------------------------------------------------------------------------
+public :: compute_trad_coriolis_matrix_code
+contains
+
+!> @brief Compute the Coriolis operator to apply the rotation vector Omega to
+!> the wind fields.
+!!
+!! @param[in] cell Identifying number of cell.
+!! @param[in] nlayers Number of layers.
+!! @param[in] ncell_3d ncell*ndf
+!! @param[in,out] matrix Local stencil or Coriolis operator.
+!! @param[in] chi_1 1st coordinate field
+!! @param[in] chi_2 2nd coordinate field
+!! @param[in] chi_3 3rd coordinate field
+!! @param[in] panel_id A field giving the ID for mesh panels.
+!! @param[in] omega Planet angular velocity
+!! @param[in] f_lat F-plane latitude
+!! @param[in] ndf Degrees of freedom per cell.
+!! @param[in] basis Vector basis functions evaluated at quadrature points.
+!! @param[in] ndf_chi Number of degrees of freedom per cell for chi
+!! @param[in] undf_chi Number of unique degrees of freedom for chi
+!! @param[in] map_chi Dofmap for the cell at the base of the column for chi
+!! @param[in] chi_basis Basis functions for Wchi evaluated at
+!! gaussian quadrature points
+!! @param[in] chi_diff_basis Differential of the Wchi basis functions
+!! evaluated at gaussian quadrature point
+!! @param[in] ndf_pid Number of degrees of freedom per cell for panel_id
+!! @param[in] undf_pid Number of unique degrees of freedom for panel_id
+!! @param[in] map_pid Dofmap for the cell at the base of the column for panel_id
+!! @param[in] nqp_h Number of horizontal quadrature points.
+!! @param[in] nqp_v Number of vertical quadrature points.
+!! @param[in] wqp_h Horizontal quadrature weights.
+!! @param[in] wqp_v Vertical quadrature weights.
+subroutine compute_trad_coriolis_matrix_code(cell, nlayers, ncell_3d, &
+ matrix, &
+ chi_1, chi_2, chi_3, &
+ panel_id, &
+ omega, f_lat, &
+ ndf, basis, &
+ ndf_chi, undf_chi, &
+ map_chi, &
+ basis_chi, diff_basis_chi, &
+ ndf_pid, undf_pid, map_pid, &
+ nqp_h, nqp_v, wqp_h, wqp_v)
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: cell, ndf, ndf_pid, ndf_chi
+ integer(kind=i_def), intent(in) :: undf_pid, undf_chi
+ integer(kind=i_def), intent(in) :: nqp_h, nqp_v
+ integer(kind=i_def), intent(in) :: nlayers
+ integer(kind=i_def), intent(in) :: ncell_3d
+
+ real(kind=r_def), dimension(3,ndf,nqp_h,nqp_v), intent(in) :: basis
+
+ integer(kind=i_def), intent(in) :: map_chi(ndf_chi)
+ integer(kind=i_def), intent(in) :: map_pid(ndf_pid)
+ real(kind=r_def), intent(inout) :: matrix(ncell_3d,ndf,ndf)
+ real(kind=r_def), intent(in) :: basis_chi(1,ndf_chi,nqp_h,nqp_v)
+ real(kind=r_def), intent(in) :: diff_basis_chi(3,ndf_chi,nqp_h,nqp_v)
+ real(kind=r_def), intent(in) :: chi_1(undf_chi)
+ real(kind=r_def), intent(in) :: chi_2(undf_chi)
+ real(kind=r_def), intent(in) :: chi_3(undf_chi)
+ real(kind=r_def), intent(in) :: panel_id(undf_pid)
+ real(kind=r_def), intent(in) :: wqp_h(nqp_h)
+ real(kind=r_def), intent(in) :: wqp_v(nqp_v)
+ real(kind=r_def), intent(in) :: omega
+ real(kind=r_def), intent(in) :: f_lat
+
+ ! Internal variables
+ integer(kind=i_def) :: df, df2, k, ik
+ integer(kind=i_def) :: qp1, qp2
+
+ real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
+ real(kind=r_def), dimension(nqp_h,nqp_v) :: dj
+ real(kind=r_def), dimension(3,3,nqp_h,nqp_v) :: jac
+ real(kind=r_def), dimension(3,nqp_h,nqp_v) :: rotation_vector
+ real(kind=r_def), dimension(3,nqp_h,nqp_v) :: vert_vec
+ real(kind=r_def), dimension(3) :: omega_cross_u
+ real(kind=r_def), dimension(3) :: jac_u, jac_v
+
+ integer(kind=i_def) :: ipanel
+
+ ipanel = int(panel_id(map_pid(1)), i_def)
+
+ ! Loop over layers: Start from 1 as in this loop k is not an offset
+ do k = 1, nlayers
+
+ ! Indirect the chi coord field here
+ do df = 1, ndf_chi
+ chi_1_e(df) = chi_1(map_chi(df) + k - 1)
+ chi_2_e(df) = chi_2(map_chi(df) + k - 1)
+ chi_3_e(df) = chi_3(map_chi(df) + k - 1)
+ end do
+
+ ! Calculate planet's rotation vector
+ if ( geometry == geometry_spherical ) then
+ call rotation_vector_sphere(ndf_chi, nqp_h, nqp_v, chi_1_e, chi_2_e, &
+ chi_3_e, ipanel, basis_chi, rotation_vector)
+ call vert_vector_sphere(ndf_chi, nqp_h, nqp_v, chi_1_e, chi_2_e, &
+ chi_3_e, ipanel, basis_chi, vert_vec)
+ else
+ call rotation_vector_fplane(nqp_h, nqp_v, omega, f_lat, rotation_vector)
+ do qp2 = 1, nqp_v
+ do qp1 = 1, nqp_h
+ vert_vec(:,qp1,qp2) = (/ 0.0_r_def, 0.0_r_def, 1.0_r_def /)
+ end do
+ end do
+ end if
+
+ ! Calculate the Jacobian and its determinant
+ call coordinate_jacobian(ndf_chi, nqp_h, nqp_v, &
+ chi_1_e, chi_2_e, chi_3_e, ipanel, &
+ basis_chi, diff_basis_chi, jac, dj)
+
+ ! To convert from reference space to physical space:
+ ! Let the volume element be dV and the unit volume be dV_hat, then:
+ ! u = J u_hat/detJ, v = J v_hat/detJ, and dV = detJ dV_hat,
+ ! where hats denote the reference element space, J is the Jacobian, and
+ ! detJ its determinant. This gives
+ ! \int (J w_hat) . (2 Omega x J u_hat) / detJ dV_hat
+ ! Note that some of the detJ factors cancel
+ ik = k + (cell-1)*nlayers
+ matrix(ik,:,:) = 0.0_r_def
+ do qp2 = 1, nqp_v
+ do qp1 = 1, nqp_h
+ do df2 = 1, ndf
+ jac_u = matmul(jac(:,:,qp1,qp2),basis(:,df2,qp1,qp2))
+
+ ! Extract only horizontal part of u
+ jac_u = jac_u &
+ - vert_vec(:,qp1,qp2) * dot_product(vert_vec(:,qp1,qp2), jac_u)
+
+ omega_cross_u = wqp_h(qp1) * wqp_v(qp2) &
+ * cross_product(rotation_vector(:,qp1,qp2), jac_u) / dj(qp1,qp2)
+
+ ! Extract only horizontal part of omega_cross_u
+ omega_cross_u = omega_cross_u &
+ - vert_vec(:,qp1,qp2) * dot_product(vert_vec(:,qp1,qp2), omega_cross_u)
+
+ do df = 1, ndf
+ jac_v = matmul(jac(:,:,qp1,qp2), basis(:,df,qp1,qp2))
+ matrix(ik,df,df2) = matrix(ik,df,df2) - dot_product(jac_v, omega_cross_u)
+ end do
+ end do
+ end do
+ end do
+ end do ! end of k loop
+
+end subroutine compute_trad_coriolis_matrix_code
+
+end module compute_trad_coriolis_matrix_kernel_mod
diff --git a/science/gungho/source/kernel/core_dynamics/kessler_kernel_mod.F90 b/science/gungho/source/kernel/core_dynamics/kessler_kernel_mod.F90
new file mode 100644
index 000000000..f7c212011
--- /dev/null
+++ b/science/gungho/source/kernel/core_dynamics/kessler_kernel_mod.F90
@@ -0,0 +1,258 @@
+!-----------------------------------------------------------------------------
+! (c) Crown copyright 2022 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 Performs a simple condensation/evaporation scheme with latent heating
+
+!> @details Given the atmospheric temperature and pressure, this kernel computes
+!! the saturation mixing ratio of water vapour. Any excess vapour is
+!! condensed to cloud liquid, while any cloud liquid in an unsaturated
+!! environment is evaporated to water vapour. The potential temperature
+!! is adjusted to capture the effects of the latent heat release or
+!! absorption associated with this phase change.
+!! Note: this only works with the lowest order spaces
+
+module kessler_kernel_mod
+
+ use argument_mod, only: arg_type, &
+ GH_FIELD, GH_WRITE, GH_READ, &
+ CELL_COLUMN, GH_REAL, GH_SCALAR, &
+ ANY_DISCONTINUOUS_SPACE_3
+ use constants_mod, only: r_def, i_def
+ use driver_water_constants_mod, only: Lv => latent_heat_h2o_condensation
+ use fs_continuity_mod, only: Wtheta
+ use kernel_mod, only: kernel_type
+ use physics_common_mod, only: qsaturation
+ use planet_config_mod, only: recip_epsilon, kappa, cp, Rd, p_zero
+
+ implicit none
+
+ !-----------------------------------------------------------------------------
+ ! Public types
+ !-----------------------------------------------------------------------------
+ !> The type declaration for the kernel. Contains the metadata needed by the Psy layer
+ type, public, extends(kernel_type) :: kessler_kernel_type
+ private
+ type(arg_type) :: meta_args(9) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, WTHETA), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, WTHETA), &
+ arg_type(GH_FIELD*6, GH_REAL, GH_WRITE, WTHETA), &
+ arg_type(GH_FIELD*6, 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, ANY_DISCONTINUOUS_SPACE_3), &
+ arg_type(GH_SCALAR, GH_REAL, GH_READ) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+
+ contains
+ procedure, nopass :: kessler_code
+ end type
+
+ !-----------------------------------------------------------------------------
+ ! Contained functions/subroutines
+ !-----------------------------------------------------------------------------
+ public :: kessler_code
+
+contains
+
+ !> @brief Performs a simple condensation/evaporation scheme with latent heating
+ !! @param[in] nlayers Integer the number of layers
+ !! @param[in,out] theta_inc Potential temperature increment
+ !! @param[in] theta_n Potential temperature in
+ !! @param[in,out] mr_v_inc Water vapour mixing ratio
+ !! @param[in,out] mr_cl_inc Liquid cloud mixing ratio
+ !! @param[in,out] mr_r_inc Rain mixing ratio
+ !! @param[in,out] mr_s_inc Snow mixing ratio
+ !! @param[in,out] mr_g_inc Graupel mixing ratio
+ !! @param[in,out] mr_ci_inc Ice cloud mixing ratio
+ !! @param[in] mr_v_n Water vapour mixing ratio
+ !! @param[in] mr_cl_n Liquid cloud mixing ratio
+ !! @param[in] mr_r_n Rain mixing ratio
+ !! @param[in] mr_s_n Snow mixing ratio
+ !! @param[in] mr_g_n Graupel mixing ratio
+ !! @param[in] mr_ci_n Ice cloud mixing ratio
+ !! @param[in] exner Exner pressure at Wtheta points
+ !! @param[in] ndf_wtheta Number of DoFs per cell for Wtheta
+ !! @param[in] undf_wtheta Universal number of DoFs for wtheta
+ !! @param[in] map_wtheta Integers mapping DoFs to columns for Wtheta
+ subroutine kessler_code(nlayers, theta_inc, theta, &
+ mr_v_inc, mr_cl_inc, mr_r_inc, &
+ mr_s_inc, mr_g_inc, mr_ci_inc, &
+ mr_v, mr_cl, mr_r, &
+ mr_s, mr_g, mr_ci, exner, &
+ rho, height_at_wth, rain2d, &
+ dt, ndf_wtheta, undf_wtheta, map_wtheta, &
+ ndf_w3_2d, undf_w3_2d, map_w3_2d)
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers, ndf_wtheta, undf_wtheta
+ integer(kind=i_def), intent(in) :: ndf_w3_2d, undf_w3_2d
+ integer(kind=i_def), dimension(ndf_wtheta), intent(in) :: map_wtheta
+ integer(kind=i_def), dimension(ndf_w3_2d), intent(in) :: map_w3_2d
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: theta_inc
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: theta
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: exner
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: rho
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: height_at_wth
+ real(kind=r_def), dimension(undf_w3_2d), intent(inout) :: rain2d
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: mr_v_inc, mr_cl_inc
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: mr_r_inc, mr_ci_inc
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: mr_s_inc, mr_g_inc
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: mr_v, mr_cl
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: mr_r, mr_ci
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: mr_s, mr_g
+ real(kind=r_def), intent(in) :: dt
+
+ ! Internal variables
+ integer(kind=i_def) :: df, wt_idx
+ real(kind=r_def) :: theta_np1(nlayers+1), mr_v_np1(nlayers+1)
+ real(kind=r_def) :: mr_cl_np1(nlayers+1), mr_r_np1(nlayers+1)
+ real(kind=r_def) :: temperature, pressure(nlayers+1)
+ real(kind=r_def) :: sed(nlayers+1), rain2d_acc
+ real(kind=r_def) :: mr_sat, dm_v, dm_r, Vr(nlayers+1), Rv
+
+ real(kind=r_def) :: Er, C, dt0
+ real(kind=r_def) :: time_counter
+ real(kind=r_def), parameter :: k1 = 0.001_r_def
+ real(kind=r_def), parameter :: k2 = 2.2_r_def
+ real(kind=r_def), parameter :: a = 0.001_r_def
+ real(kind=r_def), parameter :: rhoqr = 1000.0_r_def
+
+ ! Set max number of iterations for converging scheme
+ wt_idx = map_wtheta(1)
+
+ Rv = Rd * recip_epsilon
+
+ ! Initialise arrays
+ pressure(:) = p_zero * exner(wt_idx : wt_idx+nlayers) ** (1.0_r_def/kappa)
+ theta_np1(:) = theta(wt_idx : wt_idx+nlayers)
+ mr_v_np1(:) = mr_v(wt_idx : wt_idx+nlayers)
+ mr_cl_np1(:) = mr_cl(wt_idx : wt_idx+nlayers)
+ mr_r_np1(:) = mr_r(wt_idx : wt_idx+nlayers)
+
+ ! ======================================================================== !
+ ! Calculate terminal velocity
+ ! ======================================================================== !
+
+ do df = 1, nlayers + 1
+ Vr(df) = 36.34_r_def * SQRT(rho(wt_idx) / rho(wt_idx + df - 1)) &
+ * (MAX(0.0_r_def, mr_r_np1(df)) * 0.001_r_def * rho(wt_idx + df - 1)) ** 0.1346_r_def
+ end do
+
+ ! ======================================================================== !
+ ! Compute maximum time step
+ ! ======================================================================== !
+
+ dt0 = dt
+ do df = 1, nlayers - 1
+ ! NB: Original test for Vr /= 0 numerically unstable
+ if (abs(Vr(df)) > 1.0E-12_r_def) then
+ dt0 = min(dt0, 0.8_r_def*(height_at_wth(wt_idx+df) - height_at_wth(wt_idx+df-1)) / Vr(df))
+ end if
+ end do
+
+ ! ======================================================================== !
+ ! Loop through microphysics iteration
+ ! ======================================================================== !
+
+ ! time counter keeps track of the elapsed time during the subcycling process
+ time_counter = 0.0_r_def
+
+ ! initialize time-weighted accumulated precipitation
+ rain2d_acc = 0.0_r_def
+ ! Subcycle through the Kessler moisture processes,
+ ! time loop ends when the physics time step is reached (within a margin of 1e-5 s)
+ do while ( abs(dt - time_counter) > 1.0E-5_r_def )
+
+ ! Precipitation rate (m_water/s) over the subcycled time step
+ rain2d(map_w3_2d(1)) = rho(wt_idx) * mr_r_np1(wt_idx) * Vr(1) / rhoqr
+
+ ! accumulate the preciptation rate over the subcycled time steps
+ ! (weighted with the subcycled time step), unit is m_water
+ rain2d_acc = rain2d_acc + rain2d(map_w3_2d(1)) * dt0
+
+ ! Mass-weighted sedimentation term using upstream differencing
+ do df = 1, nlayers
+ sed(df) = dt0 * &
+ ((rho(wt_idx+df) * mr_r_np1(df+1) * Vr(df+1)) &
+ - (rho(wt_idx+df-1) * mr_r_np1(df) * Vr(df))) &
+ / (rho(wt_idx+df-1) * (height_at_wth(wt_idx+df) - height_at_wth(wt_idx+df-1)))
+ end do
+ sed(nlayers+1) = -dt0 * mr_r_np1(wt_idx+nlayers) * Vr(nlayers+1) &
+ / (0.5_r_def * (height_at_wth(wt_idx+nlayers)-height_at_wth(wt_idx+nlayers-1)))
+
+ do df = 1, nlayers+1
+
+ ! Autoconversion and collection rates following Klemp and Wilhelmson (1978), Eqs. (2.13a,b)
+ ! the collection process is handled with a semi-implicit time stepping approach
+ dm_r = mr_cl_np1(df) &
+ - (mr_cl_np1(df) - dt0 * MAX(k1 * (mr_cl_np1(df) - a), 0.0_r_def)) &
+ / (1.0_r_def + dt0 * k2 * mr_r_np1(df)**0.875_r_def)
+ mr_cl_np1(df) = MAX(mr_cl_np1(df) - dm_r, 0.0_r_def)
+ mr_r_np1(df) = MAX(mr_r_np1(df) + dm_r + sed(df), 0.0_r_def)
+
+ ! Teten's formula
+ temperature = theta_np1(df) * exner(wt_idx+df-1)
+ ! This function takes pressure in mbar so divide by 100
+ mr_sat = qsaturation(temperature, 0.01_r_def*pressure(df))
+
+ ! Determine difference to saturation amount for vapour
+ dm_v = (mr_v_np1(df) - mr_sat) / &
+ (1.0_r_def + (mr_sat * 4093.0_r_def * Lv / cp) / &
+ (temperature - 36.0_r_def)**2)
+
+ ! Evaporation of rain
+ C = 1.6_r_Def + 124.9_r_def * (rho(wt_idx+df-1)*mr_r_np1(df))**0.2046_r_def
+ Er = (1.0_r_def - mr_v_np1(df) / mr_sat)*C*(rho(wt_idx+df-1)*mr_r_np1(df))**0.525_r_def &
+ / (rho(wt_idx+df-1) * 5.4e5_r_def + 2.55e6_r_def / (pressure(df)*mr_sat))
+ Er = MIN(dt0*Er, MAX(-dm_v - mr_cl_np1(df), mr_r_np1(df)))
+
+ ! Clip to prevent negative cloud forming
+ if (dm_v < 0.0_r_def) then
+ dm_v = max(dm_v, -mr_cl_np1(df))
+ end if
+
+ ! Update fields
+ mr_v_np1(df) = MAX(mr_v_np1(df) - dm_v + Er, 0.0_r_def)
+ mr_cl_np1(df) = mr_cl_np1(df) + dm_v
+ mr_r_np1(df) = MAX(mr_r_np1(df) - Er, 0.0_r_def)
+ theta_np1(df) = theta_np1(df) + Lv / (cp * exner(wt_idx+df-1)) * (dm_v - Er)
+ end do ! Loop over DoFs
+
+ ! Compute the elapsed time
+ time_counter = time_counter + dt0
+
+ ! Recalculate liquid water terminal velocity (m/s)
+ do df = 1, nlayers + 1
+ Vr(df) = 36.34_r_def * SQRT(rho(wt_idx) / rho(wt_idx + df - 1)) &
+ * (MAX(0.0_r_def, mr_r_np1(df)) * 0.001_r_def * rho(wt_idx + df - 1)) ** 0.1346_r_def
+ end do
+
+ ! recompute the time step
+ dt0 = max(dt - time_counter, 0.0_r_def)
+ do df = 1, nlayers - 1
+ ! NB: Original test for Vr /= 0 numerically unstable
+ if (abs(Vr(df)) > 1.0E-12_r_def) then
+ dt0 = min(dt0, 0.8_r_def*(height_at_wth(wt_idx+df) - height_at_wth(wt_idx+df-1)) / Vr(df))
+ end if
+ end do
+
+ end do ! Microphysics time loop
+
+ ! ========================================================================== !
+ ! Calculate increments
+ ! ========================================================================== !
+
+ theta_inc(wt_idx : wt_idx + nlayers) = theta_np1(:) - theta(wt_idx : wt_idx + nlayers)
+ mr_v_inc(wt_idx : wt_idx + nlayers) = mr_v_np1(:) - mr_v(wt_idx : wt_idx + nlayers)
+ mr_cl_inc(wt_idx : wt_idx + nlayers) = mr_cl_np1(:) - mr_cl(wt_idx : wt_idx + nlayers)
+ mr_r_inc(wt_idx : wt_idx + nlayers) = mr_r_np1(:) - mr_r(wt_idx : wt_idx + nlayers)
+
+ end subroutine kessler_code
+
+end module kessler_kernel_mod
diff --git a/science/gungho/source/kernel/diagnostics/dbz_kernel_mod.F90 b/science/gungho/source/kernel/diagnostics/dbz_kernel_mod.F90
new file mode 100644
index 000000000..7f66028bb
--- /dev/null
+++ b/science/gungho/source/kernel/diagnostics/dbz_kernel_mod.F90
@@ -0,0 +1,122 @@
+!-----------------------------------------------------------------------------
+! (C) Crown copyright 2021 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 Computes the relative humidity from prognostic variables.
+!>
+!> @details Computes the relative humidity field, at Wtheta points from the
+!> potential temperature, the Exner pressure and the mixing ratio of
+!> water vapour.
+!>
+module dbz_kernel_mod
+
+ use argument_mod, only : arg_type, GH_SCALAR, &
+ GH_FIELD, GH_REAL, &
+ GH_WRITE, GH_READ, &
+ CELL_COLUMN
+ use constants_mod, only : r_def, i_def
+ use fs_continuity_mod, only : Wtheta
+ use kernel_mod, only : kernel_type
+ use physics_common_mod, only : qsaturation
+
+ implicit none
+
+ private
+
+ !---------------------------------------------------------------------------
+ ! Public types
+ !---------------------------------------------------------------------------
+ !> The type declaration for the kernel. Contains the metadata needed by the
+ !> Psy layer.
+ !>
+ type, public, extends(kernel_type) :: dbz_kernel_type
+ private
+ type(arg_type) :: meta_args(10) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, 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, Wtheta), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), &
+ arg_type(GH_SCALAR, GH_REAL, GH_READ), &
+ arg_type(GH_SCALAR, GH_REAL, GH_READ), &
+ arg_type(GH_SCALAR, GH_REAL, GH_READ) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+ contains
+ procedure, nopass :: dbz_code
+ end type
+
+!-----------------------------------------------------------------------------
+! Contained functions/subroutines
+!-----------------------------------------------------------------------------
+public :: dbz_code
+
+contains
+
+subroutine dbz_code( nlayers, &
+ dbz, &
+ exner_at_wt, &
+ theta, &
+ rho_at_wt, &
+ mr_v, &
+ mr_cl, &
+ mr_r, &
+ recip_epsilon, &
+ kappa, &
+ p_zero, &
+ ndf_wt, &
+ undf_wt, &
+ map_wt )
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers
+ integer(kind=i_def), intent(in) :: undf_wt, ndf_wt
+
+ integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt
+
+ real(kind=r_def), dimension(undf_wt), intent(inout) :: dbz
+ real(kind=r_def), dimension(undf_wt), intent(in) :: theta
+ real(kind=r_def), dimension(undf_wt), intent(in) :: exner_at_wt
+ real(kind=r_def), dimension(undf_wt), intent(in) :: rho_at_wt
+ real(kind=r_def), dimension(undf_wt), intent(in) :: mr_v, mr_cl, mr_r
+ real(kind=r_def), intent(in) :: recip_epsilon
+ real(kind=r_def), intent(in) :: kappa
+ real(kind=r_def), intent(in) :: p_zero
+
+ ! Internal variables
+ integer(kind=i_def) :: df, min_col_index, max_col_index
+ real(kind=r_def) :: temperature, pressure, mr_sat, rho_sqrt, rel_hum
+ real(kind=r_def) :: rho_wet, vel_r, z_e
+
+ real(kind=r_def), parameter :: vconr = 2503.23638966667_r_def
+ real(kind=r_def), parameter :: normr = 25132741228.7183_r_def
+
+ ! Find minimum and maximum DoF numberings for this column
+ min_col_index = minval(map_wt)
+ max_col_index = maxval(map_wt) + nlayers - 1
+
+ ! Directly loop over all DoFs in the column
+ do df = min_col_index, max_col_index
+
+ temperature = theta(df) * exner_at_wt(df)
+ pressure = p_zero * exner_at_wt(df) ** (1.0_r_def/kappa)
+ ! Pressure for saturation curve is needed in mbar
+ mr_sat = qsaturation(temperature, 0.01_r_def*pressure)
+ rel_hum = mr_v(df) / mr_sat * ( 1.0_r_def + mr_sat * recip_epsilon ) / &
+ ( 1.0_r_def + mr_v(df) * recip_epsilon )
+ rho_sqrt = SQRT(MIN(10.0_r_def, rho_at_wt(min_col_index) / rho_at_wt(df)))
+ rho_wet = rho_at_wt(df) * MAX(1.0e-12_r_def, mr_r(df) + DIM(mr_cl(df), 1.0e-3_r_def))
+ vel_r = MAX(1.0e-3_r_def, vconr * rho_sqrt * EXP(0.2_r_def*LOG(rho_wet / normr)))
+ z_e = 200.0_r_def * EXP(1.6_r_def * LOG(3.6e6_r_def * (rho_wet / 1.0e3_r_def) * vel_r))
+ dbz(df) = 10.0_r_def * LOG10(MAX(z_e, 0.01_r_def))
+ end do
+
+end subroutine dbz_code
+
+end module dbz_kernel_mod
diff --git a/science/gungho/source/kernel/diagnostics/theta_v_kernel_mod.F90 b/science/gungho/source/kernel/diagnostics/theta_v_kernel_mod.F90
new file mode 100644
index 000000000..3bede6e4a
--- /dev/null
+++ b/science/gungho/source/kernel/diagnostics/theta_v_kernel_mod.F90
@@ -0,0 +1,85 @@
+!-----------------------------------------------------------------------------
+! (C) Crown copyright 2021 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 Computes theta_v, the virtual potential temperature.
+module theta_v_kernel_mod
+
+ use argument_mod, only : arg_type, GH_SCALAR, &
+ GH_FIELD, GH_REAL, &
+ GH_WRITE, GH_READ, &
+ 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
+ !---------------------------------------------------------------------------
+ !> The type declaration for the kernel. Contains the metadata needed by the
+ !> Psy layer.
+ !>
+ type, public, extends(kernel_type) :: theta_v_kernel_type
+ private
+ type(arg_type) :: meta_args(4) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, Wtheta), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, 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 :: theta_v_code
+ end type
+
+!-----------------------------------------------------------------------------
+! Contained functions/subroutines
+!-----------------------------------------------------------------------------
+public :: theta_v_code
+
+contains
+
+!> @brief Computes theta_v, the virtual potential temperature.
+subroutine theta_v_code( nlayers, &
+ theta_v, &
+ theta, &
+ mr_v, &
+ recip_epsilon, &
+ ndf_wt, &
+ undf_wt, &
+ map_wt )
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers
+ integer(kind=i_def), intent(in) :: undf_wt, ndf_wt
+
+ integer(kind=i_def), dimension(ndf_wt), intent(in) :: map_wt
+
+ real(kind=r_def), dimension(undf_wt), intent(inout) :: theta_v
+ real(kind=r_def), dimension(undf_wt), intent(in) :: theta
+ real(kind=r_def), dimension(undf_wt), intent(in) :: mr_v
+ real(kind=r_def), intent(in) :: recip_epsilon
+
+ ! Internal variables
+ integer(kind=i_def) :: df, min_col_index, max_col_index
+
+ ! Find minimum and maximum DoF numberings for this column
+ min_col_index = minval(map_wt)
+ max_col_index = maxval(map_wt) + nlayers - 1
+
+ ! Directly loop over all DoFs in the column
+ do df = min_col_index, max_col_index
+ theta_v(df) = theta(df) * (1.0_r_def + mr_v(df) * recip_epsilon) / (1.0_r_def + mr_v(df))
+ end do
+
+end subroutine theta_v_code
+
+end module theta_v_kernel_mod
diff --git a/science/gungho/source/kernel/external_forcing/rayleigh_friction_kernel_mod.F90 b/science/gungho/source/kernel/external_forcing/rayleigh_friction_kernel_mod.F90
new file mode 100644
index 000000000..04ebb524b
--- /dev/null
+++ b/science/gungho/source/kernel/external_forcing/rayleigh_friction_kernel_mod.F90
@@ -0,0 +1,133 @@
+!-----------------------------------------------------------------------------
+! (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 Rayleigh Friction
+!>
+module rayleigh_friction_kernel_mod
+
+ use argument_mod, only: arg_type, &
+ GH_FIELD, GH_REAL, &
+ GH_READ, GH_INC, &
+ GH_SCALAR, &
+ CELL_COLUMN
+ use constants_mod, only: r_def, i_def, PI
+ use fs_continuity_mod, only: W2, Wtheta, W3
+ use kernel_mod, only: kernel_type
+
+ implicit none
+
+ private
+
+ !---------------------------------------------------------------------------
+ ! Public types
+ !---------------------------------------------------------------------------
+ !> The type declaration for the kernel. Contains the metadata needed by the
+ !> Psy layer.
+ !>
+ type, public, extends(kernel_type) :: rayleigh_friction_kernel_type
+ private
+ type(arg_type) :: meta_args(11) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_INC, W2), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, W2), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, W3), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), &
+ 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) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+ contains
+ procedure, nopass :: rayleigh_friction_code
+ end type
+
+ !---------------------------------------------------------------------------
+ ! Contained functions/subroutines
+ !---------------------------------------------------------------------------
+ public :: rayleigh_friction_code
+
+contains
+
+!> @brief The subroutine which is called directly by the psy layer
+!! @param[in] nlayers Integer the number of layers
+!! @param[in,out] du u increment data
+!! @param[in] u u data
+!! @param[in] w2_rmultiplicity Reciprocal of multiplicity for W2
+!! @param[in] exner The exner pressure
+!! @param[in] exner_in_wth The exner pressure in Wtheta
+!! @param[in] kappa Ratio of Rd and cp
+!! @param[in] dt The model timestep length
+!! @param[in] ndf_w2 The number of degrees of freedom per cell for W2
+!! @param[in] undf_w2 The number of unique degrees of freedom for W2
+!! @param[in] map_w2 Integer array holding the dofmap for the cell at the
+!> base of the column for W2
+!! @param[in] ndf_w3 The number of degrees of freedom per cell for W3
+!! @param[in] undf_w3 The number of unique degrees of freedom for W3
+!! @param[in] map_w3 Integer array holding the dofmap for the cell at the
+!> base of the column for W3
+!! @param[in] ndf_wth The number of degrees of freedom per cell for Wtheta
+!! @param[in] undf_wth The number of unique degrees of freedom for Wtheta
+!! @param[in] map_wth Integer array holding the dofmap for the cell at the
+!> base of the column for Wtheta
+subroutine rayleigh_friction_code( nlayers, &
+ du, u, u_ref, &
+ w2_rmultiplicity, &
+ exner, exner_in_wth, &
+ kappa, p_zero, p_cutoff, tau, dt, &
+ ndf_w2, undf_w2, map_w2, &
+ ndf_w3, undf_w3, map_w3, &
+ ndf_wth, undf_wth, map_wth &
+ )
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers
+
+ integer(kind=i_def), intent(in) :: ndf_wth, undf_wth
+ integer(kind=i_def), intent(in) :: ndf_w2, undf_w2
+ integer(kind=i_def), intent(in) :: ndf_w3, undf_w3
+
+ real(kind=r_def), dimension(undf_w2), intent(inout) :: du
+ real(kind=r_def), dimension(undf_w2), intent(in) :: u, u_ref, w2_rmultiplicity
+ real(kind=r_def), dimension(undf_wth), intent(in) :: exner_in_wth
+ real(kind=r_def), dimension(undf_w3), intent(in) :: exner
+ real(kind=r_def), intent(in) :: kappa, p_zero, p_cutoff, tau
+ real(kind=r_def), intent(in) :: dt
+
+ integer(kind=i_def), dimension(ndf_w2), intent(in) :: map_w2
+ integer(kind=i_def), dimension(ndf_wth), intent(in) :: map_wth
+ integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3
+
+ ! Internal variables
+ integer(kind=i_def) :: k, df, w2_idx
+ real(kind=r_def) :: kR(nlayers)
+ real(kind=r_def) :: exner_top, p_top
+ real(kind=r_def) :: p(nlayers), p_switch(nlayers)
+
+ exner_top = exner_in_wth(map_wth(1)+nlayers)
+ p_top = exner_top ** (1.0_r_def / kappa) * p_zero
+ p(:) = exner(map_w3(1) : map_w3(1)+nlayers-1) ** (1.0_r_def / kappa) * p_zero
+ ! Should only kick in when pressure is lower than p_cutoff
+ p_switch(:) = 0.5_r_def + 0.5_r_def*SIGN(1.0_r_def, p_cutoff - p(:))
+ kR(:) = p_switch / tau * (SIN(0.5_r_def*PI*(LOG(p_cutoff/p(:)) / LOG(p_cutoff/p_top))))**2
+
+ do df = 1, 4
+ do k = 1, nlayers
+ w2_idx = map_w2(df) + k - 1
+ ! Backward Euler discretisation
+ du(w2_idx) = du(w2_idx) + w2_rmultiplicity(w2_idx) * ( &
+ (1.0_r_def / (1.0_r_def + kR(k) * dt) - 1.0_r_def) * u(w2_idx) &
+ + kR(k) * dt / (1.0_r_def + kR(k) * dt) * u_ref(w2_idx) &
+ )
+ end do
+ end do
+
+end subroutine rayleigh_friction_code
+
+end module rayleigh_friction_kernel_mod
diff --git a/science/gungho/source/kernel/initialisation/hydrostatic_exner_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hydrostatic_exner_kernel_mod.F90
index 2b5206f60..2934e0172 100644
--- a/science/gungho/source/kernel/initialisation/hydrostatic_exner_kernel_mod.F90
+++ b/science/gungho/source/kernel/initialisation/hydrostatic_exner_kernel_mod.F90
@@ -147,10 +147,11 @@ subroutine hydrostatic_exner_code(nlayers, exner, theta, &
real(kind=r_def) :: theta_moist
real(kind=r_def) :: g_local(0:nlayers-1)
real(kind=r_def) :: coords(3), xyz(3)
- real(kind=r_def) :: exner_start
+ real(kind=r_def) :: exner_start, surface_height
real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
ipanel = int(panel_id(map_pid(1)), i_def)
+ surface_height = height_wt(map_wt(1))
if (init_exner_bt) then
layers_offset = 0
@@ -177,7 +178,7 @@ subroutine hydrostatic_exner_code(nlayers, exner, theta, &
call chi2xyz(coords(1), coords(2), coords(3), ipanel, xyz(1), xyz(2), xyz(3))
! Exner at the model surface or top
- exner_start = analytic_pressure( xyz, test, 0.0_r_def)
+ exner_start = analytic_pressure( xyz, test, 0.0_r_def, surface_height)
! Set local value of gravity at each vertical level
if (shallow) then
diff --git a/science/gungho/source/kernel/initialisation/hydrostatic_exner_squall_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/hydrostatic_exner_squall_kernel_mod.F90
new file mode 100644
index 000000000..df9fd50da
--- /dev/null
+++ b/science/gungho/source/kernel/initialisation/hydrostatic_exner_squall_kernel_mod.F90
@@ -0,0 +1,232 @@
+!-----------------------------------------------------------------------------
+! Copyright (c) 2018, Met Office, on behalf of HMSO and Queen's Printer
+! For further details please refer to the file LICENCE.original which you
+! should have received as part of this distribution.
+!-----------------------------------------------------------------------------
+!
+!-------------------------------------------------------------------------------
+
+!> @brief Computes Exner pressure distribution via hydrostatic balance equation (lowest order only)
+
+module hydrostatic_exner_squall_kernel_mod
+
+use argument_mod, only : arg_type, func_type, &
+ GH_FIELD, GH_REAL, &
+ GH_SCALAR, &
+ GH_READ, GH_WRITE, &
+ ANY_SPACE_9, GH_BASIS, &
+ ANY_DISCONTINUOUS_SPACE_3, &
+ CELL_COLUMN, GH_EVALUATOR
+use constants_mod, only : r_def, i_def
+use idealised_config_mod, only : test
+use kernel_mod, only : kernel_type
+use fs_continuity_mod, only : Wtheta, W3
+use formulation_config_mod, only : init_exner_bt, shallow
+
+implicit none
+
+private
+
+!-------------------------------------------------------------------------------
+! Public types
+!-------------------------------------------------------------------------------
+!> The type declaration for the kernel. Contains the metadata needed by the Psy layer
+type, public, extends(kernel_type) :: hydrostatic_exner_squall_kernel_type
+ private
+ type(arg_type) :: meta_args(12) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, W3), &
+ 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, Wtheta), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, W3), &
+ arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &
+ 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) &
+ /)
+
+ type(func_type) :: meta_funcs(1) = (/ &
+ func_type(ANY_SPACE_9, GH_BASIS) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+ integer :: gh_shape = GH_EVALUATOR
+ integer :: gh_evaluator_targets(1) = (/ Wtheta /)
+contains
+ procedure, nopass :: hydrostatic_exner_squall_code
+end type
+
+!-------------------------------------------------------------------------------
+! Contained functions/subroutines
+!-------------------------------------------------------------------------------
+public :: hydrostatic_exner_squall_code
+
+contains
+
+!> @brief Computes hydrostatic Exner function
+!> @param[in] nlayers Number of layers
+!> @param[in,out] exner Exner pressure field
+!> @param[in] theta Potential temperature field
+!> @param[in] moist_dyn_gas Moist dynamics factor in gas law
+!> @param[in] moist_dyn_tot Moist dynamics total mass factor
+!> @param[in] moist_dyn_fac Moist dynamics water factor
+!> @param[in] height_wt Height coordinate in Wtheta
+!> @param[in] height_w3 Height coordinate in W3
+!> @param[in] chi_1 First component of the chi coordinate field
+!> @param[in] chi_2 Second component of the chi coordinate field
+!> @param[in] chi_3 Third component of the chi coordinate field
+!> @param[in] panel_id A field giving the ID for mesh panels
+!> @param[in] gravity The planet gravity
+!> @param[in] p_zero Reference surface pressure
+!> @param[in] rd Gas constant for dry air
+!> @param[in] cp Specific heat of dry air at constant pressure
+!> @param[in] radius Radius of planet at surface
+!> @param[in] ndf_w3 Number of degrees of freedom per cell for W3
+!> @param[in] undf_w3 Number of unique degrees of freedom for W3
+!> @param[in] map_w3 Dofmap for the cell at the base of the column for W3
+!> @param[in] ndf_wt Number of degrees of freedom per cell for Wtheta
+!> @param[in] undf_wt Number of unique degrees of freedom for Wtheta
+!> @param[in] map_wt Dofmap for the cell at the base of the column for Wtheta
+!> @param[in] ndf_chi Number of degrees of freedom per cell for Wchi
+!> @param[in] undf_chi Number of unique degrees of freedom for Wchi
+!> @param[in] map_chi Dofmap for the cell at the base of the column for Wchi
+!> @param[in] basis_chi_on_wt Basis functions for Wchi evaluated at
+!! Wtheta nodal points
+!> @param[in] ndf_pid Number of degrees of freedom per cell for panel_id
+!> @param[in] undf_pid Number of unique degrees of freedom for panel_id
+!> @param[in] map_pid Dofmap for the cell at the base of the column for panel_id
+subroutine hydrostatic_exner_squall_code(nlayers, exner, theta, &
+ moist_dyn_gas, moist_dyn_tot, &
+ moist_dyn_fac, &
+ exner_surf, &
+ height_wt, height_w3, &
+ chi_1, chi_2, chi_3, panel_id, &
+ gravity, rd, cp, &
+ radius, &
+ ndf_w3, undf_w3, map_w3, &
+ ndf_wt, undf_wt, map_wt, &
+ ndf_chi, undf_chi, map_chi, &
+ basis_chi_on_wt, &
+ ndf_pid, undf_pid, map_pid )
+
+ use analytic_pressure_profiles_mod, only : analytic_pressure
+ use sci_chi_transform_mod, only : chi2xyz
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers
+ integer(kind=i_def), intent(in) :: ndf_w3, undf_w3
+ integer(kind=i_def), intent(in) :: ndf_wt, undf_wt
+ integer(kind=i_def), intent(in) :: ndf_chi, undf_chi
+ integer(kind=i_def), intent(in) :: ndf_pid, undf_pid
+ 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), dimension(ndf_chi), intent(in) :: map_chi
+ integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid
+ real(kind=r_def), dimension(undf_w3), intent(inout) :: exner
+ real(kind=r_def), dimension(undf_wt), intent(in) :: theta
+ real(kind=r_Def), dimension(undf_wt), intent(in) :: exner_surf
+ 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) :: height_wt
+ real(kind=r_def), dimension(undf_w3), intent(in) :: height_w3
+ real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3
+ real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id
+ real(kind=r_def), dimension(1,ndf_chi,ndf_wt), intent(in) :: basis_chi_on_wt
+ real(kind=r_def), intent(in) :: gravity
+ real(kind=r_def), intent(in) :: rd
+ real(kind=r_def), intent(in) :: cp
+ real(kind=r_def), intent(in) :: radius
+
+ ! Internal variables
+ integer(kind=i_def) :: k, dfc, layers_offset, wt_dof, ipanel
+ real(kind=r_def) :: dz
+ real(kind=r_def) :: theta_moist
+ real(kind=r_def) :: g_local(0:nlayers-1)
+ real(kind=r_def) :: coords(3), xyz(3)
+ real(kind=r_def) :: exner_start, surface_height
+ real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
+
+ ipanel = int(panel_id(map_pid(1)), i_def)
+ surface_height = height_wt(map_wt(1))
+
+ if (init_exner_bt) then
+ layers_offset = 0
+ wt_dof = 1
+ exner_start = exner_surf(map_wt(1))
+ else
+ layers_offset = nlayers - 1
+ wt_dof = 2
+ exner_start = exner_surf(map_wt(1)+nlayers)
+ end if
+
+ do dfc = 1, ndf_chi
+ chi_1_e(dfc) = chi_1( map_chi(dfc) + layers_offset )
+ chi_2_e(dfc) = chi_2( map_chi(dfc) + layers_offset )
+ chi_3_e(dfc) = chi_3( map_chi(dfc) + layers_offset )
+ end do
+
+ ! Horizontal coordinates of cell bottom or top
+ coords(:) = 0.0_r_def
+ do dfc = 1, ndf_chi
+ coords(1) = coords(1) + chi_1_e(dfc)*basis_chi_on_wt(1,dfc,wt_dof)
+ coords(2) = coords(2) + chi_2_e(dfc)*basis_chi_on_wt(1,dfc,wt_dof)
+ coords(3) = coords(3) + chi_3_e(dfc)*basis_chi_on_wt(1,dfc,wt_dof)
+ end do
+
+ call chi2xyz(coords(1), coords(2), coords(3), ipanel, xyz(1), xyz(2), xyz(3))
+
+ ! Set local value of gravity at each vertical level
+ if (shallow) then
+ do k = 0, nlayers-1
+ g_local(k) = gravity
+ end do
+ else
+ do k = 0, nlayers-1
+ g_local(k) = gravity * (radius / (radius + height_wt(map_wt(1)+k))) ** 2
+ end do
+ end if
+
+ if (init_exner_bt) then
+
+ ! Bottom-up initialization
+ ! Exner at the bottom level
+ dz = height_w3(map_w3(1))-height_wt(map_wt(1))
+ theta_moist = moist_dyn_gas(map_wt(1)) * theta(map_wt(1)) / &
+ moist_dyn_tot(map_wt(1))
+ exner(map_w3(1)) = exner_start - g_local(0) * dz / (cp * theta_moist)
+
+ ! Exner on other levels
+ do k = 1, nlayers-1
+ dz = height_w3(map_w3(1)+k)-height_w3(map_w3(1)+k-1)
+ theta_moist = moist_dyn_gas(map_wt(1)+k) * theta(map_wt(1)+k) / &
+ moist_dyn_tot(map_wt(1)+k)
+ exner(map_w3(1)+k) = exner(map_w3(1)+k-1) - g_local(k) * dz / (cp * theta_moist)
+ end do
+
+ else
+
+ ! Top-down initialization
+ ! Exner at the top level
+ dz = height_wt(map_wt(1)+nlayers) - height_w3(map_w3(1)+nlayers-1)
+ theta_moist = moist_dyn_gas(map_wt(1)+nlayers) * theta(map_wt(1)+nlayers) / &
+ moist_dyn_tot(map_wt(1)+nlayers)
+ exner(map_w3(1)+nlayers-1) = exner_start + g_local(nlayers-1) * dz / (cp * theta_moist)
+
+ ! Exner on other levels
+ do k = nlayers-1, 1, -1
+ dz = height_w3(map_w3(1)+k)-height_w3(map_w3(1)+k-1)
+ theta_moist = moist_dyn_gas(map_wt(1)+k) * theta(map_wt(1)+k) / &
+ moist_dyn_tot(map_wt(1)+k)
+ exner(map_w3(1)+k-1) = exner(map_w3(1)+k) + g_local(k-1) * dz / (cp * theta_moist)
+ end do
+
+ end if
+
+end subroutine hydrostatic_exner_squall_code
+
+end module hydrostatic_exner_squall_kernel_mod
diff --git a/science/gungho/source/kernel/initialisation/initial_exner_sample_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/initial_exner_sample_kernel_mod.F90
index 08f6e076c..3e4612353 100644
--- a/science/gungho/source/kernel/initialisation/initial_exner_sample_kernel_mod.F90
+++ b/science/gungho/source/kernel/initialisation/initial_exner_sample_kernel_mod.F90
@@ -15,7 +15,7 @@ module initial_exner_sample_kernel_mod
GH_BASIS, CELL_COLUMN, &
GH_EVALUATOR
use constants_mod, only : r_def, i_def
- use fs_continuity_mod, only : W3
+ use fs_continuity_mod, only : W3, Wtheta
use idealised_config_mod, only : test
use kernel_mod, only : kernel_type
@@ -31,8 +31,9 @@ module initial_exner_sample_kernel_mod
!>
type, public, extends(kernel_type) :: initial_exner_sample_kernel_type
private
- type(arg_type) :: meta_args(4) = (/ &
+ type(arg_type) :: meta_args(5) = (/ &
arg_type(GH_FIELD, GH_REAL, GH_WRITE, W3), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), &
arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &
arg_type(GH_SCALAR, GH_REAL, GH_READ) &
@@ -75,10 +76,12 @@ module initial_exner_sample_kernel_mod
!! for panel_id
subroutine initial_exner_sample_code(nlayers, &
exner, &
+ height_wth, &
chi_1, chi_2, chi_3, &
panel_id, &
current_time, &
ndf_w3, undf_w3, map_w3, &
+ ndf_wth, undf_wth, map_wth, &
ndf_chi, undf_chi, map_chi, &
basis_chi_on_w3, &
ndf_pid, undf_pid, map_pid )
@@ -90,16 +93,18 @@ subroutine initial_exner_sample_code(nlayers, &
! Arguments
integer(kind=i_def), intent(in) :: nlayers
- integer(kind=i_def), intent(in) :: ndf_w3
+ integer(kind=i_def), intent(in) :: ndf_w3, ndf_wth
integer(kind=i_def), intent(in) :: ndf_chi
integer(kind=i_def), intent(in) :: ndf_pid
- integer(kind=i_def), intent(in) :: undf_w3
+ integer(kind=i_def), intent(in) :: undf_w3, undf_wth
integer(kind=i_def), intent(in) :: undf_chi
integer(kind=i_def), intent(in) :: undf_pid
integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3
+ integer(kind=i_def), dimension(ndf_wth), intent(in) :: map_wth
integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi
integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_pid
real(kind=r_def), dimension(undf_w3), intent(inout) :: exner
+ real(kind=r_def), dimension(undf_wth), intent(in) :: height_wth
real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1
real(kind=r_def), dimension(undf_chi), intent(in) :: chi_2
real(kind=r_def), dimension(undf_chi), intent(in) :: chi_3
@@ -109,11 +114,13 @@ subroutine initial_exner_sample_code(nlayers, &
! Internal variables
real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
- real(kind=r_def) :: coords(3), xyz(3)
+ real(kind=r_def) :: coords(3), xyz(3), surface_height
integer(kind=i_def) :: df1, df, k, ipanel
ipanel = int(panel_id(map_pid(1)), i_def)
+ surface_height = height_wth(map_wth(1))
+
! Compute the RHS & LHS integrated over one cell and solve
do k = 0, nlayers-1
do df1 = 1, ndf_chi
@@ -133,7 +140,7 @@ subroutine initial_exner_sample_code(nlayers, &
call chi2xyz(coords(1), coords(2), coords(3), &
ipanel, xyz(1), xyz(2), xyz(3))
- exner(map_w3(df) + k) = analytic_pressure(xyz, test, current_time)
+ exner(map_w3(df) + k) = analytic_pressure(xyz, test, current_time, surface_height)
end do
end do
diff --git a/science/gungho/source/kernel/initialisation/initial_rel_humidity_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/initial_rel_humidity_kernel_mod.F90
new file mode 100644
index 000000000..c346db308
--- /dev/null
+++ b/science/gungho/source/kernel/initialisation/initial_rel_humidity_kernel_mod.F90
@@ -0,0 +1,116 @@
+!-----------------------------------------------------------------------------
+! Copyright (c) 2017, Met Office, on behalf of HMSO and Queen's Printer
+! For further details please refer to the file LICENCE.original which you
+! should have received as part of this distribution.
+!-----------------------------------------------------------------------------
+!
+!-------------------------------------------------------------------------------
+
+!> @brief Kernel computes the initial mixing ratio fields
+
+!> @details The kernel computes initial mixing ratio fields for mr in the same
+!> space as that of theta
+
+module initial_rel_humidity_kernel_mod
+
+ use argument_mod, only: arg_type, func_type, &
+ GH_FIELD, GH_REAL, &
+ GH_SCALAR, GH_BASIS, &
+ GH_WRITE, GH_READ, &
+ ANY_SPACE_9, &
+ ANY_DISCONTINUOUS_SPACE_3, &
+ GH_EVALUATOR, CELL_COLUMN
+ use fs_continuity_mod, only: W3, Wtheta
+ use constants_mod, only: r_def, i_def
+ use kernel_mod, only: kernel_type
+ use idealised_config_mod, only: test
+
+ implicit none
+ private
+
+ !-------------------------------------------------------------------------------
+ ! Public types
+ !-------------------------------------------------------------------------------
+ !> The type declaration for the kernel. Contains the metadata needed by the Psy layer
+ type, public, extends(kernel_type) :: initial_rel_humidity_kernel_type
+ private
+ type(arg_type) :: meta_args(3) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, Wtheta), &
+ arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) &
+ /)
+ type(func_type) :: meta_funcs(1) = (/ &
+ func_type(ANY_SPACE_9, GH_BASIS) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+ integer :: gh_shape = GH_EVALUATOR
+ contains
+ procedure, nopass :: initial_rel_humidity_code
+ end type
+
+ !-------------------------------------------------------------------------------
+ ! Contained functions/subroutines
+ !-------------------------------------------------------------------------------
+ public :: initial_rel_humidity_code
+contains
+
+ subroutine initial_rel_humidity_code(nlayers, &
+ rel_hum, &
+ chi_1, chi_2, chi_3, &
+ panel_id, &
+ ndf_wtheta, undf_wtheta, map_wtheta, &
+ ndf_chi, undf_chi, map_chi, chi_basis, &
+ ndf_pid, undf_pid, map_pid )
+
+ use analytic_moisture_profiles_mod, only : analytic_rel_humidity
+ use sci_chi_transform_mod, only : chi2xyz
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers, ndf_wtheta, undf_wtheta
+ integer(kind=i_def), intent(in) :: ndf_chi, undf_chi
+ integer(kind=i_def), intent(in) :: ndf_pid, undf_pid
+ integer(kind=i_def), dimension(ndf_wtheta), intent(in) :: map_wtheta
+ integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi
+ integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid
+
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: rel_hum
+ real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3
+ real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id
+
+ real(kind=r_def), dimension(1,ndf_chi,ndf_wtheta), intent(in) :: chi_basis
+
+ ! Internal variables
+ integer(kind=i_def) :: k, df, dfc, ipanel
+ real(kind=r_def) :: coords(3), xyz(3)
+ real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
+
+ ipanel = int(panel_id(map_pid(1)), i_def)
+
+
+ do k = 0, nlayers-1
+ do dfc = 1, ndf_chi
+ chi_1_e(dfc) = chi_1( map_chi(dfc) + k)
+ chi_2_e(dfc) = chi_2( map_chi(dfc) + k)
+ chi_3_e(dfc) = chi_3( map_chi(dfc) + k)
+ end do
+
+ do df = 1, ndf_wtheta
+ coords(:) = 0.0_r_def
+ do dfc = 1, ndf_chi
+ coords(1) = coords(1) + chi_1_e(dfc)*chi_basis(1,dfc,df)
+ coords(2) = coords(2) + chi_2_e(dfc)*chi_basis(1,dfc,df)
+ coords(3) = coords(3) + chi_3_e(dfc)*chi_basis(1,dfc,df)
+ end do
+
+ call chi2xyz(coords(1), coords(2), coords(3), &
+ ipanel, xyz(1), xyz(2), xyz(3))
+ rel_hum(map_wtheta(df) + k) = analytic_rel_humidity(xyz, test)
+
+ end do
+ end do
+
+ end subroutine initial_rel_humidity_code
+
+end module initial_rel_humidity_kernel_mod
diff --git a/science/gungho/source/kernel/initialisation/initial_theta_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/initial_theta_kernel_mod.F90
index 55dfe8df7..b7d4aa267 100644
--- a/science/gungho/source/kernel/initialisation/initial_theta_kernel_mod.F90
+++ b/science/gungho/source/kernel/initialisation/initial_theta_kernel_mod.F90
@@ -34,8 +34,9 @@ module initial_theta_kernel_mod
!> The type declaration for the kernel. Contains the metadata needed by the Psy layer
type, public, extends(kernel_type) :: initial_theta_kernel_type
private
- type(arg_type) :: meta_args(3) = (/ &
+ type(arg_type) :: meta_args(4) = (/ &
arg_type(GH_FIELD, GH_REAL, GH_WRITE, Wtheta), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), &
arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) &
/)
@@ -73,7 +74,7 @@ module initial_theta_kernel_mod
!! @param[in] undf_pid Number of unique degrees of freedom for panel_id
!! @param[in] map_pid Dofmap for the cell at the base of the column for panel_id
subroutine initial_theta_code(nlayers, &
- theta, &
+ theta, height_wth, &
chi_1, chi_2, chi_3, &
panel_id, &
ndf_wtheta, undf_wtheta, map_wtheta, &
@@ -94,15 +95,17 @@ subroutine initial_theta_code(nlayers, &
integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid
real(kind=r_def), dimension(undf_wtheta), intent(inout) :: theta
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: height_wth
real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3
real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id
real(kind=r_def), dimension(1,ndf_chi,ndf_wtheta), intent(in) :: chi_basis
! Internal variables
real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
- real(kind=r_def) :: coords(3), xyz(3)
+ real(kind=r_def) :: coords(3), xyz(3), surface_height
integer(kind=i_def) :: df, dfc, k, ipanel
+ surface_height = height_wth(map_wtheta(1))
ipanel = int(panel_id(map_pid(1)), i_def)
! Compute the pointwise theta profile
@@ -124,7 +127,7 @@ subroutine initial_theta_code(nlayers, &
call chi2xyz(coords(1), coords(2), coords(3), &
ipanel, xyz(1), xyz(2), xyz(3))
- theta(map_wtheta(df) + k) = analytic_temperature(xyz, test)
+ theta(map_wtheta(df) + k) = analytic_temperature(xyz, test, surface_height)
end do
end do
diff --git a/science/gungho/source/kernel/initialisation/initial_theta_pert_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/initial_theta_pert_kernel_mod.F90
new file mode 100644
index 000000000..a2b4e4343
--- /dev/null
+++ b/science/gungho/source/kernel/initialisation/initial_theta_pert_kernel_mod.F90
@@ -0,0 +1,130 @@
+!-----------------------------------------------------------------------------
+! Copyright (c) 2017, Met Office, on behalf of HMSO and Queen's Printer
+! For further details please refer to the file LICENCE.original which you
+! should have received as part of this distribution.
+!-----------------------------------------------------------------------------
+!
+!-------------------------------------------------------------------------------
+
+!> @brief Computes the initial theta perturbation field
+module initial_theta_pert_kernel_mod
+
+ use argument_mod, only: arg_type, func_type, &
+ GH_FIELD, GH_REAL, &
+ GH_WRITE, GH_READ, &
+ ANY_SPACE_9, GH_BASIS, &
+ ANY_DISCONTINUOUS_SPACE_3, &
+ CELL_COLUMN, GH_EVALUATOR
+ use constants_mod, only: r_def, i_def
+ use fs_continuity_mod, only: Wtheta
+ use kernel_mod, only: kernel_type
+ use idealised_config_mod, only: test
+
+ implicit none
+
+ private
+
+ !-------------------------------------------------------------------------------
+ ! Public types
+ !-------------------------------------------------------------------------------
+ !> The type declaration for the kernel. Contains the metadata needed by the Psy layer
+ type, public, extends(kernel_type) :: initial_theta_pert_kernel_type
+ private
+ type(arg_type) :: meta_args(3) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, Wtheta), &
+ arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3) &
+ /)
+ type(func_type) :: meta_funcs(1) = (/ &
+ func_type(ANY_SPACE_9, GH_BASIS) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+ integer :: gh_shape = GH_EVALUATOR
+ contains
+ procedure, nopass :: initial_theta_pert_code
+ end type
+
+ !-------------------------------------------------------------------------------
+ ! Contained functions/subroutines
+ !-------------------------------------------------------------------------------
+ public :: initial_theta_pert_code
+
+contains
+
+ !> @brief Computes the initial theta perturbation field
+ !! @param[in] nlayers Number of layers
+ !! @param[in,out] theta Potential temperature
+ !! @param[in] chi_1 First component of the chi coordinate field
+ !! @param[in] chi_2 Second component of the chi coordinate field
+ !! @param[in] chi_3 Third component of the chi coordinate field
+ !! @param[in] panel_id A field giving the ID for mesh panels
+ !! @param[in] ndf_wtheta Number of degrees of freedom per cell for wtheta
+ !! @param[in] undf_wtheta Number of total degrees of freedom for wtheta
+ !! @param[in] map_wtheta Dofmap for the cell at the base of the column
+ !! @param[in] ndf_chi Number of degrees of freedom per cell for chi
+ !! @param[in] undf_chi Number of total degrees of freedom for chi
+ !! @param[in] map_chi Dofmap for the cell at the base of the column
+ !! @param[in] chi_basis Basis functions evaluated at Wtheta points
+ !! @param[in] ndf_pid Number of degrees of freedom per cell for panel_id
+ !! @param[in] undf_pid Number of unique degrees of freedom for panel_id
+ !! @param[in] map_pid Dofmap for the cell at the base of the column for panel_id
+ subroutine initial_theta_pert_code(nlayers, &
+ theta, &
+ chi_1, chi_2, chi_3, &
+ panel_id, &
+ ndf_wtheta, undf_wtheta, map_wtheta, &
+ ndf_chi, undf_chi, map_chi, chi_basis, &
+ ndf_pid, undf_pid, map_pid )
+
+ use analytic_temperature_profiles_mod, only : analytic_theta_pert
+ use sci_chi_transform_mod, only : chi2xyz
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers, ndf_wtheta, ndf_chi, ndf_pid
+ integer(kind=i_def), intent(in) :: undf_wtheta, undf_chi, undf_pid
+
+ integer(kind=i_def), dimension(ndf_wtheta), intent(in) :: map_wtheta
+ integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi
+ integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid
+
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: theta
+ real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3
+ real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id
+ real(kind=r_def), dimension(1,ndf_chi,ndf_wtheta), intent(in) :: chi_basis
+
+ ! Internal variables
+ real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
+ real(kind=r_def) :: coords(3), xyz(3)
+ integer(kind=i_def) :: df, dfc, k, ipanel
+
+ ipanel = int(panel_id(map_pid(1)), i_def)
+
+ ! Compute the pointwise theta profile
+
+ do k = 0, nlayers-1
+ do dfc = 1, ndf_chi
+ chi_1_e(dfc) = chi_1( map_chi(dfc) + k)
+ chi_2_e(dfc) = chi_2( map_chi(dfc) + k)
+ chi_3_e(dfc) = chi_3( map_chi(dfc) + k)
+ end do
+
+ do df = 1, ndf_wtheta
+ coords(:) = 0.0_r_def
+ do dfc = 1, ndf_chi
+ coords(1) = coords(1) + chi_1_e(dfc)*chi_basis(1,dfc,df)
+ coords(2) = coords(2) + chi_2_e(dfc)*chi_basis(1,dfc,df)
+ coords(3) = coords(3) + chi_3_e(dfc)*chi_basis(1,dfc,df)
+ end do
+
+ call chi2xyz(coords(1), coords(2), coords(3), &
+ ipanel, xyz(1), xyz(2), xyz(3))
+ theta(map_wtheta(df) + k) = analytic_theta_pert(xyz, test)
+
+ end do
+ end do
+
+ end subroutine initial_theta_pert_code
+
+end module initial_theta_pert_kernel_mod
diff --git a/science/gungho/source/kernel/initialisation/initial_theta_v_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/initial_theta_v_kernel_mod.F90
new file mode 100644
index 000000000..e073b2fb6
--- /dev/null
+++ b/science/gungho/source/kernel/initialisation/initial_theta_v_kernel_mod.F90
@@ -0,0 +1,278 @@
+!-----------------------------------------------------------------------------
+! Copyright (c) 2017, Met Office, on behalf of HMSO and Queen's Printer
+! For further details please refer to the file LICENCE.original which you
+! should have received as part of this distribution.
+!-----------------------------------------------------------------------------
+!
+!-------------------------------------------------------------------------------
+
+!> @brief Computes the initial theta perturbation field
+module initial_theta_v_kernel_mod
+
+ use argument_mod, only: arg_type, func_type, &
+ GH_FIELD, GH_REAL, &
+ GH_WRITE, GH_READ, &
+ ANY_SPACE_9, GH_BASIS, &
+ ANY_DISCONTINUOUS_SPACE_3, &
+ CELL_COLUMN, GH_EVALUATOR, &
+ GH_INTEGER, GH_SCALAR
+ use constants_mod, only: r_def, i_def, l_def, PI
+ use fs_continuity_mod, only: Wtheta
+ use kernel_mod, only: kernel_type
+ use idealised_config_mod, only: test_squall_line, test_supercell
+
+ implicit none
+
+ private
+
+ !-------------------------------------------------------------------------------
+ ! Public types
+ !-------------------------------------------------------------------------------
+ !> The type declaration for the kernel. Contains the metadata needed by the Psy layer
+ type, public, extends(kernel_type) :: initial_theta_v_kernel_type
+ private
+ type(arg_type) :: meta_args(8) = (/ &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, Wtheta), &
+ arg_type(GH_FIELD, GH_REAL, GH_WRITE, 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*3, GH_REAL, GH_READ, ANY_SPACE_9), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &
+ arg_type(GH_SCALAR, GH_INTEGER, GH_READ) &
+ /)
+ type(func_type) :: meta_funcs(1) = (/ &
+ func_type(ANY_SPACE_9, GH_BASIS) &
+ /)
+ integer :: operates_on = CELL_COLUMN
+ integer :: gh_shape = GH_EVALUATOR
+ contains
+ procedure, nopass :: initial_theta_v_code
+ end type
+
+ !-------------------------------------------------------------------------------
+ ! Contained functions/subroutines
+ !-------------------------------------------------------------------------------
+ public :: initial_theta_v_code
+
+contains
+
+ subroutine initial_theta_v_code(nlayers, &
+ theta_v, exner_wt, &
+ theta_eq, exner_eq, &
+ height_wth, &
+ chi_1, chi_2, chi_3, &
+ panel_id, test, &
+ ndf_wtheta, undf_wtheta, map_wtheta, &
+ ndf_chi, undf_chi, map_chi, chi_basis, &
+ ndf_pid, undf_pid, map_pid )
+
+ use analytic_temperature_profiles_mod, only : integrate_theta_v, integrate_exner_surf
+ use sci_chi_transform_mod, only : chi2llr
+
+ implicit none
+
+ ! Arguments
+ integer(kind=i_def), intent(in) :: nlayers, ndf_wtheta, ndf_chi, ndf_pid
+ integer(kind=i_def), intent(in) :: undf_wtheta, undf_chi, undf_pid, test
+
+ integer(kind=i_def), dimension(ndf_wtheta), intent(in) :: map_wtheta
+ integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi
+ integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid
+
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: theta_v
+ real(kind=r_def), dimension(undf_wtheta), intent(inout) :: exner_wt
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: theta_eq
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: exner_eq
+ real(kind=r_def), dimension(undf_wtheta), intent(in) :: height_wth
+ real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3
+ real(kind=r_def), dimension(undf_pid), intent(in) :: panel_id
+ real(kind=r_def), dimension(1,ndf_chi,ndf_wtheta), intent(in) :: chi_basis
+
+ ! Internal variables
+ real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
+ real(kind=r_def) :: coords(3)
+ integer(kind=i_def) :: dfc, k, ipanel
+ real(kind=r_def) :: lon, lat, radius, dlat
+ real(kind=r_def) :: z_km1, z_k, z_kp1
+ integer(kind=i_def) :: i, lat_index
+ logical(kind=l_def) :: lat_index_found
+
+ integer(kind=i_def), parameter :: num_iterations = 10
+ integer(kind=i_def), parameter :: num_quad_points = 1200
+
+ real(kind=r_def), dimension(nlayers+1, num_quad_points) :: theta_v_prev
+ real(kind=r_def), dimension(nlayers+1, num_quad_points) :: theta_v_tab
+ real(kind=r_def), dimension(nlayers+1, num_quad_points) :: exner_tab
+ real(kind=r_def), dimension(num_quad_points) :: lat_points
+ real(kind=r_def), dimension(nlayers+1, num_quad_points) :: dtheta_v_dz
+ real(kind=r_def), dimension(nlayers+1) :: ueq2
+ real(kind=r_def), dimension(nlayers+1) :: dueq2_dz
+
+ real(kind=r_def) :: z, zS, dzU, Us, Uc, A, B, C
+
+ ipanel = int(panel_id(map_pid(1)), i_def)
+
+ ! ======================================================================== !
+ ! Set up coordinates
+ ! ======================================================================== !
+
+ ! Get latitude ---------------------------------------------------------
+ do dfc = 1, ndf_chi
+ chi_1_e(dfc) = chi_1( map_chi(dfc) )
+ chi_2_e(dfc) = chi_2( map_chi(dfc) )
+ chi_3_e(dfc) = chi_3( map_chi(dfc) )
+ end do
+
+ coords(:) = 0.0_r_def
+ do dfc = 1, ndf_chi
+ coords(1) = coords(1) + chi_1_e(dfc)*chi_basis(1,dfc,1)
+ coords(2) = coords(2) + chi_2_e(dfc)*chi_basis(1,dfc,1)
+ coords(3) = coords(3) + chi_3_e(dfc)*chi_basis(1,dfc,1)
+ end do
+
+ call chi2llr(coords(1), coords(2), coords(3), ipanel, lon, lat, radius)
+
+ ! Make an array of latitude points for numerical integral --------------
+ if (lat >= 0.0_r_def) then
+ dlat = 0.5_r_def * PI / REAL(num_quad_points-1, r_def)
+ else
+ dlat = -0.5_r_def * PI / REAL(num_quad_points-1, r_def)
+ end if
+ lat_points(1) = 0.0_r_def
+ lat_index_found = .false.
+ do i = 2, num_quad_points
+ lat_points(i) = lat_points(i-1) + dlat
+ if (lat >= 0.0_r_def .and. lat_points(i) > lat .and. .not. lat_index_found) then
+ lat_index = i
+ lat_index_found = .true.
+ else if (lat < 0.0_r_def .and. lat_points(i) < lat .and. .not. lat_index_found) then
+ lat_index = i
+ lat_index_found = .true.
+ end if
+ end do
+ if (lat >= 0.0_r_def) then
+ lat_points(num_quad_points) = PI/2.0_r_def
+ else
+ lat_points(num_quad_points) = -PI/2.0_r_def
+ end if
+
+ ! ======================================================================== !
+ ! Calculate du^2/dz at the equator, for each level
+ ! ======================================================================== !
+ if ( test == test_squall_line ) then
+ zS = 2500.0_r_def
+ dzU = 1000.0_r_def
+ Us = 12.0_r_def
+ Uc = 5.0_r_def
+ else if ( test == test_supercell ) then
+ zS = 5000.0_r_def
+ dzU = 1000.0_r_def
+ Us = 30.0_r_def
+ Uc = 15.0_r_def
+ end if
+
+ A = 0.25_r_def * (-zs**2 + 2.0_r_def*zS*dzU - dzU**2) / (zS * dzU)
+ B = 0.5_r_def * (zS + dzU) / dzU
+ C = - 0.25_r_def * (zS / dzU)
+
+ do k = 1, nlayers+1
+ z = height_wth(map_wtheta(1)+k-1)
+ if (z < zS - dzU) then
+ ueq2(k) = (Us*z/zS - Uc)**2
+ else if (ABS(z - zS) < dzU) then
+ ueq2(k) = ((A + B*z/zS + C*(z/zS)**2)*Us - Uc)**2
+ else
+ ueq2(k) = (Us - Uc)**2
+ end if
+ end do
+
+ ! Fit quadratic polynomial through theta_v points, and take vertical
+ ! derivative at each point
+ ! Bottom point, assume linear gradient
+ dueq2_dz(1) = &
+ (ueq2(2) - ueq2(1)) / (height_wth(map_wtheta(1)+1) - height_wth(map_wtheta(1)))
+
+ do k = 2, nlayers
+ ! Calculate derivative at height_wth(map_wtheta(1)+k)
+ z_km1 = height_wth(map_wtheta(1)+k-1)
+ z_k = height_wth(map_wtheta(1)+k)
+ z_kp1 = height_wth(map_wtheta(1)+k+1)
+ dueq2_dz(k) = ( &
+ ueq2(k-1) * (z_k - z_kp1) / ((z_km1 - z_k) * (z_km1 - z_kp1)) &
+ + ueq2(k) * (2.0_r_def*z_k - z_km1 - z_kp1) &
+ / ((z_k - z_km1) * (z_k - z_kp1)) &
+ + ueq2(k+1) * (z_k - z_km1) / ((z_kp1 - z_k) * (z_kp1 - z_km1)) &
+ )
+ end do
+ ! Top point, assume linear gradient
+ dueq2_dz(nlayers+1) = &
+ (ueq2(nlayers+1) - ueq2(nlayers)) &
+ / (height_wth(map_wtheta(1)+nlayers) - height_wth(map_wtheta(1)+nlayers-1))
+
+
+ ! ======================================================================== !
+ ! Calculate theta_v
+ ! ======================================================================== !
+
+ ! Initialise theta_v_prev
+ do k = 1, nlayers+1
+ theta_v_prev(k,:) = theta_eq(map_wtheta(1)+k-1)
+ end do
+
+ do i = 1, num_iterations
+ ! Fit quadratic polynomial through theta_v points, and take vertical
+ ! derivative at each point
+ ! Bottom point, assume linear gradient
+ dtheta_v_dz(1,:) = &
+ (theta_v_prev(2,:) - theta_v_prev(1,:)) &
+ / (height_wth(map_wtheta(1)+1) - height_wth(map_wtheta(1)))
+
+ do k = 2, nlayers
+ ! Calculate derivative at height_wth(map_wtheta(1)+k)
+ z_km1 = height_wth(map_wtheta(1)+k-1)
+ z_k = height_wth(map_wtheta(1)+k)
+ z_kp1 = height_wth(map_wtheta(1)+k+1)
+ dtheta_v_dz(k,:) = ( &
+ theta_v_prev(k-1,:) * (z_k - z_kp1) / ((z_km1 - z_k) * (z_km1 - z_kp1)) &
+ + theta_v_prev(k,:) * (2.0_r_def*z_k - z_km1 - z_kp1) &
+ / ((z_k - z_km1) * (z_k - z_kp1)) &
+ + theta_v_prev(k+1,:) * (z_k - z_km1) / ((z_kp1 - z_k) * (z_kp1 - z_km1)) &
+ )
+ end do
+ ! Top point, assume linear gradient
+ dtheta_v_dz(nlayers+1,:) = &
+ (theta_v_prev(nlayers+1,:) - theta_v_prev(nlayers,:)) &
+ / (height_wth(map_wtheta(1)+nlayers) - height_wth(map_wtheta(1)+nlayers-1))
+
+ call integrate_theta_v( &
+ theta_v_tab, theta_v_prev, dtheta_v_dz, &
+ theta_eq(map_wtheta(1) : map_wtheta(1)+nlayers), dueq2_dz, &
+ height_wth(map_wtheta(1) : map_wtheta(1)+nlayers), &
+ lat_points, nlayers+1, num_quad_points, test &
+ )
+
+ theta_v_prev(:,:) = theta_v_tab(:,:)
+ end do
+
+ call integrate_exner_surf( &
+ exner_tab, theta_v_tab, &
+ exner_eq(map_wtheta(1) : map_wtheta(1)+nlayers), &
+ height_wth(map_wtheta(1) : map_wtheta(1)+nlayers), &
+ lat_points, nlayers+1, num_quad_points, test &
+ )
+
+ ! Determine final value
+ do k = 0, nlayers
+ theta_v(map_wtheta(1)+k) = theta_v_tab(k+1,lat_index) &
+ + (theta_v_tab(k+1,lat_index+1) - theta_v_tab(k+1,lat_index)) &
+ * (lat - lat_points(lat_index)) / (lat_points(lat_index+1) - lat_points(lat_index))
+ exner_wt(map_wtheta(1)+k) = exner_tab(k+1,lat_index) &
+ + (exner_tab(k+1,lat_index+1) - exner_tab(k+1,lat_index)) &
+ * (lat - lat_points(lat_index)) / (lat_points(lat_index+1) - lat_points(lat_index))
+ end do
+
+
+ end subroutine initial_theta_v_code
+
+end module initial_theta_v_kernel_mod
diff --git a/science/gungho/source/kernel/initialisation/set_exner_kernel_mod.F90 b/science/gungho/source/kernel/initialisation/set_exner_kernel_mod.F90
index ebbde1980..cf2f71b74 100644
--- a/science/gungho/source/kernel/initialisation/set_exner_kernel_mod.F90
+++ b/science/gungho/source/kernel/initialisation/set_exner_kernel_mod.F90
@@ -15,7 +15,7 @@ module set_exner_kernel_mod
GH_BASIS, GH_DIFF_BASIS, &
CELL_COLUMN, GH_QUADRATURE_XYoZ
use constants_mod, only : r_def, i_def
- use fs_continuity_mod, only : W3
+ use fs_continuity_mod, only : W3, Wtheta
use idealised_config_mod, only : test
use kernel_mod, only : kernel_type
@@ -31,8 +31,9 @@ module set_exner_kernel_mod
!>
type, public, extends(kernel_type) :: set_exner_kernel_type
private
- type(arg_type) :: meta_args(4) = (/ &
+ type(arg_type) :: meta_args(5) = (/ &
arg_type(GH_FIELD, GH_REAL, GH_WRITE, W3), &
+ arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta), &
arg_type(GH_FIELD*3, GH_REAL, GH_READ, ANY_SPACE_9), &
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &
arg_type(GH_SCALAR, GH_REAL, GH_READ) &
@@ -82,10 +83,12 @@ module set_exner_kernel_mod
!! @param[in] wqp_v Weights of vertical quadrature points
subroutine set_exner_code(nlayers, &
exner, &
+ height_wth, &
chi_1, chi_2, chi_3, panel_id, &
time, &
ndf_w3, undf_w3, map_w3, &
w3_basis, &
+ ndf_wth, undf_wth, map_wth, &
ndf_chi, undf_chi, map_chi, &
chi_basis, chi_diff_basis, &
ndf_pid, undf_pid, map_pid, &
@@ -103,15 +106,17 @@ subroutine set_exner_code(nlayers, &
! Arguments
integer(kind=i_def), intent(in) :: nlayers
- integer(kind=i_def), intent(in) :: ndf_w3, ndf_chi, ndf_pid
- integer(kind=i_def), intent(in) :: undf_w3, undf_chi, undf_pid
+ integer(kind=i_def), intent(in) :: ndf_w3, ndf_chi, ndf_pid, ndf_wth
+ integer(kind=i_def), intent(in) :: undf_w3, undf_chi, undf_pid, undf_wth
integer(kind=i_def), intent(in) :: nqp_h, nqp_v
integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3
+ integer(kind=i_def), dimension(ndf_wth), intent(in) :: map_wth
integer(kind=i_def), dimension(ndf_pid), intent(in) :: map_pid
integer(kind=i_def), dimension(ndf_chi), intent(in) :: map_chi
real(kind=r_def), dimension(undf_w3), intent(inout) :: exner
+ real(kind=r_def), dimension(undf_wth), intent(in) :: height_wth
real(kind=r_def), dimension(1,ndf_w3,nqp_h,nqp_v), intent(in) :: w3_basis
real(kind=r_def), intent(in) :: time
real(kind=r_def), dimension(undf_chi), intent(in) :: chi_1, chi_2, chi_3
@@ -130,11 +135,13 @@ subroutine set_exner_code(nlayers, &
real(kind=r_def), dimension(nqp_h,nqp_v) :: dj
real(kind=r_def), dimension(3,3,nqp_h,nqp_v) :: jac
real(kind=r_def), dimension(ndf_chi) :: chi_1_e, chi_2_e, chi_3_e
- real(kind=r_def) :: exner_ref, integrand
+ real(kind=r_def) :: exner_ref, integrand, surface_height
real(kind=r_def) :: xyz(3), coords(3)
ipanel = int(panel_id(map_pid(1)), i_def)
+ surface_height = height_wth(map_wth(1))
+
! Compute the RHS & LHS integrated over one cell and solve
do k = 0, nlayers-1
do df1 = 1, ndf_chi
@@ -162,7 +169,7 @@ subroutine set_exner_code(nlayers, &
call chi2xyz(coords(1), coords(2), coords(3), &
ipanel, xyz(1), xyz(2), xyz(3))
- exner_ref = analytic_pressure(xyz, test, time)
+ exner_ref = analytic_pressure(xyz, test, time, surface_height)
integrand = w3_basis(1,df1,qp1,qp2) * exner_ref * dj(qp1,qp2)
rhs_e(df1) = rhs_e(df1) + wqp_h(qp1)*wqp_v(qp2)*integrand
diff --git a/science/gungho/source/mesh/gungho_extrusion_mod.f90 b/science/gungho/source/mesh/gungho_extrusion_mod.f90
index b5ac0c105..d991fe89e 100644
--- a/science/gungho/source/mesh/gungho_extrusion_mod.f90
+++ b/science/gungho/source/mesh/gungho_extrusion_mod.f90
@@ -37,6 +37,8 @@ module gungho_extrusion_mod
method_quadratic, &
method_geometric, &
method_dcmip, &
+ method_dcmip_mountain, &
+ method_dcmip_gw, &
method_um_L38_29t_9s_40km, &
method_um_L85_50t_35s_85km, &
method_um_L70_61t_9s_40km, &
@@ -169,6 +171,34 @@ module gungho_extrusion_mod
module procedure dcmip_extrusion_constructor
end interface dcmip_extrusion_type
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !> @brief Extrudes using DCMIP mountain scheme.
+ !>
+ type, public, extends(extrusion_type) :: dcmip_mountain_extrusion_type
+ private
+ contains
+ private
+ procedure, public :: extrude => dcmip_mountain_extrude
+ end type dcmip_mountain_extrusion_type
+
+ interface dcmip_mountain_extrusion_type
+ module procedure dcmip_mountain_extrusion_constructor
+ end interface dcmip_mountain_extrusion_type
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !> @brief Extrudes using DCMIP Gravity Wave scheme.
+ !>
+ type, public, extends(extrusion_type) :: dcmip_gw_extrusion_type
+ private
+ contains
+ private
+ procedure, public :: extrude => dcmip_gw_extrude
+ end type dcmip_gw_extrusion_type
+
+ interface dcmip_gw_extrusion_type
+ module procedure dcmip_gw_extrusion_constructor
+ end interface dcmip_gw_extrusion_type
+
contains
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -709,6 +739,149 @@ function dcmip_func(eta_uni) result(eta)
end function dcmip_func
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !> @brief Creates a dcmip_mountain_extrusion_type object.
+ !>
+ !> @param[in] atmosphere_bottom Bottom of the atmosphere in meters.
+ !> @param[in] atmosphere_top Top of the atmosphere in meters.
+ !> @param[in] number_of_layers Number of layers in the atmosphere.
+ !> @param[in] extrusion_id Identifier of extrusion type.
+ !>
+ !> @return New dcmip_extrusion_type object.
+ !>
+ function dcmip_mountain_extrusion_constructor( atmosphere_bottom, &
+ atmosphere_top, &
+ number_of_layers, &
+ extrusion_id ) result(new)
+
+ implicit none
+
+ real(r_def), intent(in) :: atmosphere_bottom
+ real(r_def), intent(in) :: atmosphere_top
+ integer(i_def), intent(in) :: number_of_layers
+ integer(i_def), intent(in) :: extrusion_id
+
+ type(dcmip_mountain_extrusion_type) :: new
+
+ call new%extrusion_constructor( atmosphere_bottom, atmosphere_top, &
+ number_of_layers, extrusion_id )
+
+ end function dcmip_mountain_extrusion_constructor
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !> @brief Extrudes the mesh using the DCMIP mountain scheme.
+ !>
+ !> @param[out] eta Nondimensional vertical coordinate.
+ !>
+ subroutine dcmip_mountain_extrude( this, eta )
+
+ implicit none
+
+ class(dcmip_mountain_extrusion_type), intent(in) :: this
+ real(r_def), intent(out) :: eta(0:)
+
+ real(r_def), parameter :: dz_low = 100.0_r_def
+ real(r_def), parameter :: dz_high = 500.0_r_def
+ real(r_def), parameter :: h0 = 950.0_r_def
+ real(r_def), parameter :: zU = 6000.0_r_def
+ real(r_def), parameter :: zT = 20000.0_r_def
+ real(r_def), parameter :: gamma = 1.01679_r_def
+
+
+ integer(i_def) :: k
+
+ if (this%get_number_of_layers() /= 57)then
+ call log_event( "Extrusion DCMIP mountain reqires 57 levels", log_level_error )
+ end if
+
+ eta(0) = 0.0_r_def
+ do k = 1, this%get_number_of_layers()
+ ! Lower part of extrusion
+ if (eta(k-1)*zT < h0) then
+ eta(k) = eta(k-1) + dz_low / zT
+
+ else if (eta(k-1)*zT > zU) then
+ ! Upper part of extrusion
+ eta(k) = eta(k-1) + dz_high / zT
+
+ else
+ ! Middle part of extrusion
+ eta(k) = eta(k-1) + MIN(((eta(k-1) - eta(k-2))*zT)**gamma, dz_high) / zT
+ end if
+ end do
+ eta(this%get_number_of_layers()) = 1.0_r_def
+
+ end subroutine dcmip_mountain_extrude
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !> @brief Creates a dcmip_gw_extrusion_type object.
+ !>
+ !> @param[in] atmosphere_bottom Bottom of the atmosphere in meters.
+ !> @param[in] atmosphere_top Top of the atmosphere in meters.
+ !> @param[in] number_of_layers Number of layers in the atmosphere.
+ !> @param[in] extrusion_id Identifier of extrusion type.
+ !>
+ !> @return New dcmip_extrusion_type object.
+ !>
+ function dcmip_gw_extrusion_constructor( atmosphere_bottom, &
+ atmosphere_top, &
+ number_of_layers, &
+ extrusion_id ) result(new)
+
+ implicit none
+
+ real(r_def), intent(in) :: atmosphere_bottom
+ real(r_def), intent(in) :: atmosphere_top
+ integer(i_def), intent(in) :: number_of_layers
+ integer(i_def), intent(in) :: extrusion_id
+
+ type(dcmip_gw_extrusion_type) :: new
+
+ call new%extrusion_constructor( atmosphere_bottom, atmosphere_top, &
+ number_of_layers, extrusion_id )
+
+ end function dcmip_gw_extrusion_constructor
+
+ !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ !> @brief Extrudes the mesh using the DCMIP mountain scheme.
+ !>
+ !> @param[out] eta Nondimensional vertical coordinate.
+ !>
+ subroutine dcmip_gw_extrude( this, eta )
+
+ implicit none
+
+ class(dcmip_gw_extrusion_type), intent(in) :: this
+ real(r_def), intent(out) :: eta(0:)
+
+ real(r_def), parameter :: xi = (700.0_r_def / 21000.0_r_def)
+ real(r_def), parameter :: dz_min = 100.0_r_def
+ real(r_def), parameter :: H = 40000.0_r_def
+
+ real(r_def) :: z, dz_max, dz
+ integer(i_def) :: k
+
+ if (this%get_number_of_layers() /= 88)then
+ call log_event( "Extrusion DCMIP GW reqires 88 levels", log_level_error )
+ end if
+
+ eta(0) = 0.0_r_def
+ z = 0.0_r_def
+ dz = dz_min
+ do k = 1, this%get_number_of_layers()
+ if (H - z < REAL(this%get_number_of_layers() - k + 1, r_def)*(dz_min + xi*z)) then
+ dz_max = (H - z) / REAL(this%get_number_of_layers() - k + 1, r_def)
+ dz = dz_max
+ else
+ dz = dz_min + xi*z
+ end if
+ z = z + dz
+ eta(k) = z / H
+ end do
+ eta(this%get_number_of_layers()) = 1.0_r_def
+
+ end subroutine dcmip_gw_extrude
+
!> @brief Creates vertical mesh extrusion
!> @details Creates vertical mesh with nlayers.
!> @param method Extrusion method to apply
@@ -800,6 +973,16 @@ function create_extrusion( method, &
atmosphere_bottom, domain_height, &
number_of_layers, prime_extrusion ) )
+ case (method_dcmip_mountain)
+ allocate( extrusion, source=dcmip_mountain_extrusion_type( atmosphere_bottom, &
+ domain_height, &
+ number_of_layers, &
+ PRIME_EXTRUSION ) )
+ case (method_dcmip_gw)
+ allocate( extrusion, source=dcmip_gw_extrusion_type( atmosphere_bottom, &
+ domain_height, &
+ number_of_layers, &
+ PRIME_EXTRUSION ) )
case default
write( log_scratch_space, '(A)' ) &
trim(module_name) // ': Unrecognised extrusion method: ' //&
diff --git a/science/gungho/source/orography/dcmip_2025_orography_mod.F90 b/science/gungho/source/orography/dcmip_2025_orography_mod.F90
new file mode 100644
index 000000000..fcef6eff8
--- /dev/null
+++ b/science/gungho/source/orography/dcmip_2025_orography_mod.F90
@@ -0,0 +1,192 @@
+!-----------------------------------------------------------------------------
+! (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 Calculates DCMIP 2025 orography.
+!>
+!> @details This module contains type definition and routines to calculate
+!> analytic orography profile of DCMIP200 case mountain function from
+!> spherical polar coordinates: longitude (lambda) and latitude (phi).
+!> Reference: Ulrich et al. (2012), Section 2.0.
+!> DCMIP200 mountain parameters in (lon,lat) coordinates are:
+!> mountain_height - Height of DCMIP200 mountain function (m),
+!> radius - Radius of DCMIP200 mountain function (radian),
+!> osc_half_width - Oscillation half-width of DCMIP200 mountain (radian),
+!> lambda_centre - Longitudinal centre of DCMIP200 mountain function (radian),
+!> phi_centre - Latitudinal centre of DCMIP200 mountain function (radian).
+!-------------------------------------------------------------------------------
+module dcmip_2025_orography_mod
+
+ use constants_mod, only : r_def, i_def
+ use analytic_orography_mod, only : analytic_orography_type
+ use log_mod, only : log_event, LOG_LEVEL_ERROR
+
+ implicit none
+
+ private
+
+ type, public, extends(analytic_orography_type) :: dcmip_breaking_gw_type
+
+ private
+
+ contains
+
+ procedure, public, pass(self) :: analytic_orography => dcmip_breaking_gw_orography
+
+ end type dcmip_breaking_gw_type
+
+ type, public, extends(analytic_orography_type) :: dcmip_gap_type
+
+ private
+
+ contains
+
+ procedure, public, pass(self) :: analytic_orography => dcmip_gap_orography
+
+ end type dcmip_gap_type
+
+ type, public, extends(analytic_orography_type) :: dcmip_vortex_type
+
+ private
+
+ contains
+
+ procedure, public, pass(self) :: analytic_orography => dcmip_vortex_orography
+
+ end type dcmip_vortex_type
+
+ interface dcmip_breaking_gw_type
+ module procedure dcmip_breaking_gw_constructor
+ end interface
+
+ interface dcmip_gap_type
+ module procedure dcmip_gap_constructor
+ end interface
+
+ interface dcmip_vortex_type
+ module procedure dcmip_vortex_constructor
+ end interface
+
+contains
+
+ !=============================================================================
+ type(dcmip_breaking_gw_type) function dcmip_breaking_gw_constructor( ) result(self)
+
+ implicit none
+
+ return
+ end function dcmip_breaking_gw_constructor
+
+ !=============================================================================
+ function dcmip_breaking_gw_orography(self, chi_1, chi_2) result(chi_surf)
+
+ use constants_mod, only : PI
+
+ implicit none
+
+ ! Arguments
+ class(dcmip_breaking_gw_type), intent(in) :: self
+ real(kind=r_def), intent(in) :: chi_1, chi_2
+
+ real(kind=r_def) :: chi_surf
+
+
+ call log_event('Breaking GW orography not yet implemented', LOG_LEVEL_ERROR)
+ chi_surf = 0.0_r_def
+
+ return
+ end function dcmip_breaking_gw_orography
+
+
+ !=============================================================================
+ type(dcmip_gap_type) function dcmip_gap_constructor( ) result(self)
+
+ implicit none
+
+ return
+ end function dcmip_gap_constructor
+
+ !=============================================================================
+ function dcmip_gap_orography(self, chi_1, chi_2) result(chi_surf)
+
+ use constants_mod, only : PI
+ use planet_config_mod, only : scaling_factor, scaled_radius
+
+ implicit none
+
+ ! Arguments
+ class(dcmip_gap_type), intent(in) :: self
+ real(kind=r_def), intent(in) :: chi_1, chi_2
+
+ real(kind=r_def), parameter :: e1 = 10.0_r_def
+ real(kind=r_def), parameter :: e2 = 10.0_r_def
+ real(kind=r_def), parameter :: e3 = 10.0_r_def
+ real(kind=r_def), parameter :: h0 = 1500.0_r_def
+ real(kind=r_def), parameter :: H = 0.1_r_def*h0
+ real(kind=r_def), parameter :: lon_c = PI/2.0_r_def
+ real(kind=r_def), parameter :: lat_c = 0.0_r_def
+ real(kind=r_def), parameter :: x_lon = 800000.0_r_def
+ real(kind=r_def), parameter :: x_lat = 6000000.0_r_def
+ real(kind=r_def), parameter :: x_gap = 1000000.0_r_def
+
+ real(kind=r_def) :: d1, d2, d3
+ real(kind=r_def) :: scaled_x_lon, scaled_x_lat, scaled_x_gap
+ real(kind=r_def) :: chi_surf
+
+ scaled_x_lon = x_lon / scaling_factor
+ scaled_x_lat = x_lat / scaling_factor
+ scaled_x_gap = x_gap / scaling_factor
+
+ d1 = 0.5_r_def * scaled_x_lon / scaled_radius * LOG(h0/H)**(-1.0_r_def/e1)
+ d2 = 0.5_r_def * scaled_x_lat / scaled_radius * LOG(h0/H)**(-1.0_r_def/e2)
+ d3 = 0.5_r_def * scaled_x_gap / scaled_radius * LOG(h0/H)**(-1.0_r_def/e3)
+
+ chi_surf = h0 * EXP(-((chi_1 - lon_c)/d1)**e1 - ((chi_2 - lat_c)/d2)**e2) &
+ * (1.0_r_def - EXP(-((chi_2 - lat_c)/d3)**e3))
+
+ return
+ end function dcmip_gap_orography
+
+ !=============================================================================
+ type(dcmip_vortex_type) function dcmip_vortex_constructor( ) result(self)
+
+ implicit none
+
+ return
+ end function dcmip_vortex_constructor
+
+ !=============================================================================
+ function dcmip_vortex_orography(self, chi_1, chi_2) result(chi_surf)
+
+ use constants_mod, only : PI
+ use planet_config_mod, only : scaling_factor, scaled_radius
+
+ implicit none
+
+ ! Arguments
+ class(dcmip_vortex_type), intent(in) :: self
+ real(kind=r_def), intent(in) :: chi_1, chi_2
+
+ real(kind=r_def), parameter :: h0 = 1500.0_r_def
+ real(kind=r_def), parameter :: d = 250000.0_r_def
+ real(kind=r_def), parameter :: lon_c = PI / 2.0_r_def
+ real(kind=r_def), parameter :: lat_c = PI / 9.0_r_def
+
+ real(kind=r_def) :: chi_surf, scaled_d, r
+
+ ! Great arc distance
+ r = scaled_radius * ACOS( &
+ SIN(lat_c) * SIN(chi_2) + COS(lat_c) * COS(chi_2) * COS(chi_1 - lon_c) &
+ )
+ scaled_d = d / scaling_factor
+
+ chi_surf = h0 * EXP(-(r/scaled_d)**2)
+
+ return
+ end function dcmip_vortex_orography
+
+end module dcmip_2025_orography_mod
+
diff --git a/science/gungho/source/physics/physics_common_mod.f90 b/science/gungho/source/physics/physics_common_mod.f90
index ba90e00a4..570c50ae6 100644
--- a/science/gungho/source/physics/physics_common_mod.f90
+++ b/science/gungho/source/physics/physics_common_mod.f90
@@ -40,9 +40,11 @@ function qsaturation (T, p)
if (T > qsa3 .and. p * exp (qsa2 * (t - tk0c) / (T - qsa3)) > qsa4) then
qsaturation=qsa1/(p*exp(qsa2*(t-tk0c)/(T-qsa3))-qsa4)
+ qsaturation = (3.8_r_def/p)*exp(17.27_r_def*(T-273.0_r_def)/(T-36.0_r_def))
else
qsaturation=999.0_r_def
end if
+
end function qsaturation
end module physics_common_mod