anchor.diff:                      ! Nonlinear domain
                       domain%U(i,j,UH) = -(sqrt(gravity * max(domain%U(i+1, j, STG) - bc_values(ELV), 0.0_dp) ) - &
                                            sqrt(gravity * local_h ) ) * local_h
  -                    ! Thus reduces oscillations that can otherwise arise at the boundary
  -                    domain%U(i+1,j,UH) = 0.5_dp * (domain%U(i,j,UH) + domain%U(i+1,j,UH))
  +                    !! Thus reduces oscillations that can otherwise arise at the boundary
  +                    !domain%U(i+1,j,UH) = 0.5_dp * (domain%U(i,j,UH) + domain%U(i+1,j,UH))
                       ! Extrapolate VH
                       domain%U(i,j,VH) = domain%U(i+1,j,VH)
                   end if
  
positive.code:                     domain%U(i,j,UH) = &
                          (sqrt(gravity * max(domain%U(i-1, j, STG) - bc_values(ELV), 0.0_dp) ) - &
                           sqrt(gravity * local_h ) ) * local_h
                      ! Thus reduces oscillations that can otherwise arise at the boundary
                      domain%U(i-1,j,UH) = 0.5_dp * (domain%U(i,j,UH) + domain%U(i-1,j,UH))
                      ! Extrapolate VH
                      domain%U(i,j,VH) = domain%U(i-1,j,VH)
                  end if
  
negative.code:                     ! The nonlinear solvers need all boundary values updated
                      domain%U(i,j,VH) = -(sqrt(gravity * max(domain%U(i, j+1, STG) - bc_values(ELV), 0.0_dp) ) - &
                                           sqrt(gravity * local_h ) ) * local_h
                      !! Thus reduces oscillations that can otherwise arise at the boundary
                      !domain%U(i,j+1,VH) = 0.5_dp * (domain%U(i,j,VH) + domain%U(i,j+1,VH))
                      ! Extrapolate UH
                      domain%U(i,j,UH) = domain%U(i,j+1, UH)
                  end if
  