anchor.diff:          case('nearest')
               ! Find the 'nearest' match in parent domain. This might mean we reduce the requested size of the domain
               domain%lower_left = parent_domain%lower_left + &
  -                nint((lower_left - parent_domain%lower_left)/parent_domain_dx)*parent_domain_dx 
  +                nint((lower_left - parent_domain%lower_left)/parent_domain_dx)*parent_domain_dx
   
               ur = parent_domain%lower_left + &
                   nint((upper_right - parent_domain%lower_left)/parent_domain_dx)*parent_domain_dx
  
positive.code:                             0.5_dp * (domain%area_cell_y(j) + domain%area_cell_y(j+1)) * &
                                  domain%U(i,j,VH) * domain%U(i,j,VH) * depth_N_inv  ), &
                              force_double)
   
                      else
                          ! Co-loated grid
                          depth_C = domain%U(i,j,STG) - domain%U(i,j,ELV)
  
negative.code:         !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(domain, n_ext) REDUCTION(+:domain_volume)
          domain_volume = ZERO_dp
          !$OMP DO SCHEDULE(STATIC)
          do j = (1+n_ext), (domain%nx(2) - n_ext)
              ! Be careful about precision
              local_sum = sum(&
                  real(domain%U((1+n_ext):(domain%nx(1)-n_ext),j,STG), force_double) - &
  