Beware Earth radius in single precision

CONTRIBUTORS.md

@@ -34,4 +34,5 @@
 | thomasmelvin        | Thomas Melvin      | Met Office                       | 2026-01-15 |
 | tinyendian          | Wolfgang Hayek     | Earth Sciences New Zealand       | 2026-02-02 |
 | DanStoneMO          | Daniel Stone       | Met Office                       | 2026-02-26 |
-| ericaneininger      | Erica Neininger    | Met Office                       | 2026-03-02 |
+| ericaneininger      | Erica Neininger    | Met Office                       | 2026-03-02 |
+| iboutle             | Ian Boutle         | Met Office                       | 2026-03-04 |

interfaces/physics_schemes_interface/source/kernel/bl_exp_kernel_mod.F90

@@ -517,8 +517,9 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
          bulk_cloud_fraction, rho_wet_tq, u_p, v_p, rhcpt, theta,            &
          p_rho_levels, exner_rho_levels, tgrad_bm, mix_len_tmp,              &
          exner_theta_levels,                                                 &
-         bulk_cf_conv, qcf_conv, r_rho_levels, visc_h, visc_m, rneutml_sq,    &
-         tnuc_new
+         bulk_cf_conv, qcf_conv, visc_h, visc_m, rneutml_sq, tnuc_new
+    ! Single precision is not accurate enough for distance from centre of planet
+    real(r_um), dimension(seg_len,1,nlayers) :: r_rho_levels
 
     ! profile field on boundary layer levels
     real(r_bl), dimension(seg_len,1,bl_levels) :: fqw, ftl, rhokh, bq_gb,    &
@@ -535,7 +536,9 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
 
     ! profile fields from level 0 upwards
     real(r_bl), dimension(seg_len,1,0:nlayers) :: p_theta_levels, etadot, w, &
-         q, qcl, qcf, r_theta_levels
+         q, qcl, qcf
+    ! Single precision is not accurate enough for distance from centre of planet
+    real(r_um), dimension(seg_len,1,0:nlayers) :: r_theta_levels
 
     ! profile fields with a hard-wired 2
     real(r_bl), dimension(seg_len,1,2,bl_levels) :: rad_hr, micro_tends
@@ -546,7 +549,7 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
          bl_type_1, bl_type_2, bl_type_3, bl_type_4, bl_type_5, bl_type_6,   &
          bl_type_7, uw0, vw0, zhnl, rhostar,                                 &
          recip_l_mo_sea, flandg, t1_sd, q1_sd, qcl_inv_top,                  &
-         fb_surf, rib_gb, z0m_eff_gb, zhsc, ustargbm, cos_theta_latitude,    &
+         fb_surf, rib_gb, z0m_eff_gb, zhsc, ustargbm,                        &
          max_diff, delta_smag, tnuc_nlcl_um
     real(r_um), dimension(seg_len,1) :: surf_dep_flux, zeroes
 
@@ -698,8 +701,8 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
         ! height of rho levels from centre of planet
         r_rho_levels(i,1,k) = height_w3(map_w3(1,i) + k-1) + planet_radius
         ! height of levels above surface
-        z_rho(i,1,k) = r_rho_levels(i,1,k)-r_theta_levels(i,1,0)
-        z_theta(i,1,k) = r_theta_levels(i,1,k)-r_theta_levels(i,1,0)
+        z_rho(i,1,k) = height_w3(map_w3(1,i) + k-1) - height_wth(map_wth(1,i))
+        z_theta(i,1,k) = height_wth(map_wth(1,i) + k) - height_wth(map_wth(1,i))
         ! water vapour mixing ratio
         q(i,1,k) = m_v_n(map_wth(1,i) + k)
         ! cloud liquid mixing ratio
@@ -736,11 +739,6 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
       end do
     end if
 
-    ! Set this to 1 to account for quasi-uniform grid
-    do i = 1, seg_len
-      cos_theta_latitude(i,1) = 1.0_r_um
-    end do
-
     do i = 1, seg_len
       ! surface pressure
       p_star(i,1) = p_theta_levels(i,1,0)
