! Calculate, for each dataset in XSCALE.HKL file, its delta-CC1/2
! works also for XDS_ASCII.HKL
!
! version 0.1 Kay Diederichs 12.3.2015
! compile: ifort -check bounds -init=huge,snan,arrays -traceback -warn -fpe0 -static-intel -sox getrc.o resolution.o symops.o hash.o cc_mod.o qsortfix.o xdscc12.f90
!       or gfortran -O2 -C -Wall symops.f90 getrc.f90 resolution.f90 hash.f90 cc_mod.f90 xdscc12.f90 -static-libgfortran -static-libgcc ! for Mac
! with previous compilation like: ifort -c symops.f90 getrc.f90 resolution.f90 hash.f90 cc_mod.f90 qsortfix.f90
!
! bugfix 6.9.2018: dividing by (nref-1) should be division by nref (12 occurrences)
! 6.9.18 enforce safe ATANH() argument range . Probably a better solution would be to do no Fisher transform if there are only two reflections.
! fix 5.4.2019: instead of n/(n-1) unbiasing factor for epsilon, use 1/(1-sumsig2/sumsig**2). 
! See https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Reliability_weights 
! Sep 15,2019 add WEIGHT= to XSCALE.INP.rename_me
! add option to sort on anom CC1/2
! 9.10.20 fix floating divide by zero in line 841 und 697
! 8.12.20 fix sed command 
! 29.10.21 fix same floating overflow problem as 10.12.2020: use REAL(8) xx to store sigma**2
! 17.6.22 implement -R (intensities from reference) 
MODULE xscaleitems_mod
      INTEGER :: nset=0
      REAL resmin,resmax
      CHARACTER*132, ALLOCATABLE :: inputfile(:)
      LOGICAL :: is_xds=.FALSE.
      REAL, ALLOCATABLE :: save_a(:,:),save_d(:,:)
contains
      SUBROUTINE get_nset_reso(name,batch,nbatch)
      USE resolution_mod
      IMPLICIT NONE
      CHARACTER*(*), INTENT(IN) :: name
      REAL, INTENT(IN) :: batch
      INTEGER, INTENT(OUT) :: nbatch
      CHARACTER line*80
      INTEGER :: ier,unit,iset,ih,ik,il
      REAL :: s2,resol,dum(6),osci

 1000 FORMAT(a)
      
      nbatch=1
      OPEN(NEWUNIT=unit,FILE=name,STATUS='old',ACTION='READ',iostat=ier)
      IF (ier/=0) THEN
        PRINT*,'ier=',ier
        STOP 'could not open file'
      END IF
! Open input file; find number of datasets, and allocate array inputfile
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier,' line=',line(:LEN_TRIM(line))
          STOP 'read error in xds_file'
        END IF
        IF (line(1:7)=='! ISET=') then
          READ(line(8:),*) iset
          nset=MAX(nset,iset)
        ELSE IF (line(1:26).eq.'!OUTPUT_FILE=XDS_ASCII.HKL') then 
          is_xds=.TRUE.
        ELSE IF (line(1:21).eq.'!UNIT_CELL_CONSTANTS=') then 
          READ(line(22:),*) a,b,c,alf,bet,gam
        ELSE IF (line(1:19).eq.'!OSCILLATION_RANGE=') then 
          READ(line(20:),*) osci
          WRITE(*,'(a)') TRIM(line)
          nbatch=NINT(batch/osci)
          IF (nbatch<1) nbatch=1
          IF (nbatch>1) WRITE(*,*) nbatch,' frames will be batched' 
        ELSE IF (INDEX(line,'MERGE=TRUE')>0) THEN
          STOP 'error: XSCALE file has MERGE=TRUE - this program needs MERGE=FALSE !'
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
      ALLOCATE(inputfile(nset))
      resmin=HUGE(resmin)
      resmax=-1
! find min and max resolution
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          STOP 'missing !END_OF_DATA line at end of reflection file!'
        END IF
        IF (line(1:12) == '!END_OF_DATA') EXIT
        IF (is_xds) THEN
          READ(line,*,iostat=ier) ih,ik,il,dum(:5)
          IF (ier/=0) THEN
            WRITE(*,*)'bad line in reflection file:',TRIM(line)
            STOP 'error'
          END IF
          nset=MAX(nset,INT(dum(5)/nbatch)+1)
        ELSE
          READ(line,*,iostat=ier) ih,ik,il,dum,iset
          IF (ier/=0) THEN
            WRITE(*,*)'bad line in reflection file:',TRIM(line)
            STOP 'error'
          END IF
        END IF
        CALL getres(ih,ik,il,s2,resol)
        IF (resol<resmin) resmin=resol
        IF (resol>resmax) resmax=resol      
      END DO      
      CLOSE(unit)
      END SUBROUTINE get_nset_reso
END MODULE xscaleitems_mod

MODULE reference_mod
      REAL, ALLOCATABLE :: ref(:,:,:),refminus(:,:,:)
CONTAINS
      SUBROUTINE read_reference(reference,reslim1,reslim2,anoflag,squareflag)
      USE ISO_FORTRAN_ENV  ! to get access to iostat_end
      USE symops_mod
      USE resolution_mod
      CHARACTER (LEN=132), INTENT(IN) :: reference
      LOGICAL, INTENT(IN) :: anoflag,squareflag
      CHARACTER (LEN=132) line
      INTEGER :: unit,ier,ih,ik,il,jh,jk,jl,ihkl,iop
      INTEGER :: hmin= HUGE(hmin),kmin= HUGE(hmin),lmin= HUGE(hmin),  &
                 hmax=-HUGE(hmin),kmax=-HUGE(hmin),lmax=-HUGE(hmin),n=0
      REAL s2,resol,reslim1,reslim2,amplitude,friedel
      
      OPEN(NEWUNIT=unit,FILE=reference,STATUS='old',ACTION='READ',iostat=ier)
      IF (ier/=0) THEN
        PRINT*,'ier=',ier
        STOP 'could not open reference file'
      END IF
! find min and max indices
      DO
        READ(unit,'(a)',iostat=ier) line
        IF (ier == iostat_end) EXIT
        READ(line,*,iostat=ier) ih,ik,il,amplitude
        IF (ier /= 0) THEN
          WRITE(*,*)'ignoring: ',line(:LEN_TRIM(reference))
          CYCLE
        END IF
        CALL getres(ih,ik,il,s2,resol)
        IF (resol>reslim1.OR.resol<reslim2) CYCLE
        CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,size(op,3)/ncenop)
        hmin=MIN(hmin,jh)
        kmin=MIN(kmin,jk)
        lmin=MIN(lmin,jl)
        hmax=MAX(hmax,jh)
        kmax=MAX(kmax,jk)
        lmax=MAX(lmax,jl)
        n=n+1
      END DO
      WRITE(*,*) n,' reflections read from reference'
      WRITE(*,'(a,6i6)') 'hmin,hmax,kmin,kmax,lmin,lmax=',hmin,hmax,kmin,kmax,lmin,lmax
      REWIND(unit)
      ALLOCATE(ref(hmin:hmax,kmin:kmax,lmin:lmax))
      ref=-huge(ref)
      IF (anoflag) THEN
        ALLOCATE(refminus(hmin:hmax,kmin:kmax,lmin:lmax))
        refminus=-huge(ref)
      END IF
! fill array
      n=0
      DO
        IF (anoflag) THEN
          READ(unit,*,iostat=ier) ih,ik,il,amplitude,friedel
          IF (amplitude<0.OR.friedel<0) STOP 'Error: reference I or F is <0'
        ELSE
          READ(unit,*,iostat=ier) ih,ik,il,amplitude
          IF (amplitude<0) STOP 'Error: reference I or F is <0'
        END IF
        IF (ier == iostat_end) EXIT
        IF (ier /= 0) CYCLE
        CALL getres(ih,ik,il,s2,resol)
        IF (resol>reslim1.OR.resol<reslim2) CYCLE
        CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,size(op,3)/ncenop)
        IF (ref(jh,jk,jl) < -huge(ref)/2) THEN
          ref(jh,jk,jl)=amplitude
          IF (squareflag) ref(jh,jk,jl)=amplitude**2
          IF (anoflag) THEN
            refminus(jh,jk,jl)=friedel
            IF (squareflag) refminus(jh,jk,jl)=friedel**2
! TODO: does this work with intensity reference?
            IF (ABS(amplitude-friedel)/MAX(amplitude,friedel)>0.5) &  
            PRINT '(a,3i5,2f10.3)','large anom signal in reference',ih,ik,il,amplitude,friedel
          END IF
        ELSE IF (n<10) THEN
          n=n+1
          print*,jh,jk,jl,amplitude,' duplicate reflection in reference!'
        END IF
      END DO
      
      CLOSE(unit)
      END SUBROUTINE read_reference

ELEMENTAL FUNCTION atanh_safe(x)
! purpose: prevent ATANH() from crashing if called with 1. or -1.
REAL(KIND(1.d0)), INTENT(IN) :: x
REAL(KIND(1.d0)) :: y,atanh_safe

IF (x>1.d0-EPSILON(1.d0)) THEN
  y=1.d0-EPSILON(1.d0)
ELSE IF (x<EPSILON(1.d0)-1.d0) THEN
  y=EPSILON(1.d0)-1.d0
ELSE
  y=x
END IF
atanh_safe=ATANH(y)
END FUNCTION atanh_safe

END MODULE reference_mod


