diff --git a/CONTRIBUTORS.md b/CONTRIBUTORS.md
index 377f3de06..9a9ccabc8 100644
--- a/CONTRIBUTORS.md
+++ b/CONTRIBUTORS.md
@@ -28,3 +28,4 @@
 | mo-lucy-gordon   | Lucy Gordon            | Met Office            | 2026-03-18 |
 | shreybh1         | Shrey Bhardwaj         | Met Office            | 2026-03-26 |
 | hiker            | Joerg Henrichs         | Bureau of Meteorology | 2026-02-11 |
+| sebvdb           | Sebastiaan van de Bund | EPCC & University of Edinburgh | 2026-04-08 |
diff --git a/components/science/source/kernel/algebra/matrix_vector_kernel_mod.F90 b/components/science/source/kernel/algebra/matrix_vector_kernel_mod.F90
index c53583587..1ed386378 100644
--- a/components/science/source/kernel/algebra/matrix_vector_kernel_mod.F90
+++ b/components/science/source/kernel/algebra/matrix_vector_kernel_mod.F90
@@ -12,7 +12,7 @@ module matrix_vector_kernel_mod
                             GH_REAL, GH_READ, GH_INC, &
                             ANY_SPACE_1, ANY_SPACE_2, &
                             CELL_COLUMN
-  use constants_mod, only : i_def, r_single, r_double
+  use constants_mod, only : i_def, r_single, r_double, BLOCK_SIZE
   use kernel_mod,    only : kernel_type
 
   implicit none
@@ -87,15 +87,20 @@ subroutine matrix_vector_code_r_single(cell,              &
     real(kind=r_single), dimension(ncell_3d,ndf1,ndf2), intent(in)    :: matrix
 
     ! Internal variables
-    integer(kind=i_def) :: df, ik, df2, i1, i2, nl
+    integer(kind=i_def) :: df, ik, df2, i1, i2, nl, k, kk, kend
 
     nl = nlayers-1
-    ik = (cell-1)*nlayers +1
-    do df2 = 1, ndf2
-      i2 = map2(df2)
-      do df = 1, ndf1
-        i1 = map1(df)
-        lhs(i1:i1+nl) = lhs(i1:i1+nl) + matrix(ik:ik+nl,df,df2)*x(i2:i2+nl)
+    do k = 0, nl, BLOCK_SIZE
+      ik = (cell-1)*nlayers + 1
+      kend = min(BLOCK_SIZE-1, nl-k)
+      do df2 = 1, ndf2
+        i2 = map2(df2)
+        do df = 1, ndf1
+          i1 = map1(df)
+          do kk = 0, kend
+            lhs(i1+kk+k) = lhs(i1+kk+k) + matrix(ik+kk+k,df,df2)*x(i2+kk+k)
+          end do
+        end do
       end do
     end do
 
@@ -125,15 +130,21 @@ subroutine matrix_vector_code_r_double(cell,              &
     real(kind=r_double), dimension(undf1),              intent(inout) :: lhs
     real(kind=r_double), dimension(ncell_3d,ndf1,ndf2), intent(in)    :: matrix
 
-    integer(kind=i_def) :: df, ik, df2, i1, i2, nl
+    ! Internal variables
+    integer(kind=i_def) :: df, ik, df2, i1, i2, nl, k, kk, kend
 
     nl = nlayers-1
-    ik = (cell-1)*nlayers +1
-    do df2 = 1, ndf2
-      i2 = map2(df2)
-      do df = 1, ndf1
-        i1 = map1(df)
-        lhs(i1:i1+nl) = lhs(i1:i1+nl) + matrix(ik:ik+nl,df,df2)*x(i2:i2+nl)
+    do k = 0, nl, BLOCK_SIZE
+      ik = (cell-1)*nlayers + 1
+      kend = min(BLOCK_SIZE-1, nl-k)
+      do df2 = 1, ndf2
+        i2 = map2(df2)
+        do df = 1, ndf1
+          i1 = map1(df)
+          do kk = 0, kend
+            lhs(i1+kk+k) = lhs(i1+kk+k) + matrix(ik+kk+k,df,df2)*x(i2+kk+k)
+          end do
+        end do
       end do
     end do
 
diff --git a/infrastructure/source/utilities/constants_mod.F90 b/infrastructure/source/utilities/constants_mod.F90
index 85653641b..8affddd3d 100644
--- a/infrastructure/source/utilities/constants_mod.F90
+++ b/infrastructure/source/utilities/constants_mod.F90
@@ -33,7 +33,7 @@ module constants_mod
             LARGE_REAL_NEGATIVE, LARGE_DP_NEGATIVE, LARGE_REAL_POSITIVE, &
             PI, degrees_to_radians, radians_to_degrees,                  &
             cache_block, PRECISION_REAL, PRECISION_R_SOLVER,             &
-            PRECISION_R_TRAN, EPS_R_TRAN, default_halo_depth
+            PRECISION_R_TRAN, EPS_R_TRAN, default_halo_depth, BLOCK_SIZE
 
   ! Define default application-defined kinds for all intrinsic data types
 
@@ -197,4 +197,9 @@ module constants_mod
   integer(i_def), parameter :: default_halo_depth = 1 !< Default halo depth for fields and operators
   !> @}
 
+  !> @name Vertical loop block size
+  !> @{
+  integer(i_def), parameter :: BLOCK_SIZE = 1024 !< Vertical OpenMP loop block size
+  !> @}
+
 end module constants_mod