Improve moisture conservation diagnostics

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 |
+| ericaneininger      | Erica Neininger    | Met Office                       | 2026-03-02 |
+| iboutle             | Ian Boutle         | Met Office                       | 2026-03-10 |

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 <Mx,Mx>

applications/lfric_atm/metadata/field_def_diags.xml

@@ -418,8 +418,8 @@
   <field id="microphysics__ls_rain" name="ls_rain" long_name="surface_microphysical_rainfall_rate" unit="kg m-2 s-1" domain_ref="face" />
   <field id="microphysics__ls_snow" name="ls_snow" long_name="surface_microphysical_snowfall_rate" unit="kg m-2 s-1" domain_ref="face" />
   <field id="microphysics__ls_graup" name="ls_graup" long_name="surface_microphysical_graupelfall_rate" unit="kg m-2 s-1" domain_ref="face" />
-  <field id="microphysics__total_fpr" name="total_fpr" long_name="surface_microphysical_total_frozen_precipitation_rate" unit="kg m-2 s-2" domain_ref="face"> microphysics__ls_snow + microphysics__ls_graup </field>
-  <field id="microphysics__ls_prec" name="ls_prec" long_name="surface_microphysical_precipitation_rate" unit="kg m-2 s-1" domain_ref="face"> microphysics__ls_rain + microphysics__ls_snow + microphysics__ls_graup </field>
+  <field id="microphysics__snow_no_graup" name="snow_no_graup" long_name="surface_microphysical_snowfall_rate_excluding_graupel" unit="kg m-2 s-2" domain_ref="face"> microphysics__ls_snow - microphysics__ls_graup </field>
+  <field id="microphysics__ls_prec" name="ls_prec" long_name="surface_microphysical_precipitation_rate" unit="kg m-2 s-1" domain_ref="face"> microphysics__ls_rain + microphysics__ls_snow </field>
   <field id="microphysics__ls_rain_acc" name="ls_rain_acc" long_name="surface_microphysical_rainfall_amount" unit="kg m-2" domain_ref="face"> $DT * microphysics__ls_rain </field>
   <field id="microphysics__ls_snow_acc" name="ls_snow_acc" long_name="surface_microphysical_snowfall_amount" unit="kg m-2" domain_ref="face"> $DT * microphysics__ls_snow </field>
   <field id="microphysics__ls_graup_acc" name="ls_graup_acc" long_name="surface_microphysical_graupelfall_amount" unit="kg m-2" domain_ref="face"> $DT * microphysics__ls_graup </field>

rose-stem/app/lfric_atm/opt/rose-app-cons_diag.conf

@@ -0,0 +1,2 @@
+[namelist:io]
+write_conservation_diag=.true.

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,

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 )

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

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' )

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' )
 

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
     !--------------------------------------------------------------------

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

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)

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