@@ -812,7 +810,7 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
 
     ! Calculate vertical differences
     do i = 1, seg_len
-      dzl_charney(i,1,1) = 2.0_r_bl * (r_theta_levels(i,1,1) - r_theta_levels(i,1,0))
+      dzl_charney(i,1,1) = 2.0_r_bl * z_theta(i,1,1)
       do k = 2, bl_levels
         dzl_charney(i,1,k) = dz_wth(map_wth(1,i) + k)
         rdz(i,1,k) = 1.0_r_bl/dz_wth(map_wth(1,i) + k-1)
@@ -856,7 +854,7 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
     !     IN model dimensions.
           , bl_levels, p_rho_levels, p_theta_levels(1,1,1)              &
           , exner_rho_levels, rho_wet, rho_wet_tq, z_theta, z_rho       &
-          , r_theta_levels, r_rho_levels                                &
+          , r_theta_levels                                              &
     !     IN Model switches
           , l_extra_call, no_cumulus                                    &
     !     IN cloud data
@@ -880,7 +878,7 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
     call bdy_expl2 (                                                           &
     ! IN values defining vertical grid of model atmosphere :
       bl_levels,p_theta_levels,land_field,land_index,                          &
-      r_theta_levels, r_rho_levels, cos_theta_latitude,                        &
+      r_theta_levels,                                                          &
     ! IN U, V and W momentum fields.
       u_p,v_p, u_0_p, v_0_p,                                                   &
     ! IN from other part of explicit boundary layer code
@@ -934,9 +932,9 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
       end do
       do k = 2, bl_levels-1
         do i = 1, seg_len
-          weight1 = r_theta_levels(i,1,k) - r_theta_levels(i,1,k-1)
-          weight2 = r_theta_levels(i,1,k) - r_rho_levels(i,1,k)
-          weight3 = r_rho_levels(i,1,k)   - r_theta_levels(i,1,k-1)
+          weight1 = z_theta(i,1,k) - z_theta(i,1,k-1)
+          weight2 = z_theta(i,1,k) - z_rho(i,1,k)
+          weight3 = z_rho(i,1,k)   - z_theta(i,1,k-1)
           rhokm_mix(i,1,k) = (weight3/weight1) * rhokm(i,1,k+1) &
                            + (weight2/weight1) * rhokm(i,1,k)
           ! Scale exchange coefficients by 1/dz factor, as is done for
@@ -947,8 +945,8 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
       end do
       k = bl_levels
       do i = 1, seg_len
-        weight1 = r_theta_levels(i,1,k) - r_theta_levels(i,1,k-1)
-        weight2 = r_theta_levels(i,1,k) - r_rho_levels(i,1,k)
+        weight1 = z_theta(i,1,k) - z_theta(i,1,k-1)
+        weight2 = z_theta(i,1,k) - z_rho(i,1,k)
         ! Assume rhokm(BL_LEVELS+1) is zero
         rhokm_mix(i,1,k) = (weight2/weight1) * rhokm(i,1,k)
         ! Scale exchange coefficients by 1/dz factor, as is done for
@@ -965,7 +963,7 @@ subroutine bl_exp_code(nlayers, seg_len,                      &
 
       call  tr_mix (                                                           &
            ! IN fields
-           real(r_theta_levels,r_um), real(r_rho_levels,r_um), pdims,          &
+           r_theta_levels, r_rho_levels, pdims,                                &
            bl_levels, alpha_cd,                                                &
            real(rhokm_mix(1:seg_len,1:1,2:bl_levels),r_um),                    &
            real(rhokm_mix(1:seg_len,1:1,1),r_um),                              &

interfaces/physics_schemes_interface/source/kernel/bl_imp2_kernel_mod.F90

@@ -375,7 +375,7 @@ subroutine bl_imp2_code(nlayers, seg_len,                   &
          cff_earliest, qt_force, tl_force, t_inc_pc2, q_inc_pc2, qcl_inc_pc2,&
          bcf_inc_pc2, cfl_inc_pc2, sskew, svar_turb, svar_bm, qcf_total,     &
          ri_bm, tgrad_in, mix_len_in, tau_dec_in, tau_hom_in, tau_mph_in,    &
-         wvar_in
+         wvar_in, z_rho
     real(r_bl), dimension(seg_len,1,nlayers) ::                              &
          r_rho_levels, rho_wet_tq
 
@@ -387,24 +387,16 @@ subroutine bl_imp2_code(nlayers, seg_len,                   &
 
     ! profile field on boundary layer levels
     real(r_bl), dimension(seg_len,1,bl_levels) :: fqw, ftl, rhokh,           &
-         bq_gb, bt_gb, dtrdz_charney_grid, rdz_charney_grid, rhokh_mix, qw,  &
-         tl, dqw, dtl, fqw_star, ftl_star
-
-    ! profile fields on u/v points and all levels
-    real(r_bl), dimension(seg_len,1,nlayers) :: r_u, r_v
+         bq_gb, bt_gb, dtrdz_charney_grid, rdz_charney_grid, qw,             &
+         tl, dqw, dtl, fqw_star, ftl_star, ct_ctq, dqw_nt, dtl_nt
 
     ! profile fields on u/v points and BL levels
-    real(r_bl), dimension(seg_len,1,bl_levels) :: taux, tauy,                &
-         rhokm_u, rhokm_v, dissip_u, dissip_v, taux_star, tauy_star, cq_cm_u,&
-         cq_cm_v, ct_ctq, dqw_nt, dtl_nt, du_star, dv_star
-
-    ! profile fields from level 2 upwards
-    real(r_bl), dimension(seg_len,1,2:bl_levels) :: rdz_u, rdz_v
+    real(r_bl), dimension(seg_len,1,bl_levels) ::                            &
+         dissip_u, dissip_v
 
     ! profile fields from level 0 upwards
     real(r_um), dimension(seg_len,1,0:nlayers) ::                            &
          p_theta_levels, p_rho_minus_one
-    real(r_bl), dimension(seg_len,1,0:nlayers) :: r_theta_levels
 
     ! single level real fields
     real(r_um), dimension(seg_len,1) ::                                      &
@@ -535,29 +527,26 @@ subroutine bl_imp2_code(nlayers, seg_len,                   &
          ! IN levels, switches
          bl_levels,  l_correct,                                              &
          ! IN data :
-         gamma1, gamma2, rhokm_u, rhokm_v,                                   &
-         rdz_charney_grid, r_rho_levels, dtrdz_charney_grid,rdz_u,rdz_v,     &
-         ct_ctq,cq_cm_u,cq_cm_v,dqw_nt,dtl_nt,                               &
+         gamma1, gamma2,                                                     &
+         rdz_charney_grid, r_rho_levels, dtrdz_charney_grid,                 &
+         ct_ctq,dqw_nt,dtl_nt,                                               &
          ! INOUT data :
-         qw,tl,fqw,ftl,taux,tauy,fqw_star,ftl_star,taux_star,tauy_star,      &
-         r_u,r_v,du_star,dv_star,dqw,dtl,rhokh,bl_diag,                      &
+         qw,tl,fqw,ftl,fqw_star,ftl_star,                                    &
+         dqw,dtl,rhokh,bl_diag,                                              &
          ! OUT data
-         t_latest,q_latest,rhokh_mix                                         &
+         t_latest,q_latest                                                   &
          )
 
     if (loop == 2) then
 
       do i = 1, seg_len
-        do k = 0, nlayers
-          ! height of theta levels from centre of planet
-          r_theta_levels(i,1,k) = height_wth(map_wth(1,i) + k) + planet_radius
-        end do
         zh(i,1) = zh_2d(map_2d(1,i))
       end do
       do k = 1, nlayers
         do i = 1, seg_len
           ! height of levels above surface
-          z_theta(i,1,k) = r_theta_levels(i,1,k)-r_theta_levels(i,1,0)
+          z_theta(i,1,k) = height_wth(map_wth(1,i) + k) - height_wth(map_wth(1,i))
+          z_rho(i,1,k) = height_w3(map_w3(1,i) + k-1) - height_wth(map_wth(1,i))
         end do
       end do
 
@@ -587,9 +576,9 @@ subroutine bl_imp2_code(nlayers, seg_len,                   &
           nblyr(i,1) = 1
           k = 1
 
-          weight1 = r_rho_levels(i,1,k+1) - r_theta_levels(i,1,0)
-          weight2 = r_theta_levels(i,1,k) - r_theta_levels(i,1,0)
-          weight3 = r_rho_levels(i,1,k+1) - r_theta_levels(i,1,k)
+          weight1 = z_rho(i,1,k+1)
+          weight2 = z_theta(i,1,k)
+          weight3 = z_rho(i,1,k+1) - z_theta(i,1,k)
           ftl_m = weight2 * ftl(i,1,k+1) + weight3 * ftl(i,1,k)
           fqw_m = weight2 * fqw(i,1,k+1) + weight3 * fqw(i,1,k)
           f_buoy_m = g*( bt_gb(i,1,k)*(ftl_m/cp) +                             &
@@ -601,15 +590,14 @@ subroutine bl_imp2_code(nlayers, seg_len,                   &
 
           ! Save level 1 heating increment for redistribution over
           ! boundary layer
-          fric_heating_blyr(i,1) = fric_heating_inc *                          &
-                       (r_rho_levels(i,1,2)-r_theta_levels(i,1,0))
+          fric_heating_blyr(i,1) = fric_heating_inc * z_rho(i,1,2)
         end do
 
         do k = 2, bl_levels-1
           do i = 1, seg_len
-            weight1 = r_rho_levels(i,1,k+1) - r_rho_levels(i,1,k)
-            weight2 = r_theta_levels(i,1,k) - r_rho_levels(i,1,k)
-            weight3 = r_rho_levels(i,1,k+1) - r_theta_levels(i,1,k)
+            weight1 = z_rho(i,1,k+1) - z_rho(i,1,k)
+            weight2 = z_theta(i,1,k) - z_rho(i,1,k)
+            weight3 = z_rho(i,1,k+1) - z_theta(i,1,k)
             ftl_m = weight2 * ftl(i,1,k+1) + weight3 * ftl(i,1,k)
             fqw_m = weight2 * fqw(i,1,k+1) + weight3 * fqw(i,1,k)
 
@@ -628,7 +616,7 @@ subroutine bl_imp2_code(nlayers, seg_len,                   &
               nblyr(i,1) = k
               fric_heating_blyr(i,1) = fric_heating_blyr(i,1) +                &
                                        fric_heating_incv(i,1) *                &
-                                  (r_rho_levels(i,1,k+1)-r_rho_levels(i,1,k))
+                                  (z_rho(i,1,k+1)-z_rho(i,1,k))
             else
               t_latest(i,1,k) = t_latest(i,1,k) + fric_heating_incv(i,1)
               if (bl_diag%l_dtfric) then
@@ -644,8 +632,7 @@ subroutine bl_imp2_code(nlayers, seg_len,                   &
         ! for lack of BL mixing of these increments
         !-----------------------------------------------------------------------
         do i = 1, seg_len
-          z_blyr = r_rho_levels(i,1,nblyr(i,1)+1)                              &
-                 - r_theta_levels(i,1,0)
+          z_blyr = z_rho(i,1,nblyr(i,1)+1)
           fric_heating_blyr(i,1) = fric_heating_blyr(i,1) / z_blyr
 
           do k = 1, nblyr(i,1)

interfaces/physics_schemes_interface/source/kernel/bl_imp_kernel_mod.F90

@@ -223,25 +223,15 @@ subroutine bl_imp_code(nlayers, seg_len,                   &
          dtrdz_charney_grid, rdz_charney_grid, qw, tl, dqw, dtl, ct_ctq, &
          dqw_nt, dtl_nt
 
-    ! profile fields on u/v points and all levels
-    real(r_bl), dimension(seg_len,1,nlayers) :: r_u, r_v
-
-    ! profile fields on u/v points and BL levels
-    real(r_bl), dimension(seg_len,1,bl_levels) :: taux, tauy, &
-         dtrdz_u, dtrdz_v, rhokm_u, rhokm_v, cq_cm_u, cq_cm_v
-
-    ! profile fields from level 2 upwards
-    real(r_bl), dimension(seg_len,1,2:bl_levels) :: rdz_u, rdz_v
-
     ! profile fields from level 0 upwards
     real(r_bl), dimension(seg_len,1,0:nlayers) :: q, qcl, qcf, r_theta_levels
 
     ! single level real fields
     real(r_bl), dimension(seg_len,1) :: gamma1, gamma2, ctctq1_1, &
-         dqw1_1, dtl1_1, cq_cm_u_1, du_1, cq_cm_v_1, dv_1
+         dqw1_1, dtl1_1
 
     ! single level integer fields
-    integer(i_um), dimension(seg_len,1) :: k_blend_tq, k_blend_uv
+    integer(i_um), dimension(seg_len,1) :: k_blend_tq
 
     ! parameters for new BL solver
     real(r_bl) :: pnonl,p1,p2
@@ -372,15 +362,13 @@ subroutine bl_imp_code(nlayers, seg_len,                   &
          bl_levels, l_correct,                                               &
          ! IN fields
          q, qcl, qcf, q_latest, qcl_latest, qcf_latest, t, t_latest,         &
-         dtrdz_charney_grid, dtrdz_u, dtrdz_v, rhokh, rhokm_u, rhokm_v,      &
-         rdz_charney_grid, rdz_u, rdz_v, gamma1, gamma2, real(alpha_cd,r_bl),&
-         r_u, r_v, r_theta_levels, r_rho_levels, k_blend_tq,                 &
-         k_blend_uv, k_blend_uv,                                             &
+         dtrdz_charney_grid, rhokh,                                          &
+         rdz_charney_grid, gamma1, gamma2, real(alpha_cd,r_bl),              &
+         r_theta_levels, r_rho_levels, k_blend_tq,                           &
          ! INOUT fields
-         fqw, ftl, taux, tauy, r_u, r_v, dqw, dtl,                           &
+         fqw, ftl, dqw, dtl,                                                 &
          ! OUT fields
-         dqw_nt, dtl_nt, qw, tl, ct_ctq, cq_cm_u, cq_cm_v, cq_cm_u_1,        &
-         cq_cm_v_1, du_1, dv_1, dqw1_1,dtl1_1,ctctq1_1                       &
+         dqw_nt, dtl_nt, qw, tl, ct_ctq, dqw1_1, dtl1_1, ctctq1_1            &
          )
 
     do k = 1, bl_levels

rose-stem/site/meto/groups/groups_lfric_atm.cylc

@@ -27,6 +27,7 @@
         "lfric_atm_nwp_gal9-pert-C12_azspice_gnu_fast-debug-32bit",
         "lfric_atm_nwp_casim-C12_azspice_gnu_fast-debug-32bit",
         "lfric_atm_nwp_coma9-C12_azspice_gnu_fast-debug-32bit",
+        "lfric_atm_nwp_comorph_dev-C12_azspice_gnu_fast-debug-32bit",
         "lfric_atm_nwp_gal9_mol-C12_azspice_gnu_fast-debug-32bit",
     ],
     "lfric_atm_nwp_azspice": [

rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_aquaplanet-C12_azspice_gnu_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         468AD911
-Inner product checksum theta =         4FD68987
-Inner product checksum u =         6A4E39D0
-Inner product checksum mr1 =         3F080026
-Inner product checksum mr2 =         369843BE
-Inner product checksum mr3 =         3416CDCB
-Inner product checksum mr4 =         36905047
+Inner product checksum rho =         468AD914
+Inner product checksum theta =         4FD68974
+Inner product checksum u =         6A4E33B7
+Inner product checksum mr1 =         3F080A3F
+Inner product checksum mr2 =         36906C45
+Inner product checksum mr3 =         341CB09C
+Inner product checksum mr4 =         3695C395
 Inner product checksum mr5 =                0
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_clim_gal9-C12_azspice_gnu_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         46D6671E
-Inner product checksum theta =         518E9CC1
-Inner product checksum u =         6AF4B3E5
-Inner product checksum mr1 =         3FD1DD2C
-Inner product checksum mr2 =         3744637E
-Inner product checksum mr3 =         3533FB34
-Inner product checksum mr4 =         36C1322A
+Inner product checksum rho =         46D65ECE
+Inner product checksum theta =         518E98C6
+Inner product checksum u =         6AF44C3A
+Inner product checksum mr1 =         3FD207A2
+Inner product checksum mr2 =         37337A50
+Inner product checksum mr3 =         352DAE88
+Inner product checksum mr4 =         36C10ECE
 Inner product checksum mr5 =                0
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_clim_gal9_chem-C12_azspice_gnu_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         46D66653
-Inner product checksum theta =         518E9990
-Inner product checksum u =         6AF43E1A
-Inner product checksum mr1 =         3FD1D04B
-Inner product checksum mr2 =         3745B69A
-Inner product checksum mr3 =         353019A0
-Inner product checksum mr4 =         36BDF494
+Inner product checksum rho =         46D663E2
+Inner product checksum theta =         518E98C1
+Inner product checksum u =         6AF4759A
+Inner product checksum mr1 =         3FD1CB60
+Inner product checksum mr2 =         373582D9
+Inner product checksum mr3 =         3529BC4E
+Inner product checksum mr4 =         36BC7ED8
 Inner product checksum mr5 =                0
 Inner product checksum mr6 =                0

rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_casim-C12_azspice_gnu_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         46D8578E
-Inner product checksum theta =         518BD702
-Inner product checksum u =         6A8A9FF5
-Inner product checksum mr1 =         3FCD79BE
-Inner product checksum mr2 =         38309A00
-Inner product checksum mr3 =         353B5ACC
-Inner product checksum mr4 =         36D670E9
-Inner product checksum mr5 =         2E296A8A
-Inner product checksum mr6 =         3534119C
+Inner product checksum rho =         46D858A9
+Inner product checksum theta =         518BD3D8
+Inner product checksum u =         6A8B6469
+Inner product checksum mr1 =         3FCCF28A
+Inner product checksum mr2 =         3833704D
+Inner product checksum mr3 =         3563EB58
+Inner product checksum mr4 =         36B2D559
+Inner product checksum mr5 =         2F158C5F
+Inner product checksum mr6 =         354EEAFB

rose-stem/site/meto/kgos/lfric_atm/azspice/checksum_lfric_atm_nwp_coma9-C12_azspice_gnu_fast-debug-32bit.txt

@@ -1,9 +1,9 @@
-Inner product checksum rho =         46D8451F
-Inner product checksum theta =         518BD5CF
-Inner product checksum u =         6A871C80
-Inner product checksum mr1 =         3FCF8CE1
-Inner product checksum mr2 =         37E92CD1
-Inner product checksum mr3 =         378E0B27
-Inner product checksum mr4 =         378C60AF
-Inner product checksum mr5 =         369846E2
+Inner product checksum rho =         46D84584
+Inner product checksum theta =         518BD574
+Inner product checksum u =         6A87BC4E
+Inner product checksum mr1 =         3FCFD56F
+Inner product checksum mr2 =         37F1D70A
+Inner product checksum mr3 =         3798DB1D
+Inner product checksum mr4 =         379054C3
+Inner product checksum mr5 =         369EAFD8
 Inner product checksum mr6 =                0