! 
! Input is a file in XSCALE format 
! version 0.1 Kay Diederichs 19.10.2015
!
! ifort -check all,noarg_temp_created -traceback -warn -fpe0 -init=snan,arrays \
! -static-intel xscale_maxcc12.f90 symops.o  cc_mod.o resolution.o
! or 
! gfortran -O2 -C -Wall cc12_maxcc12.f90 symops.o cc_mod.o resolution.o
! with previous compilation like: ifort -c symops.f90
!

MODULE reference_mod
      REAL, ALLOCATABLE :: ref(:,:,:)
CONTAINS
      SUBROUTINE read_reference(reference,reslim1,reslim2)
      USE ISO_FORTRAN_ENV  ! to get access to iostat_end
      USE symops_mod,ONLY: op,ncenop,asuput
      USE resolution_mod, ONLY: getres
      CHARACTER (LEN=132), INTENT(IN) :: reference
      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
      
      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: ',TRIM(line)
          CYCLE
        END IF
        CALL getres(ih,ik,il,s2,resol)
        IF (resol>reslim2.OR.resol<reslim1) 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=-1.
! fill array
      n=0
      DO
        READ(unit,*,iostat=ier) ih,ik,il,amplitude
        IF (ier == iostat_end) EXIT
        IF (ier /= 0) CYCLE
        CALL getres(ih,ik,il,s2,resol)
        IF (resol>reslim2.OR.resol<reslim1) CYCLE
        CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,size(op,3)/ncenop)
        IF (ref(jh,jk,jl) < 0) THEN
          ref(jh,jk,jl)=amplitude**2
        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
END MODULE reference_mod

MODULE xscaleitems_mod
      INTEGER nset
      REAL resmin,resmax
      CHARACTER*132, ALLOCATABLE :: inputfile(:)
contains
      SUBROUTINE get_nset_reso(name)
      USE resolution_mod
      IMPLICIT NONE
      CHARACTER*(*), INTENT(IN) :: name
      CHARACTER line*80
      INTEGER :: ier,unit,iset,ih,ik,il
      REAL :: s2,resol
 1000 FORMAT(a)
      
      nset=0
      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=',TRIM(line)
          STOP 'read error in xscale_file'
        END IF
        IF (line(1:7)=='! ISET=') then
          READ(line(8:),*) iset
          nset=max(nset,iset)
        ELSE IF (line(1:21).eq.'!UNIT_CELL_CONSTANTS=') then 
          READ(line(22:),*) a,b,c,alf,bet,gam
        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
          PRINT*,'ier=',ier
          PRINT*,'line=',TRIM(line)
          STOP 'read error in xscale_file'
        END IF
        IF (line(1:12) == '!END_OF_DATA') EXIT
        READ(line,*)ih,ik,il
        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

SUBROUTINE xscale_file(unit,reslim1,reslim2,nbin,abcdef,reference)
      USE symops_mod, ONLY: asuput,getops,is_centric,op,ncenop,hkluncompress
      USE resolution_mod, ONLY: a,b,c,alf,bet,gam,getres,getbin,countbin
      USE cc_mod, ONLY: get_cc,put_cc,cc_type
      USE reference_mod, ONLY: ref,read_reference
      IMPLICIT NONE
     
      INTEGER, INTENT(IN) :: unit,nbin
      REAL, INTENT(IN) :: reslim1,reslim2
      LOGICAL, INTENT(IN) :: abcdef(9)
      CHARACTER (LEN=132),INTENT(IN) :: reference
