Skip to content

Some kernels that write to W2 fields aren't safe for threading #195

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
tommbendall wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Some kernels that write to W2 fields aren't safe for threading #195

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
a606b9f
Use face selector for kernels which write to W2
tommbendall Jan 28, 2026
de5774e
fix to linear app
tommbendall Jan 28, 2026
3c23b9b
sign contributors on this branch
tommbendall Jan 28, 2026
c0492e9
update KGOs
tommbendall Jan 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CONTRIBUTORS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,3 +5,4 @@
| james-bruten-mo | James Bruten | Met Office | 2025-12-09 |
| jennyhickson | Jenny Hickson | Met Office | 2025-12-10 |
| mo-marqh | mark Hedley | Met Office | 2025-12-11 |
| tommbendall | Thomas Bendall | Met Office | 2026-01-28 |
18 changes: 14 additions & 4 deletions interfaces/physics_schemes_interface/source/algorithm/bl_exp_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,9 @@ module bl_exp_alg_mod
use fs_continuity_mod, only: W2, W3, Wtheta
use sci_geometric_constants_mod, only: get_height_fv, get_da_at_w2, &
get_delta_at_wtheta, &
get_dz_at_wtheta
get_dz_at_wtheta, &
get_face_selector_ew, &
get_face_selector_ns
use physics_constants_mod, only: get_max_diff, get_rdz_fd1, &
get_rdz_w3, get_dtrdz_wth
use map_fd_to_prognostics_alg_mod, only: set_wind
Expand Down Expand Up @@ -208,6 +210,9 @@ contains
type( field_type ), pointer :: w2_rmultiplicity
type( field_type ), pointer :: dz_wth

type( integer_field_type ), pointer :: face_selector_ew
type( integer_field_type ), pointer :: face_selector_ns

! local variables
type(mesh_type), pointer :: mesh
type( field_type ) :: dtl_mphys, dmt_mphys, ngstress_w2, &
Expand Down Expand Up @@ -321,6 +326,8 @@ contains

height_wth => get_height_fv( Wtheta, mesh%get_id() )
height_w3 => get_height_fv( W3, mesh%get_id() )
face_selector_ew => get_face_selector_ew( mesh%get_id() )
face_selector_ns => get_face_selector_ns( mesh%get_id() )
! 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() )
Expand Down Expand Up @@ -428,7 +435,8 @@ contains
call invoke(bl_exp_du_kernel_type(4_i_def, tau_w2, tau_land_w2, tau_ssi_w2,&
rhokm_w2, rdz_fd1, u_physics, &
surf_interp_w2, ngstress_w2, fd_tau_w2, &
sea_current_w2_ptr) )
sea_current_w2_ptr, &
face_selector_ew, face_selector_ns) )

if (l_calc_tau_at_p) then
call dz_wth%copy_field_properties(rdz_wth)
Expand All @@ -446,11 +454,13 @@ contains
bl_exp_du_kernel_type(1_i_def, taux, taux_land, taux_ssi, &
rhokm_bl, rdz_wth, u_in_w3, &
surf_interp, ngstress_bl, fd_taux, &
sea_u_current_ptr), &
sea_u_current_ptr, &
face_selector_ew, face_selector_ns), &
bl_exp_du_kernel_type(1_i_def, tauy, tauy_land, tauy_ssi, &
rhokm_bl, rdz_wth, v_in_w3, &
surf_interp, ngstress_bl, fd_tauy, &
sea_v_current_ptr) )
sea_v_current_ptr, &
face_selector_ew, face_selector_ns) )
end if

if ( subroutine_timers ) call timer('bl_exp_alg')
Expand Down
13 changes: 11 additions & 2 deletions interfaces/physics_schemes_interface/source/algorithm/bl_imp_alg_mod.x90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,9 @@ module bl_imp_alg_mod
use integer_field_mod, only: integer_field_type
use field_collection_mod, only: field_collection_type
use fs_continuity_mod, only: W1, W2, W3, Wtheta
use sci_geometric_constants_mod, only: get_height_fv, get_da_at_w2
use sci_geometric_constants_mod, only: get_height_fv, get_da_at_w2, &
get_face_selector_ew, &
get_face_selector_ns
use physics_constants_mod, only: get_rdz_fd1, get_dtrdz_fd2, &
get_rdz_w3
use mr_indices_mod, only: nummr, imr_v, imr_cl, imr_ci, &
Expand Down Expand Up @@ -198,6 +200,9 @@ contains
type( field_type ), pointer :: dtrdz => null()
type( field_type ), pointer :: rdz_w3 => null()