SUBROUTINE reflection_file(unit,reslim1,reslim2,nbin,abcdef,reference,nbatch,anoflag)
      USE symops_mod
      USE resolution_mod
      USE rc_mod
      USE hash_mod
      USE reference_mod
      USE cc_mod
      USE xscaleitems_mod, ONLY: is_xds,nset,save_a,save_d
      IMPLICIT NONE
      
      CHARACTER (LEN=132),INTENT(IN) :: reference
      INTEGER, INTENT(IN) :: unit,nbin,nbatch
      REAL, INTENT(IN) :: reslim1,reslim2
      LOGICAL, INTENT(IN) :: abcdef(10),anoflag
      INTEGER, PARAMETER  :: maxref=100000   ! only needed for the hash function
      INTEGER :: ih,ik,il,oldhkl=0,jh,jk,jl,ihkl=0,iop,i,ispgr,ncheck,ier, &
                 iset,nref,kset,nobs,nxa,ibin,misfits=0,ntot(nbin),ncontrib=0,m
      REAL :: s2,resol,fsq,sigma,dum(4),  &
              vvr(nbin),varr(nbin),wwr(nbin),s2save,resmin=HUGE(resmin),resmax=0., &
              ref_int,ref_minus
      REAL(KIND(1.d0)) :: xx,yy,vv,ww,var,totsum(nbin),totsum2(nbin), &
                          totvar(nbin),xxr(nbin),yyr(nbin)
      REAL(KIND(1.d0)) :: sumsig2,sumxap,sumxap2,sumxam,sumxam2,t1,t2,sumsigp,sumsigm, &
                          sumsigp2,sumsigm2,sumxa,sumxa2,sumsig
      CHARACTER (LEN=80) :: line,fmt,fmti,fmtr
      LOGICAL :: last=.FALSE.,centric
      INTEGER :: nanotot(nbin),nxap,nxam
      REAL(KIND(xx)) :: anosum(nbin),anosum2(nbin),anovar(nbin),ano,varp,varm
      REAL(KIND(xx)), ALLOCATABLE :: rcset(:,:,:),anoset(:,:,:)
      INTEGER, ALLOCATABLE :: nrcset(:,:),nrcset2(:,:),nanoset(:,:),nanoset2(:,:), &
                              nap(:),nam(:),icontrib(:)            
! floating overflow in line 573 xap2(..... 10.12.2020 therefore kind(1.d0)
      REAL(KIND(1.d0)), ALLOCATABLE :: xap(:),xap2(:),xam(:),xam2(:),sigap(:),sigam(:),sigap2(:),sigam2(:)
                           
      TYPE(cc_type), ALLOCATABLE :: cc_only(:,:),cc_onlyano(:,:),cc_overall(:),cc_overallano(:)

      fmt ='(a,40f7.2,a,f7.2)'
      fmtr='(a,i5,40f6.1,a,f6.1)'
      fmti='(a,5x,40i6,  a,i7)'
      WRITE(fmt(4:5),'(i2.2)')nbin
      WRITE(fmtr(7:8),'(i2.2)')nbin
      WRITE(fmti(7:8),'(i2.2)')nbin
!      WRITE(*,*)fmt
 1000 FORMAT(A)
      ntot=0
      totsum=0.
      totsum2=0.
      totvar=0.
      nanotot=0
      anosum=0.
      anosum2=0.
      anovar=0.
      xxr=0.
      yyr=0.
  ! the factor of 2 makes sure that the hash table is never more than half full:
      CALL hash_init(2*maxref)
      ncheck=0
! Open input file; get spacegroup and cell constants
      DO
        READ(unit,1000,iostat=ier) line
        IF (ier/=0) THEN
          PRINT*,'ier=',ier
          PRINT*,'line=',TRIM(line)
          STOP 'read error in xds_file'
        END IF
        IF (line(1:20)=='!SPACE_GROUP_NUMBER=') then
          READ(line(21:),*) ispgr
          CALL getops(ispgr)
          WRITE(*,'(a)') TRIM(line)
        ELSE IF (line(1:21).eq.'!UNIT_CELL_CONSTANTS=') then 
          READ(line(22:),*) a,b,c,alf,bet,gam
          WRITE(*,'(a)') TRIM(line)
        ELSE IF (line(1:9) == '!ITEM_H=1') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:9) == '!ITEM_K=2') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:9) == '!ITEM_L=3') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:12) == '!ITEM_IOBS=4') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:19) == '!ITEM_SIGMA(IOBS)=5') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:13) == '!ITEM_ISET=10') THEN
          ncheck=ncheck+1
        ELSE IF (line(1:7) == '! ISET=') THEN
          READ(line(8:),*) iset
          nset=max(nset,iset)
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
      IF (ncheck <5) THEN
        STOP 'header items not found'  ! we tolerate it if ISET is not there
      END IF
      IF (is_xds.AND.nbatch>1) THEN
        WRITE(*,*) '# of batches=',nset
      ELSE IF (is_xds) THEN
        WRITE(*,*) '# of frames=',nset
      ELSE
        WRITE(*,*) '# of datasets=',nset
      END IF
      
      IF (.NOT.abcdef(9)) CALL read_reference(reference,reslim1,reslim2,anoflag,.TRUE.)
      IF (.NOT.abcdef(10)) CALL read_reference(reference,reslim1,reslim2,anoflag,.FALSE.)
      
      ALLOCATE(rcset(9,nset,nbin),anoset(9,nset,nbin),nrcset(nset,nbin), &
               nrcset2(nset,nbin),nanoset(nset,nbin),nanoset2(nset,nbin),&
               icontrib(nset),nap(nset),nam(nset),xap(nset),xam(nset),   &
               xap2(nset),xam2(nset),sigap(nset),sigam(nset),sigap2(nset),&
               sigam2(nset),save_a(nset,2),save_d(nset,2))
      save_d=0.
      ALLOCATE(cc_only(nset,nbin),cc_overall(nbin)) ! ,cc_with(nset,nbin),cc_without(nset,nbin)
! bugfix 5.12.19: remove "IF (anoflag)"
!     IF (anoflag) ALLOCATE(cc_onlyano(nset,nbin),cc_overallano(nbin))
      ALLOCATE(cc_onlyano(nset,nbin),cc_overallano(nbin))
      nobs=0
      nref=0
      rcset=0.
      anoset=0.
      nrcset=0
      nrcset2=0
      nanoset=0
      nanoset2=0
      icontrib=0
      nap=0
      nam=0
      xap=0
      xam=0
      xap2=0
      xam2=0
      sigap=0
      sigam=0
      sigap2=0
      sigam2=0
! read reflections
      DO
        READ(unit,1000) line
        IF (line(1:12) == '!END_OF_DATA') last=.TRUE. 
        IF (.NOT.last) THEN
          IF (is_xds) THEN
            READ(line,*)ih,ik,il,fsq,sigma,dum(1:3)
            iset=INT(dum(3)/nbatch)+1
          ELSE 
            READ(line,*)ih,ik,il,fsq,sigma,dum,iset
          END IF
          CALL getres(ih,ik,il,s2,resol)
          IF (resol<resmin) resmin=resol
          IF (resol>resmax) resmax=resol
          IF (resol>reslim1.OR.resol<reslim2) CYCLE
          IF (sigma<0.) THEN
            misfits = misfits + 1
            CYCLE
          END IF 
          IF (.NOT.abcdef(7)) THEN  ! switch off weighting
            sigma=1
          END IF
          CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,SIZE(op,3)/ncenop)
          nobs=nobs+1
          IF (nobs==1) THEN
            oldhkl=ihkl
            s2save=s2
          END IF
        ELSE IF (nobs==0) THEN
          EXIT ! the reflection loop
        END IF
        IF (ihkl /= oldhkl.OR.last) THEN  ! process unique reflection
          nref=nref+1                  ! update # of unique refs
          IF (nref==maxref+1) THEN
            PRINT*,'too many unique reflections - no checking for duplicates'
          ELSE IF (nref<=maxref) THEN
            IF (hash_exists(oldhkl)) THEN
              PRINT*,'duplicate unique reflection:',oldhkl
              STOP 'reflections are not properly sorted - is the file from XDS/XSCALE??'
            ELSE
              CALL hash_insert(oldhkl)
            END IF
          END IF
          ref_int=-1
          IF ((.NOT.abcdef(9)).OR.(.NOT.abcdef(10))) THEN  ! find out reference intensities
            CALL hkluncompress(oldhkl,ih,ik,il) 
            IF (ih>= LBOUND(ref,1).AND.ih<=UBOUND(ref,1).AND.&
                ik>= LBOUND(ref,2).AND.ik<=UBOUND(ref,2).AND.&
                il>= LBOUND(ref,3).AND.il<=UBOUND(ref,3)) THEN
              ref_int=ref(ih,ik,il)
              IF (anoflag) ref_minus=refminus(ih,ik,il)
            END IF
          END IF
! find resolution bin:
          CALL getbin(s2save,ibin)
          nxap=0
          nxam=0
          sumxap=0.
          sumxap2=0.
          sumxam=0.
          sumxam2=0.
          sumsigp=0.
          sumsigm=0.
          sumsigp2=0.
          sumsigm2=0.
