anchor.diff:          call co_sum(global_energy_total_on_rho)
   #endif
   
  -        write(log_output_unit, "(A)"         ) ''
  -        write(log_output_unit, "(A)"         ) '-----------'
  -        write(log_output_unit, "(A)"         ) 'Global stage range (over all domains and images): '
  -        write(log_output_unit, "(A, ES25.12E3)") '        ', global_max_stage
  -        write(log_output_unit, "(A, ES25.12E3)") '        ', global_min_stage
  -        write(log_output_unit, "(2A)"        ) 'Global speed range (over all domains and images): '
  -        write(log_output_unit, "(A, ES25.12E3)") '        ', global_max_speed
  -        write(log_output_unit, "(A, ES25.12E3)") '        ', global_min_speed
  -        write(log_output_unit, "(2A)"        ) &
  +        write(md%log_output_unit, "(A)"         ) ''
  +        write(md%log_output_unit, "(A)"         ) '-----------'
  +        write(md%log_output_unit, "(A)"         ) 'Global stage range (over all domains and images): '
  +        write(md%log_output_unit, "(A, ES25.12E3)") '        ', global_max_stage
  +        write(md%log_output_unit, "(A, ES25.12E3)") '        ', global_min_stage
  +        write(md%log_output_unit, "(2A)"        ) 'Global speed range (over all domains and images): '
  +        write(md%log_output_unit, "(A, ES25.12E3)") '        ', global_max_speed
  +        write(md%log_output_unit, "(A, ES25.12E3)") '        ', global_min_speed
  +        write(md%log_output_unit, "(2A)"        ) &
               'Global energy-potential / rho (over all domains and images), zero when stage=domain%msl_linear: '
  -        write(log_output_unit, "(A, ES25.12E3)") '        ', global_energy_potential_on_rho
  -        write(log_output_unit, "(2A)"        ) 'Global energy-kinetic / rho (over all domains and images): '
  -        write(log_output_unit, "(A, ES25.12E3)") '        ', global_energy_kinetic_on_rho
  -        write(log_output_unit, "(2A)"        ) &
  +        write(md%log_output_unit, "(A, ES25.12E3)") '        ', global_energy_potential_on_rho
  +        write(md%log_output_unit, "(2A)"        ) 'Global energy-kinetic / rho (over all domains and images): '
  +        write(md%log_output_unit, "(A, ES25.12E3)") '        ', global_energy_kinetic_on_rho
  +        write(md%log_output_unit, "(2A)"        ) &
               'Global energy-total / rho (over all domains and images): '
  -        write(log_output_unit, "(A, ES25.12E3)") '        ', global_energy_total_on_rho
  -        write(log_output_unit, "(A)") '-----------'
  +        write(md%log_output_unit, "(A, ES25.12E3)") '        ', global_energy_total_on_rho
  +        write(md%log_output_unit, "(A)") '-----------'
           call md%report_mass_conservation_statistics()
   
           ! Unstable models may produce NaN energy -- this is useful to detect.
  
positive.code:     !
      subroutine convert_priority_domain_info_to_boxes(priority_domain_index, &
          priority_domain_image, halo_width, myimage, myindex, box_metadata, &
          xs, ys, periodic_xs, periodic_ys, max_parent_dx_ratio)
  
          integer(ip), intent(in) :: priority_domain_index(:,:), &
              priority_domain_image(:,:)
  
negative.code:                                     abs(cell_area1 - cell_area2), &
                                      sqrt(cell_area1)
  
                                  write(md_log_output_unit,*) 'ii: ', ii, ' di: ', di
                                  write(md_log_output_unit,*) 'xy: ', xs(xi), ys(xi)
                                  write(md_log_output_unit,*) 'bbox: '
                                  write(md_log_output_unit,*) bbox(1,:)
                                  write(md_log_output_unit,*) bbox(2,:)
                                  write(md_log_output_unit,*) bbox(3,:)
                                  write(md_log_output_unit,*) bbox(4,:)
                                  stop
  
                              end if
  