anchor.diff:                  merge(ONE_dp, ZERO_dp, h_jph_wet_strict(xL:(xU-1)) > ZERO_dp)
   #endif
   
  +#if defined(SPHERICAL)
  +            ! Spherical metric terms related to advection. These have a form like
  +            ! (+ tan(lat)/R * uvh , -tan(lat)/R * uuh ) 
  +            ! although here we include dt, and exclude the leading sign (treated later).
  +            ! In practice these terms are often small and ignored in many papers.
  +            do i = xL, xU-1
  +
  +                ! Inverse depth at i+1/2, j
  +                inv_depth_local = ZERO_dp
  +                if(h_iph_wet_strict(i) > ZERO_dp) inv_depth_local = ONE_dp/h_iph_vec(i)
  +
  +                ! dt * tan(lat)/R * uvh at i+1/2, j
  +                dt_tanlat_uvh_on_radius_earth_iph_j(i) = dt * &
  +                    domain%tanlat_on_radius_earth(j) * &
  +                    domain%U(i,j,UH) * & ! uh
  +                    0.5_dp * (vh_iph_jmh(i) + vh_iph_jph(i)) * & ! vh
  +                    inv_depth_local ! 1/h
  +
  +                ! Inverse depth at i, j+1/2
  +                inv_depth_local = ZERO_dp
  +                if(h_jph_wet_strict(i) > ZERO_dp) inv_depth_local = ONE_dp/h_jph_vec(i)
  +
  +                ! dt * tan(lat)/R * uuh at i, j+1/2
  +                dt_tanlat_uuh_on_radius_earth_i_jph(i) = dt * &
  +                    0.5_dp * (domain%tanlat_on_radius_earth(j) + domain%tanlat_on_radius_earth(j+1)) * &
  +                    0.5_dp * (uh_i_jp1(i) + uh_i_j(i)) * & ! uh
  +                    0.5_dp * (uh_i_jp1(i) + uh_i_j(i)) * & ! uh
  +                    inv_depth_local ! 1/h
  +            end do
  +#endif
  +
               ! Compute a multiplier that applies a semi-implicit treatment of manning friction.
               !
               ! This is equivalent to adding a term "g * depth * friction_slope" to the equations.
  
positive.code:                 domain%U(i, j, VH) = domain%U(i, j, VH)  &
                      - dt_half_coriolis_jph(i) * (uh_i_j(i) + uh_i_jp1(i))
  
  #endif
  
                  ! Add the semi-implicit friction terms
  
negative.code:         uh_i_jp1_max_loop_index(xL:(xU-1)) = ZERO_dp
          dt_half_coriolis(xL:(xU-1)) = ZERO_dp
          dt_half_coriolis_jph(xL:(xU-1)) = ZERO_dp
          dt_tanlat_uvh_on_radius_earth_iph_j(xL:(xU-1)) = ZERO_dp
          dt_tanlat_uuh_on_radius_earth_i_jph(xL:(xU-1)) = ZERO_dp
  
          !
          ! Before starting the loop, get the value of VH at
  