! sum over all datasets
          DO i=1,ncontrib
            kset=icontrib(i)   ! dataset number
            nxap=nxap+nap(kset)
            nxam=nxam+nam(kset)
            sumxap  = sumxap  + xap(kset)
            sumxap2 = sumxap2 + xap2(kset)
            sumxam  = sumxam  + xam(kset)
            sumxam2 = sumxam2 + xam2(kset)
            sumsigp = sumsigp + sigap(kset)
            sumsigm = sumsigm + sigam(kset)
            sumsigp2 = sumsigp2 + sigap2(kset)
            sumsigm2 = sumsigm2 + sigam2(kset)
! CC of each data set against reference:
            IF (ref_int>0) THEN
              IF (anoflag) THEN  ! ref_int and ref_minus have to be considered. sum of them gives same CC as average
                CALL put_cc(cc_only(kset,ibin),(xap(kset)+xam(kset))/(sigap(kset)+sigam(kset)),ref_int+ref_minus)
                IF (sigap(kset)*sigam(kset)>0) &
                CALL put_cc(cc_onlyano(kset,ibin),xap(kset)/sigap(kset)-xam(kset)/sigam(kset),ref_int-ref_minus)
              ELSE
! bugfix 24.11.18 hier hat die Klammer um xap(kset)+xam(kset) gefehlt!
                CALL put_cc(cc_only(kset,ibin),(xap(kset)+xam(kset))/(sigap(kset)+sigam(kset)),ref_int)
              END IF
            END IF
          END DO
          nxa    = nxap    + nxam
          sumxa  = sumxap  + sumxam
          sumxa2 = sumxap2 + sumxam2
          sumsig = sumsigp + sumsigm
          sumsig2 = sumsigp2 + sumsigm2
! CC overall against reference:
          IF (ref_int>0) THEN
            IF (anoflag) THEN
              CALL put_cc(cc_overall(ibin),sumxa/sumsig,ref_int+ref_minus)
              IF (sumsigp*sumsigm>0) CALL put_cc(cc_overallano(ibin),sumxap/sumsigp-sumxam/sumsigm,ref_int-ref_minus)
            ELSE
              CALL put_cc(cc_overall(ibin),sumxa/sumsig,ref_int)   ! bugfix 26.4.19 had been missing /sumsig 
            END IF
          END IF
! isomorphous signal: store contributions to this unique reflection
          IF (ncontrib>1.OR.nap(icontrib(1))+nam(icontrib(1))>1) THEN  ! >2 is needed for CC12-overall, >2 for delta-CC12
! contributions to overall CC1/2
            sumxa = sumxa/sumsig
            sumxa2 = sumxa2/sumsig
            ntot(ibin)   = ntot(ibin) + 1                          ! number of contributions
            totsum(ibin) = totsum(ibin) + sumxa             ! sum of mean_intensities
            totsum2(ibin)= totsum2(ibin) + sumxa**2       ! sum of mean_intensities**2
!           var    = (sumxa2-sumxa**2)*nxa/(nxa-1.)     ! variance of epsilon. 
!            var    = (sumxa2-sumxa**2)/(1.-sumsig2/sumsig**2) ! (weighted) variance of epsilon. Fix 5.4.19
!           var    = var / (nxa/2.)                    ! variance of half-dataset
!            var    = var / (sumsig**2/sumsig2/2.)                    ! variance of half-dataset. Fix 6.4.19
            var = 2.*sumsig2*(sumxa2-sumxa**2)/(sumsig**2-sumsig2)
            totvar(ibin) = totvar(ibin) + var                      ! sum of variances of half-dataset
            DO i=1,ncontrib
              kset=icontrib(i)
              IF (nap(kset)+nam(kset)>nxa-2) THEN   ! not enough other obs for variance calc
                nrcset2(kset,ibin) = nrcset2(kset,ibin) + 1             ! number of contributions
                rcset(7,kset,ibin)= rcset(7,kset,ibin) + sumxa    ! sum of mean intensities
                rcset(8,kset,ibin)= rcset(8,kset,ibin) + sumxa**2 ! sum of mean_intensities**2
                rcset(9,kset,ibin)= rcset(9,kset,ibin) + var       ! sum of variances of half-dataset               
              ELSE
! _with_ this dataset: same contributions as above for ntot,tosum,totsum2,totvar
                nrcset(kset,ibin) = nrcset(kset,ibin) + 1             ! number of contributions
                rcset(1,kset,ibin)= rcset(1,kset,ibin) + sumxa    ! sum of mean intensities
                rcset(2,kset,ibin)= rcset(2,kset,ibin) + sumxa**2 ! sum of mean_intensities**2
                rcset(3,kset,ibin)= rcset(3,kset,ibin) + var       ! sum of variances of half-dataset
! _without_ this dataset: modify contributions
                xx=sumxa*sumsig -(xap(kset) + xam(kset))
                xx=xx/(sumsig-(sigap(kset)+sigam(kset)))     
                yy=sumxa2*sumsig-(xap2(kset) + xam2(kset))
                yy=yy/(sumsig-(sigap(kset)+sigam(kset)))
                rcset(4,kset,ibin)=rcset(4,kset,ibin) + xx  ! sum of mean_intensities
                rcset(5,kset,ibin)=rcset(5,kset,ibin) + xx**2 ! sum of mean_intensities**2
                m= nxa - ( nap(kset) + nam(kset))
!                xx = (yy-xx**2)*m/(m-1.)           ! variance of epsilon
                t1 = sumsig  - ( sigap(kset) +  sigam(kset))
                t2 = sumsig2 - (sigap2(kset) + sigam2(kset))
!                xx = (yy-xx**2)/(1.-t2/t1**2)           ! (weighted) variance of epsilon
!               xx = xx / (m/2.)                           ! variance of half-dataset
!                xx = xx / (t1**2/t2/2.)                           ! variance of half-dataset
                xx = 2.*t2*(yy-xx**2)/(t1**2-t2)
                rcset(6,kset,ibin)=rcset(6,kset,ibin) + xx     ! sum of variances of half-dataset
              END IF
            END DO
          END IF  ! isomorphous signal
          
! anomalous signal: store contributions to this unique reflection
          IF (nxap>1.AND.nxam>1) THEN  ! >2 is needed for anoCC12-overall; >2 for delta-anoCC12
! contributions to overall anoCC1/2
            nanotot(ibin)   = nanotot(ibin) + 1                    ! number of contributions
            ano = sumxap/sumsigp - sumxam/sumsigm           
            anosum(ibin) = anosum(ibin) +  ano              ! sum of anomalous signals
            anosum2(ibin)= anosum2(ibin) + ano**2           ! sum of anomalous signals**2
!           varp   = (sumxap2/sumsigp-(sumxap/sumsigp)**2)*nxap/(nxap-1.) ! variance of +epsilon
!            varp   = (sumxap2/sumsigp-(sumxap/sumsigp)**2)/(1.-sumsigp2/sumsigp**2) ! (weighted) variance of +epsilon
!           varp   = varp / (nxap/2.)                  ! variance of +half-dataset
!            varp   = varp / (sumsigp**2/sumsigp2/2.)                  ! variance of +half-dataset
            varp = 2*sumsigp2*(sumxap2/sumsigp-(sumxap/sumsigp)**2)/(sumsigp**2-sumsigp2)

!           varm   = (sumxam2/sumsigm-(sumxam/sumsigm)**2)*nxam/(nxam-1.) ! variance of -epsilon
!            varm   = (sumxam2/sumsigm-(sumxam/sumsigm)**2)/(1.-sumsigm2/sumsigm**2) ! (weighted) variance of -epsilon
!           varm   = varm / (nxam/2.)                  ! variance of -half-dataset
!            varm   = varm / (sumsigm**2/sumsigm2/2.)                  ! variance of -half-dataset
            varm = 2*sumsigm2*(sumxam2/sumsigm-(sumxam/sumsigm)**2)/(sumsigm**2-sumsigm2)
            
            anovar(ibin) = anovar(ibin) + varp + varm      ! sum of variances of half-datasets
! contributions to dataset-CC1/2s
            DO i=1,ncontrib
              kset=icontrib(i)
              IF (nap(kset)>nxap-2.OR.nam(kset)>nxam-2) THEN   ! not enough other obs for variance calc
                nanoset2(kset,ibin)= nanoset2(kset,ibin) + 1       ! number of contributions
                anoset(7,kset,ibin)= anoset(7,kset,ibin) + ano    ! sum 
                anoset(8,kset,ibin)= anoset(8,kset,ibin) + ano**2 ! sum of squares
                anoset(9,kset,ibin)= anoset(9,kset,ibin) + varp + varm ! sum of variances of half-dataset               
              ELSE
! _with_ this dataset: same contributions as above for overall anoCC1/2
                nanoset(kset,ibin) = nanoset(kset,ibin) + 1         ! number of contributions
                anoset(1,kset,ibin)= anoset(1,kset,ibin) + ano    ! sum 
                anoset(2,kset,ibin)= anoset(2,kset,ibin) + ano**2 ! sum of squares
                anoset(3,kset,ibin)= anoset(3,kset,ibin) + varp + varm ! sum of variances
! _without_ this dataset: modify contributions
                t1=sumsigp  -  sigap(kset)
                xx=(sumxap-xap(kset))/t1                            ! I+ without
                yy=(sumxam-xam(kset))/(sumsigm-sigam(kset))         ! I- without
                anoset(4,kset,ibin)=anoset(4,kset,ibin) + xx - yy        ! sum 
                anoset(5,kset,ibin)=anoset(5,kset,ibin) + (xx - yy)**2   ! sum of squares
! variance of epsilon: contribution from I+ :
!                xx = (sumxap  - xap(kset)) /(sumsigp - sigap(kset))
                yy = (sumxap2 - xap2(kset))/t1
