anchor.diff:      use cliffs_tolkova_mod, only: cliffs, setSSLim
       use extrapolation_limiting_mod, only: limited_gradient_dx_vectorized
   
  -#ifdef SPHERICAL    
  +#ifdef SPHERICAL
       ! Compile with -DSPHERICAL to get the code to run in spherical coordinates
       use spherical_mod, only: area_lonlat_rectangle, deg2rad, earth_angular_freq
       use global_mod, only: radius_earth
  
positive.code:             end do
              !$OMP END DO
              !$OMP END PARALLEL
              
              ! NOTE: We don't need a flux for elevation 
              allocate(domain%flux_NS(nx, ny+1_ip, STG:VH))
              !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(domain, ny)
              !$OMP DO SCHEDULE(STATIC)
  
negative.code:         if(present(max_dt_out)) max_dt_out = max_dt_out_local
  
          !
          ! Spatially integrate boundary fluxes. Order is N, E, S, W.
          ! Note we integrate in the INTERIOR (i.e. ignoring the outer n_ext rows/columns which can be updated by boundary-condition
          ! routines, then integrate the boundary of what remains).
          !
          n_ext = domain%exterior_cells_width
          nx = domain%nx(1)
          ny = domain%nx(2)
  
          ! Outward boundary flux over the north
          domain%boundary_flux_store(1) = sum(domain%flux_NS( (1+n_ext):(nx-n_ext), ny+1-n_ext, STG))
  