diff --git a/rose-stem/app/linear_model/opt/rose-app-ffsl.conf b/rose-stem/app/linear_model/opt/rose-app-ffsl.conf index e6b0e24c5c..4d971437f0 100644 --- a/rose-stem/app/linear_model/opt/rose-app-ffsl.conf +++ b/rose-stem/app/linear_model/opt/rose-app-ffsl.conf @@ -1,6 +1,6 @@ [namelist:transport] ffsl_vertical_order=5*1 horizontal_method=2,1,2,2,2 -log_space=.false.,.true.,.false.,.false.,.false. +log_space=5*.false. reversible=.false.,.true.,.false.,.false.,.false. -vertical_method=2,1,2,2,2 +vertical_method=2,3,2,2,2 diff --git a/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_dcmip301_ffsl-C24_azspice_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_dcmip301_ffsl-C24_azspice_gnu_fast-debug-64bit.txt index 94030fbcc8..a0b6d7e845 100644 --- a/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_dcmip301_ffsl-C24_azspice_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_dcmip301_ffsl-C24_azspice_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ -Inner product checksum rho = 3F202FBD547CFA2D -Inner product checksum theta = 403285E05B4CEE85 -Inner product checksum u = 432FFEFCD1E82680 +Inner product checksum rho = 3F2019E259693832 +Inner product checksum theta = 4031F066852F6E6E +Inner product checksum u = 432F625A30D467E1 diff --git a/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_nwp_gal9_ffsl-C12_MG_azspice_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_nwp_gal9_ffsl-C12_MG_azspice_gnu_fast-debug-64bit.txt index 2b6afa4932..63ff2c8f72 100644 --- a/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_nwp_gal9_ffsl-C12_MG_azspice_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/linear_model/azspice/checksum_linear_model_nwp_gal9_ffsl-C12_MG_azspice_gnu_fast-debug-64bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 3FB292C1A8DEAC36 -Inner product checksum theta = 4170C7E80A07EB51 -Inner product checksum u = 456DDD0B26100D01 -Inner product checksum mr1 = 3F96D21A5BF1ED5C -Inner product checksum mr2 = 3F73B1D0D7B555A1 -Inner product checksum mr3 = 3EF232F1B59B8986 -Inner product checksum mr4 = 3EF2CFD7D89D0460 +Inner product checksum rho = 3FB2CA8556C55BB7 +Inner product checksum theta = 41797208FAB6082E +Inner product checksum u = 45701C36CF132A6E +Inner product checksum mr1 = 3F96CFC5AE7C887E +Inner product checksum mr2 = 3F73B1BC0FE7F5D3 +Inner product checksum mr3 = 3EF232F5C6451F32 +Inner product checksum mr4 = 3EF2CFC507886DA2 Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_dcmip301_ffsl-C24_ex1a_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_dcmip301_ffsl-C24_ex1a_gnu_fast-debug-64bit.txt index f321aa8c6f..d7b49628d7 100644 --- a/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_dcmip301_ffsl-C24_ex1a_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_dcmip301_ffsl-C24_ex1a_gnu_fast-debug-64bit.txt @@ -1,3 +1,3 @@ -Inner product checksum rho = 3F202FBE9989F062 -Inner product checksum theta = 403285DFD9BCF1D0 -Inner product checksum u = 432FFEF63CBF50AD +Inner product checksum rho = 3F2019E2F24B8D66 +Inner product checksum theta = 4031F066379053E7 +Inner product checksum u = 432F6254B8B7054B diff --git a/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_nwp_gal9_ffsl-C12_MG_ex1a_gnu_fast-debug-64bit.txt b/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_nwp_gal9_ffsl-C12_MG_ex1a_gnu_fast-debug-64bit.txt index e73bed860f..0101df990c 100644 --- a/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_nwp_gal9_ffsl-C12_MG_ex1a_gnu_fast-debug-64bit.txt +++ b/rose-stem/site/meto/kgos/linear_model/ex1a/checksum_linear_model_nwp_gal9_ffsl-C12_MG_ex1a_gnu_fast-debug-64bit.txt @@ -1,9 +1,9 @@ -Inner product checksum rho = 3FB292B6ABD3E01E -Inner product checksum theta = 4170C7E8288381C8 -Inner product checksum u = 456DDD0B80AF63F2 -Inner product checksum mr1 = 3F96D2137D2B27E5 -Inner product checksum mr2 = 3F73B1D0D7B6CD3A -Inner product checksum mr3 = 3EF232F1B44659C0 -Inner product checksum mr4 = 3EF2CFD7E1969965 +Inner product checksum rho = 3FB2CA785066B264 +Inner product checksum theta = 41797208F3830022 +Inner product checksum u = 45701C36A8FED04A +Inner product checksum mr1 = 3F96CFBC6DB7C8DE +Inner product checksum mr2 = 3F73B1BC0FE0A525 +Inner product checksum mr3 = 3EF232F5C6381B6C +Inner product checksum mr4 = 3EF2CFC505737AA9 Inner product checksum mr5 = 0 Inner product checksum mr6 = 0 diff --git a/science/linear/source/algorithm/transport/control/tl_split_transport_mod.x90 b/science/linear/source/algorithm/transport/control/tl_split_transport_mod.x90 index 451a32d3db..328043d0cc 100644 --- a/science/linear/source/algorithm/transport/control/tl_split_transport_mod.x90 +++ b/science/linear/source/algorithm/transport/control/tl_split_transport_mod.x90 @@ -194,6 +194,7 @@ contains use tl_mol_conservative_alg_mod, only: tl_mol_conservative_alg use tl_mol_advective_alg_mod, only: tl_mol_advective_alg use tl_ffsl_control_alg_mod, only: tl_ffsl_control_alg + use tl_vertical_sl_advective_alg_mod, only: tl_vertical_sl_advective_alg implicit none @@ -303,19 +304,14 @@ contains ! Choose form of transport equation for vertical select case ( transport_metadata%get_equation_form() ) - case ( equation_form_conservative ) - call log_event( & - 'TL: Vertical SL conservative not coded yet', & - LOG_LEVEL_ERROR & - ) case ( equation_form_advective ) - call log_event( & - 'TL: Vertical SL advective not coded yet', & - LOG_LEVEL_ERROR & - ) + call tl_vertical_sl_advective_alg( & + field_np1, field_n, ls_field_n, tl_transport_controller & + ) + case default call log_event( & - 'Trying to solve unrecognised form of transport equation', & + 'TL vertical SL is only coded for advective form transport', & LOG_LEVEL_ERROR & ) end select diff --git a/science/linear/source/algorithm/transport/sl/tl_vertical_sl_advective_alg_mod.x90 b/science/linear/source/algorithm/transport/sl/tl_vertical_sl_advective_alg_mod.x90 new file mode 100644 index 0000000000..a4b6867cf7 --- /dev/null +++ b/science/linear/source/algorithm/transport/sl/tl_vertical_sl_advective_alg_mod.x90 @@ -0,0 +1,158 @@ +!----------------------------------------------------------------------------- +! (c) Crown copyright 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 An algorithm for performing 1D vertical semi-Lagrangian advective TL transport. +!> @details The algorithm performs a 1D vertical semi-Lagrangian advective +!! transport of a field for the linear model. It computes the perturbation field +!! at the departure point using linear or cubic interpolation, and the +!! derivative of the ls field interpolation at the departure point. + +module tl_vertical_sl_advective_alg_mod + + ! Constants and types + use constants_mod, only: r_tran, i_def, l_def + use integer_field_mod, only: integer_field_type + use log_mod, only: log_event, LOG_LEVEL_ERROR + use mesh_mod, only: mesh_type + use r_tran_field_mod, only: r_tran_field_type + use timing_mod, only: start_timing, stop_timing, & + tik, LPROF + + ! Transport control + use tl_transport_controller_mod, only: tl_transport_controller_type + use transport_controller_mod, only: transport_controller_type + use transport_counter_mod, only: transport_counter_type + use transport_metadata_mod, only: transport_metadata_type + use wind_precomputations_alg_mod, only: wind_precomputations_type + + ! Algorithms and kernels + use tl_end_of_transport_step_alg_mod, only: tl_end_of_advective_step_alg + use tl_vertical_cubic_sl_kernel_mod, only: tl_vertical_cubic_sl_kernel_type + + ! Configs + use transport_config_mod, only: vertical_sl_order, & + vertical_sl_order_cubic, & + vertical_sl_order_cubic_hermite + + implicit none + + private + + public :: tl_vertical_sl_advective_alg + +contains + + !----------------------------------------------------------------------------- + !> @brief An algorithm to perform 1D vertical advection in the linear model + !! by interpolating data to departure points (SL-advection) and + !! computing the ls field gradient. + !> @param[in,out] field_np1 ACTIVE Perturbation field at the end of + !! the transport step + !> @param[in] field ACTIVE Perturbation field at the start + !! of the transport step + !> @param[in] ls_field PASSIVE Linearisation state field at the + !! start of the transport step + !> @param[in,out] tl_transport_controller + !! Object controlling the transport + subroutine tl_vertical_sl_advective_alg( field_np1, field_n, ls_field_n, & + tl_transport_controller ) + + implicit none + + ! Arguments + type(r_tran_field_type), intent(inout) :: field_np1 + type(r_tran_field_type), intent(in) :: field_n + type(r_tran_field_type), intent(in) :: ls_field_n + type(tl_transport_controller_type), intent(inout) :: tl_transport_controller + + ! Mesh and transport objects + type(mesh_type), pointer :: mesh + type(wind_precomputations_type), pointer :: ls_wind_precomputations + type(wind_precomputations_type), pointer :: pert_wind_precomputations + type(transport_controller_type), pointer :: ls_transport_controller + type(transport_controller_type), pointer :: pert_transport_controller + type(transport_counter_type), pointer :: transport_counter + type(transport_metadata_type), pointer :: transport_metadata + integer(kind=i_def) :: sl_order + integer(kind=i_def) :: mesh_id + integer(kind=i_def) :: space + integer(kind=i_def) :: step + integer(kind=i_def) :: splitting + logical(kind=l_def) :: reversibility + integer(tik) :: id + + ! Interpolation coefficients + type(r_tran_field_type), pointer :: interp_coeffs(:) + type(integer_field_type), pointer :: interp_indices(:) + type(r_tran_field_type), pointer :: dep_dist_pert + + if ( LPROF ) call start_timing( id, 'transport.tl_sl_vertical' ) + + ! Get transport objects + ls_transport_controller => tl_transport_controller%get_ls_wind_pert_rho_controller() + pert_transport_controller => tl_transport_controller%get_pert_wind_ls_rho_controller() + ls_wind_precomputations => ls_transport_controller%get_wind_precomputations() + pert_wind_precomputations => pert_transport_controller%get_wind_precomputations() + transport_counter => pert_transport_controller%get_transport_counter() + transport_metadata => pert_transport_controller%get_transport_metadata() + + ! Get mesh and transport options + mesh => field_n%get_mesh() + mesh_id = mesh%get_id() + space = field_n%which_function_space() + reversibility = transport_metadata%get_reversible() + step = transport_counter%get_split_step_of_substep_counter() + splitting = transport_metadata%get_splitting() + + select case ( vertical_sl_order ) + ! Note TL only set up for cubic coefficients + case ( vertical_sl_order_cubic ) + sl_order = vertical_sl_order + if ( reversibility ) then + sl_order = vertical_sl_order_cubic_hermite + end if + case default + call log_event('TL vertical semi-Lagrangian only compatible with cubic coefficients', LOG_LEVEL_ERROR) + end select + + ! Set field_np1 = field_n and compute TL transport as + ! field_np1 = field_n_D - dt u grad ls_field + call invoke( setval_X(field_np1, field_n) ) + + ! Get the SL interpolation coefficients and indices + interp_coeffs => ls_wind_precomputations%get_vert_sl_coeff( & + mesh_id, space, sl_order, splitting, step & + ) + interp_indices => ls_wind_precomputations%get_vert_sl_index( & + mesh_id, space, sl_order, splitting, step & + ) + ! Require departure distance from perturbation wind + dep_dist_pert => pert_wind_precomputations%get_dep_dist_z( & + mesh_id, splitting, step & + ) + + ! Cubic (Lagrange or Hermite) interpolation of field + call invoke( tl_vertical_cubic_sl_kernel_type( field_np1, & + ls_field_n, & + dep_dist_pert, & + interp_coeffs(1), & + interp_coeffs(2), & + interp_coeffs(3), & + interp_coeffs(4), & + interp_indices(1), & + interp_indices(2), & + interp_indices(3), & + interp_indices(4) ) ) + + ! End of step: if necessary enforce min val and overwrite in blending zone + call tl_end_of_advective_step_alg( & + field_np1, field_n, transport_counter, transport_metadata & + ) + + if ( LPROF ) call stop_timing( id, 'transport.tl_sl_vertical' ) + + end subroutine tl_vertical_sl_advective_alg + +end module tl_vertical_sl_advective_alg_mod diff --git a/science/linear/source/kernel/transport/sl/tl_vertical_cubic_sl_kernel_mod.F90 b/science/linear/source/kernel/transport/sl/tl_vertical_cubic_sl_kernel_mod.F90 new file mode 100644 index 0000000000..bd2ea77ffb --- /dev/null +++ b/science/linear/source/kernel/transport/sl/tl_vertical_cubic_sl_kernel_mod.F90 @@ -0,0 +1,185 @@ +!----------------------------------------------------------------------------- +! (c) Crown copyright 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 Kernel to compute the vertical cubic semi-Lagragian advection of a field +!! in the vertical direction for the linear model. +!> @Details The 1D vertical advective transport equation for a W3/Wtheta variable +!! is solved using a cubic semi-Lagragian advection scheme. There are two +!! parts to the TL advection equation, and the update has the form +!! f^{n+1} = f^{n} - ls_u dt grad f - u_pert dt grad ls_f +!! The first two terms on the RHS are solved using the SL scheme, and the +!! third term is computed using the gradient of the ls field in the +!! departure cell. + +module tl_vertical_cubic_sl_kernel_mod + + use argument_mod, only : arg_type, & + GH_FIELD, GH_SCALAR, & + GH_REAL, GH_INTEGER, & + GH_READWRITE, GH_READ, & + CELL_COLUMN, GH_LOGICAL, & + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2 + use fs_continuity_mod, only : W2v + use constants_mod, only : r_tran, i_def, l_def, EPS_R_TRAN + 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) :: tl_vertical_cubic_sl_kernel_type + private + type(arg_type) :: meta_args(11) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, ANY_DISCONTINUOUS_SPACE_1), & ! field + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & ! ls_field + arg_type(GH_FIELD, GH_REAL, GH_READ, W2v), & ! dep_dist_pert + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & ! cubic_coef + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & ! cubic_coef + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & ! cubic_coef + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & ! cubic_coef + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & ! cubic_indices + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & ! cubic_indices + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & ! cubic_indices + arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_2) & ! cubic_indices + /) + integer :: operates_on = CELL_COLUMN + contains + procedure, nopass :: tl_vertical_cubic_sl_code + end type + + !------------------------------------------------------------------------------- + ! Contained functions/subroutines + !------------------------------------------------------------------------------- + public :: tl_vertical_cubic_sl_code + + contains + + !------------------------------------------------------------------------------- + !> @details This kernel interpolates the field to the departure point + !! using 1d-Cubic-Lagrange interpolation. It then adds the + !! contribution from the gradient of the ls field in the + !! departure cell. + !> @param[in] nlayers The number of layers + !> @param[in,out] field The perturbation field to be advected + !> @param[in] ls_field The ls field + !> @param[in] dep_dist_pert The perturbation wind departure point + !> @param[in] cubic_coef The cubic interpolation coefficients (1-4) + !> @param[in] cubic_indices The cubic interpolation indices (1-4) + !> @param[in] ndf_wf Num Dofs per cell for the field + !> @param[in] undf_wf Num Dofs in this partition for the field + !> @param[in] map_wf Dofmap for the field + !> @param[in] ndf_w2 Num Dofs per cell for dep_dist_pert + !> @param[in] undf_w2 Num Dofs in this partition for dep_dist_pert + !> @param[in] map_w2 Dofmap for dep_dist_pert + !> @param[in] ndf_wc Num Dofs per cell for the coefficients + !> @param[in] undf_wc Num Dofs per cell in this partition + !! for the coefficients + !> @param[in] map_wc Dofmap for the coefficients + !------------------------------------------------------------------------------- + subroutine tl_vertical_cubic_sl_code( nlayers, & + field, & + ls_field, & + dep_dist_pert, & + cubic_coef_1, & + cubic_coef_2, & + cubic_coef_3, & + cubic_coef_4, & + cubic_indices_1, & + cubic_indices_2, & + cubic_indices_3, & + cubic_indices_4, & + ndf_wf, undf_wf, map_wf, & + ndf_w2, undf_w2, map_w2, & + ndf_wc, undf_wc, map_wc ) + + implicit none + + ! Arguments + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ndf_wf + integer(kind=i_def), intent(in) :: undf_wf + integer(kind=i_def), intent(in) :: ndf_w2 + integer(kind=i_def), intent(in) :: undf_w2 + integer(kind=i_def), intent(in) :: ndf_wc + integer(kind=i_def), intent(in) :: undf_wc + integer(kind=i_def), intent(in) :: map_wf(ndf_wf) + integer(kind=i_def), intent(in) :: map_w2(ndf_w2) + integer(kind=i_def), intent(in) :: map_wc(ndf_wc) + real(kind=r_tran), intent(inout) :: field(undf_wf) + real(kind=r_tran), intent(in) :: ls_field(undf_wf) + real(kind=r_tran), intent(in) :: dep_dist_pert(undf_w2) + real(kind=r_tran), intent(in) :: cubic_coef_1(undf_wc) + real(kind=r_tran), intent(in) :: cubic_coef_2(undf_wc) + real(kind=r_tran), intent(in) :: cubic_coef_3(undf_wc) + real(kind=r_tran), intent(in) :: cubic_coef_4(undf_wc) + integer(kind=i_def), intent(in) :: cubic_indices_1(undf_wc) + integer(kind=i_def), intent(in) :: cubic_indices_2(undf_wc) + integer(kind=i_def), intent(in) :: cubic_indices_3(undf_wc) + integer(kind=i_def), intent(in) :: cubic_indices_4(undf_wc) + + ! Local arrays + real(kind=r_tran) :: field_local(nlayers+ndf_wf-1,4) + real(kind=r_tran) :: ls_field_local(nlayers+ndf_wf-1,2) + real(kind=r_tran) :: field_dep(nlayers+ndf_wf-1) + real(kind=r_tran) :: grad_ls_field(nlayers+ndf_wf-1) + real(kind=r_tran) :: pert_dist(nlayers+ndf_wf-1) + + ! Indices + integer(kind=i_def) :: k, nl, wf_idx, wc_idx, w2_idx + + ! nl = nlayers for w3 + ! = nlayers+1 for wtheta + nl = nlayers + ndf_wf - 1 + wf_idx = map_wf(1) + w2_idx = map_w2(1) + wc_idx = map_wc(1) + + ! Create local arrays + do k = 1, nl + field_local(k,1) = field(wf_idx + cubic_indices_1(wc_idx+k-1) - 1) + field_local(k,2) = field(wf_idx + cubic_indices_2(wc_idx+k-1) - 1) + field_local(k,3) = field(wf_idx + cubic_indices_3(wc_idx+k-1) - 1) + field_local(k,4) = field(wf_idx + cubic_indices_4(wc_idx+k-1) - 1) + ! Require indices 2 and 3 for ls_field as these lie around + ! the departure point + ls_field_local(k,1) = ls_field(wf_idx + cubic_indices_2(wc_idx+k-1) - 1) + ls_field_local(k,2) = ls_field(wf_idx + cubic_indices_3(wc_idx+k-1) - 1) + end do + + ! Interpolate field + field_dep(:) = ( & + cubic_coef_1(wc_idx : wc_idx+nl-1)*field_local(:,1) & + + cubic_coef_2(wc_idx : wc_idx+nl-1)*field_local(:,2) & + + cubic_coef_3(wc_idx : wc_idx+nl-1)*field_local(:,3) & + + cubic_coef_4(wc_idx : wc_idx+nl-1)*field_local(:,4) & + ) + + ! Compute gradient of ls_field in departure cell + grad_ls_field(:) = ls_field_local(:,2)-ls_field_local(:,1) + + ! Compute pert dist based on whether this is W3 or W3theta field + if (ndf_wf == 1) then + ! W3 field so require pert dist averaged to W3 point + pert_dist(1:nlayers) = ( dep_dist_pert(w2_idx : w2_idx + nlayers - 1) & + + dep_dist_pert(w2_idx + 1 : w2_idx + nlayers) ) / 2.0_r_tran + else + ! Wtheta field so can use pert dist at W2v point + pert_dist(1:nl) = dep_dist_pert(w2_idx : w2_idx + nl - 1) + end if + + ! Put answer back from local array into global field including + ! the gradient of the ls_field + field(wf_idx : wf_idx+nl-1) = field_dep(:) - grad_ls_field(:)*pert_dist(:) + + end subroutine tl_vertical_cubic_sl_code + +end module tl_vertical_cubic_sl_kernel_mod \ No newline at end of file diff --git a/science/linear/unit-test/kernel/transport/sl/tl_vertical_cubic_sl_kernel_mod_test.pf b/science/linear/unit-test/kernel/transport/sl/tl_vertical_cubic_sl_kernel_mod_test.pf new file mode 100644 index 0000000000..581a133f6b --- /dev/null +++ b/science/linear/unit-test/kernel/transport/sl/tl_vertical_cubic_sl_kernel_mod_test.pf @@ -0,0 +1,100 @@ +!----------------------------------------------------------------------------- +! (c) Crown copyright 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 TL vertical cubic semi Lagrangian kernel. +!> +module tl_vertical_cubic_sl_kernel_mod_test + + use constants_mod, only : i_def, l_def, r_tran + + implicit none + +contains + + !------------------------------------------------------------------ + + @test + subroutine tl_vertical_cubic_sl_kernel_test( ) + + use funit + use, intrinsic :: iso_fortran_env, only: real64 + use tl_vertical_cubic_sl_kernel_mod, only: tl_vertical_cubic_sl_code + + implicit none + + real(r_tran), parameter :: tol = 1.0e-12_r_tran + real(r_tran) :: use_tol + integer(i_def), parameter :: nlayers = 5 + integer(i_def), parameter :: ndf_w3 = 1 + integer(i_def), parameter :: undf_w3 = nlayers + integer(i_def), parameter :: ndf_w2 = 2 + integer(i_def), parameter :: undf_w2 = nlayers+1 + integer(i_def), parameter :: ndf_wc = 1 + integer(i_def), parameter :: undf_wc = nlayers + integer(i_def), dimension(ndf_w3) :: map_w3 + integer(i_def), dimension(ndf_w2) :: map_w2 + integer(i_def), dimension(ndf_wc) :: map_wc + real(r_tran), dimension(undf_w3) :: field_in, field, ls_field, field_out + real(r_tran), dimension(undf_w2) :: dep_dist_pert + real(r_tran), dimension(undf_wc) :: cubic_coeff_1, cubic_coeff_2, & + cubic_coeff_3, cubic_coeff_4 + integer(i_def), dimension(undf_wc) :: cubic_index_1, cubic_index_2, & + cubic_index_3, cubic_index_4 + integer(i_def) :: i + + ! Create the dof map for any number of layers + map_w3(:) = 1 + map_wc(:) = 1 + map_w2(:) = (/ 1, 2 /) + + ! Set the index + cubic_index_1(:) = (/ 1_i_def, 1_i_def, 1_i_def, 2_i_def, 3_i_def /) + cubic_index_2(:) = (/ 1_i_def, 1_i_def, 2_i_def, 3_i_def, 4_i_def /) + cubic_index_3(:) = (/ 1_i_def, 2_i_def, 3_i_def, 4_i_def, 5_i_def /) + cubic_index_4(:) = (/ 2_i_def, 3_i_def, 4_i_def, 5_i_def, 5_i_def /) + + ! Set up field, ls_field, and dep_dist_pert + field_in(:) = (/1.0_r_tran, 2.0_r_tran, 3.0_r_tran, 4.0_r_tran, 5.0_r_tran /) + ls_field(:) = (/1.0_r_tran, 2.0_r_tran, 3.0_r_tran, 2.0_r_tran, 1.0_r_tran /) + dep_dist_pert(:) = (/0.0_r_tran, 0.5_r_tran, -0.5_r_tran, & + 0.25_r_tran, 1.0_r_tran, 0.0_r_tran /) + + ! Test with variable coefficients + cubic_coeff_1(:) = (/ 1.0_r_tran, 0.25_r_tran, 0.0_r_tran, 0.0_r_tran, 0.0_r_tran /) + cubic_coeff_2(:) = (/ 0.0_r_tran, 0.25_r_tran, 0.5_r_tran, 0.5_r_tran, 0.0_r_tran /) + cubic_coeff_3(:) = (/ 0.0_r_tran, 0.25_r_tran, 0.5_r_tran, 0.5_r_tran, 0.0_r_tran /) + cubic_coeff_4(:) = (/ 0.0_r_tran, 0.25_r_tran, 0.0_r_tran, 0.0_r_tran, 1.0_r_tran /) + + field(:) = field_in(:) + ! field_out is the sum of the coefficients multiplied by the field at the + ! corresponding index, subtract dep_dist_pert (averaged to cell centre) + ! multiplied by the gradient of the ls_field in the departure cell + field_out(:) = (/1.0_r_tran, 1.75_r_tran, 2.625_r_tran, 4.125_r_tran, 5.5_r_tran /) + + call tl_vertical_cubic_sl_code( nlayers, & + field, ls_field, dep_dist_pert, & + cubic_coeff_1, cubic_coeff_2, & + cubic_coeff_3, cubic_coeff_4, & + cubic_index_1, cubic_index_2, & + cubic_index_3, cubic_index_4, & + ndf_w3, undf_w3, map_w3, & + ndf_w2, undf_w2, map_w2, & + ndf_wc, undf_wc, map_wc ) + + if ( r_tran == real64 ) then + use_tol = tol + else + use_tol = 10.0_r_tran*spacing(maxval(field_in(:))) + end if + + do i = 1,nlayers + @assertEqual(field_out(i), field(i) , use_tol) + end do + + end subroutine tl_vertical_cubic_sl_kernel_test + +end module tl_vertical_cubic_sl_kernel_mod_test