!                m= nxap - nap(kset)                
!                xx = (yy-xx**2)*m/(m-1.) / (m/2.) ! (m/2.) to get half-dataset variance
                t2=sumsigp2 - sigap2(kset)
!                xx = (yy-xx**2)/(1.-t2/t1**2) / (m/2.) ! (m/2.) to get (weighted) half-dataset variance
!                xx = (yy-xx**2)/(1.-t2/t1**2) / (t1**2/t2/2.) ! t1**2/t2/2. to get (weighted) half-dataset variance
                xx = 2*t2*(yy-xx**2)/(t1**2-t2)  ! simplified, and less risk of divide-by-zero
                anoset(6,kset,ibin)=anoset(6,kset,ibin) + xx ! add variance of I+ half-dataset
! variance of epsilon: contribution from I- :
! FIXME: repeated calculation: yy above is here called xx
                t1=sumsigm  -  sigam(kset)
                xx = (sumxam  - xam(kset)) /t1
                yy = (sumxam2 - xam2(kset))/t1
!                m= nxam - nam(kset)                
!                xx = (yy-xx**2)*m/(m-1.) / (m/2.) ! (m/2.) to get half-dataset variance
                t2=sumsigm2 - sigam2(kset)
!                xx = (yy-xx**2)/(1.-t2/t1**2) / (m/2.) ! (m/2.) to get (weighted) half-dataset variance
!                xx = (yy-xx**2)/(1.-t2/t1**2) / (t1**2/t2/2.) ! t1**2/t2/2. to get (weighted) half-dataset variance
                xx = 2*t2*(yy-xx**2)/(t1**2-t2)  ! simplified
                anoset(6,kset,ibin)=anoset(6,kset,ibin) + xx ! add variance of I- half-dataset
              END IF
            END DO
          END IF  ! anomalous signal
          
          IF (last) EXIT  ! the big loop reading the XSCALE file    
          DO i=1,ncontrib   ! clear for next unique ref
            kset=icontrib(i)
            nap(kset)=0
            xap(kset)=0.
            xap2(kset)=0.
            nam(kset)=0
            xam(kset)=0.
            xam2(kset)=0.
            sigap(kset)=0.
            sigam(kset)=0.
            sigap2(kset)=0.
            sigam2(kset)=0.
          END DO
          ncontrib=0
          oldhkl=ihkl
          s2save=s2
        END IF
! store this observation
        IF (ncontrib==0) THEN   ! happens once for each unique ref
          centric=is_centric(jh,jk,jl,op,size(op,3)/ncenop)
        END IF
        IF (nap(iset)+nam(iset)==0) THEN  ! new iset
          ncontrib=ncontrib+1
          icontrib(ncontrib)=iset
        END IF
        xx=sigma**2
        IF (iop>0.OR.centric) THEN  ! this is taken as I+
          xap(iset) = xap(iset) + fsq/xx
          xap2(iset)= xap2(iset) + fsq**2/xx
          nap(iset) = nap(iset) + 1
          sigap(iset) = sigap(iset) + 1/xx
          sigap2(iset) = sigap2(iset) + 1/xx**2
        ELSE                        !  I-
          xam(iset) = xam(iset) + fsq/xx
          xam2(iset)= xam2(iset) + fsq**2/xx
          nam(iset) = nam(iset) + 1
          sigam(iset) = sigam(iset) + 1/xx
          sigam2(iset) = sigam2(iset) + 1/xx**2
        END IF
      END DO   ! the big loop reading the XSCALE file 

! output of statistics

      WRITE(*,*) '# obs (w/o misfits), unique, misfits =',nobs,nref,misfits
      WRITE(*,*) 'max and min resolution of data in file =',resmax,resmin
      WRITE(*,*) 'data between',reslim1,' and',reslim2,'A will be used'
      WRITE(*,'(i3,a)') nbin,' resolution shells (for lines starting with b,c,e,f,r,s):'
!      WRITE(*,*) 'resolution shells (for lines starting with b,c,e,f,r,s):'
      CALL getlim(vvr)
      WRITE(*,'(40f7.3)') vvr
      WRITE(*,*)
! overall CC1/2
      IF (SUM(ntot)>1) THEN
        xx = (SUM(totsum2)-SUM(totsum)**2/SUM(ntot)) / (SUM(ntot)-1)  
        var = SUM(totvar)/SUM(ntot)  ! bugfix 6.9.2018; was: (SUM(ntot)-1)
        xx = (xx - 0.5*var)/(xx + 0.5*var)
      ELSE
        xx = 0
      END IF
      WRITE(*,'(a,f10.3,a,i8,a)') ' overall CC1/2:',100*xx,' nref=',SUM(ntot), &
                 ' (but the overall CC1/2 is meaningless!)'
      
      WHERE(ntot>1)
        vvr = (totsum2-totsum**2/ntot) / (ntot-1)
        varr = totvar/ntot         ! bugfix 6.9.2018; was: (ntot-1)
        vvr = (vvr - 0.5*varr)/(vvr + 0.5*varr)
      ELSE WHERE
        vvr=0
      END WHERE
      vv=DOT_PRODUCT(vvr,ntot)/MAX(1,SUM(ntot))
      WRITE(*,'(a,f10.3,a)') ' <CC1/2>:',100*vv,' (frequency-weighted average of CC1/2 in resolution shells)'
      WRITE(*,'(a)') ' CC1/2 in resolution shells:'
      WRITE(*,'(40f7.1)') 100*vvr
      WRITE(*,'(a)') ' CC* in resolution shells:'
      WRITE(*,'(40f7.1)') 100*SQRT((2*MAX(0.,vvr))/(1+vvr))
      WRITE(*,'(a)') ' frequency, i.e. number of unique reflections in resolution shells:'     
      WRITE(*,'(40i7)') ntot
      IF (.NOT.abcdef(9)) THEN
        IF (anoflag) THEN
          WRITE(*,'(a)') ' CC Iobs vs Iref (Iref calculated from reference Fref+ and Fref-), in resolution shells:'
          WRITE(*,'(40f7.1)') 100*get_cc(cc_overall(:))
          WRITE(*,'(a)') ' # of Fref reference reflections in resolution shells (with -s option: only acentric refs!):'
          WRITE(*,'(40i7)') cc_overall(:)%n
        ELSE
          WRITE(*,'(a)') ' CC Iobs vs Iref (Fref from reference), in resolution shells:'
          WRITE(*,'(40f7.1)') 100*get_cc(cc_overall(:))
          WRITE(*,'(a)') ' # of Fref reference reflections used, in resolution shells:'
          WRITE(*,'(40i7)') cc_overall(:)%n
        END IF
      END IF
      IF (.NOT.abcdef(10)) THEN
        IF (anoflag) THEN
          WRITE(*,'(a)') ' CC Iobs vs Iref (Iref from reference Fref+ or Fref-), in resolution shells:'
          WRITE(*,'(40f7.1)') 100*get_cc(cc_overall(:))
          WRITE(*,'(a)') ' # of Iref reference reflections in resolution shells (with -s option: only acentric refs!):'
          WRITE(*,'(40i7)') cc_overall(:)%n
        ELSE
          WRITE(*,'(a)') ' CC Iobs vs Iref (Iref from reference), in resolution shells:'
          WRITE(*,'(40f7.1)') 100*get_cc(cc_overall(:))
          WRITE(*,'(a)') ' # of Iref reference reflections used, in resolution shells:'
          WRITE(*,'(40i7)') cc_overall(:)%n
        END IF
      END IF

      WRITE(*,*)
      IF ((.NOT.abcdef(9)).OR.(.NOT.abcdef(10))) WRITE(*,*)'lines starting with r or R or s or # refer to the CC(Iobs,Iref)'
      WRITE(*,*) 'headings for lines starting with a,b,c:'

      IF (abcdef(1)) THEN
        IF (is_xds.AND.nbatch==1) THEN
          WRITE(*,'(a)') 'a:  <CC1/2> of each frame:'
          WRITE(*,'(a)') 'a:   reflections of this frame only       reflections of all frames'
          WRITE(*,'(a)') 'a: frm   nref    with without   delta     nref    with without   delta'
        ELSE IF (is_xds) THEN
          WRITE(*,'(a)') 'a:  <CC1/2> of each batch (batch # is given as highest frame #):'
          WRITE(*,'(a)') 'a:   reflections of this batch only        reflections of all batches'
          WRITE(*,'(a)') 'a: batch nref    with without   delta     nref    with without   delta'
        ELSE
          WRITE(*,'(a)') 'a:  <CC1/2> of each dataset:'
          WRITE(*,'(a)') 'a:   reflections of this dataset only      reflections of all datasets'
          WRITE(*,'(a)') 'a: set   nref    with without   delta     nref    with without   delta'
        END IF
      END IF
      IF (abcdef(2)) WRITE(*,'(a)') 'b: delta-CC1/2 in resolution shells'
      IF (abcdef(3)) WRITE(*,'(a)') 'c: # reflections for delta-CC1/2'     

