diff --git a/rose-stem/app/gungho_model/rose-app.conf b/rose-stem/app/gungho_model/rose-app.conf index f1cb5c754..a3915eb35 100644 --- a/rose-stem/app/gungho_model/rose-app.conf +++ b/rose-stem/app/gungho_model/rose-app.conf @@ -783,6 +783,8 @@ si_tolerance=1.0e-1 split_w=.true. [namelist:mixing] +!!conservative_diffusion=.true. +!!density_weighted=.true. !!leonard_kl=1.0 leonard_term=.false. !!method='blending' diff --git a/rose-stem/app/lfric_atm/rose-app.conf b/rose-stem/app/lfric_atm/rose-app.conf index bd907a551..5cfa112b8 100644 --- a/rose-stem/app/lfric_atm/rose-app.conf +++ b/rose-stem/app/lfric_atm/rose-app.conf @@ -806,6 +806,8 @@ si_tolerance=1.0e-1 split_w=.true. [namelist:mixing] +!!conservative_diffusion=.true. +!!density_weighted=.true. !!leonard_kl=2.0 leonard_term=.false. !!method='blend_1dbl_fa' diff --git a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf index 4380951b7..ea70f93a7 100644 --- a/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf +++ b/science/gungho/rose-meta/lfric-gungho/HEAD/rose-meta.conf @@ -3534,6 +3534,36 @@ help=Currently available schemes are: ns=namelist/Science/Dynamics/Mixing sort-key=Section-A05 +[namelist:mixing=conservative_diffusion] +compulsory=true +description=Use conservative form of tracer diffusion operator +help=When true, the diffusion operator for tracers is implemented in + =conservative form, i.e. the fluxes are computed and then the divergence of + =the fluxes is taken to obtain the tendency. + = + =When false, a non-conservative form is used, i.e. the Laplacian of the + =tracer field is computed directly to obtain the tendency. +!kind=default +ns=namelist/Science/Dynamics/Mixing +sort-key=Panel-A08 +trigger=namelist:mixing=density_weighted: .true.; +type=logical + +[namelist:mixing=density_weighted] +compulsory=true +description=Whether to include the dry density in the tracer diffusion operator +help=When true, the conservative diffusion operator for tracers includes the + =density: i.e. the operator is (1/rho)*div(rho*nu*grad(tracer)), where + =nu is the diffusivity. + = + =When false, the operator is then simply div(nu*grad(tracer)). + = + =If conservative_diffusion is false, then this setting has no effect. +!kind=default +ns=namelist/Science/Dynamics/Mixing +sort-key=Panel-A09 +type=logical + [namelist:mixing=leonard_kl] compulsory=true description=Leonard term parameter @@ -3563,6 +3593,14 @@ sort-key=Panel-A06 trigger=namelist:mixing=leonard_kl: .true. ; type=logical +[namelist:mixing=max_diff_factor] +compulsory=true +description=Maximum factor for capping Smagorinsky diffusivity coefficient +help=Typical value is 1.0 +!kind=default +sort-key=Panel-A10 +type=real + [namelist:mixing=method] compulsory=true description=Smagorinsky subgrid mixing scheme option @@ -3625,6 +3663,8 @@ sort-key=Panel-A03 trigger=namelist:mixing=method: .true. ; =namelist:mixing=mix_factor: .true. ; =namelist:mixing=smag_l_calc: .true. ; + =namelist:mixing=max_diff_factor: .true.; + =namelist:mixing=conservative_diffusion: .true.; =namelist:physics=smagorinsky_placement: .true. ; type=logical diff --git a/science/gungho/rose-meta/lfric-gungho/versions.py b/science/gungho/rose-meta/lfric-gungho/versions.py index 01798ad2b..6f8337d1c 100644 --- a/science/gungho/rose-meta/lfric-gungho/versions.py +++ b/science/gungho/rose-meta/lfric-gungho/versions.py @@ -31,3 +31,21 @@ def upgrade(self, config, meta_config=None): # Add settings return config, self.reports """ + +class vn31_t378(MacroUpgrade): + """Upgrade macro for ticket #378 by Thomas Bendall.""" + + BEFORE_TAG = "vn3.1" + AFTER_TAG = "vn3.1_t378" + + def upgrade(self, config, meta_config=None): + self.add_setting( + config, ["namelist:mixing", "conservative_diffusion"], ".false." + ) + self.add_setting( + config, ["namelist:mixing", "density_weighted"], ".true." + ) + self.add_setting( + config, ["namelist:mixing", "max_diff_factor"], "1.0" + ) + return config, self.reports diff --git a/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 b/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 index e24f3487d..6de7abcdd 100644 --- a/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 +++ b/science/gungho/source/algorithm/diffusion/mixing_alg_mod.x90 @@ -22,8 +22,9 @@ contains !> @param[in] derived_fields Group of derived fields !> @param[in] rho Density !> @param[in] dt The model timestep length + !> @param[in] model_clock The model clock subroutine mixing_alg(mr, theta, u, derived_fields, & - rho, dt) + rho, dt, model_clock) use constants_mod, only: i_def,r_def use log_mod, only: log_event, & @@ -39,6 +40,7 @@ contains use momentum_viscosity_kernel_mod, only: momentum_viscosity_kernel_type use physics_config_mod, only: smagorinsky_placement, & smagorinsky_placement_end + use model_clock_mod, only: model_clock_type implicit none @@ -51,6 +53,7 @@ contains ! Timestepping information real( r_def ), intent(in) :: dt + class( model_clock_type ), intent(in) :: model_clock ! Increment fields type( field_type ) :: du, dtheta @@ -68,7 +71,7 @@ contains call du%initialise( u%get_function_space() ) call invoke( setval_c(dtheta, 0.0_r_def), & setval_c(du, 0.0_r_def) ) - call smagorinsky_alg(dtheta, du, mr, theta, u, derived_fields, rho, dt) + call smagorinsky_alg(dtheta, du, mr, theta, u, derived_fields, rho, dt, model_clock) call invoke( inc_X_plus_Y( theta, dtheta ), & inc_X_plus_Y( u, du ) ) diff --git a/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 b/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 index c66bf5883..5ea5f14f6 100644 --- a/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 +++ b/science/gungho/source/algorithm/diffusion/smagorinsky_alg_mod.x90 @@ -4,42 +4,67 @@ ! under which the code may be used. !------------------------------------------------------------------------------- -!>@brief Contains code for Smagorinsky diffusion calculation +!> @brief Contains code for Smagorinsky diffusion calculation module smagorinsky_alg_mod - use constants_mod, only: i_def, r_def - use field_mod, only: field_type - use field_collection_mod, only: field_collection_type - use mr_indices_mod, only: nummr, imr_v, imr_cl - - use formulation_config_mod, only: moisture_formulation, & - moisture_formulation_dry - use section_choice_config_mod, only: boundary_layer, & - boundary_layer_um - use mixing_config_mod, only: mix_factor, smag_l_calc, & - smag_l_calc_use_geo - - use tracer_smagorinsky_diff_kernel_mod, only: tracer_smagorinsky_diff_kernel_type - use momentum_smagorinsky_kernel_mod, only: momentum_smagorinsky_kernel_type - - use log_mod, only: log_event, & - LOG_LEVEL_INFO + ! Infrastructure + use constants_mod, only: i_def, r_def, l_def + use extrusion_mod, only: SHIFTED + use field_mod, only: field_type + use field_collection_mod, only: field_collection_type + use fs_continuity_mod, only: W1, W2, W2H, W3 + use function_space_mod, only: function_space_type + use function_space_collection_mod, & + only: function_space_collection + use log_mod, only: log_event, LOG_LEVEL_DEBUG use mesh_mod, only: mesh_type - use io_config_mod, only: write_conservation_diag + use mesh_collection_mod, only: mesh_collection + use model_clock_mod, only: model_clock_type + use mr_indices_mod, only: nummr, imr_v, imr_cl use timing_mod, only: start_timing, stop_timing, tik, LPROF - use moisture_conservation_alg_mod, & - only: moisture_conservation_alg - use sci_geometric_constants_mod, & - only: get_height_fv, get_panel_id, & - get_delta_at_wtheta, get_dx_at_w2, & + + ! Pointers to objects + use sci_fem_constants_mod, only: get_rmultiplicity_fv + use sci_geometric_constants_mod, & + only: get_height_fv, get_panel_id, & + get_delta_at_wtheta, get_dx_at_w2, & + get_dA_at_w2, get_detj_at_w3_fv, & get_dz_at_wtheta - use fs_continuity_mod, only: W1, W2 + use physics_constants_mod, only: get_max_diff + + ! Algorithms + use moisture_conservation_alg_mod, & + only: moisture_conservation_alg + + ! Kernels + use cap_diffusion_coeff_kernel_mod, & + only: cap_diffusion_coeff_kernel_type + use diffusion_flux_kernel_mod, only: diffusion_flux_kernel_type + use fv_divergence_2d_kernel_mod, & + only: fv_divergence_2d_kernel_type + use momentum_smagorinsky_kernel_mod, & + only: momentum_smagorinsky_kernel_type + use tracer_smagorinsky_diff_kernel_mod, & + only: tracer_smagorinsky_diff_kernel_type + use sci_w3_to_w2_average_kernel_mod, & + only: w3_to_w2_average_kernel_type + + ! Config + use formulation_config_mod, only: moisture_formulation, & + moisture_formulation_dry + use io_config_mod, only: subroutine_timers, & + write_conservation_diag + use mixing_config_mod, only: mix_factor, conservative_diffusion, & + smag_l_calc, smag_l_calc_use_geo, & + density_weighted + use section_choice_config_mod, only: boundary_layer, boundary_layer_um implicit none private public :: smagorinsky_alg + public :: smagorinsky_tracer_alg contains @@ -55,53 +80,64 @@ contains !> @param[in] derived_fields Group of derived fields !> @param[in] rho Density !> @param[in] dt The model timestep length -subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, & - derived_fields, rho, dt ) +subroutine smagorinsky_alg( & + dtheta_io, du_io, mr, theta, u, derived_fields, rho, dt, model_clock & +) implicit none ! Prognostic fields - type( field_type ), intent(inout) :: du_io, dtheta_io - type( field_type ), intent(inout) :: mr ( nummr ) - - type( field_collection_type ), intent(in) :: derived_fields - - type( field_type ), intent(in) :: rho, theta, u - real( r_def ), intent(in) :: dt + type(field_type), intent(inout) :: du_io, dtheta_io + type(field_type), intent(inout) :: mr (nummr) + type(field_collection_type), intent(in) :: derived_fields + type(field_type), intent(in) :: rho, theta, u + real(kind=r_def), intent(in) :: dt + class(model_clock_type), intent(in) :: model_clock ! Increments on state fields - type( field_type ) :: du, dtheta - - type( field_type ), pointer :: shear => null() - type( field_type ), pointer :: visc_h - type( field_type ), pointer :: visc_m - type( field_type ), pointer :: height_w2 => null() - type( field_type ), pointer :: height_w1 => null() - type( field_type ), pointer :: delta => null() - type( field_type ), pointer :: dz_wth => null() - type( field_type ), pointer :: dx_at_w2 => null() - type( field_type ), pointer :: panel_id => null() - type( field_type ) :: dmr_v, dmr_cl - - type( mesh_type ), pointer :: mesh => null() + type(field_type) :: du, dtheta + type(field_type) :: dmr_v, dmr_cl + + integer(kind=i_def) :: shifted_mesh_id + type(field_type), pointer :: shear + type(field_type), pointer :: visc_h + type(field_type), pointer :: visc_m + type(field_type), pointer :: height_w2 + type(field_type), pointer :: height_w1 + type(field_type), pointer :: dz_wth + type(field_type), pointer :: rho_in_wth + type(field_type) :: rho_in_w2h + type(field_type), target :: delta_geo + type(field_type), pointer :: delta + type(field_type), pointer :: dx_at_w2 + type(field_type), pointer :: panel_id + type(field_type), pointer :: rmult + type(field_type), target :: visc_w2h + type(field_type), pointer :: visc_h_ptr + type(field_type), pointer :: max_diff + type(mesh_type), pointer :: mesh + type(mesh_type), pointer :: shifted_mesh + + type(function_space_type), pointer :: w2h_sh_fs + + type(field_type), target :: rho_one ! Stencil depth for the Smagorinsky diffusion kernel integer(kind=i_def), parameter :: smag_stencil_depth = 1 + real(kind=r_def), parameter :: one_third = 1.0_r_def / 3.0_r_def - logical :: use_moisture - real( r_def ) :: one_third - integer(tik) :: id + logical(kind=l_def) :: use_moisture + integer(kind=tik) :: id if ( LPROF ) call start_timing( id, 'smagorinsky' ) - call log_event( 'Applying Smagorinsky mixing', LOG_LEVEL_INFO ) + call log_event( 'Applying Smagorinsky mixing', LOG_LEVEL_DEBUG ) use_moisture = ( moisture_formulation /= moisture_formulation_dry ) mesh => theta%get_mesh() call du%initialise( u%get_function_space() ) - call invoke(setval_c(du,0.0_r_def)) call dtheta%initialise( theta%get_function_space() ) call derived_fields%get_field('visc_h', visc_h) call derived_fields%get_field('visc_m', visc_m) @@ -112,54 +148,28 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, & dx_at_w2 => get_dx_at_w2(mesh%get_id()) panel_id => get_panel_id(mesh%get_id()) - if ( boundary_layer == boundary_layer_um ) then + call invoke( setval_c(du, 0.0_r_def) ) + + ! ======================================================================== ! + ! Get mixing coefficients + ! ======================================================================== ! + + ! With UM BL scheme -------------------------------------------------------- + if (boundary_layer == boundary_layer_um) then ! Use stability-dependent diffusion coefficient from UM BL scheme ! This will be the blended BL-Smag coefficient if mixing_option='blending' ! Vertical mixing is done by the BL scheme - call log_event( 'Using stability-dependent Smagorinsky coefficient', LOG_LEVEL_INFO ) - - ! Potential temperature: - call invoke( tracer_smagorinsky_diff_kernel_type( dtheta, & - theta, & - smag_stencil_depth, & - visc_h, & - dx_at_w2 ) ) - - if ( use_moisture ) then - ! Water vapour and cloud liquid: - call dtheta%copy_field_properties(dmr_v) - call dtheta%copy_field_properties(dmr_cl) - call invoke( tracer_smagorinsky_diff_kernel_type( dmr_v, & - mr(imr_v), & - smag_stencil_depth, & - visc_h, & - dx_at_w2 ), & - tracer_smagorinsky_diff_kernel_type( dmr_cl, & - mr(imr_cl), & - smag_stencil_depth, & - visc_h, & - dx_at_w2 ) ) - end if + call log_event( & + 'Using stability-dependent Smagorinsky coefficient', LOG_LEVEL_DEBUG & + ) + + else ! Without UM BL scheme ------------------------------------------------ - ! Momentum: - call invoke(momentum_smagorinsky_kernel_type( du, & - u, & - smag_stencil_depth, & - dx_at_w2, & - smag_stencil_depth, & - height_w2, & - height_w1, & - visc_m, & - smag_stencil_depth, & - panel_id, & - smag_stencil_depth) ) - - else ! Without UM BL scheme - - call log_event( 'Using pure Smagorinsky coefficient', LOG_LEVEL_INFO ) + call log_event( 'Using pure Smagorinsky coefficient', LOG_LEVEL_DEBUG ) ! Calculate visc_h and visc_m as (mix_factor * delta)**2 * shear call derived_fields%get_field('shear', shear) + ! delta is calculated in sci_calc_delta_at_wtheta_kernel_mod.F90 ! in lfric core as the mininum of the cell horizontal edge lengths delta => get_delta_at_wtheta(mesh%get_id()) @@ -167,72 +177,206 @@ subroutine smagorinsky_alg(dtheta_io, du_io, mr, theta, u, & if (smag_l_calc == smag_l_calc_use_geo) then ! Use geometric mean mixing length dz_wth => get_dz_at_wtheta(mesh%get_id()) - - one_third = 1.0_r_def / 3.0_r_def - call invoke( inc_X_powint_n( delta, 2 ), & - inc_X_times_Y( delta, dz_wth ), & - inc_X_powreal_a( delta, one_third ) ) - end if - ! Potential temperature: - call invoke( a_times_X( visc_h, mix_factor, delta ), & - inc_X_powint_n( visc_h, 2 ), & - inc_X_times_Y( visc_h, shear ), & - tracer_smagorinsky_diff_kernel_type( dtheta, & - theta, & - smag_stencil_depth, & - visc_h, & - dx_at_w2 ) ) - - if (use_moisture) then - ! Water vapour and cloud liquid: - call dtheta%copy_field_properties(dmr_v) - call dtheta%copy_field_properties(dmr_cl) - call invoke( tracer_smagorinsky_diff_kernel_type( dmr_v, & - mr(imr_v), & - smag_stencil_depth, & - visc_h, & - dx_at_w2 ), & - tracer_smagorinsky_diff_kernel_type( dmr_cl, & - mr(imr_cl), & - smag_stencil_depth, & - visc_h, & - dx_at_w2 ) ) + call delta_geo%initialise(delta%get_function_space()) + + call invoke( & + setval_X(delta_geo, delta), & + inc_X_powint_n(delta_geo, 2), & + inc_X_times_Y(delta_geo, dz_wth), & + inc_X_powreal_a(delta_geo, one_third) & + ) + delta => delta_geo end if - ! Momentum: - call invoke( a_times_X( visc_m, mix_factor, delta ), & - inc_X_powint_n( visc_m, 2 ), & - inc_X_times_Y( visc_m, shear ), & - momentum_smagorinsky_kernel_type( du, & - u, & - smag_stencil_depth, & - dx_at_w2, & - smag_stencil_depth, & - height_w2, & - height_w1, & - visc_m, & - smag_stencil_depth, & - panel_id, & - smag_stencil_depth) ) + call invoke( & + a_times_X(visc_h, mix_factor, delta), & + inc_X_powint_n(visc_h, 2), & + inc_X_times_Y(visc_h, shear), & + a_times_X(visc_m, mix_factor, delta), & + inc_X_powint_n(visc_m, 2), & + inc_X_times_Y(visc_m, shear) & + ) end if ! With or without BL scheme + ! Cap viscosities to avoid instabilities + max_diff => get_max_diff(mesh%get_id(), model_clock) + call invoke( & + cap_diffusion_coeff_kernel_type(visc_h, max_diff), & + cap_diffusion_coeff_kernel_type(visc_m, max_diff) & + ) + + ! ======================================================================== ! + ! Convert fields to W2H, if needed + ! ======================================================================== ! + + if (conservative_diffusion) then + shifted_mesh => mesh_collection%get_mesh(mesh, SHIFTED) + shifted_mesh_id = shifted_mesh%get_id() + w2h_sh_fs => function_space_collection%get_fs(shifted_mesh, 0, 0, W2H) + rmult => get_rmultiplicity_fv(W2H, shifted_mesh_id) + + call visc_w2h%initialise(w2h_sh_fs) + call rho_in_w2h%initialise(w2h_sh_fs) + call derived_fields%get_field('rho_in_wth', rho_in_wth) + + ! Get density at W2H points + if (density_weighted) then + call invoke( & + setval_c(rho_in_w2h, 0.0_r_def), & + w3_to_w2_average_kernel_type(rho_in_w2h, rho_in_wth, rmult) & + ) + else + ! Set rho to 1 everywhere if not using density weighting + call rho_one%initialise(rho_in_wth%get_function_space()) + call invoke( & + setval_c(rho_one, 1.0_r_def), & + setval_c(rho_in_w2h, 1.0_r_def) & + ) + rho_in_wth => rho_one + end if + + ! Obtain viscosity at W2H points + call invoke( & + setval_c(visc_w2h, 0.0_r_def), & + w3_to_w2_average_kernel_type(visc_w2h, visc_h, rmult) & + ) + visc_h_ptr => visc_w2h + else + ! Directly use visc_h + visc_h_ptr => visc_h + end if + + ! ======================================================================== ! + ! Perform Smagorinsky diffusion on tracers and momentum + ! ======================================================================== ! + ! Potential temperature: + call smagorinsky_tracer_alg( & + dtheta, theta, visc_h_ptr, rho_in_w2h, rho_in_wth & + ) + if (use_moisture) then - call invoke( inc_X_plus_bY( mr(imr_v), dt, dmr_v ), & - inc_X_plus_bY( mr(imr_cl), dt, dmr_cl ) ) + ! Water vapour and cloud liquid: + call smagorinsky_tracer_alg( & + dmr_v, mr(imr_v), visc_h_ptr, rho_in_w2h, rho_in_wth & + ) + call smagorinsky_tracer_alg( & + dmr_cl, mr(imr_cl), visc_h_ptr, rho_in_w2h, rho_in_wth & + ) end if - ! Apply increments from Smagorinsky mixing - call invoke( inc_X_plus_bY( dtheta_io, dt, dtheta ), & - inc_X_plus_bY( du_io, dt, du ) ) + ! Momentum: + call invoke( & + momentum_smagorinsky_kernel_type( & + du, u, smag_stencil_depth, dx_at_w2, smag_stencil_depth, & + height_w2, height_w1, visc_m, smag_stencil_depth, & + panel_id, smag_stencil_depth & + ) & + ) + + ! ======================================================================== ! + ! Apply increments + ! ======================================================================== ! + if (use_moisture) then + call invoke( & + inc_X_plus_bY(mr(imr_v), dt, dmr_v), & + inc_X_plus_bY(mr(imr_cl), dt, dmr_cl) & + ) + end if - if (write_conservation_diag .and. use_moisture) & - call moisture_conservation_alg( rho, mr, 'After mixing' ) + call invoke( & + inc_X_plus_bY(dtheta_io, dt, dtheta), & + inc_X_plus_bY(du_io, dt, du) & + ) - nullify( height_w1, height_w2 ) + if (write_conservation_diag .and. use_moisture) then + call moisture_conservation_alg(rho, mr, 'After mixing') + end if if ( LPROF ) call stop_timing( id, 'smagorinsky' ) end subroutine smagorinsky_alg + !> @details Performs Smagorinsky diffusion on a Wtheta tracer field + !> @param[in,out] increment The diffusion increment for the tracer + !> @param[in] tracer The tracer field to be diffused + !> @param[in] visc_h The horizontal viscosity field + !> @param[in] rho_in_w2h The density field at W2H points. Unused if not + !! using conservative diffusion. + !> @param[in] rho_in_wth The density field at Wtheta points. Unused if + !! not using conservative diffusion. + subroutine smagorinsky_tracer_alg( & + increment, tracer, visc_h, rho_in_w2h, rho_in_wth & + ) + + implicit none + + ! Arguments + type(field_type), intent(inout) :: increment + type(field_type), intent(in) :: tracer + type(field_type), intent(in) :: visc_h + type(field_type), intent(in) :: rho_in_w2h + type(field_type), intent(in) :: rho_in_wth + + ! Internal variables + integer(kind=i_def) :: shifted_mesh_id + type(field_type) :: tracer_shifted + type(field_type) :: increment_shifted + type(field_type) :: flux + type(field_type), pointer :: dx_at_w2 + type(field_type), pointer :: dA, volume + type(field_type), pointer :: rmult + type(mesh_type), pointer :: mesh, shifted_mesh + type(function_space_type), pointer :: w3_sh_fs, w2h_sh_fs + + ! Stencil depth for the Smagorinsky diffusion kernel + integer(kind=i_def), parameter :: smag_stencil_depth = 1 + + ! Get mesh and required fields + mesh => tracer%get_mesh() + + call tracer%copy_field_properties(increment) + + if (conservative_diffusion) then + ! Get shifted mesh and function spaces + shifted_mesh => mesh_collection%get_mesh(mesh, SHIFTED) + shifted_mesh_id = shifted_mesh%get_id() + w3_sh_fs => function_space_collection%get_fs(shifted_mesh, 0, 0, W3) + w2h_sh_fs => function_space_collection%get_fs(shifted_mesh, 0, 0, W2H) + dA => get_dA_at_w2(shifted_mesh_id) + volume => get_detj_at_w3_fv(shifted_mesh_id) + rmult => get_rmultiplicity_fv(W2H, shifted_mesh_id) + dx_at_w2 => get_dx_at_w2(shifted_mesh_id) + + call tracer_shifted%initialise(w3_sh_fs) + call increment_shifted%initialise(w3_sh_fs) + call flux%initialise(w2h_sh_fs) + + call invoke( & + ! Map tracer to shifted mesh + setval_X(tracer_shifted, tracer), & + setval_c(flux, 0.0_r_def), & + ! Compute flux on shifted mesh + diffusion_flux_kernel_type( & + flux, tracer_shifted, visc_h, rho_in_w2h, dx_at_w2, dA, rmult & + ), & + ! Take divergence to get increment + fv_divergence_2d_kernel_type(increment_shifted, flux, volume), & + ! Map increment back from shifted mesh and normalise by rho + X_divideby_Y(increment, increment_shifted, rho_in_wth) & + ) + + else + dx_at_w2 => get_dx_at_w2(mesh%get_id()) + + ! Non-conservative Smagorinsky diffusion, using Laplacian directly + call invoke( & + tracer_smagorinsky_diff_kernel_type( & + increment, tracer, smag_stencil_depth, visc_h, dx_at_w2 & + ) & + ) + end if + + end subroutine smagorinsky_tracer_alg + end module smagorinsky_alg_mod 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 851921548..67df16781 100644 --- a/science/gungho/source/algorithm/physics/fast_physics_alg_mod.X90 +++ b/science/gungho/source/algorithm/physics/fast_physics_alg_mod.X90 @@ -392,7 +392,7 @@ contains call derived_fields%get_field('theta_star', theta_star) call derived_fields%get_field('u_star', u_star) call smagorinsky_alg(dtheta, du, mr, theta_star, u_star, & - derived_fields, rho, dt) + derived_fields, rho, dt, clock) end if call print_field_stats_alg(dtheta, LOG_LEVEL_INFO, printmin=-20.0_r_def, & diff --git a/science/gungho/source/algorithm/runtime_constants/physics_constants_mod.X90 b/science/gungho/source/algorithm/runtime_constants/physics_constants_mod.X90 index 94a95dcf1..463a633f6 100644 --- a/science/gungho/source/algorithm/runtime_constants/physics_constants_mod.X90 +++ b/science/gungho/source/algorithm/runtime_constants/physics_constants_mod.X90 @@ -91,6 +91,7 @@ contains use extrusion_config_mod, only: planet_radius use field_parent_mod, only: write_interface use lfric_xios_write_mod, only: write_field_generic + use mixing_config_mod, only: max_diff_factor implicit none @@ -155,7 +156,7 @@ contains ! Calculate maximum diffusion coefficient allowed in this run for ! stability based on UMDP 28 equation 32: max_diff = dx^2/(8*dt) - diffusion_const = 0.125_r_def / dt_stored + diffusion_const = 0.125_r_def * max_diff_factor / dt_stored call r_squared_w2%initialise( vector_space = w2_k0_fs ) call r_squared_wth%initialise( vector_space = wtheta_k0_fs ) 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..bd9ba01d8 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 @@ -928,7 +928,7 @@ contains call smagorinsky_alg(self%dtheta, self%du, mr, self%state(igh_t), & self%state(igh_u), & derived_fields, self%state(igh_d), & - cast_dt) + cast_dt, model_clock) end if if (use_wavedynamics) then @@ -1084,7 +1084,7 @@ contains !-------------------------------------------------------------------- call mixing_alg(mr, self%state(igh_t), & self%state(igh_u), derived_fields, & - self%state(igh_d), cast_dt ) + self%state(igh_d), cast_dt, model_clock) ! ---------------------------------------------------------------------- ! Call cloud scheme to generate cloud and latent heating after pressure diff --git a/science/gungho/source/kernel/diffusion/cap_diffusion_coeff_kernel_mod.F90 b/science/gungho/source/kernel/diffusion/cap_diffusion_coeff_kernel_mod.F90 new file mode 100644 index 000000000..23a371cc6 --- /dev/null +++ b/science/gungho/source/kernel/diffusion/cap_diffusion_coeff_kernel_mod.F90 @@ -0,0 +1,77 @@ +!------------------------------------------------------------------------------- +! (c) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!------------------------------------------------------------------------------- +! +!> @brief Caps the diffusion coefficient +module cap_diffusion_coeff_kernel_mod + + use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_READWRITE, 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) :: cap_diffusion_coeff_kernel_type + private + type(arg_type) :: meta_args(2) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, Wtheta), & + arg_type(GH_FIELD, GH_REAL, GH_READ, Wtheta) & + /) + integer :: operates_on = CELL_COLUMN + contains + procedure, nopass :: cap_diffusion_coeff_code + end type + + !----------------------------------------------------------------------------- + ! Contained functions/subroutines + !----------------------------------------------------------------------------- + public :: cap_diffusion_coeff_code + +contains + + !> @brief Caps the diffusion coefficient to ensure stability + !> @param[in] nlayers The number of layers + !> @param[in,out] visc Diffusion coefficient to be capped + !> @param[in] cap Field in Wtheta to cap to + !> @param[in] ndf_wt Number of degrees of freedom per cell for Wtheta + !> @param[in] undf_wt Number of partition degrees of freedom for Wtheta + !> @param[in] map_wt Dofmap for Wtheta + subroutine cap_diffusion_coeff_code( & + nlayers, & + visc, cap, & + ndf_wt, undf_wt, map_wt & + ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ndf_wt, undf_wt + integer(kind=i_def), intent(in) :: map_wt(ndf_wt) + real(kind=r_def), intent(inout) :: visc(undf_wt) + real(kind=r_def), intent(in) :: cap(undf_wt) + + integer(kind=i_def) :: b_idx, t_idx + + b_idx = map_wt(1) + t_idx = map_wt(1)+nlayers + + visc(b_idx:t_idx) = MIN(visc(b_idx:t_idx), cap(b_idx:t_idx)) + + end subroutine cap_diffusion_coeff_code + +end module cap_diffusion_coeff_kernel_mod diff --git a/science/gungho/source/kernel/diffusion/diffusion_flux_kernel_mod.F90 b/science/gungho/source/kernel/diffusion/diffusion_flux_kernel_mod.F90 new file mode 100644 index 000000000..bceda6eaf --- /dev/null +++ b/science/gungho/source/kernel/diffusion/diffusion_flux_kernel_mod.F90 @@ -0,0 +1,133 @@ +!------------------------------------------------------------------------------- +! (c) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!------------------------------------------------------------------------------- +! +!> @brief Computes the horizontal flux for the tracer diffusion equation +!> @details Computes the flux at W2H points for the diffusion of a tracer, +!! using a linear reconstruction. For tracer m, flux F, viscosity nu, +!! and density rho, the flux is given by: +!! F = **grad(m) +!! where < > indicates averaging to W2H points, and grad(m) is given +!! by a difference divided by the grid spacing. +module diffusion_flux_kernel_mod + + use argument_mod, only : arg_type, & + GH_FIELD, GH_REAL, & + GH_INC, GH_READ, & + CELL_COLUMN + use constants_mod, only : r_def, i_def, EPS + use fs_continuity_mod, only : W2H, W3, W2 + 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) :: diffusion_flux_kernel_type + private + type(arg_type) :: meta_args(7) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_INC, W2H), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W3), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W2H), & + arg_type(GH_FIELD, GH_REAL, GH_READ, W2H), & + 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, W2H) & + /) + integer :: operates_on = CELL_COLUMN + contains + procedure, nopass :: diffusion_flux_code + end type + + !----------------------------------------------------------------------------- + ! Contained functions/subroutines + !----------------------------------------------------------------------------- + public :: diffusion_flux_code + +contains + + !> @brief Computes the horizontal flux for the tracer diffusion equation + !> @param[in] nlayers The number of layers (in the shifted mesh) + !> @param[in,out] flux Flux field in shifted W2H to compute + !> @param[in] tracer Tracer field in shifted W3 space + !> @param[in] visc Viscosity field in shifted W2H space + !> @param[in] rho Density field in shifted W2H space + !> @param[in] dx_at_w2 Distance between cell centres at W2 points + !> @param[in] dA_at_w2 Area of faces at W2 points + !> @param[in] rmult Reciprocal of the multiplicity at W2H points + !> @param[in] ndf_w2h Number of degrees of freedom per cell for W2h + !> @param[in] undf_w2h Number of partition degrees of freedom for W2h + !> @param[in] map_w2h Dofmap for W2h + !> @param[in] ndf_w3 Number of degrees of freedom per cell for W3 + !> @param[in] undf_w3 Number of partition degrees of freedom for W3 + !> @param[in] map_w3 Dofmap for W3 + !> @param[in] ndf_w2 Number of degrees of freedom per cell for W2 + !> @param[in] undf_w2 Number of partition degrees of freedom for W2 + !> @param[in] map_w2 Dofmap for W2 + subroutine diffusion_flux_code( & + nlayers, & + flux, tracer, visc, rho, dx_at_w2, dA_at_w2, rmult, & + ndf_w2h, undf_w2h, map_w2h, & + ndf_w3, undf_w3, map_w3, & + ndf_w2, undf_w2, map_w2 & + ) + + 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) :: map_w3(ndf_w3) + integer(kind=i_def), intent(in) :: ndf_w2h, undf_w2h + integer(kind=i_def), intent(in) :: map_w2h(ndf_w2h) + integer(kind=i_def), intent(in) :: ndf_w2, undf_w2 + integer(kind=i_def), intent(in) :: map_w2(ndf_w2) + real(kind=r_def), intent(inout) :: flux(undf_w2h) + real(kind=r_def), intent(in) :: tracer(undf_w3) + real(kind=r_def), intent(in) :: visc(undf_w2h) + real(kind=r_def), intent(in) :: rho(undf_w2h) + real(kind=r_def), intent(in) :: dx_at_w2(undf_w2) + real(kind=r_def), intent(in) :: dA_at_w2(undf_w2) + real(kind=r_def), intent(in) :: rmult(undf_w2h) + + ! Encodes the direction of the gradient at the W2H points (W, S, E, N) + ! Note that the y-direction is flipped + integer(kind=i_def), parameter :: grad_sign(4) = (/ 1, -1, -1, 1 /) + + integer(kind=i_def) :: df + integer(kind=i_def) :: w2h_b_idx, w2h_t_idx + integer(kind=i_def) :: w2_b_idx, w2_t_idx + integer(kind=i_def) :: w3_b_idx, w3_t_idx + + w3_b_idx = map_w3(1)+1 + w3_t_idx = map_w3(1)+nlayers-2 + + do df = 1, ndf_w2h + ! Extract array indices for bottom and top values to use in column + ! Don't compute any values for top and bottom layers + w2h_b_idx = map_w2h(df)+1 + w2h_t_idx = map_w2h(df)+nlayers-2 + w2_b_idx = map_w2(df)+1 + w2_t_idx = map_w2(df)+nlayers-2 + + ! Only compute calculation when not at domain boundary + if (rmult(w2h_b_idx) < 1.0_r_def - EPS) then + flux(w2h_b_idx:w2h_t_idx) = flux(w2h_b_idx:w2h_t_idx) & + + visc(w2h_b_idx:w2h_t_idx) * rho(w2h_b_idx:w2h_t_idx) & + * dA_at_w2(w2_b_idx:w2_t_idx) / dx_at_w2(w2_b_idx:w2_t_idx) & + * grad_sign(df) * tracer(w3_b_idx:w3_t_idx) + end if + end do + + + end subroutine diffusion_flux_code + +end module diffusion_flux_kernel_mod diff --git a/science/gungho/source/kernel/transport/ffsl/fv_divergence_2d_kernel_mod.F90 b/science/gungho/source/kernel/transport/ffsl/fv_divergence_2d_kernel_mod.F90 index dfe5d1011..f5193bb55 100644 --- a/science/gungho/source/kernel/transport/ffsl/fv_divergence_2d_kernel_mod.F90 +++ b/science/gungho/source/kernel/transport/ffsl/fv_divergence_2d_kernel_mod.F90 @@ -17,7 +17,7 @@ module fv_divergence_2d_kernel_mod GH_FIELD, GH_REAL, & GH_WRITE, GH_READ, & CELL_COLUMN - use constants_mod, only : r_tran, i_def + use constants_mod, only : r_single, r_double, i_def use fs_continuity_mod, only : W2H, W3 use kernel_mod, only : kernel_type @@ -39,8 +39,6 @@ module fv_divergence_2d_kernel_mod arg_type(GH_FIELD, GH_REAL, GH_READ, W3) & ! detj_w3 /) integer :: operates_on = CELL_COLUMN - contains - procedure, nopass :: fv_divergence_2d_code end type !--------------------------------------------------------------------------- @@ -48,6 +46,12 @@ module fv_divergence_2d_kernel_mod !--------------------------------------------------------------------------- public :: fv_divergence_2d_code + interface fv_divergence_2d_code + module procedure & + fv_divergence_2d_r_single, & + fv_divergence_2d_r_double + end interface + contains !> @brief Computes the finite-volume divergence in the 2d horizontal directions. @@ -61,16 +65,65 @@ module fv_divergence_2d_kernel_mod !> @param[in] ndf_w2h Number of degrees of freedom for W2h per cell !> @param[in] undf_w2h Number of unique degrees of freedom for W2h !> @param[in] map_w2h Dofmap for W2h - subroutine fv_divergence_2d_code( nlayers, & - divergence, & - mass_flux, & - detj_w3, & - ndf_w3, & - undf_w3, & - map_w3, & - ndf_w2h, & - undf_w2h, & - map_w2h ) + subroutine fv_divergence_2d_r_single( nlayers, & + divergence, & + mass_flux, & + detj_w3, & + ndf_w3, & + undf_w3, & + map_w3, & + ndf_w2h, & + undf_w2h, & + map_w2h ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ndf_w3 + integer(kind=i_def), intent(in) :: undf_w3 + integer(kind=i_def), intent(in) :: ndf_w2h + integer(kind=i_def), intent(in) :: undf_w2h + integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 + integer(kind=i_def), dimension(ndf_w2h), intent(in) :: map_w2h + real(kind=r_single), dimension(undf_w3), intent(inout) :: divergence + real(kind=r_single), dimension(undf_w2h), intent(in) :: mass_flux + real(kind=r_single), dimension(undf_w3), intent(in) :: detj_w3 + + integer(kind=i_def) :: nl, w3_idx, E_idx, W_idx, N_idx, S_idx + + ! This is based on the lowest order W2 dof map + ! + ! ---4--- + ! | | + ! 1 3 horizontal + ! | | + ! ---2--- + + w3_idx = map_w3(1) + W_idx = map_w2h(1) + S_idx = map_w2h(2) + E_idx = map_w2h(3) + N_idx = map_w2h(4) + nl = nlayers - 1 + + divergence(w3_idx : w3_idx+nl) = ( & + mass_flux(E_idx : E_idx+nl) - mass_flux(W_idx : W_idx+nl) & + + mass_flux(S_idx : S_idx+nl) - mass_flux(N_idx : N_idx+nl) & + ) / detj_w3(w3_idx : w3_idx+nl) + + end subroutine fv_divergence_2d_r_single + + subroutine fv_divergence_2d_r_double( nlayers, & + divergence, & + mass_flux, & + detj_w3, & + ndf_w3, & + undf_w3, & + map_w3, & + ndf_w2h, & + undf_w2h, & + map_w2h ) implicit none @@ -82,9 +135,9 @@ subroutine fv_divergence_2d_code( nlayers, & integer(kind=i_def), intent(in) :: undf_w2h integer(kind=i_def), dimension(ndf_w3), intent(in) :: map_w3 integer(kind=i_def), dimension(ndf_w2h), intent(in) :: map_w2h - real(kind=r_tran), dimension(undf_w3), intent(inout) :: divergence - real(kind=r_tran), dimension(undf_w2h), intent(in) :: mass_flux - real(kind=r_tran), dimension(undf_w3), intent(in) :: detj_w3 + real(kind=r_double), dimension(undf_w3), intent(inout) :: divergence + real(kind=r_double), dimension(undf_w2h), intent(in) :: mass_flux + real(kind=r_double), dimension(undf_w3), intent(in) :: detj_w3 integer(kind=i_def) :: nl, w3_idx, E_idx, W_idx, N_idx, S_idx @@ -108,6 +161,6 @@ subroutine fv_divergence_2d_code( nlayers, & + mass_flux(S_idx : S_idx+nl) - mass_flux(N_idx : N_idx+nl) & ) / detj_w3(w3_idx : w3_idx+nl) - end subroutine fv_divergence_2d_code + end subroutine fv_divergence_2d_r_double end module fv_divergence_2d_kernel_mod diff --git a/science/gungho/unit-test/kernel/diffusion/cap_diffusion_coeff_kernel_mod_test.pf b/science/gungho/unit-test/kernel/diffusion/cap_diffusion_coeff_kernel_mod_test.pf new file mode 100644 index 000000000..b6d7901ac --- /dev/null +++ b/science/gungho/unit-test/kernel/diffusion/cap_diffusion_coeff_kernel_mod_test.pf @@ -0,0 +1,62 @@ +!------------------------------------------------------------------------------- +! (c) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!------------------------------------------------------------------------------- + +module cap_diffusion_coeff_kernel_mod_test + + use constants_mod, only : i_def, r_def + use funit + + implicit none + + private + public :: test_all + + @TestCase + type, extends(TestCase), public :: cap_diffusion_coeff_test_type + private + contains + procedure test_all + end type cap_diffusion_coeff_test_type + +contains + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @Test + subroutine test_all( this ) + + use cap_diffusion_coeff_kernel_mod, only : cap_diffusion_coeff_code + + implicit none + + class(cap_diffusion_coeff_test_type), intent(inout) :: this + + integer(kind=i_def), parameter :: nlayers = 3 + integer(kind=i_def), parameter :: ndf_wt = 2 + integer(kind=i_def), parameter :: undf_wt = nlayers+1 + real(kind=r_def), parameter :: tol = 1.0e-6_r_def + + integer(kind=i_def) :: map_wt(ndf_wt) + real(kind=r_def) :: visc(undf_wt) + real(kind=r_def) :: cap(undf_wt) + real(kind=r_def) :: answer(undf_wt) + + map_wt = (/ 1, 2 /) + ! Choose a mix of numbers for visc above and below the cap values + visc = (/ 0.1_r_def, 0.5_r_def, 1.0_r_def, 0.7_r_def /) + cap = (/ 0.2_r_def, 0.4_r_def, 0.6_r_def, 0.8_r_def /) + answer = (/ 0.1_r_def, 0.4_r_def, 0.6_r_def, 0.7_r_def /) + + call cap_diffusion_coeff_code( & + nlayers, & + visc, cap, & + ndf_wt, undf_wt, map_wt & + ) + + @assertEqual(answer, visc, tol) + + end subroutine test_all + +end module cap_diffusion_coeff_kernel_mod_test diff --git a/science/gungho/unit-test/kernel/diffusion/diffusion_flux_kernel_mod_test.pf b/science/gungho/unit-test/kernel/diffusion/diffusion_flux_kernel_mod_test.pf new file mode 100644 index 000000000..c9bfc289a --- /dev/null +++ b/science/gungho/unit-test/kernel/diffusion/diffusion_flux_kernel_mod_test.pf @@ -0,0 +1,147 @@ +!----------------------------------------------------------------------------- +! (c) Crown copyright 2026 Met Office. All rights reserved. +! The file LICENCE, distributed with this code, contains details of the terms +! under which the code may be used. +!----------------------------------------------------------------------------- + +!> Test the diffusion flux kernel. +!> +module diffusion_flux_kernel_mod_test + + use constants_mod, only : i_def, r_def + use funit + + implicit none + + private + public :: diffusion_flux_kernel_test_type, test_all + + @TestCase + type, extends(TestCase) :: diffusion_flux_kernel_test_type + private + contains + procedure test_all + end type diffusion_flux_kernel_test_type + +contains + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + @Test + subroutine test_all( this ) + + use diffusion_flux_kernel_mod, only : diffusion_flux_code + + implicit none + + class(diffusion_flux_kernel_test_type), intent(inout) :: this + + integer(kind=i_def) :: nlayers + integer(kind=i_def) :: ndf_w3, undf_w3 + integer(kind=i_def), allocatable :: map_w3(:) + integer(kind=i_def) :: ndf_w2h, undf_w2h + integer(kind=i_def), allocatable :: map_w2h(:) + integer(kind=i_def) :: ndf_w2, undf_w2 + integer(kind=i_def), allocatable :: map_w2(:) + real(kind=r_def), allocatable :: flux(:) + real(kind=r_def), allocatable :: flux_answer(:) + real(kind=r_def), allocatable :: rmultiplicity(:) + real(kind=r_def), allocatable :: tracer(:) + real(kind=r_def), allocatable :: visc(:) + real(kind=r_def), allocatable :: rho(:) + real(kind=r_def), allocatable :: dx_at_w2(:) + real(kind=r_def), allocatable :: dA_at_w2(:) + + real(kind=r_def), parameter :: tol = 1.0e-6_r_def + + ! Sizes of arrays + nlayers = 4 + ndf_w3 = 1 + undf_w3 = nlayers * ndf_w3 + ndf_w2h = 4 + undf_w2h = nlayers * ndf_w2h + ndf_w2 = 6 + undf_w2 = undf_w2h + nlayers + 1 + + ! Allocate arrays + allocate(map_w3(ndf_w3)) + allocate(map_w2h(ndf_w2h)) + allocate(map_w2(ndf_w2)) + allocate(flux(undf_w2h)) + allocate(flux_answer(undf_w2h)) + allocate(rmultiplicity(undf_w2h)) + allocate(tracer(undf_w3)) + allocate(visc(undf_w2h)) + allocate(rho(undf_w2h)) + allocate(dx_at_w2(undf_w2)) + allocate(dA_at_w2(undf_w2)) + + ! Set DoF maps + map_w3(1) = 1 + map_w2h(:) = (/ 1, 1+nlayers, 1+2*nlayers, 1+3*nlayers /) + map_w2(:) = (/ & + 1, 1+nlayers, 1+2*nlayers, 1+3*nlayers, 1+4*nlayers, 2+4*nlayers & + /) + + ! Set initial values + rmultiplicity(:) = 0.5_r_def + tracer(:) = (/ 10.0_r_def, 20.0_r_def, 30.0_r_def, 40.0_r_def /) + visc(:) = (/ & + 1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def, & + 1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def, & + 1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def, & + 1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def & + /) + rho(:) = (/ & + 0.1_r_def, 0.2_r_def, 0.3_r_def, 0.4_r_def, & + 0.1_r_def, 0.2_r_def, 0.3_r_def, 0.4_r_def, & + 0.1_r_def, 0.2_r_def, 0.3_r_def, 0.4_r_def, & + 0.1_r_def, 0.2_r_def, 0.3_r_def, 0.4_r_def & + /) + dx_at_w2(:) = (/ & + 100.0_r_def, 200.0_r_def, 300.0_r_def, 400.0_r_def, & + 100.0_r_def, 200.0_r_def, 300.0_r_def, 400.0_r_def, & + 100.0_r_def, 200.0_r_def, 300.0_r_def, 400.0_r_def, & + 100.0_r_def, 200.0_r_def, 300.0_r_def, 400.0_r_def, & + 1.0_r_def, 1.0_r_def, 1.0_r_def, 1.0_r_def & + /) + dA_at_w2(:) = (/ & + 1000.0_r_def, 2000.0_r_def, 3000.0_r_def, 4000.0_r_def, & + 1000.0_r_def, 2000.0_r_def, 3000.0_r_def, 4000.0_r_def, & + 1000.0_r_def, 2000.0_r_def, 3000.0_r_def, 4000.0_r_def, & + 1000.0_r_def, 2000.0_r_def, 3000.0_r_def, 4000.0_r_def, & + 1.0_r_def, 1.0_r_def, 1.0_r_def, 1.0_r_def & + /) + flux(:) = 0.0_r_def + + ! Expected answer + flux_answer(:) = (/ & + 0.0_r_def, 80.0_r_def, 270.0_r_def, 0.0_r_def, & + 0.0_r_def, -80.0_r_def, -270.0_r_def, 0.0_r_def, & + 0.0_r_def, -80.0_r_def, -270.0_r_def, 0.0_r_def, & + 0.0_r_def, 80.0_r_def, 270.0_r_def, 0.0_r_def & + /) + + call diffusion_flux_code( & + nlayers, & + flux, tracer, visc, rho, dx_at_w2, dA_at_w2, rmultiplicity, & + ndf_w2h, undf_w2h, map_w2h, & + ndf_w3, undf_w3, map_w3, & + ndf_w2, undf_w2, map_w2 & + ) + + @assertEqual(flux, flux_answer, tol) + + deallocate(map_w3) + deallocate(map_w2h) + deallocate(map_w2) + deallocate(flux) + deallocate(flux_answer) + deallocate(rmultiplicity) + deallocate(tracer) + deallocate(visc) + deallocate(rho) + deallocate(dx_at_w2) + deallocate(dA_at_w2) + + end subroutine test_all + +end module diffusion_flux_kernel_mod_test diff --git a/science/gungho/unit-test/kernel/diffusion/leonard_term_kl_kernel_mod_test.pf b/science/gungho/unit-test/kernel/diffusion/leonard_term_kl_kernel_mod_test.pf index b11b8e28a..b4b7d4481 100644 --- a/science/gungho/unit-test/kernel/diffusion/leonard_term_kl_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/diffusion/leonard_term_kl_kernel_mod_test.pf @@ -52,7 +52,10 @@ contains method = method_3d_smag, & smag_l_calc = smag_l_calc_UseDx, & leonard_term = .false., & - leonard_kl = 4.0_r_def ) + leonard_kl = 4.0_r_def, & + conservative_diffusion = .false., & + density_weighted = .false., & + max_diff_factor = 1.0_r_def ) call feign_finite_element_config( & cellshape=cellshape_quadrilateral, & diff --git a/science/gungho/unit-test/kernel/diffusion/tracer_smagorinsky_diff_kernel_mod_test.pf b/science/gungho/unit-test/kernel/diffusion/tracer_smagorinsky_diff_kernel_mod_test.pf index 445276032..cf0804bca 100644 --- a/science/gungho/unit-test/kernel/diffusion/tracer_smagorinsky_diff_kernel_mod_test.pf +++ b/science/gungho/unit-test/kernel/diffusion/tracer_smagorinsky_diff_kernel_mod_test.pf @@ -53,7 +53,10 @@ contains method = method_3d_smag, & smag_l_calc = smag_l_calc_UseDx, & leonard_term = .false., & - leonard_kl = 1.0_r_def ) + leonard_kl = 1.0_r_def, & + conservative_diffusion = .false., & + density_weighted = .false., & + max_diff_factor = 1.0_r_def ) call feign_finite_element_config( & cellshape=cellshape_quadrilateral, & 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..3877812f7 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 @@ -914,7 +914,7 @@ contains !-------------------------------------------------------------------- call mixing_alg(mr, state(igh_t), & state(igh_u), derived_fields, & - state(igh_d), cast_dt ) + state(igh_d), cast_dt, modeldb%clock ) ! Update derived variables for time level n+1 if (use_moisture) then