! local variables
      INTEGER :: ih,ik,il,oldhkl,jh,jk,jl,ihkl,iop,i,j,ispgr,ncheck,minset, &
                 ier,iset,maxset=0,nunique,nobs,ibin,firstframe=0,lastframe=0
      INTEGER, ALLOCATABLE :: uniquecount(:,:),nap(:),nam(:),niso(:,:),nano(:,:),shellcount(:)
      REAL, ALLOCATABLE :: xap(:),xam(:),xap2(:),xam2(:),variso(:,:),varano(:,:), &
                           sigap(:),sigam(:),sumvar(:,:)
      REAL(8), ALLOCATABLE :: sumiso(:,:),sumiso2(:,:),sumano(:,:),sumano2(:,:)
      REAL s2,s2save,resol,fsq,sigma,dum4(4),xx,yy,ref_int
      CHARACTER :: line*80,format*20
      LOGICAL :: first=.TRUE.,last=.FALSE.,centric,xscale=.TRUE.
      TYPE(cc_type), ALLOCATABLE :: cc_ref(:,:)

 1000 FORMAT(A)
 
      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   ! XSCALE.HKL
          ncheck=ncheck+1
        ELSE IF (line(1:10) == '!ITEM_ZD=8') THEN     ! XDS_ASCII.HKL
          ncheck=ncheck+1
	ELSE IF (line(1:7) == '! ISET=') THEN
	  READ(line(8:),*) iset
	  maxset=max(maxset,iset)
! FIXME: miniset sollte auch berechnet werden, und alle Schleifen von
! miniset bis maxiset gehen statt von 1 bis maxiset. miniset>1 passiert,
! wenn XSCALE mit mehr als 1 OUTPUT_FILE laeuft
	ELSE IF (line(1:12) == '!DATA_RANGE=') THEN  ! XDS_ASCII.HKL !!
	  READ(line(13:),*) firstframe,lastframe
        ELSE IF (line(1:14) == '!END_OF_HEADER') THEN
          EXIT
        END IF
      END DO
! XSCALE: h,k,l,i,sigma,iset
! XDS_ASCII: h,k,l,i,sigma,zd
      IF (ncheck < 6) stop 'header items not found'
      IF (maxset==0) THEN   ! single dataset; use ZD instead of ISET
        xscale=.FALSE.
!! find max ZD
!        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:12) == '!END_OF_DATA') EXIT
!          READ(line,*)ih,ik,il,fsq,sigma,dum4(1:3)
!          maxset=MAX(maxset,1+INT(dum4(3)))    ! ITEM_ZD
!        END DO      
!        REWIND(unit)
!        DO
!          READ(unit,1000,iostat=ier) line
!          IF (line(1:14) == '!END_OF_HEADER') EXIT
!        END DO
        IF (firstframe*lastframe==0) THEN
          STOP 'did not find DATA_RANGE information in XDS_ASCII.HKL'
        END IF
        maxset=lastframe
        PRINT*,'number of frames=',maxset
      ELSE
        PRINT*,'number of datasets=',maxset
      END IF

      IF (.NOT.abcdef(9)) CALL read_reference(reference,reslim1,reslim2)

      ALLOCATE(cc_ref(nbin,maxset),xap(maxset),xap2(maxset),xam(maxset),xam2(maxset), &
               sumano(nbin,maxset),sumano2(nbin,maxset),sumiso(nbin,maxset), &
               sumiso2(nbin,maxset),varano(nbin,maxset),nap(maxset),nam(maxset),&
               sigap(maxset),sigam(maxset),uniquecount(nbin,maxset), &
               niso(nbin,maxset),nano(nbin,maxset),variso(nbin,maxset), &
               shellcount(nbin),sumvar(nbin,maxset))
! initialize
      xap=0;xap2=0;xam=0;xam2=0;sigap=0;sigam=0;nap=0;nam=0;sumvar=0
      sumiso=0;sumiso2=0;variso=0;sumano=0;sumano2=0;varano=0;niso=0;nano=0
      uniquecount=0
      oldhkl=0
      nunique=0     ! total # uniques
      nobs=0     ! total # observations
      minset=HUGE(minset)