!      
! for each dataset: calculate and print delta-CC1/2
!  
! we use CC12 = var_tau / (var_tau + var_eps) , the formula in supplement of Science 2012 paper
! isomorphous signal
      DO i=1,nset
        IF (.NOT.abcdef(9).AND.SUM(cc_only(i,:)%n)>0) THEN
          WRITE(*,fmtr) 'r',i*nbatch,100*get_cc(cc_only(i,:)),' avg:', &
            DOT_PRODUCT(100*get_cc(cc_only(i,:)),cc_only(i,:)%n)/SUM(cc_only(i,:)%n)
          WRITE(*,fmti) '#',cc_only(i,:)%n,' sum',SUM(cc_only(i,:)%n)
        END IF
        IF (.NOT.abcdef(10).AND.SUM(cc_only(i,:)%n)>0) THEN
          WRITE(*,fmtr) 'R',i*nbatch,100*get_cc(cc_only(i,:)),' avg:', &
            DOT_PRODUCT(100*get_cc(cc_only(i,:)),cc_only(i,:)%n)/SUM(cc_only(i,:)%n)
          WRITE(*,fmti) '#',cc_only(i,:)%n,' sum',SUM(cc_only(i,:)%n)
        END IF
        IF (SUM(nrcset(i,:))<2) CYCLE
! mode 1: only reflections present in dataset contribute to CC12:
        WHERE (nrcset(i,:)>1)  ! 30.7.15 "floating invalid" wenn >0 
!          xxr =(rcset(2,i,:)-rcset(1,i,:)**2/nrcset(i,:)) ! / (nrcset(i,:)-1)  ! var_tau +0.5*var_eps, with
!	   varr=rcset(3,i,:) !/(nrcset(i,:)-1)     ! half-dataset variance of epsilon, with
! bugfix 6.9.18:
          xxr =(rcset(2,i,:)-rcset(1,i,:)**2/nrcset(i,:)) / (nrcset(i,:)-1)  ! var_tau +0.5*var_eps, with
          varr=rcset(3,i,:)/nrcset(i,:)     ! half-dataset variance of epsilon, with
          xxr = (xxr - 0.5*varr)/(xxr + 0.5*varr)
!          yyr =(rcset(5,i,:)-rcset(4,i,:)**2/nrcset(i,:)) !/ (nrcset(i,:)-1) ! var_tau +0.5*var_eps, without
!	   varr=rcset(6,i,:) !/(nrcset(i,:)-1)   ! var_eps= half-dataset variance of epsilon, without
! bugfix 6.9.18:
          yyr =(rcset(5,i,:)-rcset(4,i,:)**2/nrcset(i,:)) / (nrcset(i,:)-1) ! var_tau +0.5*var_eps, without
          varr=rcset(6,i,:)/nrcset(i,:)   ! var_eps= half-dataset variance of epsilon, without
          yyr = (yyr - 0.5*varr)/(yyr + 0.5*varr)
        ELSE WHERE
          xxr=0.
          yyr=0.
        END WHERE

! mode 2: all reflections contribute to CC12:
! remove contributions from uniques where this dataset contributed more than n-2 observations
        nrcset2(i,:)=ntot - nrcset2(i,:)

! CC12 without this dataset
        rcset(4,i,:)=rcset(4,i,:) - rcset(1,i,:) + totsum(:) -rcset(7,i,:)
        rcset(5,i,:)=rcset(5,i,:) - rcset(2,i,:) + totsum2(:)-rcset(8,i,:)
        rcset(6,i,:)=rcset(6,i,:) - rcset(3,i,:) + totvar(:) -rcset(9,i,:)   
        WHERE (nrcset2(i,:)>1)   ! bugfix 9.10.20 was: >0 
!          wwr=(rcset(5,i,:)-rcset(4,i,:)**2/nrcset2(i,:)) ! / (nrcset2(i,:)-1)  ! var_tau +0.5*var_eps, with
!          varr=rcset(6,i,:) !/(nrcset2(i,:)-1)     ! var_eps= half-dataset variance of epsilon, with
! bugfix 6.9.18:
          wwr=(rcset(5,i,:)-rcset(4,i,:)**2/nrcset2(i,:)) / (nrcset2(i,:)-1)  ! var_tau +0.5*var_eps, with
          varr=rcset(6,i,:)/nrcset2(i,:)     ! var_eps= half-dataset variance of epsilon, with
          wwr = (wwr - 0.5*varr)/(wwr + 0.5*varr)

! CC12 with this dataset:
          rcset(1,i,:)=totsum(:) - rcset(7,i,:)
          rcset(2,i,:)=totsum2(:)- rcset(8,i,:)
          rcset(3,i,:)=totvar(:) - rcset(9,i,:)
!          vvr=(rcset(2,i,:)-rcset(1,i,:)**2/nrcset2(i,:)) !/ (nrcset2(i,:)-1)  ! var_tau +0.5*var_eps, with
!          varr=rcset(3,i,:) !/(nrcset2(i,:)-1)     ! half-dataset variance of epsilon, with
! bugfix 6.9.18:
          vvr=(rcset(2,i,:)-rcset(1,i,:)**2/nrcset2(i,:)) / (nrcset2(i,:)-1)  ! var_tau +0.5*var_eps, with
          varr=rcset(3,i,:)/nrcset2(i,:)     ! half-dataset variance of epsilon, with
          vvr = (vvr - 0.5*varr)/(vvr + 0.5*varr)
        ELSE WHERE
          vvr=0.
          wwr=0.
        END WHERE
        
        xx=DOT_PRODUCT(xxr,nrcset(i,:))/SUM(nrcset(i,:))
        yy=DOT_PRODUCT(yyr,nrcset(i,:))/SUM(nrcset(i,:))
        vv=DOT_PRODUCT(vvr,nrcset2(i,:))/SUM(nrcset2(i,:))
        ww=DOT_PRODUCT(wwr,nrcset2(i,:))/SUM(nrcset2(i,:))
        IF (.NOT.abcdef(8)) THEN  ! Fisher transform switched off
!          IF (abcdef(1)) WRITE(*,'(a,i6,i7,3f8.3,i9,3f8.3)') 'a',i*nbatch,                        &
          IF (abcdef(1)) WRITE(*,'(a,i6,i7,2f8.3,f9.4,i9,2f8.3,f9.4)') 'a',i*nbatch,  &  ! 10.6 changed to 9.4 2.7.22 KD
            SUM(nrcset(i,:)),100*xx,100*yy,100*(xx-yy),  &
            SUM(nrcset2(i,:)),100*vv,100*ww,100*(vv-ww)
          save_a(i,1)=100*(xx-yy)
          save_a(i,2)=100*(vv-ww)
          IF (abcdef(2)) WRITE(*,'(a,40(1x,f0.3))') 'b',100*(xxr-yyr)
        ELSE
          IF (abcdef(1)) THEN
            WRITE(*,'(a,i7,i7,3f8.3,i9,3f8.3)') 'a',i*nbatch,                        &
              SUM(nrcset(i,:)),100*xx,100*yy,100*TANH(ATANH_safe(xx)-ATANH_safe(yy)),  &
              SUM(nrcset2(i,:)),100*vv,100*ww,100*TANH(ATANH_safe(vv)-ATANH_safe(ww))
            save_a(i,1)=100*TANH(ATANH_safe(xx)-ATANH_safe(yy))
            save_a(i,2)=100*TANH(ATANH_safe(vv)-ATANH_safe(ww))
          END IF
          IF (abcdef(2)) &
            WRITE(*,'(a,40(1x,f0.3))') 'b',100*TANH(ATANH_safe(xxr)-ATANH_safe(yyr))  
        END IF 
        IF (abcdef(3)) WRITE(*,'(a,40i7)') 'c',nrcset(i,:)
      END DO
      
      WRITE(*,*)'-------------------------------------------------------------'
      WRITE(*,*)      
! overall anoCC1/2
      IF (SUM(nanotot)>1) THEN
        xx = (SUM(anosum2)-SUM(anosum)**2/SUM(nanotot)) / (SUM(nanotot)-1)
        var = SUM(anovar)/SUM(nanotot)  ! bugfix 6.9.2018, was: (SUM(nanotot)-1)
        xx = (xx - 0.5*var)/(xx + 0.5*var)
      ELSE
        xx = 0
      END IF
      WRITE(*,'(a,f10.3,a,i8,a)') 'overall CC1/2ano:',100*xx,' nref=',SUM(nanotot),&
              ' (but the overall value is meaningless!)'
      WHERE(nanotot>1)
        vvr = (anosum2-anosum**2/nanotot) / (nanotot-1)
        varr = anovar/nanotot         ! bugfix 6.9.2018, was: (nanotot-1) 
        vvr = (vvr - 0.5*varr)/(vvr + 0.5*varr)
      ELSE WHERE
        vvr = 0
      END WHERE
      vv=DOT_PRODUCT(vvr,nanotot)/MAX(1,SUM(nanotot))  ! bugfix 9.5.19 avoid divide-by-zero
      WRITE(*,'(a,f10.3,a)') '<CC1/2ano> :',100*vv,' (frequency-weighted average of CC1/2ano in resolution shells)'
      WRITE(*,'(a)') 'CC1/2ano in resolution shells:'
      WRITE(*,'(40f7.1)') 100*vvr
      WRITE(*,'(a)') 'frequency i.e. number of unique Friedel pairs:'
      WRITE(*,'(40i7)') nanotot
      IF (.NOT.abcdef(9)) THEN
        IF (anoflag) THEN
          WRITE(*,'(a)') ' CC (Iobs+ - Iobs-) vs (Iref+ - Iref-) (Iref+/- is squared Fref+/-), in shells:'
          WRITE(*,'(40f7.1)') 100*get_cc(cc_overallano(:))
          WRITE(*,'(a)') ' # of reference reflections in resolution shells (with -s option: only acentric refs!):'
          WRITE(*,'(40i7)') cc_overallano(:)%n
        END IF
      END IF
      WRITE(*,*)
      IF (abcdef(4)) THEN
        WRITE(*,*)'headings for lines starting with d,e,f:'
        IF (is_xds.AND.nbatch==1) THEN
          WRITE(*,'(a)')'d: <CC1/2ano> of each frame:'
          WRITE(*,'(a)')'d:   reflections of this frame only        reflections of all frames'
          WRITE(*,'(a)')'d: frm   nref    with without   delta     nref    with without   delta'
        ELSE IF (is_xds) THEN
          WRITE(*,'(a)')'d: <CC1/2ano> of each batch (batch # is given as highest included frame #):'
          WRITE(*,'(a)')'d:   reflections of this batch only        reflections of all batches'
          WRITE(*,'(a)')'d: batch nref    with without   delta     nref    with without   delta'
        ELSE
          WRITE(*,'(a)')'d: <CC1/2ano> of each dataset:'
          WRITE(*,'(a)')'d: reflections of this dataset only      reflections of all datasets'
          WRITE(*,'(a)')'d: set   nref    with without   delta     nref    with without   delta'
        END IF          
      END IF
      IF (abcdef(5)) WRITE(*,'(a)') 'e: delta-CC1/2ano in resolution shells'
      IF (abcdef(6)) WRITE(*,'(a)') 'f: # reflections for delta-CC1/2ano'     
      
