anchor.diff:                          end if
                       end do
                   end do
  +                !$OMP END DO
  +
               else if(.not. domain%is_staggered_grid) then
                   !
                   ! Compute energy statistics using non-staggered grid
                   !
   
                   ! Integrate over the model interior
  +                !$OMP DO
                   do j = 1 + ecw, (domain%nx(2) - ecw)
                       do i = (1+ecw), (domain%nx(1)-ecw)
                           depth = domain%U(i,j,STG) - domain%U(i,j,ELV)
  
positive.code:                         end if 
                      end do
                  end do
  
              end if
  
              ! In the case we use a "not-truely-linear" variant of the linear solver, the energy
              ! calculations do not really make sense
  
negative.code:         ! Optionally report integrated energy
          if(report_energy_statistics) then
  
              !$OMP PARALLEL DEFAULT(PRIVATE) SHARED(domain, dry_depth_threshold) REDUCTION(+:energy_potential), &
              !$OMP REDUCTION(+:energy_kinetic)
  
              ecw = domain%exterior_cells_width
              energy_potential = 0.0_force_double
              energy_kinetic = 0.0_force_double
  
  
              if(domain%is_staggered_grid .and. domain%linear_solver_is_truely_linear) then
                  !
  