! read reflections
      DO
        READ(unit,1000) line
        last = line(1:12) == '!END_OF_DATA'
        IF (.NOT.last) THEN
          IF (xscale) THEN
            READ(line,*)ih,ik,il,fsq,sigma,dum4,iset
          ELSE                   ! XDS_ASCII.HKL
            READ(line,*)ih,ik,il,fsq,sigma,dum4(1:3)
            iset=1+INT(dum4(3))    ! ITEM_ZD
          END IF
          IF (sigma<0.) CYCLE
          IF (.NOT.abcdef(7)) sigma=1.
	  CALL getres(ih,ik,il,s2,resol)
	  IF (resol>reslim2.OR.resol<reslim1) CYCLE
	  CALL asuput(ih,ik,il,jh,jk,jl,ihkl,iop,op,size(op,3)/ncenop)
          IF (nobs==0) THEN
            s2save=s2
            oldhkl=ihkl
          END IF
          nobs=nobs+1
        END IF
	IF (ihkl /= oldhkl.OR.last) THEN ! process the _previous_ unique reflection
	  nunique=nunique+1
          ref_int=-1
          IF (.NOT.abcdef(9)) THEN  ! look at reference intensities
            CALL hkluncompress(oldhkl,ih,ik,il) ! ih,ik,il needed for ref(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)) ref_int=ref(ih,ik,il)
          END IF
! find resolution bin
	  CALL getbin(s2save,ibin)
! store contributions to CC1/2
!          DO j=1,maxset   ! faster: start the loop at MIN(iset-values of this refl)
!!$omp parallel do private(j,xx,yy) schedule(guided) ! bringt nichts - zuwenig
          DO j=minset,maxset
            IF (nap(j)+nam(j)>0) THEN
              uniquecount(ibin,j)=uniquecount(ibin,j)+1
              xx=(xap(j)+xam(j))/(sigap(j)+sigam(j))    ! weighted mean intensity
              IF (ref_int>0) CALL put_cc(cc_ref(ibin,j),xx,ref_int)
! isomorphous signal: contribution to CC1/2
              IF (nap(j)+nam(j)>1) THEN
                sumiso(ibin,j) = sumiso(ibin,j) + xx
                sumiso2(ibin,j)= sumiso2(ibin,j)+ DPROD(xx,xx)
                yy=(xap2(j)+xam2(j))/(sigap(j)+sigam(j)) ! weighted mean intensity**2          
	        xx    = (yy-xx**2)/(nap(j)+nam(j)-1)   ! variance 
	        variso(ibin,j) = variso(ibin,j) + xx   ! sum of variances of half-dataset
                niso(ibin,j) = niso(ibin,j) + 1
                sumvar(ibin,j) = sumvar(ibin,j) + 1/(sigap(j)+sigam(j)) 
! anomalous signal: contribution to CC1/2ano
                IF (nap(j)>1.AND.nam(j)>1) THEN
                  xap(j) = xap(j)/sigap(j)
                  xam(j) = xam(j)/sigam(j)
                  xap2(j) = xap2(j)/sigap(j)
                  xam2(j) = xam2(j)/sigam(j)
                  xx = xap(j)-xam(j)
                  sumano(ibin,j)  = sumano(ibin,j)  + xx
                  sumano2(ibin,j) = sumano2(ibin,j) + DPROD(xx,xx)
                  xx = (xap2(j)-xap(j)**2)/(nap(j)-1) ! variance of +half-dataset
                  yy = (xam2(j)-xam(j)**2)/(nam(j)-1) ! variance of -half-dataset
                  varano(ibin,j) = varano(ibin,j) + xx + yy ! sum of variances of half-datasets
                  nano(ibin,j) = nano(ibin,j) + 1
                END IF
              END IF
              xap(j)=0
              xam(j)=0
              xap2(j)=0
              xam2(j)=0
              nap(j)=0
              nam(j)=0
              sigam(j)=0
              sigap(j)=0
            END IF
	  END DO
          IF (last) EXIT  ! the big loop reading the XSCALE file    
          first=.TRUE.
          oldhkl=ihkl
          s2save=s2
          minset=HUGE(minset)
        END IF