! anomalous signal:
      DO i=1,nset
        IF (((.NOT.abcdef(9)).OR.(.NOT.abcdef(10))).AND.SUM(cc_onlyano(i,:)%n)>0) THEN
!          WRITE(*,fmt) 'r',100*get_cc(cc_only(i,:)),' avg:', &
!            DOT_PRODUCT(100*get_cc(cc_only(i,:)),cc_only(i,:)%n)/SUM(cc_only(i,:)%n)
          IF (anoflag) THEN 
            WRITE(*,fmtr) 's',i*nbatch,100*get_cc(cc_onlyano(i,:)),' avg:', &
            DOT_PRODUCT(100*get_cc(cc_onlyano(i,:)),cc_onlyano(i,:)%n)/MAX(1,SUM(cc_onlyano(i,:)%n))
            WRITE(*,fmti) '#',cc_onlyano(i,:)%n,' sum',SUM(cc_onlyano(i,:)%n)
          END IF
        END IF
        IF (SUM(nanoset(i,:))<2) CYCLE
! mode 1: only reflections present in dataset contribute to anoCC12:
        WHERE (nanoset(i,:)>1)  ! 30.7.15 "floating divide by zero" wenn >0 
!        xxr=(anoset(2,i,:)-anoset(1,i,:)**2/nanoset(i,:)) !/ (nanoset(i,:)-1)  ! var_tau +0.5*var_eps, with
!	varr=anoset(3,i,:) !/(nanoset(i,:)-1)     ! half-dataset variance of epsilon, with
! bugfix 6.9.18:
        xxr=(anoset(2,i,:)-anoset(1,i,:)**2/nanoset(i,:)) / (nanoset(i,:)-1)  ! var_tau +0.5*var_eps, with
        varr=anoset(3,i,:)/nanoset(i,:)     ! half-dataset variance of epsilon, with
        xxr = (xxr - 0.5*varr)/(xxr + 0.5*varr)

!        yyr=(anoset(5,i,:)-anoset(4,i,:)**2/nanoset(i,:)) !/ (nanoset(i,:)-1) ! var_tau +0.5*var_eps, without
!	varr=anoset(6,i,:) ! /(nanoset(i,:)-1)   ! var_eps= half-dataset variance of epsilon, without
! bugfix 6.9.18:
        yyr=(anoset(5,i,:)-anoset(4,i,:)**2/nanoset(i,:)) / (nanoset(i,:)-1) ! var_tau +0.5*var_eps, without
        varr=anoset(6,i,:)/nanoset(i,:)  ! var_eps= half-dataset variance of epsilon, without
        yyr = (yyr - 0.5*varr)/(yyr + 0.5*varr)
        ELSE WHERE
        xxr=0.
        yyr=0.
        END WHERE
        
! mode 2: all reflections contribute to anoCC12:
! remove contributions from uniques where this dataset contributed more than n-2 observations
        nanoset2(i,:)=nanotot - nanoset2(i,:)

! anoCC12 without this dataset
        anoset(4,i,:)=anoset(4,i,:) - anoset(1,i,:) + anosum -anoset(7,i,:)
        anoset(5,i,:)=anoset(5,i,:) - anoset(2,i,:) + anosum2-anoset(8,i,:)
        anoset(6,i,:)=anoset(6,i,:) - anoset(3,i,:) + anovar -anoset(9,i,:)   
        WHERE (nanoset2(i,:)>1)   ! bugfix 9.10.20 . was : >0 
!        wwr=(anoset(5,i,:)-anoset(4,i,:)**2/nanoset2(i,:)) !/ (nanoset2(i,:)-1)  ! var_tau +0.5*var_eps, with
!        varr=anoset(6,i,:) !/(nanoset2(i,:)-1)     ! var_eps= half-dataset variance of epsilon, with
! bugfix 6.9.18:
        wwr=(anoset(5,i,:)-anoset(4,i,:)**2/nanoset2(i,:)) / (nanoset2(i,:)-1)  ! var_tau +0.5*var_eps, with
        varr=anoset(6,i,:)/nanoset2(i,:)     ! var_eps= half-dataset variance of epsilon, with
        wwr = (wwr - 0.5*varr)/(wwr + 0.5*varr)

! anoCC12 with this dataset:
        anoset(1,i,:)=anosum - anoset(7,i,:)
        anoset(2,i,:)=anosum2- anoset(8,i,:)
        anoset(3,i,:)=anovar - anoset(9,i,:)
!        vvr=(anoset(2,i,:)-anoset(1,i,:)**2/nanoset2(i,:)) !/ (nanoset2(i,:)-1)  ! var_tau +0.5*var_eps, with
!        varr=anoset(3,i,:) !/(nanoset2(i,:)-1)     ! half-dataset variance of epsilon, with
! bugfix 6.9.18:
        vvr=(anoset(2,i,:)-anoset(1,i,:)**2/nanoset2(i,:)) / (nanoset2(i,:)-1)  ! var_tau +0.5*var_eps, with
        varr=anoset(3,i,:)/nanoset2(i,:)     ! half-dataset variance of epsilon, with
        vvr = (vvr - 0.5*varr)/(vvr + 0.5*varr)
        ELSE WHERE
        vvr=0.
        wwr=0.
        END WHERE
        
        xx=DOT_PRODUCT(xxr,nanoset(i,:))/SUM(nanoset(i,:))    
        yy=DOT_PRODUCT(yyr,nanoset(i,:))/SUM(nanoset(i,:))
        vv=DOT_PRODUCT(vvr,nanoset2(i,:))/SUM(nanoset2(i,:))   
        ww=DOT_PRODUCT(wwr,nanoset2(i,:))/SUM(nanoset2(i,:))  
        IF (.NOT.abcdef(8)) THEN  ! Fisher transform switched off
!          IF (abcdef(4)) WRITE(*,'(a,i6,i7,3f8.3,i9,3f8.3)') 'd',nbatch*i, &
          IF (abcdef(4)) WRITE(*,'(a,i6,i7,2f8.3,f9.4,i9,2f8.3,f9.4)') 'd',nbatch*i, &  ! 2.7.22 KD
            SUM(nanoset(i,:)),100*xx,100*yy,100*(xx-yy),  &
            SUM(nanoset2(i,:)),100*vv,100*ww,100*(vv-ww)
          save_d(i,1)=100*(xx-yy)
          save_d(i,2)=100*(vv-ww)
          IF (abcdef(5)) WRITE(*,'(a,40(1x,f0.3))') 'e',100*(xxr-yyr)
        ELSE
          IF (abcdef(4)) WRITE(*,'(a,i6,i7,3f8.3,i9,3f8.3)') 'd',nbatch*i, &
            SUM(nanoset(i,:)),100*xx,100*yy,100*TANH(ATANH_safe(xx)-ATANH_safe(yy)),  &
            SUM(nanoset2(i,:)),100*vv,100*ww,100*TANH(ATANH_safe(vv)-ATANH_safe(ww))
          save_d(i,1)=100*TANH(ATANH_safe(xx)-ATANH_safe(yy))
          save_d(i,2)=100*TANH(ATANH_safe(vv)-ATANH_safe(ww))
          IF (abcdef(5)) WRITE(*,'(a,40(1x,f0.3))') 'e',100*TANH(ATANH_safe(xxr)-ATANH_safe(yyr))
        END IF
        IF (abcdef(6)) WRITE(*,'(a,40i7)') 'f',nanoset(i,:)        
      END DO
!    
END SUBROUTINE reflection_file

