diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md index 00c6e6edb..c2ff84a0f 100644 --- a/CONTRIBUTORS.md +++ b/CONTRIBUTORS.md @@ -34,4 +34,5 @@ | thomasmelvin | Thomas Melvin | Met Office | 2026-01-15 | | tinyendian | Wolfgang Hayek | Earth Sciences New Zealand | 2026-02-02 | | DanStoneMO | Daniel Stone | Met Office | 2026-02-26 | -| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | \ No newline at end of file +| ericaneininger | Erica Neininger | Met Office | 2026-03-02 | +| iboutle | Ian Boutle | Met Office | 2026-03-10 | \ No newline at end of file diff --git a/applications/adjoint_tests/source/algorithm/solver/adjt_semi_implicit_solver_step_alg_mod.x90 b/applications/adjoint_tests/source/algorithm/solver/adjt_semi_implicit_solver_step_alg_mod.x90 index 3e61588ff..75fc3b6e9 100644 --- a/applications/adjoint_tests/source/algorithm/solver/adjt_semi_implicit_solver_step_alg_mod.x90 +++ b/applications/adjoint_tests/source/algorithm/solver/adjt_semi_implicit_solver_step_alg_mod.x90 @@ -266,7 +266,7 @@ contains ! -------------------------------------------------------------------------- call semi_implicit_solver_alg_step( state, rhs, moist_dyn_gas_law, mr, & - .false., first_iteration ) + first_iteration ) ! -------------------------------------------------------------------------- ! Compute diff --git a/applications/lfric_atm/metadata/field_def_diags.xml b/applications/lfric_atm/metadata/field_def_diags.xml index a17d884f9..8774940a2 100644 --- a/applications/lfric_atm/metadata/field_def_diags.xml +++ b/applications/lfric_atm/metadata/field_def_diags.xml @@ -418,8 +418,8 @@ - microphysics__ls_snow + microphysics__ls_graup - microphysics__ls_rain + microphysics__ls_snow + microphysics__ls_graup + microphysics__ls_snow - microphysics__ls_graup + microphysics__ls_rain + microphysics__ls_snow $DT * microphysics__ls_rain $DT * microphysics__ls_snow $DT * microphysics__ls_graup diff --git a/rose-stem/app/lfric_atm/opt/rose-app-cons_diag.conf b/rose-stem/app/lfric_atm/opt/rose-app-cons_diag.conf new file mode 100644 index 000000000..47475a468 --- /dev/null +++ b/rose-stem/app/lfric_atm/opt/rose-app-cons_diag.conf @@ -0,0 +1,2 @@ +[namelist:io] +write_conservation_diag=.true. diff --git a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc index a49293d58..f2b08ee4f 100644 --- a/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc +++ b/rose-stem/site/common/lfric_atm/tasks_lfric_atm.cylc @@ -421,7 +421,7 @@ {% elif task_ns.conf_name == "nwp_gal9_debug-C12" %} {% do task_dict.update({ - "opt_confs": ["um_dump","no_diags"], + "opt_confs": ["um_dump","no_diags","cons_diag"], "resolution": "C12", "DT": 1800, "tsteps": 4, diff --git a/science/adjoint/source/algorithm/solver/adj_semi_implicit_solver_alg_mod.x90 b/science/adjoint/source/algorithm/solver/adj_semi_implicit_solver_alg_mod.x90 index 3a1f0ee48..f0078ffe9 100644 --- a/science/adjoint/source/algorithm/solver/adj_semi_implicit_solver_alg_mod.x90 +++ b/science/adjoint/source/algorithm/solver/adj_semi_implicit_solver_alg_mod.x90 @@ -306,10 +306,8 @@ contains !> @param[in] moist_dyn_gas_law Gas law component of moist dynamics !! factors !> @param[in] mr Mixing ratio array - !> @param[in] write_moisture_diag Flag to control output of moisture - !! conservation diagnostics !> @param[in] first_iteration Flag for first inner iteration - subroutine step( self, state, rhs, moist_dyn_gas_law, mr, write_moisture_diag, first_iteration ) + subroutine step( self, state, rhs, moist_dyn_gas_law, mr, first_iteration ) use solver_constants_mod, only: get_normalisation_r_solver, get_im3_div_r_solver, get_normalisation use mr_indices_mod, only: nummr @@ -333,7 +331,7 @@ contains type(field_type), dimension(bundle_size), intent(inout) :: state, rhs type(field_type), intent(in) :: moist_dyn_gas_law type(field_type), dimension(nummr), intent(in) :: mr - logical(kind=l_def), intent(in) :: write_moisture_diag, first_iteration + logical(kind=l_def), intent(in) :: first_iteration ! Local variables type(field_type), pointer :: t_normalisation, u_normalisation @@ -359,13 +357,6 @@ contains mesh_id = state(igh_p)%get_mesh_id() - if (write_moisture_diag) then - call log_event( "adj_semi_implicit_solver_type%step: write_moisture_diag not coded", LOG_LEVEL_ERROR ) - ! Best guess for the prognostic fields state can be calculated as: - ! state_best_guess = rhs_n - rhs_np1 + state + rhs_adv + rhs_phys - ! where fortunately all of the above except "state" is already in "rhs_np1" - end if - ! Normalise theta & u residual t_normalisation => get_normalisation( Wtheta, mesh_id ) u_normalisation => get_normalisation( W2, mesh_id ) diff --git a/science/adjoint/source/algorithm/timestepping/atl_si_timestep_alg_mod.x90 b/science/adjoint/source/algorithm/timestepping/atl_si_timestep_alg_mod.x90 index 679e30c38..13ea7113c 100644 --- a/science/adjoint/source/algorithm/timestepping/atl_si_timestep_alg_mod.x90 +++ b/science/adjoint/source/algorithm/timestepping/atl_si_timestep_alg_mod.x90 @@ -530,7 +530,7 @@ contains self%ls_rhs_np1(:,outer,inner), & self%ls_moist_dyn_itns(gas_law, next_outer, next_inner), & self%ls_mr_itns(:,next_outer,next_inner), & - .false., first_iteration=(inner==1) ) + first_iteration=(inner==1) ) ! If not already done, update factors for moist dynamics if (.not. guess_np1 .and. self%use_moisture) then call moist_dyn_factors_alg( self%ls_moist_dyn_itns(:,next_outer,next_inner), & @@ -652,7 +652,7 @@ contains if (.not. bundle_is_zero( mixed_solver_a_tol, self%state, bundle_size )) then call self%adj_semi_implicit_solver%step( self%state, self%rhs_np1, & moist_dyn(gas_law), mr, & - .false., first_iteration=(inner==1) ) + first_iteration=(inner==1) ) end if ! Accelerators for inner loop convergence diff --git a/science/gungho/source/algorithm/diagnostics/moisture_conservation_alg_mod.x90 b/science/gungho/source/algorithm/diagnostics/moisture_conservation_alg_mod.x90 index 48ab982ca..ce952f204 100644 --- a/science/gungho/source/algorithm/diagnostics/moisture_conservation_alg_mod.x90 +++ b/science/gungho/source/algorithm/diagnostics/moisture_conservation_alg_mod.x90 @@ -68,8 +68,10 @@ contains ! Fields and values type(field_type) :: rho_X_shifted type(field_type) :: rho_d_shifted + type(field_type) :: total_water real(kind=r_def) :: total_dry real(kind=r_def) :: total_water_all(nummr) + real(kind=r_def) :: total_water_sum integer(kind=i_def) :: i_mr, shifted_mesh_id integer(kind=i_def) :: element_order_h, element_order_v integer(kind=i_def) :: nqp_h_exact, nqp_v_exact @@ -106,13 +108,19 @@ contains chi_shifted => get_coordinates(shifted_mesh_id) w3_sh_fs => function_space_collection%get_fs( shifted_mesh, 0, 0, W3 ) call rho_X_shifted%initialise( vector_space = w3_sh_fs ) + call mr(1)%copy_field_properties(total_water) + call invoke(setval_c(total_water, 0.0_r_def)) call obtain_shifted_rho(rho_d_shifted, rho_d) do i_mr = 1, nummr - call invoke( X_times_Y(rho_X_shifted, mr(i_mr), rho_d_shifted) ) + call invoke( inc_X_plus_Y(total_water, mr(i_mr)), & + X_times_Y(rho_X_shifted, mr(i_mr), rho_d_shifted) ) call compute_total_mass_alg(total_water_all(i_mr), rho_X_shifted, & shifted_mesh, method=integral_method_fe) end do + call invoke( X_times_Y(rho_X_shifted, total_water, rho_d_shifted) ) + call compute_total_mass_alg(total_water_sum, rho_X_shifted, & + shifted_mesh, method=integral_method_fe) ! Output diagnostic values to log write( log_scratch_space, '(2A,E32.12)' ) & @@ -124,7 +132,9 @@ contains total_water_all(i_mr) call log_event( log_scratch_space, LOG_LEVEL_INFO ) end do - + write( log_scratch_space, '(2A,E32.12)' ) & + 'Conservation: total water, ', stage, total_water_sum + call log_event( log_scratch_space, LOG_LEVEL_INFO ) nullify( chi, panel_id, w3_fs, w3_fs, mesh, shifted_mesh, chi_shifted ) if ( LPROF ) call stop_timing( id, 'diags.moisture_conservation' ) diff --git a/science/gungho/source/algorithm/diagnostics/moisture_fluxes_alg_mod.x90 b/science/gungho/source/algorithm/diagnostics/moisture_fluxes_alg_mod.x90 index 77c216e35..92a8ac4c7 100644 --- a/science/gungho/source/algorithm/diagnostics/moisture_fluxes_alg_mod.x90 +++ b/science/gungho/source/algorithm/diagnostics/moisture_fluxes_alg_mod.x90 @@ -68,12 +68,11 @@ contains type(field_type) :: fqw_total type(field_type), pointer :: ls_rain => null() type(field_type), pointer :: ls_snow => null() - type(field_type), pointer :: ls_graup => null() type(field_type), pointer :: conv_rain => null() type(field_type), pointer :: conv_snow => null() type(field_type), pointer :: moist_flux_bl => null() real(kind=r_def) :: bl_source_mass, mp_sink_mass, conv_sink_mass - real(kind=r_def) :: ls_rain_mass, ls_snow_mass, ls_graup_mass + real(kind=r_def) :: ls_rain_mass, ls_snow_mass real(kind=r_def) :: conv_rain_mass, conv_snow_mass integer(kind=i_def) :: mesh integer(tik) :: id @@ -83,7 +82,7 @@ contains mesh = area%get_mesh_id() call microphysics_fields%get_field('ls_rain', ls_rain) call microphysics_fields%get_field('ls_snow', ls_snow) - call microphysics_fields%get_field('ls_graup', ls_graup) + ! Graupel is included in ls_snow call convection_fields%get_field('conv_rain', conv_rain) call convection_fields%get_field('conv_snow', conv_snow) call turbulence_fields%get_field('moist_flux_bl', moist_flux_bl) @@ -120,13 +119,6 @@ contains dt & ), & sum_X(ls_snow_mass, mass), & - setval_c(mass, 0.0_r_def), & - compute_bottom_mass_from_flux_kernel_type( mass, & - ls_graup, & - area, & - dt & - ), & - sum_X(ls_graup_mass, mass), & ! Compute convective preciptation sink mass setval_c(mass, 0.0_r_def), & compute_bottom_mass_from_flux_kernel_type( mass, & @@ -144,7 +136,7 @@ contains sum_X(conv_snow_mass, mass) & ) - mp_sink_mass = ls_rain_mass + ls_snow_mass + ls_graup_mass + mp_sink_mass = ls_rain_mass + ls_snow_mass conv_sink_mass = conv_rain_mass + conv_snow_mass write( log_scratch_space, '(A, E32.12)' ) & @@ -161,7 +153,7 @@ contains call log_event( log_scratch_space, LOG_LEVEL_INFO ) nullify( fs_2d ) - nullify( ls_rain, ls_snow, ls_graup, conv_rain, conv_snow, moist_flux_bl ) + nullify( ls_rain, ls_snow, conv_rain, conv_snow, moist_flux_bl ) if ( LPROF ) call stop_timing( id, 'diags.moisture_fluxes' ) 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 f6523d797..53b224f34 100644 --- a/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90 +++ b/science/gungho/source/algorithm/physics/slow_physics_alg_mod.X90 @@ -365,7 +365,8 @@ contains #endif integer(tik) :: id if ( LPROF ) call start_timing( id, 'slow_physics' ) - + if (write_conservation_diag) & + call moisture_conservation_alg( rho, mr, 'Before slow' ) !-------------------------------------------------------------------- ! Initialisation of fields and flags !-------------------------------------------------------------------- diff --git a/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 b/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 index ca02473b8..aa4804fcc 100644 --- a/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 +++ b/science/gungho/source/algorithm/solver/semi_implicit_solver_alg_mod.x90 @@ -491,20 +491,16 @@ contains !> @param[in,out] rhs Residuals !> @param[in] moist_dyn_gas_law Gas law component of moist dynamics factors !> @param[in] mr Mixing ratio array - !> @param[in] write_moisture_diag Flag to control output of moisture - !! conservation diagnostics !> @param[in] first_iteration Flag for first inner iteration subroutine semi_implicit_solver_alg_step(state, rhs, & moist_dyn_gas_law, & - mr, write_moisture_diag, & - first_iteration) + mr, first_iteration) use solver_constants_mod, only: get_normalisation_r_solver, & get_im3_div_r_solver, & get_normalisation use fs_continuity_mod, only: W2, Wtheta use mr_indices_mod, only: nummr - use moisture_conservation_alg_mod, only: moisture_conservation_alg use si_operators_alg_mod, only: get_m3_rho_star, & get_rho_at_u, & get_p2theta, & @@ -534,7 +530,6 @@ contains type( field_type ), dimension(bundle_size), intent( inout ) :: rhs type( field_type ), intent( in ) :: moist_dyn_gas_law type( field_type ), dimension(nummr), intent( in ) :: mr - logical( kind=l_def ), intent( in ) :: write_moisture_diag logical( kind=l_def), intent( in ) :: first_iteration real( kind=r_def ) :: si_err(bundle_size) @@ -546,8 +541,6 @@ contains type( field_type ) :: rhs_rdef, & inc_theta - type( field_type ) :: rho_guess - type( r_solver_field_type ), dimension(bundle_size) :: rhs_rsol type( r_solver_field_type ) :: div_u, & rhs_tmp, & @@ -580,15 +573,6 @@ contains ! Input fields are r_def fields so preliminary work uses field_types mesh_id = state(igh_p)%get_mesh_id() - if ( write_moisture_diag ) then - ! Best guess for the prognostic fields state can be calculated as - ! state_best_guess = rhs_n - rhs_np1 + state + rhs_adv + rhs_phys - ! where fortunately all of the above except "state" is already in "rhs_np1" - call rho_guess%initialise( state(igh_d)%get_function_space() ) - call invoke(X_plus_Y(rho_guess, state(igh_d), rhs(igh_d))) - call moisture_conservation_alg( rho_guess, mr, 'Before solve' ) - end if - ! Normalise theta & u residual ! @TODO #416: can these be at r_solver precision? t_normalisation => get_normalisation(Wtheta, mesh_id) @@ -740,9 +724,6 @@ contains call invoke_inc_rdefX_plus_rsolverY( state(igh_t), inc_theta_rsol ) call invoke_inc_rdefX_plus_rsolverY( state(igh_d), inc_rho_rsol ) - if ( write_moisture_diag ) & - call moisture_conservation_alg( state(igh_d), mr, 'After solve' ) - if ( LPROF ) call stop_timing( id_si, 'dynamics.solver' ) end subroutine semi_implicit_solver_alg_step diff --git a/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 b/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 index eb7dc2d09..a402adc99 100644 --- a/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 +++ b/science/gungho/source/algorithm/timestepping/semi_implicit_timestep_alg_mod.X90 @@ -111,6 +111,7 @@ module semi_implicit_timestep_alg_mod use sci_mass_matrix_solver_alg_mod, & only: mass_matrix_solver_alg use moist_dyn_factors_alg_mod, only: moist_dyn_factors_alg + use moisture_conservation_alg_mod, only: moisture_conservation_alg use update_prognostic_scalars_alg_mod, & only: update_prognostic_scalars_alg use mixing_alg_mod, only: mixing_alg @@ -617,6 +618,7 @@ contains ! Moisture field to transport type( field_type ), pointer :: mr_to_adv(:) + type( field_type ) :: rho_after_transport type( field_type ) :: dcfl_tot, dcff_tot, dbcf_tot, & dcfl_adv, dcff_adv, dbcf_adv @@ -1012,15 +1014,28 @@ contains outer == outer_iterations .and. & inner == inner_iterations .and. & self%use_moisture + if (write_moisture_diag) then + ! Calculate transported density + call self%state(igh_d)%copy_field_properties(rho_after_transport) + call invoke(X_plus_Y(rho_after_transport, & + self%state_after_slow(igh_d), & + self%rhs_adv(igh_d))) + call moisture_conservation_alg(rho_after_transport, mr, & + 'Before solve') + end if !-------------------------------------------------------------------- ! Solve semi-implicit system: A*inc = rhs, and incement state by inc !-------------------------------------------------------------------- call semi_implicit_solver_alg_step( self%state, self%rhs_np1, & moist_dyn(gas_law), & mr, & - write_moisture_diag, & first_iteration=(inner==1) ) + if (write_moisture_diag) then + call moisture_conservation_alg(self%state(igh_d), mr, & + 'After solve') + end if + ! If not already done update factors for moist dynamics if ( .not. guess_np1 .and. self%use_moisture ) & call moist_dyn_factors_alg(moist_dyn, mr) diff --git a/science/linear/source/algorithm/timestepping/tl_si_timestep_alg_mod.x90 b/science/linear/source/algorithm/timestepping/tl_si_timestep_alg_mod.x90 index d912604f0..67bdf8c8c 100644 --- a/science/linear/source/algorithm/timestepping/tl_si_timestep_alg_mod.x90 +++ b/science/linear/source/algorithm/timestepping/tl_si_timestep_alg_mod.x90 @@ -396,9 +396,6 @@ contains ! these may differ from input ! values during the spinup period - ! Don't write moisture diagnostics in linear app - logical( kind=l_def ) :: write_moisture_diag = .false. - ! Configuration type( namelist_type ), pointer :: mixed_solver_nml real( kind=r_def ) :: mixed_solver_a_tol @@ -643,7 +640,7 @@ contains ls_rhs_np1(:, outer, inner), & ls_moist_dyn_itns(gas_law, next_outer, next_inner), & ls_mr_itns(:, next_outer, next_inner), & - write_moisture_diag, first_iteration=(inner==1) ) + first_iteration=(inner==1) ) ! If not already done update factors for moist dynamics if ( .not. guess_np1 .and. use_moisture) then call moist_dyn_factors_alg( & @@ -858,7 +855,6 @@ contains call semi_implicit_solver_alg_step( state, rhs_np1, & moist_dyn(gas_law), & mr, & - write_moisture_diag, & first_iteration=(inner==1) ) end if