anchor.diff:  
           real(dp):: maxstage, minstage
           integer:: i,j, ecw, nx, ny
  -        real(dp):: dry_depth_threshold, energy_total, energy_potential, energy_kinetic
  +        real(dp):: dry_depth_threshold
  +        real(force_double) :: energy_total, energy_potential, energy_kinetic
           real(dp):: depth, depth_iplus, depth_jplus
           logical, parameter:: report_energy_statistics=.TRUE.
   
  
positive.code:                         ! Integrate over wet cells
                          if(domain%msl_linear > domain%U(i,j,ELV) + dry_depth_threshold) then
              
                              energy_potential = energy_potential + domain%area_cell_y(j) *&
                                  (domain%U(i,j,STG) - domain%msl_linear)**2
  
                              !
                              ! Kinetic energy integration needs to account for grid staggering -- cell
                              ! areas are constant for constant lat, but changing as lat changes.
                              !
                              depth_iplus = HALF_dp * (domain%U(i,j,STG) - domain%U(i,j,ELV) + &
                                  domain%U(i+1,j,STG) - domain%U(i+1,j,ELV))
                              if(depth_iplus > dry_depth_threshold) then
                                  energy_kinetic = energy_kinetic + domain%area_cell_y(j)*&
                                      (domain%U(i,j,UH)**2)/depth_iplus
                              end if
  
                              depth_jplus = HALF_dp * (domain%U(i,j,STG) - domain%U(i,j,ELV) + &
                                  domain%U(i,j+1,STG) - domain%U(i,j+1,ELV))
                              if(depth_jplus > dry_depth_threshold) then
                                  energy_kinetic = energy_kinetic + HALF_dp * (domain%area_cell_y(j) + domain%area_cell_y(j+1))*&
                                      (domain%U(i,j,VH)**2)/depth_jplus
                              end if
  
                          end if
  
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
                  !
  