anchor.diff:                      do i = (1+ecw), (domain%nx(1)-ecw)
                           depth = domain%U(i,j,STG) - domain%U(i,j,ELV)
                           if(depth > dry_depth_threshold) then
  -                            energy_potential = energy_potential + domain%area_cell_y(j) * &
  -                                (domain%U(i,j,STG) - domain%msl_linear)**2
  +                            energy_potential = energy_potential + real(domain%area_cell_y(j) * &
  +                                (domain%U(i,j,STG) - domain%msl_linear)**2, force_double)
                                   !depth * domain%U(i,j,STG) 
  -                            energy_kinetic = energy_kinetic + domain%area_cell_y(j) * &
  -                                (domain%U(i,j,UH)**2 + domain%U(i,j,VH)**2)/depth
  +                            energy_kinetic = energy_kinetic + real(domain%area_cell_y(j) * &
  +                                (domain%U(i,j,UH)**2 + domain%U(i,j,VH)**2)/depth, force_double)
                           end if 
                       end do
                   end do
  
positive.code:                 ! Do nothing
              else
                  ! Rescale energy statistics appropriately 
                  energy_potential = energy_potential * gravity * HALF_dp
                  energy_kinetic = energy_kinetic * HALF_dp
                  energy_total = energy_potential + energy_kinetic
  
                  write(domain%logfile_unit, *) 'Energy Total / rho: ', energy_total
  
negative.code:         if(report_energy_statistics) then
  
              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
                  !
  