! store this observation
        minset=MIN(iset,minset)
!        sigma=1/sigma**2
        IF (first) THEN   ! happens once for each unique ref
          centric=is_centric(jh,jk,jl,op,size(op,3)/ncenop)
          first=.FALSE.
        END IF
        IF (iop>0.OR.centric) THEN  ! this is taken as I+
          DO j=iset,maxset ! observation contributes to this ISET, and following
            xap(j) = xap(j) + fsq/sigma**2
            xap2(j)= xap2(j) + fsq**2/sigma**2
            nap(j) = nap(j) + 1
            sigap(j) = sigap(j) + 1/sigma**2
          END DO
        ELSE                        !  I-
          DO j=iset,maxset
            xam(j) = xam(j) + fsq/sigma**2
            xam2(j)= xam2(j) + fsq**2/sigma**2
            nam(j) = nam(j) + 1
            sigam(j) = sigam(j) + 1/sigma**2
          END DO
        END IF
      END DO ! outer_loop over observations
      PRINT*,'# obs (excluding misfits), # unique =',nobs,nunique
      PRINT*,' '
      PRINT*,'a = CC1/2_iso as function of resolution' !, and difference to previous'
      PRINT*,'b = number of unique reflections used for CC1/2_iso'
      PRINT*,'c = completeness (disregarding screw axes)'
      PRINT*,'d = CC1/2_ano as function of resolution' !, and difference to previous'
      PRINT*,'e = number of unique reflections used for CC1/2_ano'
      IF (.NOT.abcdef(9)) PRINT*,'f = CC against reference'
!      PRINT*,'g = CC1/2_iso expected if all datasets were isomorphous'
!      PRINT*,'h = CC1/2_iso expected taking non-isomorphism into account'
      PRINT*,' '
!      PRINT*,hmin,hmax,kmin,kmax,lmin,lmax
      CALL countbin(ispgr,shellcount)
!print*,shellcount
      IF (.NOT.abcdef(9)) THEN
        PRINT*,'ISET  a,b,c,d,e,f,g,h in resolution shells:'
      ELSE
        PRINT*,'ISET  a,b,c,d,e,g,h in resolution shells:'
      END IF
      PRINT*,' '
      DO i=1,maxset
!      print*,sumiso(:,i)
!      print*,sumvar(:,i)
!      end do
!      stop
!      DO
        DO j=1,nbin
          IF (niso(j,i)<=1) THEN
            variso(j,i) = 0
            niso(j,i) = 0
            sumvar(j,i) = 0
            CYCLE
          END IF
          sumiso2(j,i) = (sumiso2(j,i) - sumiso(j,i)**2/niso(j,i)) / (niso(j,i)-1)
          variso(j,i) = variso(j,i)/niso(j,i)
! CC1/2_iso
          variso(j,i) = (sumiso2(j,i)-variso(j,i))/(sumiso2(j,i)+variso(j,i))
! expected CC1/2_iso calculated from <I>/<sigma>: formula in Box 1 of K&D 2015
          sumvar(j,i) = 1/(1+4*niso(j,i)/sumiso(j,i)**2*sumvar(j,i))
	END DO
!        format='(i5,10f6.1,f7.2,a)'
        format='(i5,10f6.1,a)'
        WRITE(format(5:6),'(i2.2)')nbin
        WRITE(*,format) i,100*variso(:,i),' a'
        format='(5x,10i6  ,a)'
        WRITE(format(5:6),'(i2.2)')nbin
        WRITE(*,format) niso(:,i),' b'
        format='(5x,10f6.1,a)'
        WRITE(format(5:6),'(i2.2)')nbin
        WRITE(*,format) (100.*uniquecount(:,i))/shellcount,' c'
        DO j=1,nbin
          IF (nano(j,i)<=1) THEN
            varano(j,i) = 0
            nano(j,i) = 0
            CYCLE
          END IF
          sumano2(j,i) = (sumano2(j,i) - sumano(j,i)**2/nano(j,i)) / (nano(j,i)-1)
          varano(j,i) = varano(j,i)/nano(j,i)