type( integer_field_type ), pointer :: face_selector_ew
type( integer_field_type ), pointer :: face_selector_ns

! Local variables
type( field_type ) :: u_physics_star
type( field_type ) :: du_conv_w2, dw_conv
Expand Down Expand Up @@ -312,6 +317,9 @@ contains

mesh => zh%get_mesh()

face_selector_ew => get_face_selector_ew(mesh%get_id())
face_selector_ns => get_face_selector_ns(mesh%get_id())

! Set up the BL increments
call theta%copy_field_properties(dtheta_bl)

Expand Down Expand Up @@ -367,7 +375,8 @@ contains
rhokm_w2, rdz, dtrdz, wetrho_in_w2, &
u_physics, u_physics_star, &
surf_interp_w2, du_conv_w2, dw_bl, dA, &
height_w1, height_w2), &
height_w1, height_w2, &
face_selector_ew, face_selector_ns), &
! Map the wind back into dynamics quantities
inc_X_times_Y(du_bl, dA), &
inc_X_times_Y(dissip, dA) )
Expand Down
221 changes: 122 additions & 99 deletions interfaces/physics_schemes_interface/source/kernel/bl_exp_du_kernel_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,12 +13,14 @@ module bl_exp_du_kernel_mod
ANY_SPACE_1, ANY_SPACE_2, &
ANY_SPACE_3, ANY_SPACE_4, &
ANY_SPACE_5, ANY_SPACE_6, &
ANY_DISCONTINUOUS_SPACE_3, &
GH_REAL, GH_WRITE, GH_SCALAR, GH_INTEGER
use constants_mod, only: r_def, i_def, r_bl
use fs_continuity_mod, only: W1, W2
use kernel_mod, only: kernel_type
use nlsizes_namelist_mod, only: bl_levels
use jules_surface_config_mod, only: formdrag, formdrag_dist_drag
use reference_element_mod, only: N

implicit none

Expand All @@ -30,19 +32,21 @@ module bl_exp_du_kernel_mod
!> Kernel metadata type.
type, public, extends(kernel_type) :: bl_exp_du_kernel_type
private
type(arg_type) :: meta_args(11) = (/ &
arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! nfaces
arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! tau w3/w2
arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_land 2d w3/w2
arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_ssi 2d w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! rhokm wth/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_4), &! rdz wth/fd1
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! u_physics w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_5), &! surf_interp mult w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! ngstress wth/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! fd_tau w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_6) &! sea_current mix w3/w2
/)
type(arg_type) :: meta_args(13) = (/ &
arg_type(GH_SCALAR, GH_INTEGER, GH_READ), &! nfaces
arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_1), &! tau w3/w2
arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_land 2d w3/w2
arg_type(GH_FIELD, GH_REAL, GH_WRITE, ANY_SPACE_2), &! tau_ssi 2d w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! rhokm wth/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_4), &! rdz wth/fd1
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! u_physics w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_5), &! surf_interp mult w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_3), &! ngstress wth/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_1), &! fd_tau w3/w2
arg_type(GH_FIELD, GH_REAL, GH_READ, ANY_SPACE_6), &! sea_current mix w3/w2
arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3), &! face_selector_ew
arg_type(GH_FIELD, GH_INTEGER, GH_READ, ANY_DISCONTINUOUS_SPACE_3) &! face_selector_ns
/)
integer :: operates_on = CELL_COLUMN
contains
procedure, nopass :: bl_exp_du_code
Expand All @@ -69,6 +73,10 @@ module bl_exp_du_kernel_mod
!> @param[in] ngstress NG stress function
!> @param[in] fd_tau Stress from turbulent form-drag
!> @param[in] sea_current Ocean surface current
!> @param[in] face_selector_ew 2D field indicating which W/E faces
!! to loop over in this column
!> @param[in] face_selector_ns 2D field indicating which N/S faces
!! to loop over in this column
!> @param[in] ndf_half Number of DOFs per cell for half levels
!> @param[in] undf_half Number of unique DOFs for half levels
!> @param[in] map_half Dofmap for the cell at the base of the column
Expand All @@ -87,36 +95,44 @@ module bl_exp_du_kernel_mod
!> @param[in] ndf_curr Number of DOFs per cell for current space
!> @param[in] undf_curr Number of unique DOFs for current space
!> @param[in] map_curr Dofmap for the cell at the base of the column
subroutine bl_exp_du_code(nlayers, &
nfaces, &
tau, &
tau_land, &
tau_ssi, &
rhokm, &
rdz, &
u_physics, &
surf_interp, &
ngstress, &
fd_tau, &
sea_current, &
ndf_half, &
undf_half, &
map_half, &
ndf_2d, &
undf_2d, &
map_2d, &
ndf_full, &
undf_full, &
map_full, &
ndf_rdz, &
undf_rdz, &
map_rdz, &
ndf_surf, &
undf_surf, &
map_surf, &
ndf_curr, &
undf_curr, &
map_curr)
!> @param[in] ndf_w3_2d Num of DoFs for 2D W3 per cell
!> @param[in] undf_w3_2d Num of DoFs for this partition for 2D W3
!> @param[in] map_w3_2d Map for 2D W3
subroutine bl_exp_du_code(nlayers, &
nfaces, &
tau, &
tau_land, &
tau_ssi, &
rhokm, &
rdz, &
u_physics, &
surf_interp, &
ngstress, &
fd_tau, &
sea_current, &
face_selector_ew, &
face_selector_ns, &
ndf_half, &
undf_half, &
map_half, &
ndf_2d, &
undf_2d, &
map_2d, &
ndf_full, &
undf_full, &
map_full, &
ndf_rdz, &
undf_rdz, &
map_rdz, &
ndf_surf, &
undf_surf, &
map_surf, &
ndf_curr, &
undf_curr, &
map_curr, &
ndf_w3_2d, &
undf_w3_2d, &
map_w3_2d )