PROGRAM deltacc12
      USE resolution_mod
      USE reference_mod
      USE xscaleitems_mod, ONLY: resmin,resmax,get_nset_reso,is_xds
      IMPLICIT NONE
      REAL :: reslim1=-1,reslim2=-2,batch=0.
      CHARACTER (LEN=132) :: arg='',name='',reference=''
      CHARACTER :: date*8,time*10,zone*5
      INTEGER :: unit,ier,nbin=10,mode=2,i,narg,nbatch
      LOGICAL :: abcdef(10)=.TRUE.,skiparg=.FALSE.,anoflag=.FALSE.,sort_anom=.FALSE.
      INTEGER,DIMENSION(8) :: itime1
      EXTERNAL REFLECTION_FILE, WRITE_XSCALEINPRENAMEME
 
      narg=COMMAND_ARGUMENT_COUNT()
      WRITE(*,'(a)') 'xdscc12 KD 2023-01-08. Academic use only; no redistribution. -h option shows options.'
      WRITE(*,'(a)') 'Please cite Assmann, G., Brehm, W., Diederichs, K. (2016) J.Appl.Cryst. 49, 1021-1028'  
      CALL DATE_AND_TIME(DATE=date,TIME=time,ZONE=zone,VALUES=itime1)
      CALL GET_COMMAND(arg)
      WRITE(*,'(a)') ' running '''//TRIM(arg)//''' on '//date//' at '//time(1:2)//':'//time(3:4)//':'//time(5:6)//' '//zone
!      IF (itime1(1)*10000 + itime1(2)*100 + itime1(3)>20211231) STOP 'expired'
      IF (narg==1) THEN
        CALL GET_COMMAND_ARGUMENT(1, arg)
        IF (arg(1:2)=='-h') skiparg=.TRUE.
      END IF
      IF (narg==0.OR.skiparg) THEN
        WRITE(*,'(a)') 'usage: xdscc12 [-dmin <lowres>] [-dmax <highres>] [-nbin <nbin>] [-mode <1 or 2>] ' // &
                       '[-<abcdefstwz>] [-r|-R <ref>] FILE_NAME'
        WRITE(*,'(a)') 'dmin (default 999A), dmax (default 1A) and nbin (default 10) have the usual meanings.'
        WRITE(*,'(a)') 'mode can be 1 (equal volumes of resolution shells) or 2 (increasing volumes; default).'
        WRITE(*,'(a)') '  -r: <ref> is ASCII reference file with lines: h,k,l,Fref or h,k,l,Fref+,Fref-'
        WRITE(*,'(a)') '  -R: <ref> is ASCII reference file with lines: h,k,l,Iref or h,k,l,Iref+,Iref-'
        WRITE(*,'(a)') '      -r and -R allow calculation of CC of isomorphous signal with reference'
        WRITE(*,'(a)') '  -s: read two columns from reference: Fref+, Fref- or Iref+, Iref-. '
        WRITE(*,'(a)') '      this allows calculation of CC of anomalous signal with that of reference'
        WRITE(*,'(a)') '  -A: sort INPUT_FILEs in XSCALE.INP.rename_me by anomalous instead of isomorphous delta-CC1/2'
        WRITE(*,*) ' -t: total oscillation (degree) to batch fine-sliced frames into'
        WRITE(*,*) 'FILE_NAME can be XDS or XSCALE reflection file'
        WRITE(*,'(a)') ' other options can be combined (e.g. -def), and switch the following OFF:'
        WRITE(*,*) ' -a: individual isomorphous summary values'
        WRITE(*,*) ' -b: individual (Fisher-transformed) delta-CC1/2 values'
        WRITE(*,*) ' -c: individual delta-CC1/2 reflection numbers'
        WRITE(*,*) ' -d: individual anomalous summary values'
        WRITE(*,*) ' -e: individual (Fisher-transformed) delta-CC1/2ano values'
        WRITE(*,*) ' -f: individual delta-CC1/2ano reflection numbers'
        WRITE(*,*) ' -w: weighting of intensities with their sigmas'
        WRITE(*,*) ' -z: Fisher transformation of delta-CC1/2 values'
        WRITE(*,*)
        WRITE(*,'(a)') 'The program writes a commented XSCALE.INP.rename_me that is sorted on delta-CC1/2'
        STOP
      END IF
 
      DO i=1,narg
        IF (skiparg) THEN
          skiparg=.FALSE.
          CYCLE
        END IF
        CALL GET_COMMAND_ARGUMENT(i, arg)
        IF      (arg(1:5)=='-DMIN'.OR.arg(1:5)=='-dmin') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) reslim2
          ELSE IF (i<narg) THEN
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) reslim2
            skiparg=.TRUE.
          ENDIF
        ELSE IF (arg(1:5)=='-DMAX'.OR.arg(1:5)=='-dmax') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) reslim1
          ELSE IF (i<narg) THEN
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) reslim1 
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:5)=='-NBIN'.OR.arg(1:5)=='-nbin') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) nbin
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) nbin 
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:5)=='-MODE'.OR.arg(1:5)=='-mode') THEN
          IF (arg(6:6)=='=') THEN
            READ(arg(7:),*) mode
          ELSE IF (i<narg) THEN 
            CALL GET_COMMAND_ARGUMENT(i+1, arg)
            READ(arg,*) mode 
            WRITE(*,'(a,i2)')' mode for resolution shells is',mode
            skiparg=.TRUE.
          END IF
        ELSE IF (arg(1:1)=='-') THEN
          IF (INDEX(arg,'a')>0) THEN
            abcdef(1)=.FALSE.
            WRITE(*,*)'option -a found'
          END IF
          IF (INDEX(arg,'b')>0) THEN
            abcdef(2)=.FALSE.
            WRITE(*,*)'option -b found: no per-dataset (Fisher-transformed) delta-CC1/2 values'
          END IF 
          IF (INDEX(arg,'c')>0) THEN
            abcdef(3)=.FALSE.
            WRITE(*,*)'option -c found: no per-dataset delta-CC1/2 reflection numbers'
          END IF
          IF (INDEX(arg,'d')>0) THEN
            abcdef(4)=.FALSE.
            WRITE(*,*)'option -d found: no per-dataset anomalous summary values'
          END IF 
          IF (INDEX(arg,'e')>0) THEN
            abcdef(5)=.FALSE.
            WRITE(*,*)'option -e found: no per-dataset (Fisher-transformed) delta-CC1/2ano values'
          END IF 
          IF (INDEX(arg,'f')>0) THEN
            abcdef(6)=.FALSE.
            WRITE(*,*)'option -f found: no per-dataset delta-CC1/2ano reflection numbers'
          END IF 
          IF (INDEX(arg,'r')>0) THEN
            abcdef(9)=.FALSE.
            IF (i<narg) THEN 
              CALL GET_COMMAND_ARGUMENT(i+1, reference)
            ELSE
              STOP 'expected name of reference file'
            END IF
            WRITE(*,*)'option -r found: amplitude reference file is ',TRIM(reference)
            WRITE(*,*)'lines starting with r: CC of Iobs vs Iref'
            skiparg=.TRUE.
          END IF 
          IF (INDEX(arg,'R')>0) THEN
            IF (.NOT.abcdef(9)) STOP 'cannot use both -r and -R'
            abcdef(10)=.FALSE.
            IF (i<narg) THEN 
              CALL GET_COMMAND_ARGUMENT(i+1, reference)
            ELSE
              STOP 'expected name of reference file'
            END IF
            WRITE(*,*)'option -R found: intensity reference file is ',TRIM(reference)
            WRITE(*,*)'lines starting with R: CC of Iobs vs Iref'
            skiparg=.TRUE.
          END IF 
          IF (INDEX(arg,'s')>0) THEN
            anoflag=.TRUE.
            WRITE(*,'(a)')' option -s found: reference file has lines with h,k,l,Fref+,Fref- or h,k,l,Iref+,Iref-'
            WRITE(*,'(a)')' lines starting with s: CC of anomalous signal with that of reference'
          END IF 
          IF (INDEX(arg,'t')>0) THEN
            IF (arg(6:6)=='=') THEN
              READ(arg(7:),*) batch
            ELSE IF (i<narg) THEN 
              CALL GET_COMMAND_ARGUMENT(i+1, arg)
              READ(arg,*) batch 
              WRITE(*,'(a,f5.1)')' option -t found: frames batched such that total rotation >= ',batch
              skiparg=.TRUE.
            END IF
          END IF 
          IF (INDEX(arg,'w')>0) THEN
            abcdef(7)=.FALSE.
            WRITE(*,*)'option -w found: weighting switched off'
          END IF 
          IF (INDEX(arg,'z')>0) THEN
            abcdef(8)=.FALSE.
            WRITE(*,*)'option -z found: no Fisher transform of delta-CC1/2 values'
          END IF 
          IF (INDEX(arg,'A')>0) THEN
            sort_anom=.TRUE.
            WRITE(*,*)'option -A found: sort INPUT_FILEs by anomalous delta-CC1/2 values'
          END IF 
        ELSE  ! file name
          name=arg
        END IF
      END DO 
      IF (.NOT.abcdef(9).AND.(.NOT.anoflag)) WRITE(*,'(a)') &
        ' no option -s found, therefore reference file expected to have h,k,l,Fref (both centric & acentric)'
      IF (.NOT.abcdef(10).AND.(.NOT.anoflag)) WRITE(*,'(a)') &
        ' no option -s found, therefore reference file expected to have h,k,l,Iref (both centric & acentric)'
      IF (abcdef(7)) WRITE(*,*) 'no option -w found, therefore statistics are weighted by sigma values'
      IF (abcdef(8)) WRITE(*,*) 'no option -z found, therefore delta-CC1/2 values are Fisher-transformed'
      WRITE(*,*)
      WRITE(*,'(a,a)') ' reflection file is ',TRIM(name)
      CALL get_nset_reso(name,batch,nbatch)   
      IF (reslim1<0) reslim1=resmin
      IF (reslim2<0) reslim2=resmax
      CALL setbin(nbin,max(reslim1,reslim2),min(reslim1,reslim2),mode)
      OPEN(NEWUNIT=unit,FILE=name,STATUS='old',ACTION='READ',iostat=ier)
      IF (ier/=0) THEN
        PRINT*,'ier=',ier
        STOP 'could not open file'
      END IF
      CALL reflection_file(unit,max(reslim1,reslim2),min(reslim1,reslim2),nbin,&
           abcdef,reference,nbatch,anoflag)
      CLOSE(unit)
      IF (.NOT.is_xds) CALL write_xscaleinprenameme(name,abcdef,sort_anom)
      STOP 'normal termination'
