From 0fb29a0e84e81d8650cdbaafee14a2cc3385b4e8 Mon Sep 17 00:00:00 2001 From: Thomas Bendall <14180399+tommbendall@users.noreply.github.com> Date: Mon, 24 Aug 2026 09:33:37 +0100 Subject: [PATCH 1/6] bring over changes to implement substepping for tracers --- .../gungho_transport_control_alg_mod.X90 | 3 +- .../control/transport_controller_mod.x90 | 69 ++++++++++++++++--- 2 files changed, 60 insertions(+), 12 deletions(-) diff --git a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 index f799cef47f..4656203b1c 100644 --- a/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 +++ b/science/gungho/source/algorithm/transport/control/gungho_transport_control_alg_mod.X90 @@ -244,6 +244,7 @@ contains ! Internal variables logical(kind=l_def) :: do_moisture_diagnostics logical(kind=l_def) :: cheap_update_step + logical(kind=l_def) :: safe_ref ! Temporary fields or pointers type(field_type) :: fields_np1(bundle_size) @@ -346,7 +347,7 @@ contains ! Check negative reference fields at this point if (check_any_eqn_consistent()) then - call transport_controller%check_reference_field() + safe_ref = transport_controller%check_reference_field(fail_on_neg=.true.) end if ! ------------------------------------------------------------------------ ! diff --git a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 index fe2eaeeb55..6965ba45d9 100644 --- a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 +++ b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 @@ -557,7 +557,10 @@ contains !> @brief Calculates the reference fields to be used in tracer transport, and !! if any of these are going to be negative, may perform the Watkins !! algorithm to adjust how the fluxes are distributed between steps. - subroutine check_reference_field(self) + !> @param[in] fail_on_neg If true, the model will stop if a reference field + !! is negative. If false, the model will continue. + !> @return True if the ref field is safe to use, false if not + function check_reference_field(self, fail_on_neg) result(ref_field_safe) use transport_config_mod, only: ffsl_unity_3d, & adjust_vhv_wind, & @@ -572,6 +575,7 @@ contains implicit none class(transport_controller_type), target, intent(inout) :: self + logical(kind=l_def), intent(in) :: fail_on_neg type(flux_precomputations_type), pointer :: flux_precomputations @@ -688,22 +692,34 @@ contains call self%flux_precomputations%initialise_step(step=(3+3*(i-1)), ref_flux=step_3_flux) end do + ref_field_safe = .true. else ! We're about to fail... log the Lipschitz numbers + ref_field_safe = .false. + call lipschitz_diag_alg( & self%config, & self%reference_splitting, self%transporting_wind, & self%dt_substep, LOG_LEVEL_WARNING & ) - if (adjust_tracer_equation) then + if (.not. fail_on_neg .and. self%num_substeps < 3) then + ! Time to attempt substepping... + call log_event( & + 'The Watkins algorithm failed to keep all reference fields ' // & + 'positive for the transport of tracers. Now try substepping...', & + LOG_LEVEL_WARNING & + ) + + else if (adjust_tracer_equation) then ! Use advective form for tracer transport instead of consistent form call log_event( & 'Reference field for tracer transport is negative, so ' // & 'making all tracers advective for this step. This will ' // & - 'make tracer transport non-conservative. ', LOG_LEVEL_INFO & + 'make tracer transport non-conservative. ', LOG_LEVEL_WARNING & ) self%make_tracers_advective = .true. + else ! Fail gracefully call log_event( & @@ -734,7 +750,7 @@ contains end if end if - end subroutine check_reference_field + end function check_reference_field ! ============================================================================ ! ! SETTERS @@ -1027,6 +1043,8 @@ contains type(r_tran_field_type), intent(in) :: ref_flux type(mesh_type), pointer :: mesh + logical(kind=l_def) :: safe_ref + integer(kind=i_def) :: i, start_substeps if (self%flux_precomputations%is_initialised()) then call log_event( & @@ -1043,13 +1061,42 @@ contains ! coarser mesh. In either case, the fluxes calculated during the transport ! of the density are not the same as those that will be used to transport ! the tracers. We set the correct flux here. - call self%flux_precomputations%initialise( & - self%config, & - mesh, self%reference_splitting, self%num_substeps, & - self%transporting_wind, self%dt_substep, & - flux=ref_flux, ref_field=self%ref_field_rtran & - ) - call self%check_reference_field() + + ! Loop through substeps, as we may need to increase the number to keep + ! the reference fields positive + safe_ref = .false. + start_substeps = self%num_substeps + do i = start_substeps, 3 ! Force max substeps of 3 + call self%flux_precomputations%initialise( & + self%config, & + mesh, self%reference_splitting, self%num_substeps, & + self%transporting_wind, self%dt_substep, & + flux=ref_flux, ref_field=self%ref_field_rtran & + ) + safe_ref = self%check_reference_field(fail_on_neg=.false.) + if (.not. safe_ref) then + ! Watkins algorithm failed to keep reference field positive + ! This is likely to be BDF tracer transport near start of run + ! Try enforcing substepping for tracer transport + self%dt_substep = ( & + self%dt_substep * self%num_substeps & + / real(self%num_substeps + 1, r_tran) & + ) + self%num_substeps = self%num_substeps + 1 + + ! Clear all precomputations in case they were set + call self%ffsl_precomputations%finalise() + call self%flux_precomputations%finalise() + call self%wind_precomputations%finalise() + + write(log_scratch_space, '(A,I4)') & + 'Transport: Reference field for tracer transport cannot stay ' // & + 'positive. Increasing number of substeps to ', self%num_substeps + call log_event(log_scratch_space, LOG_LEVEL_WARNING) + else + EXIT + end if + end do end subroutine initialise_flux_precomputations From 8143cc1cc469aaa5e73a985cd02e345e9027e8ab Mon Sep 17 00:00:00 2001 From: Thomas Bendall <14180399+tommbendall@users.noreply.github.com> Date: Tue, 25 Aug 2026 09:16:18 +0100 Subject: [PATCH 2/6] only substep tracers if namelist says so! --- .../transport/control/transport_controller_mod.x90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 index 6965ba45d9..75015144e8 100644 --- a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 +++ b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 @@ -1044,7 +1044,7 @@ contains type(mesh_type), pointer :: mesh logical(kind=l_def) :: safe_ref - integer(kind=i_def) :: i, start_substeps + integer(kind=i_def) :: i, start_substeps, end_substeps if (self%flux_precomputations%is_initialised()) then call log_event( & @@ -1066,7 +1066,13 @@ contains ! the reference fields positive safe_ref = .false. start_substeps = self%num_substeps - do i = start_substeps, 3 ! Force max substeps of 3 + if (substep_transport == substep_transport_adaptive) then + end_substeps = 3 ! Force max substeps of 3 + else + end_substeps = 1 ! No substepping + end if + + do i = start_substeps, end_substeps call self%flux_precomputations%initialise( & self%config, & mesh, self%reference_splitting, self%num_substeps, & From 351c5693f9fd9fb6b04a1a97797ddb72c840c1ca Mon Sep 17 00:00:00 2001 From: Thomas Bendall <14180399+tommbendall@users.noreply.github.com> Date: Tue, 25 Aug 2026 10:08:12 +0100 Subject: [PATCH 3/6] tighten up substepping criteria --- .../control/transport_controller_mod.x90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 index 75015144e8..74515aa3de 100644 --- a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 +++ b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 @@ -703,7 +703,9 @@ contains self%dt_substep, LOG_LEVEL_WARNING & ) - if (.not. fail_on_neg .and. self%num_substeps < 3) then + if (.not. fail_on_neg & + .and. substep_transport == substep_transport_adaptive & + .and. self%num_substeps < 3) then ! Time to attempt substepping... call log_event( & 'The Watkins algorithm failed to keep all reference fields ' // & @@ -1037,6 +1039,9 @@ contains !> @param[in] ref_flux The reference flux field subroutine initialise_flux_precomputations(self, ref_flux) + use transport_config_mod, only: substep_transport, & + substep_transport_adaptive + implicit none class(transport_controller_type), target, intent(inout) :: self @@ -1079,7 +1084,13 @@ contains self%transporting_wind, self%dt_substep, & flux=ref_flux, ref_field=self%ref_field_rtran & ) - safe_ref = self%check_reference_field(fail_on_neg=.false.) + if (i == end_substeps) then + ! Fail if we have reached max substeps + safe_ref = self%check_reference_field(fail_on_neg=.true.) + else + safe_ref = self%check_reference_field(fail_on_neg=.false.) + end if + if (.not. safe_ref) then ! Watkins algorithm failed to keep reference field positive ! This is likely to be BDF tracer transport near start of run From eba60b0add18564b66dc5471d8eac8eaed948fb6 Mon Sep 17 00:00:00 2001 From: Thomas Bendall <14180399+tommbendall@users.noreply.github.com> Date: Thu, 24 Sep 2026 20:23:29 +0100 Subject: [PATCH 4/6] fixes: ensure advective transport kicks in correctly, clear shared stencil extent --- .../algorithm/transport/control/transport_controller_mod.x90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 index 74515aa3de..91d158304e 100644 --- a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 +++ b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 @@ -721,6 +721,7 @@ contains 'make tracer transport non-conservative. ', LOG_LEVEL_WARNING & ) self%make_tracers_advective = .true. + ref_field_safe = .true. else ! Fail gracefully @@ -1105,6 +1106,7 @@ contains call self%ffsl_precomputations%finalise() call self%flux_precomputations%finalise() call self%wind_precomputations%finalise() + self%dep_stencil_extent_computed = .false. write(log_scratch_space, '(A,I4)') & 'Transport: Reference field for tracer transport cannot stay ' // & From ba71df5efbf94fd75a5ad8b8aad4e37c97b0566a Mon Sep 17 00:00:00 2001 From: Thomas Bendall <14180399+tommbendall@users.noreply.github.com> Date: Thu, 24 Sep 2026 20:33:14 +0100 Subject: [PATCH 5/6] fixes: use correct dt for share stencil extent, take substepping into account for Watkins --- .../control/transport_controller_mod.x90 | 10 +++-- .../transport/common/watkins_kernel_mod.F90 | 42 ++++++++++++++----- .../common/watkins_kernel_mod_test.pf | 2 + 3 files changed, 41 insertions(+), 13 deletions(-) diff --git a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 index 91d158304e..f0d6f6b7ec 100644 --- a/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 +++ b/science/gungho/source/algorithm/transport/control/transport_controller_mod.x90 @@ -396,7 +396,8 @@ contains watkins_fail_field, & wind_flux, & dt_one, & - detj_at_w3), & + detj_at_w3, & + 1), & ! Determine the adjusted wind for the final step X_minus_Y(step_3_wind_flux, vert_wind_flux, step_1_wind_flux) ) @@ -672,7 +673,8 @@ contains watkins_fail_field, & ref_flux_one_step, & dt_one, & - ref_mass), & + ref_mass, & + self%num_substeps), & ! Determine the adjusted wind for the final step X_minus_Y(step_3_flux, vert_flux, step_1_flux) ) @@ -883,14 +885,16 @@ contains type(flux_precomputations_type), pointer :: ffsl_precomputations type(mesh_type), pointer :: mesh type(r_tran_field_type) :: wind_flux + real(kind=r_tran) :: dt_whole_step ! If object is not created, then create it here if (.not. self%ffsl_precomputations%is_initialised()) then mesh => self%transporting_wind%get_mesh() ! Wind needs multiplying by dt + dt_whole_step = self%dt_substep * real(self%num_substeps, r_tran) call self%transporting_wind%copy_field_properties(wind_flux) - call invoke( a_times_X(wind_flux, self%dt_substep, self%transporting_wind) ) + call invoke( a_times_X(wind_flux, dt_whole_step, self%transporting_wind) ) ! Since the object has not already been created, we initialise it with the ! 3D wind, rather than a wind flux for each split step diff --git a/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 b/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 index 2aacc73f3c..a2c95a9494 100644 --- a/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/common/watkins_kernel_mod.F90 @@ -31,13 +31,14 @@ module watkins_kernel_mod !> The type declaration for the kernel. Contains the metadata needed by the Psy layer type, public, extends(kernel_type) :: watkins_kernel_type private - type(arg_type) :: meta_args(6) = (/ & + type(arg_type) :: meta_args(7) = (/ & arg_type(GH_FIELD, GH_REAL, GH_WRITE, W2v), & ! first_v_wind arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & ! lipschitz_max_field arg_type(GH_FIELD, GH_INTEGER, GH_WRITE, ANY_DISCONTINUOUS_SPACE_2), & ! watkins_failures arg_type(GH_FIELD, GH_REAL, GH_READ, W2), & ! wind arg_type(GH_SCALAR, GH_REAL, GH_READ), & ! dt - arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & ! detj + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & ! detj + arg_type(GH_SCALAR, GH_INTEGER, GH_READ) & ! num_substeps /) integer :: operates_on = CELL_COLUMN contains @@ -60,6 +61,7 @@ module watkins_kernel_mod !> @param[in] wind 3D wind !> @param[in] dt Transport time step !> @param[in] detj Det(J) at W3: the volume of cells +!> @param[in] num_substeps Maximum number of transport substeps !> @param[in] ndf_w2v Num of DoFs per cell for W2V !> @param[in] undf_w2v Num of W2V DoFs for this partition !> @param[in] map_w2v Map of lowest-cell W2V DoFs @@ -79,6 +81,7 @@ subroutine watkins_code( nlayers, & wind, & dt, & detj, & + num_substeps, & ndf_w2v, & undf_w2v, & map_w2v, & @@ -102,6 +105,7 @@ subroutine watkins_code( nlayers, & integer(kind=i_def), intent(in) :: map_w2v(ndf_w2v) integer(kind=i_def), intent(in) :: map_w3(ndf_w3) integer(kind=i_def), intent(in) :: map_w3_2d(ndf_w3_2d) + integer(kind=i_def), intent(in) :: num_substeps real(kind=r_tran), intent(in) :: dt real(kind=r_tran), intent(in) :: detj(undf_w3) real(kind=r_tran), intent(in) :: wind(undf_w2) @@ -115,12 +119,14 @@ subroutine watkins_code( nlayers, & real(kind=r_tran) :: lip_x_kp1, lip_y_kp1, lip_hori_kp1, lip_z_half_kp1 real(kind=r_tran) :: lip_lim_km1, lip_lim, lip_lim_kp1, slack_km1, remainder real(kind=r_tran) :: max_lip_init, max_lip_new + real(kind=r_tran) :: threshold, threshold_kp1 + real(kind=r_tran) :: lip_3d_kp1, lip_3d, lip_3d_km1 ! The threshold sets the maximum allowed Lipschitz number under this scheme, ! while the tolerance is used for checking whether the algorithm has failed, ! and this is slacker than machine precision to avoid machine precision errors ! being detected as failures. - real(kind=r_tran), parameter :: threshold = 0.9_r_tran + real(kind=r_tran), parameter :: max_threshold = 0.9_r_tran real(kind=r_tran), parameter :: tolerance = 1.0E-3_r_tran ! Set failures to zero @@ -154,8 +160,10 @@ subroutine watkins_code( nlayers, & lip_z_half_kp1 = dt*(first_v_wind(map_w2v(2)) - first_v_wind(map_w2v(1)))/detj(map_w3(1)) lip_hori_kp1 = MAX(lip_x_kp1, lip_y_kp1, lip_x_kp1 + lip_y_kp1) lip_lim_kp1 = MAX(lip_z_half_kp1, lip_z_half_kp1 + lip_hori_kp1) - max_lip_init = lip_lim_kp1 + lip_3d_kp1 = MAX(0.0_r_tran, lip_x_kp1 + lip_y_kp1 + 2.0_r_tran*lip_z_half_kp1) + max_lip_init = lip_lim_kp1 + REAL(num_substeps - 1, r_tran)*lip_3d_kp1 max_lip_new = 0.0_r_tran + threshold_kp1 = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d_kp1 do k = 0, nlayers - 1 @@ -167,6 +175,8 @@ subroutine watkins_code( nlayers, & if (k > 0) then lip_hori_km1 = lip_hori + lip_3d_km1 = lip_3d + threshold = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d_km1 lip_z_half_km1 = dt*(first_v_wind(map_w2v(2)+k-1) - first_v_wind(map_w2v(1)+k-1))/detj(map_w3(1)+k-1) lip_lim_km1 = MAX(lip_z_half_km1, lip_z_half_km1 + lip_hori_km1) @@ -176,8 +186,12 @@ subroutine watkins_code( nlayers, & ! Lipschitz numbers for this cell ------------------------------------------ lip_hori = lip_hori_kp1 + lip_3d = lip_3d_kp1 lip_z_half = dt*(first_v_wind(map_w2v(2)+k) - first_v_wind(map_w2v(1)+k))/detj(map_w3(1)+k) lip_lim = MAX(lip_z_half, lip_z_half + lip_hori) + ! Anticipate subsequent substeps, so build 3D Lipschitz number into threshold + ! Each substep will add lip_3d to the Lipschitz number + threshold = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d ! Lipschitz numbers for cell above ----------------------------------------- if (k < nlayers - 1) then @@ -186,7 +200,14 @@ subroutine watkins_code( nlayers, & lip_z_half_kp1 = dt*(first_v_wind(map_w2v(2)+k+1) - first_v_wind(map_w2v(1)+k+1))/detj(map_w3(1)+k+1) lip_hori_kp1 = MAX(lip_x_kp1, lip_y_kp1, lip_x_kp1 + lip_y_kp1) lip_lim_kp1 = MAX(lip_z_half_kp1, lip_z_half_kp1 + lip_hori_kp1) - max_lip_init = MAX(max_lip_init, lip_lim_kp1) + lip_3d_kp1 = MAX(0.0_r_tran, & + lip_x_kp1 + lip_y_kp1 & + + dt*(wind(map_w2(T)+k+1) - wind(map_w2(B)+k+1))/detj(map_w3(1)+k+1) & + ) + ! For tracking if max Lipschitz number has increased + max_lip_init = MAX(max_lip_init, & + lip_lim_kp1 + REAL(num_substeps-1, r_tran)*lip_3d_kp1) + threshold_kp1 = max_threshold - REAL(num_substeps - 1, r_tran)*lip_3d_kp1 end if ! ------------------------------------------------------------------------ ! @@ -201,7 +222,7 @@ subroutine watkins_code( nlayers, & if (slack_km1 > (lip_lim - threshold) * detj(map_w3(1)+k) / dt) then ! Check whether cell above might need some of this slack too - if (lip_lim_kp1 > threshold .AND. k < nlayers - 1) then + if (lip_lim_kp1 > threshold_kp1 .AND. k < nlayers - 1) then ! TODO: don't need to use all slack here! Arbitrarily using half first_v_wind(map_w2v(1)+k) = first_v_wind(map_w2v(1)+k) & + (lip_lim - threshold) * detj(map_w3(1)+k) / dt & @@ -229,7 +250,8 @@ subroutine watkins_code( nlayers, & ! Calculate Lipschitz number from cell below (which now won't change) if (k > 0) then lip_z_half_km1 = dt*(first_v_wind(map_w2v(2)+k-1) - first_v_wind(map_w2v(1)+k-1))/detj(map_w3(1)+k-1) - lip_lim_km1 = MAX(lip_z_half_km1, lip_z_half_km1 + lip_hori_km1) + lip_lim_km1 = MAX(lip_z_half_km1, lip_z_half_km1 + lip_hori_km1) & + + REAL(num_substeps - 1, r_tran)*lip_3d_km1 max_lip_new = MAX(max_lip_new, lip_lim_km1) end if @@ -248,14 +270,14 @@ subroutine watkins_code( nlayers, & k = nlayers - 1 lip_z_half = dt*(first_v_wind(map_w2v(2)+k) - first_v_wind(map_w2v(1)+k))/detj(map_w3(1)+k) lip_lim = MAX(lip_z_half, lip_z_half + lip_hori) - max_lip_new = MAX(max_lip_new, lip_lim) + max_lip_new = MAX(max_lip_new, lip_lim + REAL(num_substeps - 1, r_tran)*lip_3d) ! -------------------------------------------------------------------------- ! ! Check if algorithm has failed, and if so take original wind ! -------------------------------------------------------------------------- ! if (max_lip_new > max_lip_init + tolerance .AND. & - max_lip_new > threshold + tolerance) then + max_lip_new > max_threshold + tolerance) then watkins_failures(map_w3_2d(1)) = 1_i_def ! Set the bottom value @@ -269,7 +291,7 @@ subroutine watkins_code( nlayers, & ! Set the top values first_v_wind(map_w2v(1)+nlayers) = 0.0_r_tran - else if (max_lip_new > threshold + tolerance) then + else if (max_lip_new > max_threshold + tolerance) then watkins_failures(map_w3_2d(1)) = 1_i_def end if diff --git a/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf b/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf index 88956d7366..296fc5f928 100644 --- a/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/transport/common/watkins_kernel_mod_test.pf @@ -160,6 +160,7 @@ contains wind, & dt, & detj_at_w3, & + 1, & ndf_w2v, & undf_w2v, & map_w2v, & @@ -261,6 +262,7 @@ contains wind, & dt, & detj_at_w3, & + 1, & ndf_w2v, & undf_w2v, & map_w2v, & From 5000f7f8f422150a0411350e1b69b99527f77de3 Mon Sep 17 00:00:00 2001 From: Thomas Bendall <14180399+tommbendall@users.noreply.github.com> Date: Mon, 28 Sep 2026 21:51:47 +0100 Subject: [PATCH 6/6] protect total_ref_flux so that it is correct when set --- .../common/flux_precomputations_mod.x90 | 27 ++++++++++++++++--- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 b/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 index eada1ca879..5de81674e7 100644 --- a/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 +++ b/science/gungho/source/algorithm/transport/common/flux_precomputations_mod.x90 @@ -107,8 +107,9 @@ module flux_precomputations_alg_mod integer(kind=i_def) :: max_stencil_extent logical(kind=l_def) :: max_stencil_extent_computed = .false. - ! Total reference flux (over whole transport step), needed for cheap update - type(r_tran_field_type) :: total_ref_flux + ! Total reference flux (over whole transport step), needed when there is + ! a separate tracer transport, e.g. cheap update, coarse tracers, TR-BDF2 + type(r_tran_field_type), allocatable :: total_ref_flux ! Wind and dt for whole transport step, needed for logging type(r_tran_field_type), pointer :: transporting_wind => null() @@ -294,6 +295,7 @@ contains self%to_initialise_by_step = present(flux) self%reset_unity = .false. if (present(reset_unity)) self%reset_unity = reset_unity + allocate(self%total_ref_flux) ! Determine settings for extended mesh remapping of density if (panel_edge_treatment == panel_edge_treatment_extended_mesh) then @@ -603,6 +605,7 @@ contains nullify(self%transporting_wind) self%is_initialised_flag = .false. self%max_stencil_extent_computed = .false. + if ( allocated( self%total_ref_flux ) ) deallocate( self%total_ref_flux ) if ( allocated( self%mesh_ids ) ) deallocate( self%mesh_ids ) if ( allocated( self%ref_field ) ) deallocate( self%ref_field ) if ( allocated( self%ref_mass ) ) deallocate( self%ref_mass ) @@ -1758,7 +1761,8 @@ contains !! reference field over the *whole* transport step. This is needed when !! using the "cheap update" transport formulation, as the total flux !! needs storing from previous outer loops for the transport of tracers - !! in the final outer loop. + !! in the final outer loop. It is also needed for coarse transport of + !! tracers, and the transport of tracers in TR-BDF2. !> @return The total reference flux function get_total_ref_flux(self) result(total_ref_flux) @@ -1768,6 +1772,13 @@ contains type(r_tran_field_type), pointer :: total_ref_flux character(len=str_def) :: field_name + if (.not. allocated(self%total_ref_flux)) then + call log_event( & + 'flux_precomp%get_total_ref_flux: self%total_ref_flux not allocated', & + LOG_LEVEL_ERROR & + ) + end if + total_ref_flux => self%total_ref_flux ! Check field has been initialised @@ -1780,7 +1791,8 @@ contains !! field over the *whole* transport step. This is needed when using the !! "cheap update" transport formulation, as the total flux needs !! storing from previous outer loops for the transport of tracers in - !! the final outer loop. + !! the final outer loop. It is also needed for coarse transport of + !! tracers, and the transport of tracers in TR-BDF2. !! As we may be substepping, the field is incremented rather than set !> @param[in] total_ref_flux The total reference flux to set subroutine set_total_ref_flux(self, total_ref_flux) @@ -1790,6 +1802,13 @@ contains class(flux_precomputations_type), target, intent(inout) :: self type(r_tran_field_type), intent(in) :: total_ref_flux + if (.not. allocated(self%total_ref_flux)) then + call log_event( & + 'flux_precomp%set_total_ref_flux: self%total_ref_flux not allocated', & + LOG_LEVEL_ERROR & + ) + end if + if (.not. self%total_ref_flux%is_initialised()) then call self%total_ref_flux%initialise(total_ref_flux%get_function_space()) call invoke( setval_c(self%total_ref_flux, 0.0_r_tran) )