!---------------------------------------
! UM modules containing switches or global constants
Expand All @@ -141,6 +157,8 @@ subroutine bl_exp_du_code(nlayers, &
integer(kind=i_def), intent(in) :: map_full(ndf_full)
integer(kind=i_def), intent(in) :: ndf_surf, undf_surf
integer(kind=i_def), intent(in) :: map_surf(ndf_surf)
integer(kind=i_def), intent(in) :: ndf_w3_2d, undf_w3_2d
integer(kind=i_def), intent(in) :: map_w3_2d(ndf_w3_2d)

real(kind=r_def), dimension(undf_half), intent(inout) :: tau
real(kind=r_def), dimension(undf_2d), intent(inout) :: tau_land, tau_ssi
Expand All @@ -152,70 +170,75 @@ subroutine bl_exp_du_code(nlayers, &
real(kind=r_def), dimension(undf_surf), intent(in) :: surf_interp
real(kind=r_def), dimension(undf_curr), intent(in) :: sea_current

integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ew
integer(kind=i_def), dimension(undf_w3_2d), intent(in) :: face_selector_ns

! Internal variables
integer(kind=i_def) :: df, k
integer(kind=i_def) :: df, k, j, total_faces
real(kind=r_bl) :: rhokm_land, rhokm_ssi, fland, flandfac, fseafac
real(kind=r_bl) :: zh_nonloc(1)
real(kind=r_bl), dimension(0:bl_levels-1) :: tau_grad, tau_non_grad, u_sp, &
rdz_sp, rhokm_sp, ngstress_sp, tau_sp, fd_tau_sp

total_faces = MIN( &
nfaces, face_selector_ew(map_w3_2d(1)) + face_selector_ns(map_w3_2d(1)) &
)

! loop over all faces of the cell
do df = 1,nfaces

! Only calculate face if it's not already been done
if (tau(map_half(df)) == 0.0_r_def) then

!================================================================
! In the UM this happens in bdy_expl3
!================================================================
fland = surf_interp(map_surf(df) + 0)
rhokm_land = surf_interp(map_surf(df) + 1)
rhokm_ssi = surf_interp(map_surf(df) + 2)
flandfac = surf_interp(map_surf(df) + 3)
fseafac = surf_interp(map_surf(df) + 4)
zh_nonloc = surf_interp(map_surf(df) + 5)

do k = 0, bl_levels-1
u_sp(k) = u_physics(map_half(df)+k)
rdz_sp(k) = rdz(map_rdz(df)+k)
rhokm_sp(k) = rhokm(map_full(df)+k)
ngstress_sp(k) = ngstress(map_full(df)+k)
fd_tau_sp(k) = fd_tau(map_half(df)+k)
tau_sp(k) = tau(map_half(df)+k)
end do

tau_land(map_2d(df)) = rhokm_land * u_sp(0) * flandfac

! If sea surface current (sea_current) has not been obtained
! from coupling fields or from an ancillary file then it will simply
! contain uniform zeros.
tau_ssi(map_2d(df)) = rhokm_ssi * &
(u_sp(0) - sea_current(map_curr(df))) * fseafac
tau_sp(0) = fland * tau_land(map_2d(df)) &
+ (1.0_r_bl - fland) * tau_ssi(map_2d(df))

if (formdrag == formdrag_dist_drag) then
if (fland > 0.0_r_bl) then
tau_land(map_2d(df)) = tau_land(map_2d(df)) &
+ fd_tau(map_half(df)) / fland
end if
do j = 1, total_faces
df = j
if (j == 3 .and. face_selector_ns(map_w3_2d(1)) == 2 &
.and. face_selector_ew(map_w3_2d(1)) == 1) df = N

!================================================================
! In the UM this happens in bdy_expl3
!================================================================
fland = surf_interp(map_surf(df) + 0)
rhokm_land = surf_interp(map_surf(df) + 1)
rhokm_ssi = surf_interp(map_surf(df) + 2)
flandfac = surf_interp(map_surf(df) + 3)
fseafac = surf_interp(map_surf(df) + 4)
zh_nonloc = surf_interp(map_surf(df) + 5)

do k = 0, bl_levels-1
u_sp(k) = u_physics(map_half(df)+k)
rdz_sp(k) = rdz(map_rdz(df)+k)
rhokm_sp(k) = rhokm(map_full(df)+k)
ngstress_sp(k) = ngstress(map_full(df)+k)
fd_tau_sp(k) = fd_tau(map_half(df)+k)
tau_sp(k) = tau(map_half(df)+k)
end do

tau_land(map_2d(df)) = rhokm_land * u_sp(0) * flandfac

! If sea surface current (sea_current) has not been obtained
! from coupling fields or from an ancillary file then it will simply
! contain uniform zeros.
tau_ssi(map_2d(df)) = rhokm_ssi * &
(u_sp(0) - sea_current(map_curr(df))) * fseafac
tau_sp(0) = fland * tau_land(map_2d(df)) &
+ (1.0_r_bl - fland) * tau_ssi(map_2d(df))

if (formdrag == formdrag_dist_drag) then
if (fland > 0.0_r_bl) then
tau_land(map_2d(df)) = tau_land(map_2d(df)) &
+ fd_tau(map_half(df)) / fland
end if

call ex_flux_uv(pdims, pdims, pdims, bl_levels, &
u_sp(0:bl_levels-1), &
zh_nonloc, &
rdz_sp(1:bl_levels-1), &
rhokm_sp(0:bl_levels-1), &
ngstress_sp(1:bl_levels-1), &
fd_tau_sp(0:bl_levels-1), &
tau_sp(0:bl_levels-1), &
tau_grad(0:bl_levels-1), tau_non_grad(0:bl_levels-1))

do k = 0, bl_levels-1
tau(map_half(df)+k) = tau_sp(k)
end do

end if ! this face needs calculating
end if

call ex_flux_uv(pdims, pdims, pdims, bl_levels, &
u_sp(0:bl_levels-1), &
zh_nonloc, &
rdz_sp(1:bl_levels-1), &
rhokm_sp(0:bl_levels-1), &
ngstress_sp(1:bl_levels-1), &
fd_tau_sp(0:bl_levels-1), &
tau_sp(0:bl_levels-1), &
tau_grad(0:bl_levels-1), tau_non_grad(0:bl_levels-1))

do k = 0, bl_levels-1
tau(map_half(df)+k) = tau_sp(k)
end do

end do ! loop over df

Expand Down
Loading
Loading