END PROGRAM deltacc12

SUBROUTINE write_xscaleinprenameme(name,abcdef,sort_anom)
      USE xscaleitems_mod, ONLY: save_a,save_d
      USE qsortfix_m
      CHARACTER, INTENT(IN) :: name*132
      LOGICAL, INTENT(IN) :: abcdef(10),sort_anom 
! local variables
      EXTERNAL select
      CHARACTER line*180
      CHARACTER(LEN=180), ALLOCATABLE :: inputfile(:)
      INTEGER i,nset,inunit,outunit
      REAL, ALLOCATABLE :: temp(:)
      INTEGER, ALLOCATABLE :: inx(:)
      REAL xmed,save_a_max !,weight
      
      OPEN(NEWUNIT=inunit,FILE=name,STATUS='old',ACTION='READ')
      OPEN(NEWUNIT=outunit,FILE='XSCALE.INP.rename_me')
      i=0
      nset=SIZE(save_a,DIM=1)
      ALLOCATE(inputfile(nset),temp(nset),inx(nset))
      IF (sort_anom) THEN
        temp=save_d(:,1)
      ELSE
        temp=save_a(:,1)
      END IF
      DO
        READ(inunit,'(a)') line
        IF (INDEX(line,'INPUT_FILE=')>0) THEN
          i=i+1
          inputfile(i)=line(INDEX(line,'INPUT_FILE='):)
        ELSE IF (INDEX(line,'!SPACE_GROUP_NUMBER=')==1) THEN
          WRITE(outunit,'(a)') TRIM(line(2:))
        ELSE IF (INDEX(line,'!UNIT_CELL_CONSTANTS=')==1) THEN
          WRITE(outunit,'(a)') TRIM(line(2:))
        ELSE IF (INDEX(line,'!END_OF_HEADER')==1) THEN
          EXIT
        END IF
      END DO
      IF (i/=nset) STOP 'i /= nset'
      WRITE(outunit,'(a)') 'OUTPUT_FILE= '//TRIM(name)
      WRITE(outunit,'(a)') 'PRINT_CORRELATIONS= FALSE'
      WRITE(outunit,'(a)') 'SAVE_CORRECTION_IMAGES= FALSE'
      WRITE(outunit,'(a)') 'FRIEDEL''S_LAW= FALSE'
      inx=(/(i,i=1,nset)/)
      CALL myqsort(nset,temp,inx)
      save_a_max=MAXVAL(save_a(:,1))
      WRITE(*,*)
      WRITE(*,*) 'best delta-CC1/2_only=',save_a_max

      i=(nset+1)/2
      
! median of 6th col of ^a lines
      temp=save_a(:,1)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'median of delta-cc1/2 ("only" i.e. 6th col of "a" lines) =',xmed
      WRITE(outunit,'(a,f10.3)') '! median of delta-cc1/2 "only" values=',xmed

! median absolute deviation of 6th col of ^a lines
      temp=ABS(temp-xmed)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'noise= (MAD, median absolute deviation) from this median =',xmed
      WRITE(outunit,'(a,f10.3)') '! noise (MAD) of these values=',xmed

! median of 10th col of ^a lines
      temp=save_a(:,2)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'median of delta-cc1/2 ("all" i.e. 10th col of "a" lines) =',xmed
      WRITE(outunit,'(a,f10.3)') '! median of delta-cc1/2 "all" values=',xmed

! median absolute deviation of 10th col of ^a lines
      temp=ABS(temp-xmed)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'noise= (MAD, median absolute deviation) from this median =',xmed
      WRITE(outunit,'(a,f10.3)') '! noise (MAD) of these values=',xmed

! same for anomalous:

! median of 6th col of ^d lines
      temp=save_d(:,1)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'median of delta-cc1/2-ano ("only" i.e. 6th col of "d" lines) =',xmed
      WRITE(outunit,'(a,f10.3)') '! median of delta-cc1/2-ano "only" values=',xmed

! median absolute deviation of 6th col of ^d lines
      temp=ABS(temp-xmed)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'noise= (MAD, median absolute deviation) from this median =',xmed
      WRITE(outunit,'(a,f10.3)') '! noise (MAD) of these values=',xmed

! median of 10th col of ^d lines
      temp=save_d(:,2)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'median of delta-cc1/2-ano ("all" i.e. 10th col of "d" lines) =',xmed
      WRITE(outunit,'(a,f10.3)') '! median of delta-cc1/2-ano "all" values=',xmed

! median absolute deviation of 10th col of ^d lines
      temp=ABS(temp-xmed)
      CALL select(temp,nset,i)
      xmed=temp(i)
      IF (MOD(nset,2)==0) xmed=(xmed + MINVAL(temp(i+1:nset)))/2.
      WRITE(*,*) 'noise= (MAD, median absolute deviation) from this median =',xmed
      WRITE(outunit,'(a,f10.3)') '! noise (MAD) of these values=',xmed

      WRITE(*,*)
      IF (sort_anom) THEN
        WRITE(*,'(a)') 'Wrote a commented XSCALE.INP.rename_me that is sorted on'// &
                     ' delta-CC1/2-anom "only"'
      ELSE
        WRITE(*,'(a)') 'Wrote a commented XSCALE.INP.rename_me that is sorted on'// &
                     ' delta-CC1/2 "only"'
      END IF
      WRITE(*,'(a)') 'You may edit that file, or e.g. add a line after each '// &
           'INPUT_FILE line with sed ''/INPUT_FILE/a INCLUDE_RESOLUTION_RANGE=99 3''' 
!      WRITE(*,'(a)') 'Consider trying the WEIGHT values; these are '// &
!           'exp(ln(2)*(deltacc12_only/deltacc12_only_max-1))'
      WRITE(outunit,'(a)') '! input files sorted by deltacc12_only (highest first):'
      DO i=nset,1,-1
        IF (abcdef(4)) THEN
          WRITE(outunit,'(a,2f10.4,a,2f10.4)') '! deltacc12 only / all:',save_a(inx(i),:), &
                                             ' deltacc12-ano only /all:',save_d(inx(i),:)
        ELSE
          WRITE(outunit,'(a,2f10.4)') '! deltacc12 only / all:',save_a(inx(i),:)
        END IF
        WRITE(outunit,'(a)') TRIM(inputfile(inx(i)))
!        weight = LOG(2.) * (save_a(inx(i),1)/save_a_max-1)**2
!        IF (weight>LOG(1000.)) weight=LOG(1000.)
!        WRITE(outunit,'(a,f10.3)') '!WEIGHT=',EXP(weight) 
      END DO
      CLOSE(inunit)
      CLOSE(outunit)
END SUBROUTINE write_xscaleinprenameme     

SUBROUTINE select(arr,n,k)
! given k, return arr partially sorted such that the k-th element in sort order
! is at position k of arr
! mathematically: arr(i) <= arr(k) for i<k, and arr(i) >= arr(k) for i>k 
! adapted from Numerical Recipes in Fortran90
        IMPLICIT NONE

!       INTERFACE 
!         SUBROUTINE SWAP(A,B,MASK)
!           REAL(KIND=4), INTENT(INOUT) :: A
!           REAL(KIND=4), INTENT(INOUT) :: B
!           LOGICAL(KIND=4) ,OPTIONAL, INTENT(IN) :: MASK
!         END SUBROUTINE SWAP
!       END INTERFACE 

        INTEGER, INTENT(IN) :: k
        REAL, INTENT(INOUT) :: arr(n)
        INTEGER :: i,r,j,l,n
        REAL :: a

        l=1
        r=n
        do
          if (r-l <= 1) then
            if (r-l == 1) call swap(arr(l),arr(r),arr(l)>arr(r))
            RETURN
          else
! use median of (left,middle,right) element as pivot
            i=(l+r)/2 
            call swap(arr(i),arr(l+1))
            call swap(arr(l),arr(r),arr(l)>arr(r))
            call swap(arr(l+1),arr(r),arr(l+1)>arr(r))
            call swap(arr(l),arr(l+1),arr(l)>arr(l+1))
            i=l+1
            j=r
            a=arr(l+1)
            do
              do
                i=i+1
                if (arr(i) >= a) exit
              end do
              do
                j=j-1
                if (arr(j) <= a) exit
              end do
              if (j < i) exit
              call swap(arr(i),arr(j))
            end do
            arr(l+1)=arr(j)
            arr(j)=a
            if (j >= k) r=j-1
            if (j <= k) l=i
          end if
        end do
contains
        SUBROUTINE swap(a,b,mask)
        REAL, INTENT(INOUT) :: a,b
        LOGICAL, INTENT(IN), OPTIONAL :: mask
        REAL :: swp

        IF (PRESENT(mask)) THEN
          IF (.NOT.mask) RETURN
        END IF
        swp=a
        a=b
        b=swp
        END SUBROUTINE swap
	
END subroutine select
 
