diff --git a/applications/lfric2lfric/source/kernel/lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod.F90 b/applications/lfric2lfric/source/kernel/lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod.F90 new file mode 100644 index 000000000..be7e17aa7 --- /dev/null +++ b/applications/lfric2lfric/source/kernel/lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod.F90 @@ -0,0 +1,196 @@ +!----------------------------------------------------------------------------- +! (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 SHARKS Perform the conservative weighted injection-prolongation operation +!! of a coarse grid scalar to a fine grid face field +!> @details SHARKS Prolong the coarse grid correction into all cells on the fine grid +!! that are contained in that coarse grid cell. This is an injection +!! on the "mass" field, so that the prolongation is conservative. +!! This kernel only works for the lowest-order face spaces. +!! This is a copy of prolong_scalar_weighted_kernel_mod, with a +!! modification to take into account multidata values in grid points + +module lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod + +use constants_mod, only: i_def, r_double, r_single +use kernel_mod, only: kernel_type +use argument_mod, only: arg_type, CELL_COLUMN, & + GH_FIELD, GH_INTEGER, GH_REAL, & + GH_READ, GH_READWRITE, GH_SCALAR, W3,& + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2 + +implicit none + +private + +type, public, extends(kernel_type) :: lfric2lfric_vertical_linear_interp_linear_extrap_kernel_type + private + type(arg_type) :: meta_args(6) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2) & + /) + integer :: operates_on = CELL_COLUMN +end type lfric2lfric_vertical_linear_interp_linear_extrap_kernel_type + +public :: lfric2lfric_vertical_linear_interp_linear_extrap_kernel_code + + ! Generic interface for real32 and real64 types + interface lfric2lfric_vertical_linear_interp_linear_extrap_kernel_code + module procedure & + lfric2lfric_vertical_linear_interp_linear_extrap_code_r_single !, & + !lfric2lfric_vertical_linear_interp_linear_extrap_code_r_double + end interface + +contains + + !> @brief Performs the weighted injection-prolongation for scalar fields + !> @param[in] nlayers Number of layers in a model column (in the destination field) + !> @param[in,out] destination_field The destination field to interpolate onto + !> @param[in] source_field Source field to be interpolated + !> @param[in] source_layers Number of layers in the source field + !> @param[in] ncell Number of cells in the partition + !! of both grids + !> @param[in] source_heights Heights of vertical layers in the source field + !> @param[in] dest_heights Heights of vertical layers in the destination field + !> @param[in] ndf_dest Num of DoFs per cell on the destination grid + !> @param[in] undf_dest Total num of DoFs on the destination grid + !! for this mesh partition + !> @param[in] map_dest DoFmap of cells on the destination grid + !> @param[in] ndf_source Num of DoFs per cell on the source grid + !> @param[in] undf_source Total num of DoFs on the source + !! grid for this mesh partition + !> @param[in] map_source DoFmap of cells on the source grid + + ! R_SINGLE PRECISION + ! ================== + subroutine lfric2lfric_vertical_linear_interp_linear_extrap_code_r_single( & + nlayers, & + destination_field, & + source_field, & + source_layers, & + ncell, & + dest_heights, & + source_heights, & + ndf_dest, & + undf_dest, & + map_dest, & + ndf_source, & + undf_source, & + map_source) + + implicit none + + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ncell + integer(kind=i_def), intent(in) :: source_layers + integer(kind=i_def), intent(in) :: ndf_dest, ndf_source + integer(kind=i_def), intent(in) :: undf_dest, undf_source + integer(kind=i_def), intent(in), dimension(ndf_dest) :: map_dest(ndf_dest) + integer(kind=i_def), intent(in), dimension(ndf_source) :: map_source(ndf_source) + real(kind=r_single), intent(inout), dimension(ndf_dest) :: destination_field(undf_dest) + real(kind=r_single), intent(in), dimension(ndf_source) :: source_field(undf_source) + real(kind=r_single), intent(in), dimension(ndf_dest) :: dest_heights(undf_dest) + real(kind=r_single), intent(in), dimension(ndf_source) :: source_heights(undf_source) + + integer(kind=i_def) :: multidata, df, k, m, top_df, level_below(nlayers) + + ! Assume lowest order W3 or Wtheta space + df = 1 + ! Loop is 0 -> nlayers-1 for W3 fields, but 0 -> nlayers for Wtheta fields + dest_top_df = nlayers - 2 + ndf_dest + source_top_df = source_layers - 2 + ndf_source + ! Number of multidata values per grid cell + multidata = undf_dest/((dest_top_df+1)*ncell) - 1 + + do kk= 0, dest_top_df + do k= 0, source_top_df + if ( (source_heights(k) > dest_heights(kk)) .AND. & + (level_below(kk) == source_layers) ) THEN + level_below(kk) = k-1 + ! potential optimisation: start from level_below(kk-1) + end if + end do + end do + + do m = 0, multidata + do kk=0, dest_top_df + ! EXTRAPOLATION METHOD - ! Linear extrapolation at top and bottom + + ! IF ( desired_r(j) > r_at_data(j,data_levels) ) THEN + if (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) + > source_heights(DATALEVELS)) then + + ! If requested level is above top of model, do linear + ! extrapolation using data on top and second top levels. + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = source_field(DATALEVELS) + + + data_out(j) = data_in(j,data_levels) + (desired_r(j) & + - r_at_data(j,data_levels)) * (data_in(j,data_levels) & + - data_in(j,data_levels-1))/(r_at_data(j,data_levels) & + - r_at_data(j,data_levels-1)) + + ! ELSE IF (desired_r(j) == r_at_data(j,data_levels) ) THEN + else if (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) + == source_heights(DATALEVELS)) then + + ! data_out(j) = data_in(j,data_levels) + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = source_field(DATALEVELS) + + ! IF ( desired_r(j) < r_at_data(j,1) ) THEN + else if (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) + < source_heights(map_source(df) + m*(source_top_df+1) + level_below(1))) then + + ! If requested level is below bottom of model, do linear + ! extrapolation using data on first and second levels. + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = source_field(map_source(df) + m*(source_top_df+1) + 1) + + data_out(j) = data_in(j,1) + (desired_r(j) & + - r_at_data(j,1)) * (data_in(j,1) & + - data_in(j,2))/(r_at_data(j,1) - r_at_data(j,2)) + + ! ELSE IF (desired_r(j) == r_at_data(j,1) ) THEN + else if (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) + == source_heights(map_source(df) + m*(source_top_df+1) + level_below(1))) then + + ! data_out(j) = data_in(j,1) + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = source_field(map_source(df) + m*(source_top_df+1) + 1) + + else + + ! Linearly interpolate + + ! dk(kk) = ( (dh(kk) - sh(lb(kk))) * sf(lb(kk)+1) - (dh(kk) - sh(lb(kk)+1)) * sf(lb(kk)) ) + ! / (sh(lb(kk)+1) - sh(lb(kk))) + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = & + ( (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) & + - source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)) ) & + * source_field(map_source(df) + m*(source_top_df+1) + level_below(kk)+1) & + - (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) & + - source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)+1)) & + * source_field(map_source(df) + m*(source_top_df+1) + level_below(kk)) ) & + / ( source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)+1) & + - source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)) ) + end if + end do + end do + + end subroutine lfric2lfric_vertical_linear_interp_linear_extrap_code_r_single + + ! R_DOUBLE PRECISION + ! ================== + !subroutine lfric2lfric_vertical_linear_interp_linear_extrap_code_r_double & + +!!!!!!! SHARKS COPY R SINGLE TO HERE, change real fields from r_single to r_double + + ! end subroutine lfric2lfric_vertical_linear_interp_linear_extrap_code_r_double + +end module lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod diff --git a/applications/lfric2lfric/source/kernel/lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod.F90 b/applications/lfric2lfric/source/kernel/lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod.F90 new file mode 100644 index 000000000..c0b6d060d --- /dev/null +++ b/applications/lfric2lfric/source/kernel/lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod.F90 @@ -0,0 +1,173 @@ +!----------------------------------------------------------------------------- +! (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 SHARKS Perform the conservative weighted injection-prolongation operation +!! of a coarse grid scalar to a fine grid face field +!> @details SHARKS Prolong the coarse grid correction into all cells on the fine grid +!! that are contained in that coarse grid cell. This is an injection +!! on the "mass" field, so that the prolongation is conservative. +!! This kernel only works for the lowest-order face spaces. +!! This is a copy of prolong_scalar_weighted_kernel_mod, with a +!! modification to take into account multidata values in grid points + +module lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod + +use constants_mod, only: i_def, r_double, r_single +use kernel_mod, only: kernel_type +use argument_mod, only: arg_type, CELL_COLUMN, & + GH_FIELD, GH_INTEGER, GH_REAL, & + GH_READ, GH_READWRITE, GH_SCALAR, W3,& + ANY_DISCONTINUOUS_SPACE_1, & + ANY_DISCONTINUOUS_SPACE_2 + +implicit none + +private + +type, public, extends(kernel_type) :: lfric2lfric_vertical_linear_interp_no_extrap_kernel_type + private + type(arg_type) :: meta_args(6) = (/ & + arg_type(GH_FIELD, GH_REAL, GH_READWRITE, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_SCALAR, GH_INTEGER, GH_READ), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_1), & + arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_DISCONTINUOUS_SPACE_2) & + /) + integer :: operates_on = CELL_COLUMN +end type lfric2lfric_vertical_linear_interp_no_extrap_kernel_type + +public :: lfric2lfric_vertical_linear_interp_no_extrap_kernel_code + + ! Generic interface for real32 and real64 types + interface lfric2lfric_vertical_linear_interp_no_extrap_kernel_code + module procedure & + lfric2lfric_vertical_linear_interp_no_extrap_code_r_single !, & + !lfric2lfric_vertical_linear_interp_no_extrap_code_r_double + end interface + +contains + + !> @brief Performs the weighted injection-prolongation for scalar fields + !> @param[in] nlayers Number of layers in a model column (in the destination field) + !> @param[in,out] destination_field The destination field to interpolate onto + !> @param[in] source_field Source field to be interpolated + !> @param[in] source_layers Number of layers in the source field + !> @param[in] ncell Number of cells in the partition + !! of both grids + !> @param[in] source_heights Heights of vertical layers in the source field + !> @param[in] dest_heights Heights of vertical layers in the destination field + !> @param[in] ndf_dest Num of DoFs per cell on the destination grid + !> @param[in] undf_dest Total num of DoFs on the destination grid + !! for this mesh partition + !> @param[in] map_dest DoFmap of cells on the destination grid + !> @param[in] ndf_source Num of DoFs per cell on the source grid + !> @param[in] undf_source Total num of DoFs on the source + !! grid for this mesh partition + !> @param[in] map_source DoFmap of cells on the source grid + + ! R_SINGLE PRECISION + ! ================== + subroutine lfric2lfric_vertical_linear_interp_no_extrap_code_r_single( & + nlayers, & + destination_field, & + source_field, & + source_layers, & + ncell, & + dest_heights, & + source_heights, & + ndf_dest, & + undf_dest, & + map_dest, & + ndf_source, & + undf_source, & + map_source) + + implicit none + + integer(kind=i_def), intent(in) :: nlayers + integer(kind=i_def), intent(in) :: ncell + integer(kind=i_def), intent(in) :: source_layers + integer(kind=i_def), intent(in) :: ndf_dest, ndf_source + integer(kind=i_def), intent(in) :: undf_dest, undf_source + integer(kind=i_def), intent(in), dimension(ndf_dest) :: map_dest(ndf_dest) + integer(kind=i_def), intent(in), dimension(ndf_source) :: map_source(ndf_source) + real(kind=r_single), intent(inout), dimension(ndf_dest) :: destination_field(undf_dest) + real(kind=r_single), intent(in), dimension(ndf_source) :: source_field(undf_source) + real(kind=r_single), intent(in), dimension(ndf_dest) :: dest_heights(undf_dest) + real(kind=r_single), intent(in), dimension(ndf_source) :: source_heights(undf_source) + + integer(kind=i_def) :: multidata, df, k, m, top_df, level_below(nlayers) + + ! Assume lowest order W3 or Wtheta space + df = 1 + ! Loop is 0 -> nlayers-1 for W3 fields, but 0 -> nlayers for Wtheta fields + dest_top_df = nlayers - 2 + ndf_dest + source_top_df = source_layers - 2 + ndf_source + ! Number of multidata values per grid cell + multidata = undf_dest/((dest_top_df+1)*ncell) - 1 + + do kk= 0, dest_top_df + do k= 0, source_top_df + if ( (source_heights(k) > dest_heights(kk)) .AND. & + (level_below(kk) == source_layers) ) THEN + level_below(kk) = k-1 + ! potential optimisation: start from level_below(kk-1) + end if + end do + end do + + do m = 0, multidata + do kk=0, dest_top_df + ! EXTRAPOLATION METHOD - ! No linear extrapolation at top or bottom + + ! IF ( desired_r(j) >= r_at_data(j,data_levels) ) THEN + ! data_out(j) = data_in(j,data_levels) + if (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) + >= source_heights(DATALEVELS)) then + + ! Top: Set to top input data + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = source_field(DATALEVELS) + + + ! ELSE IF ( desired_r(j) <= r_at_data(j,1) ) THEN + ! data_out(j) = data_in(j,1) + else if (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) + <= source_heights(map_source(df) + m*(source_top_df+1) + level_below(1))) then + + ! Bottom: Set to bottom input data + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = source_field(map_source(df) + m*(source_top_df+1) + 1) + + else + + ! Linearly interpolate + + ! dk(kk) = ( (dh(kk) - sh(lb(kk))) * sf(lb(kk)+1) - (dh(kk) - sh(lb(kk)+1)) * sf(lb(kk)) ) + ! / (sh(lb(kk)+1) - sh(lb(kk))) + destination_field(map_dest(df) + m*(dest_top_df+1) + kk) = & + ( (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) & + - source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)) ) & + * source_field(map_source(df) + m*(source_top_df+1) + level_below(kk)+1) & + - (dest_heights(map_dest(df) + m*(dest_top_df+1) + kk) & + - source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)+1)) & + * source_field(map_source(df) + m*(source_top_df+1) + level_below(kk)) ) & + / ( source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)+1) & + - source_heights(map_source(df) + m*(source_top_df+1) + level_below(kk)) ) + end if + end do + end do + + end subroutine lfric2lfric_vertical_linear_interp_no_extrap_code_r_single + + ! R_DOUBLE PRECISION + ! ================== + !subroutine lfric2lfric_vertical_linear_interp_no_extrap_code_r_double & + +!!!!!!! SHARKS COPY R SINGLE TO HERE, change real fields from r_single to r_double + + ! end subroutine lfric2lfric_vertical_linear_interp_no_extrap_code_r_double + +end module lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod diff --git a/applications/lfric2lfric/unit-test/kernel/lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod_test.pf b/applications/lfric2lfric/unit-test/kernel/lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod_test.pf new file mode 100644 index 000000000..e7e9c1bb4 --- /dev/null +++ b/applications/lfric2lfric/unit-test/kernel/lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod_test.pf @@ -0,0 +1,239 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2025 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 vertical linear interpolation kernel +module lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod_test + + use constants_mod, only : r_def, i_def + use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w3_m3x3_q3x3x3_size + use get_unit_test_m3x3_dofmap_mod, only : get_w3_m3x3_dofmap + + use funit + + implicit none + + private + public:: test_all + + contains + +@Test + subroutine test_all() + + use lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod, only: lfric2lfric_vertical_linear_interp_linear_extrap_kernel_code + + implicit none + + integer(kind=i_def) :: nlayers + integer(kind=i_def) :: ncell + integer(kind=i_def) :: dim_space, dim_space_diff, nqp_h, nqp_v + integer(kind=i_def), parameter :: source_layers + integer(kind=i_def), parameter :: ndf_dest, ndf_source + integer(kind=i_def), parameter :: undf_dest, undf_source + integer(kind=i_def), allocatable :: map_dest(:,:) + integer(kind=i_def), allocatable :: map_source(:,:) + real(kind=r_single), allocatable :: destination_field(:) + real(kind=r_single), allocatable :: source_field(:) + real(kind=r_single), allocatable :: dest_heights(:) + real(kind=r_single), allocatable :: source_heights(:) + real(r_def), parameter :: tol = 1.0e-12_r_def + + nlayers = 3_i_def + source_layers = 4_i_def + call get_w3_m3x3_q3x3x3_size( ndf_source, undf_source, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + source_layers ) + + call get_w3_m3x3_dofmap(map_source, 4_i_def) + + call get_w3_m3x3_q3x3x3_size( ndf_dest, undf_dest, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + nlayers ) + call get_w3_m3x3_dofmap(map_dest, 5_i_def) + + ! allocate + allocate(destination_field(undf_dest)) + allocate(source_field(undf_source)) + allocate(dest_heights(undf_dest)) + allocate(source_heights(undf_source)) + + source_field = (/& + (1.1_r_def, 2.3_r_def, 3.2_r_def, 4.1_r_def), + (1.1_r_def, 2.2_r_def, 3.1_r_def, 4.4_r_def), + (0.7_r_def, 2.3_r_def, 2.9_r_def, 4.0_r_def), + (0.5_r_def, 1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def, 4.5_r_def), + (1.4_r_def, 2.0_r_def, 3.4_r_def, 4.2_r_def), + (1.0_r_def, 2.4_r_def, 3.3_r_def, 4.4_r_def), + (1.1_r_def, 2.1_r_def, 3.3_r_def, 4.1_r_def) + /) + expected_field = (/& + (1.7_r_def, 2.75_r_def, 3.65_r_def), + (1.65_r_def, 2.65_r_def, 3.75_r_def), + (1.5_r_def, 2.6_r_def, 3.45_r_def), + (1.0_r_def, 2.0_r_def, 3.0_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (2.0_r_def, 3.0_r_def, 4.0_r_def), + (1.7_r_def, 2.7_r_def, 3.8_r_def), + (1.7_r_def, 2.85_r_def, 3.85_r_def), + (1.6_r_def, 2.7_r_def, 3.7_r_def) + /) + source_heights = (/& + (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), + (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), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def) + /) + dest_heights = (/& + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def) + /) + + call lfric2lfric_vertical_linear_interp_linear_extrap_code_r_single( & + nlayers, & + destination_field, & + source_field, & + source_layers, & + ncell, & + dest_heights, & + source_heights, & + ndf_dest, & + undf_dest, & + map_dest, & + ndf_source, & + undf_source, & + map_source) + + @assertEqual(expected_field, destination_field, tol) + + end subroutine test_all + +@Test + subroutine test_all_with_extrap() + + use lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod, only: lfric2lfric_vertical_linear_interp_linear_extrap_kernel_code + + implicit none + + integer(kind=i_def) :: nlayers + integer(kind=i_def) :: ncell + integer(kind=i_def) :: dim_space, dim_space_diff, nqp_h, nqp_v + integer(kind=i_def), parameter :: source_layers + integer(kind=i_def), parameter :: ndf_dest, ndf_source + integer(kind=i_def), parameter :: undf_dest, undf_source + integer(kind=i_def), allocatable :: map_dest(:,:) + integer(kind=i_def), allocatable :: map_source(:,:) + real(kind=r_single), allocatable :: destination_field(:) + real(kind=r_single), allocatable :: source_field(:) + real(kind=r_single), allocatable :: dest_heights(:) + real(kind=r_single), allocatable :: source_heights(:) + real(r_def), parameter :: tol = 1.0e-12_r_def + + nlayers = 3_i_def + source_layers = 4_i_def + call get_w3_m3x3_q3x3x3_size( ndf_source, undf_source, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + source_layers ) + + call get_w3_m3x3_dofmap(map_source, 4_i_def) + + call get_w3_m3x3_q3x3x3_size( ndf_dest, undf_dest, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + nlayers ) + call get_w3_m3x3_dofmap(map_dest, 5_i_def) + + ! allocate + allocate(destination_field(undf_dest)) + allocate(source_field(undf_source)) + allocate(dest_heights(undf_dest)) + allocate(source_heights(undf_source)) + + source_field = (/& + (1.1_r_def, 2.3_r_def, 3.2_r_def, 4.1_r_def), + (1.1_r_def, 2.2_r_def, 3.1_r_def, 4.4_r_def), + (0.7_r_def, 2.3_r_def, 2.9_r_def, 4.0_r_def), + (0.5_r_def, 1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def, 4.5_r_def), + (1.4_r_def, 2.0_r_def, 3.4_r_def, 4.2_r_def), + (1.0_r_def, 2.4_r_def, 3.3_r_def, 4.4_r_def), + (1.1_r_def, 2.1_r_def, 3.3_r_def, 4.1_r_def) + /) + expected_field = (/& + (0.5_r_def, 1.7_r_def, 2.75_r_def, 3.65_r_def, 4.55_r_def), + (0.55_r_def, 1.65_r_def, 2.65_r_def, 3.75_r_def, 5.05_r_def), + (-0.1_r_def, 1.5_r_def, 2.6_r_def, 3.45_r_def, 4.55_r_def), + (0.0_r_def, 1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def), + (0.5_r_def, 1.5_r_def, 2.5_r_def, 3.5_r_def, 4.5_r_def), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def, 5.0_r_def), + (1.1_r_def, 1.7_r_def, 2.7_r_def, 3.8_r_def, 4.6_r_def), + (0.3_r_def, 1.7_r_def, 2.85_r_def, 3.85_r_def, 4.95_r_def), + (0.6_r_def, 1.6_r_def, 2.7_r_def, 3.7_r_def, 4.5_r_def) + /) + source_heights = (/& + (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), + (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), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def) + /) + dest_heights = (/& + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def) + /) + + call lfric2lfric_vertical_linear_interp_linear_extrap_code_r_single( & + nlayers, & + destination_field, & + source_field, & + source_layers, & + ncell, & + dest_heights, & + source_heights, & + ndf_dest, & + undf_dest, & + map_dest, & + ndf_source, & + undf_source, & + map_source) + + @assertEqual(expected_field, destination_field, tol) + + end subroutine test_all_with_extrap + +end module diff --git a/applications/lfric2lfric/unit-test/kernel/lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod_test.pf b/applications/lfric2lfric/unit-test/kernel/lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod_test.pf new file mode 100644 index 000000000..b98e13294 --- /dev/null +++ b/applications/lfric2lfric/unit-test/kernel/lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod_test.pf @@ -0,0 +1,239 @@ +!----------------------------------------------------------------------------- +! (C) Crown copyright 2025 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 vertical linear interpolation kernel +module lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod_test + + use constants_mod, only : r_def, i_def + use get_unit_test_m3x3_q3x3x3_sizes_mod, only : get_w3_m3x3_q3x3x3_size + use get_unit_test_m3x3_dofmap_mod, only : get_w3_m3x3_dofmap + + use funit + + implicit none + + private + public:: test_all + + contains + +@Test + subroutine test_all() + + use lfric2lfric_vertical_linear_interp_no_extrap_kernel_mod, only: lfric2lfric_vertical_linear_interp_no_extrap_kernel_code + + implicit none + + integer(kind=i_def) :: nlayers + integer(kind=i_def) :: ncell + integer(kind=i_def) :: dim_space, dim_space_diff, nqp_h, nqp_v + integer(kind=i_def), parameter :: source_layers + integer(kind=i_def), parameter :: ndf_dest, ndf_source + integer(kind=i_def), parameter :: undf_dest, undf_source + integer(kind=i_def), allocatable :: map_dest(:,:) + integer(kind=i_def), allocatable :: map_source(:,:) + real(kind=r_single), allocatable :: destination_field(:) + real(kind=r_single), allocatable :: source_field(:) + real(kind=r_single), allocatable :: dest_heights(:) + real(kind=r_single), allocatable :: source_heights(:) + real(r_def), parameter :: tol = 1.0e-12_r_def + + nlayers = 3_i_def + source_layers = 4_i_def + call get_w3_m3x3_q3x3x3_size( ndf_source, undf_source, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + source_layers ) + + call get_w3_m3x3_dofmap(map_source, 4_i_def) + + call get_w3_m3x3_q3x3x3_size( ndf_dest, undf_dest, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + nlayers ) + call get_w3_m3x3_dofmap(map_dest, 5_i_def) + + ! allocate + allocate(destination_field(undf_dest)) + allocate(source_field(undf_source)) + allocate(dest_heights(undf_dest)) + allocate(source_heights(undf_source)) + + source_field = (/& + (1.1_r_def, 2.3_r_def, 3.2_r_def, 4.1_r_def), + (1.1_r_def, 2.2_r_def, 3.1_r_def, 4.4_r_def), + (0.7_r_def, 2.3_r_def, 2.9_r_def, 4.0_r_def), + (0.5_r_def, 1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def, 4.5_r_def), + (1.4_r_def, 2.0_r_def, 3.4_r_def, 4.2_r_def), + (1.0_r_def, 2.4_r_def, 3.3_r_def, 4.4_r_def), + (1.1_r_def, 2.1_r_def, 3.3_r_def, 4.1_r_def) + /) + expected_field = (/& + (1.7_r_def, 2.75_r_def, 3.65_r_def), + (1.65_r_def, 2.65_r_def, 3.75_r_def), + (1.5_r_def, 2.6_r_def, 3.45_r_def), + (1.0_r_def, 2.0_r_def, 3.0_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (2.0_r_def, 3.0_r_def, 4.0_r_def), + (1.7_r_def, 2.7_r_def, 3.8_r_def), + (1.7_r_def, 2.85_r_def, 3.85_r_def), + (1.6_r_def, 2.7_r_def, 3.7_r_def) + /) + source_heights = (/& + (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), + (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), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def) + /) + dest_heights = (/& + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def) + /) + + call lfric2lfric_vertical_linear_interp_no_extrap_code_r_single( & + nlayers, & + destination_field, & + source_field, & + source_layers, & + ncell, & + dest_heights, & + source_heights, & + ndf_dest, & + undf_dest, & + map_dest, & + ndf_source, & + undf_source, & + map_source) + + @assertEqual(expected_field, destination_field, tol) + + end subroutine test_all + +@Test + subroutine test_all_with_extrap() + + use lfric2lfric_vertical_linear_interp_linear_extrap_kernel_mod, only: lfric2lfric_vertical_linear_interp_linear_extrap_kernel_code + + implicit none + + integer(kind=i_def) :: nlayers + integer(kind=i_def) :: ncell + integer(kind=i_def) :: dim_space, dim_space_diff, nqp_h, nqp_v + integer(kind=i_def), parameter :: source_layers + integer(kind=i_def), parameter :: ndf_dest, ndf_source + integer(kind=i_def), parameter :: undf_dest, undf_source + integer(kind=i_def), allocatable :: map_dest(:,:) + integer(kind=i_def), allocatable :: map_source(:,:) + real(kind=r_single), allocatable :: destination_field(:) + real(kind=r_single), allocatable :: source_field(:) + real(kind=r_single), allocatable :: dest_heights(:) + real(kind=r_single), allocatable :: source_heights(:) + real(r_def), parameter :: tol = 1.0e-12_r_def + + nlayers = 3_i_def + source_layers = 4_i_def + call get_w3_m3x3_q3x3x3_size( ndf_source, undf_source, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + source_layers ) + + call get_w3_m3x3_dofmap(map_source, 4_i_def) + + call get_w3_m3x3_q3x3x3_size( ndf_dest, undf_dest, & + ncell, & + dim_space, dim_space_diff, & + nqp_h, nqp_v, & + nlayers ) + call get_w3_m3x3_dofmap(map_dest, 5_i_def) + + ! allocate + allocate(destination_field(undf_dest)) + allocate(source_field(undf_source)) + allocate(dest_heights(undf_dest)) + allocate(source_heights(undf_source)) + + source_field = (/& + (1.1_r_def, 2.3_r_def, 3.2_r_def, 4.1_r_def), + (1.1_r_def, 2.2_r_def, 3.1_r_def, 4.4_r_def), + (0.7_r_def, 2.3_r_def, 2.9_r_def, 4.0_r_def), + (0.5_r_def, 1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def, 4.5_r_def), + (1.4_r_def, 2.0_r_def, 3.4_r_def, 4.2_r_def), + (1.0_r_def, 2.4_r_def, 3.3_r_def, 4.4_r_def), + (1.1_r_def, 2.1_r_def, 3.3_r_def, 4.1_r_def) + /) + expected_field = (/& + (1.1_r_def, 1.7_r_def, 2.75_r_def, 3.65_r_def, 4.1_r_def), + (1.1_r_def, 1.65_r_def, 2.65_r_def, 3.75_r_def, 4.4_r_def), + (0.7_r_def, 1.5_r_def, 2.6_r_def, 3.45_r_def, 4.0_r_def), + (0.5_r_def, 1.0_r_def, 2.0_r_def, 3.0_r_def, 3.5_r_def), + (1.0_r_def, 1.5_r_def, 2.5_r_def, 3.5_r_def, 4.0_r_def), + (1.5_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def, 4.5_r_def), + (1.4_r_def, 1.7_r_def, 2.7_r_def, 3.8_r_def, 4.2_r_def), + (1.0_r_def, 1.7_r_def, 2.85_r_def, 3.85_r_def, 4.4_r_def), + (1.1_r_def, 1.6_r_def, 2.7_r_def, 3.7_r_def, 4.1_r_def) + /) + source_heights = (/& + (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), + (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), + (1.0_r_def, 2.0_r_def, 3.0_r_def, 4.0_r_def) + /) + dest_heights = (/& + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def), + (1.5_r_def, 2.5_r_def, 3.5_r_def) + /) + + call lfric2lfric_vertical_linear_interp_linear_extrap_code_r_single( & + nlayers, & + destination_field, & + source_field, & + source_layers, & + ncell, & + dest_heights, & + source_heights, & + ndf_dest, & + undf_dest, & + map_dest, & + ndf_source, & + undf_source, & + map_source) + + @assertEqual(expected_field, destination_field, tol) + + end subroutine test_all_with_extrap + +end module