! CC1/2_ano
          varano(j,i) = (sumano2(j,i)-varano(j,i))/(sumano2(j,i)+varano(j,i))
	END DO
!        format='(5x,10f6.1,f7.2,a)'
        format='(5x,10f6.1,a)'
        WRITE(format(5:6),'(i2.2)')nbin
        WRITE(*,format) 100*varano(:,i),' d'
        format='(5x,10i6  ,a)'
        WRITE(format(5:6),'(i2.2)')nbin  
        WRITE(*,format) nano(:,i),' e'
        format='(5x,10f6.1,a)'
        WRITE(format(5:6),'(i2.2)')nbin        
        IF (.NOT.abcdef(9)) WRITE(*,format) 100*get_cc(cc_ref(:,i)),' f'
        format='(5x,10f6.1,a)'
        WRITE(format(5:6),'(i2.2)')nbin
!        WRITE(*,format) 100*sumvar(:,i),' g'
      END DO
!    
END SUBROUTINE xscale_file

PROGRAM xscale_maxcc12
      USE resolution_mod, ONLY: getlim,setbin
      USE xscaleitems_mod, ONLY: get_nset_reso,resmin,resmax
      IMPLICIT NONE
      REAL :: reslim1=0.,reslim2=0.
      REAL, ALLOCATABLE :: vvr(:)
      CHARACTER (LEN=132) :: arg='',name='',reference=''
      INTEGER :: unit,ier,nbin=10,mode=2,i,narg
!      INTEGER,DIMENSION(8) :: itime1
      LOGICAL :: abcdef(9)=.TRUE.,skiparg=.FALSE.

      narg=COMMAND_ARGUMENT_COUNT()
      WRITE(*,'(a)') 'xscale_maxcc12 KD 2019-02-06. -h option shows options'
      WRITE(*,'(a)') 'Academic use only; no redistribution.'
!      CALL DATE_AND_TIME(VALUES=itime1)
!      IF (itime1(1)*10000 + itime1(2)*100 + itime1(3)>20181231) 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: xscale_maxcc12 -dmin <lowres> -dmax <highres> -nbin <nbin> -mode <1 or 2> -<rx> XSCALE_FILE_NAME'
        WRITE(*,'(a)') 'dmax, dmin and nbin (default 10) have the usual meanings.'
        WRITE(*,'(a)') 'mode can be 1 (equal volumes of resolution shells) or 2 (increasing volumes of resolution shells; default).'
        WRITE(*,*) '  -x: no weighting of intensities with their sigmas'
        WRITE(*,*) '  -r: next parameter: ASCII reference file with h,k,l,Fcalc'
        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,'x')>0) THEN
            abcdef(7)=.FALSE.
            WRITE(*,*)'option -x found: weighting switched off'
          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: reference file is ', &
                            TRIM(reference)
            skiparg=.TRUE.
          END IF 
        ELSE  ! file name
          name=arg
        END IF
      END DO 
      ALLOCATE(vvr(nbin))
      IF (abcdef(7)) WRITE(*,*) 'statistics weighted by sigma values'
      WRITE(*,*)
      WRITE(*,'(a,a)') 'XSCALE file is ',TRIM(name)
      CALL get_nset_reso(name)
      IF (reslim1==0.) reslim1=resmin
      IF (reslim2==0.) reslim2=resmax
      PRINT*,'reslim1,reslim2=',reslim1,reslim2
      CALL setbin(nbin,max(reslim1,reslim2),min(reslim1,reslim2),mode)
      WRITE(*,*) 'resolution shells:'
      CALL getlim(vvr)
      PRINT '(5x,40f6.3)',vvr
      
      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 xscale_file(unit,reslim1,reslim2,nbin,abcdef,reference)
      CLOSE(unit)
END PROGRAM xscale_maxcc12
