anchor.diff:                  domain%U(i, j, VH) = domain%U(i, j, VH)  &
                       - dt_half_coriolis_jph(i) * (uh_i_j(i) + uh_i_jp1(i))
   
  +#endif
  +
  +#if defined(SPHERICAL)
  +                ! Spherical metric terms related to advection
  +
  +                ! duh/dt = tan(lat)/earth_radius * uvh
  +                domain%U(i,j,UH) = domain%U(i,j,UH) + dt_tanlat_uvh_on_radius_earth_iph_j(i)
  +
  +                ! dvh/dt = -tan(lat)/earth_radius * uuh
  +                domain%U(i,j,VH) = domain%U(i,j,VH) - dt_tanlat_uuh_on_radius_earth_i_jph(i)
   #endif
   
                   ! Add the semi-implicit friction terms
  
positive.code:         real(dp):: fr_limit_uh, fr_limit_vh
          ! d(stage)/dy; depth_{i, j+1/2}, depth_{j, i+1/2, j}
          real(dp):: dw_j(domain%nx(1)), h_jph_vec(domain%nx(1)), h_iph_vec(domain%nx(1))
          real(dp) :: h_iph_wet_strict(domain%nx(1)), h_jph_wet_strict(domain%nx(1))
          real(dp) :: h_iph_wet(domain%nx(1)), h_jph_wet(domain%nx(1))
  
  
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
  