anchor.diff:              end do density
           end do temperature
           
  +        ! comparison for this composition
           write(output_unit,'(a,f7.3)') 'rms(|diff|) = ',norm2(diff(:,:,k))/sqrt(real(Nrho*5,dp))
           loc = maxloc(abs(diff(:,:,k)))
           write(output_unit,'(a,f7.3,2(a,f6.2))') 'max(|diff|) = ',maxval(abs(diff(:,:,k))), &
  
positive.code:         &   '; lg(T) = ',lgT(loc(2))
      end do composition
  
      write(output_unit,'(/,a)') 'overall differences'
      write(output_unit,'(a,f7.3)') 'rms(|diff|) = ',norm2(diff)/sqrt(real(Nrho*NT*NX,dp))
      locz = maxloc(abs(diff))
  
negative.code:         call compute_composition_moments(2,chem_ids,Y,ionic,Xsum, &
              & ncharged, charged_ids, Yion, exclude_neutrons=.TRUE.)
          
          ! table header
          write(output_unit,'(/,2("X(",a,") = ",f3.1,tr4))') &
          &   (trim(name(j)),X(j),j=1,2)
          write(output_unit,'(2(a,f6.2,tr4))') '<Z> = ',ionic%Z,'<A> = ',ionic% A
  