      PROGRAM REFINE
************************************************************************
*         P R O F I L E   R E F I N E M E N T   P R O G R A M          *
*======================================================================*
*                       Version  3-JAN-1989                            *
*                                                                      *
************************************************************************
*************************  H I S T O R Y  ******************************
************************************************************************
*                  H. M. Rietveld - December 1970.                     *
************************************************************************
*    Modified December 1972 for anisotropic temperature factors by     *
*    A. W. Hewat, UKAERE Harwell, Didcot, Berkshire, England.          *
************************************************************************
*    Modified for Voigtian shape function on 1983 by M. Nurmela,       *
*    University of Helsinki.                                           *
************************************************************************
*    Modified June 1986 for x-ray data by M. Nurmela,                  *
*    Imatran Voima Oy, Helsinki                                        *
************************************************************************
*    Modified July 1986 by M. Nurmela, Imatran Voima Oy, Helsinki      *
*    The concept 'number of overlapping reflections' is obsolete.      *
************************************************************************
*    Modified July 1987 by M. Nurmela, Imatran Voima Oy, Helsinki      *
*    Texture correction is changed according to M.Jarvinen.
*    Preparation program attached as a subroutine.                     *
************************************************************************
*    Modified january 1989 by M. Nurmela, Imatran Voima Oy, Helsinki   *
*    A reference to the theory of texture correction is given. No      *
*    changes in code between versions 8-OCT-1987 and 3-JAN-1989.       *
************************************************************************
*                         VARIABLE DIMENSIONS                          *
*                         -------------------                          *
*           B=MAXIMUM NUMBER OF BACKGROUND POSITIONS                   *
*           R=MAXIMUM NUMBER OF EQUIVALENT POSITIONS + 1               *
*           K=MAXIMUM NUMBER OF SCATTERING LENGTHS                     *
*           F=MAXIMUM NUMBER OF FORMFACTORS                            *
*           N=MAXIMUM NUMBER OF ATOMS                                  *
*           M=MAXIMUM NUMBER OF MAGNETIC ROTATION MATRICES             *
*           P=MAXIMUM NUMBER OF LEAST SQUARES PARAMETERS               *
*           C=MAXIMUM NUMBER OF CONSTRAINT FUNCTIONS                   *
*           T=MAXIMUM NUMBER OF TERMS IN CONSTRAINT FUNCTION           *
*           I=MAXIMUM NUMBER OF REFLECTIONS                            *
*           L=MAXIMUM NUMBER OF EXCLUDED REGIONS                       *
*           X=MAXIMUM NUMBER OF PROFILE INTENSITIES                    *
*          PC=P+C                                                      *
*                                   *****                              *
*                      Arrays with variable dimensions                 *
*                      -------------------------------                 *
*          ICODE(I)            IH(3,I)             IHIGH(L)            *
*          IL(C,T)             INDEX(I)            IORD(X)             *
*          IQ(C,T)             JCODE(I)            JL(C,T)             *
*          JQ(C,T)             LL(I)               LN(C)               *
*          LOW(L)              M(I)                MEQ(N)              *
*          MTYP(N)             NTYP(N)             ATEXT(N)            *
*          BAK(B)              BRAGG(I)            C(C,PC)             *
*          CF(C)               COEF(C,T)           COSAR(R,N)          *
*          COSTHE(I)           COSPHI(I)           DERIV(PC)           *
*          DERSTO(I,PC)        DF(C)               DISPIM(K)           *
*          DISPRE(K)           DUMP(I)             EQUIV(4,3,R)        *
*          F(21,F)             FAC(I)              FESD(I)             *
*          FF(I,F)             FMAG(I)             FNUC(I)             *
*          FOBS(I)             H(R,3)              HKL(3,I)            *
*          HW(I)               PLOR(I)             POS(B)              *
*          RM(9,M,R)           RSA(N,3,3)          RSB(N,3,3)          *
*          S(21,F)             SA(N)               SB(N)               *
*          SCAT(K)             SINAR(R,N)          SJJ(I)              *
*          SMM(PC,PC)          SNEX(N)             SOMEGA(I)           *
*          SZ(6,N)             T(I)                TEMP(N)             *
*          TEXCOR(I)           TL(I)               TR(R)               *
*          V(PC)               WEIGHT(X)                               *
************************************************************************
*                            ******************                        *
*                            * PRESENT VALUES *                        *
*                            *      R =   24  *                        *
*                            *      K =   10  *                        *
*                            *      F =   10  *                        *
*                            *      N =   30  *                        *
*                            *      M =    2  *                        *
*                            *      P =   47  *                        *
*                            *      C =    4  *                        *
*                            *      PC=   51  *                        *
*                            *      I =  999  *                        *
*                            *      T =    9  *                        *
*                            *      B =  100  *                        *
*                            *      L =   10  *                        *
*                            *      X =15000  *                        *
*                            ******************                        *
*               *****************************************              *
*               *   UNIT  3 = Raw measurement data      *              *
*               *   UNIT  5 = Standard input            *              *
*               *   UNIT  6 = Standard output           *              *
*               *   UNIT  7 = Disk file for reflection  *              *
*               *             output                    *              *
*               *   UNIT  8 = Disk file for profile     *              *
*               *             output                    *              *
*               *   UNIT 13 = Output file for observed  *              *
*               *             and calculated profile    *              *
*               *   UNIT 18 = Output file for observed  *              *
*               *             and calculated structure  *              *
*               *             factors.                  *              *
*               *****************************************              *
*                                                                      *
************************************************************************
*             PARAMETER NAME CONVENTIONS IN THIS PROGRAM               *
*                                                                      *
*      1.          Instrument dependent parameters                     *
*      --          -------------------------------                     *
*                                                                      *
*               PINSTR(1)  = Asymmetry parameter                       *
*               PINSTR(2)  = Overall scale factor                      *
*               PINSTR(3)  = Gaussian halfwidth parameter T            *
*               PINSTR(4)  = Gaussian halfwidth parameter U            *
*               PINSTR(5)  = Gaussian halfwidth parameter V            *
*               PINSTR(6)  = Gaussian halfwidth parameter W            *
*               PINSTR(7)  = Lorentzian halfwidth parameter X          *
*               PINSTR(8)  = Lorentzian halfwidth parameter Y          *
*               PINSTR(9)  = Zeropoint parameter Z                     *
*               PINSTR(10) = Zeropoint parameter S                     *
*----------------------------------------------------------------------*
*                                                                      *
*      2.          Structure dependent parameters                      *
*      --          ------------------------------                      *
*                                                                      *
*               PCRYST(1)  = Cell constant A in the equation below     *
*                  .                    .                              *
*                  .                    .                              *
*                  .                    .                              *
*               PCRYST(6)  = Cell constant F in the equation below     *
*               PCRYST(7)  = Overall temperature factor                *
*               PCRYST(8)  = Texture parameter #1                      *
*               PCRYST(9)  = Texture parameter #2                      *
*                  .                    .                              *
*                  .                    .                              *
*                  .                    .                              *
*               PCRYST(22) = Texture parameter #15                     *
*                                                                      *
* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -*
*        1/d**2 = A*h**2 + B*k**2 + C*l**2 + D*k*l + E*h*l + F*h*k     *
*                                                                      *
*----------------------------------------------------------------------*
*                                                                      *
*      3.                Atom parameters                               *
*      --                ---------------                               *
*                                                                      *
*               PATOM(1)   = Fractional atomic coordinate x            *
*               PATOM(2)   = Fractional atomic coordinate y            *
*               PATOM(3)   = Fractional atomic coordinate z            *
*               PATOM(4)   = Isotropic atomic temperature factor B     *
*               PATOM(5)   = Occupation number                         *
*               PATOM(6)   = Magnetic vector component K(x)            *
*               PATOM(7)   = Magnetic vector component K(y)            *
*               PATOM(8)   = Magnetic vector component K(z)            *
*               PATOM(9)   = Anisotropic temperature factor B(11)      *
*               PATOM(10)  = Anisotropic temperature factor B(22)      *
*               PATOM(11)  = Anisotropic temperature factor B(33)      *
*               PATOM(12)  = Anisotropic temperature factor B(12)      *
*               PATOM(13)  = Anisotropic temperature factor B(13)      *
*               PATOM(14)  = Anisotropic temperature factor B(23)      *
*                                                                      *
************************************************************************
*                                                                      *
      INCLUDE 'COMMONS.FOR'
*
*     Read the data and check for proper initialization.
*
      CALL INPUT
      ICYCLE=0
      FORCED=.FALSE.
    1 ICYCLE=ICYCLE+1
      IF(ICYCLE.GT.NCYCLE) GO TO 2
*
*     Take a new iteration step.
*
      CALL ITER
      GO TO 1
*
*     After the final cycle the results are written on disk files
*     if certain options are set.
*
    2 CALL EXPUT
      END
*     ******************
      FUNCTION APPROX(Y)
*     ******************
*
*     This function gives the value of expression
*
*     (1 + A*Y + B*Y**2)/(1 + C*Y + D*Y**2) .
*
      DATA A,B,C,D/0.903964515,0.769954808,1.36421616,1.13619508/
*
      APPROX=(1 + (A + B*Y)*Y)/(1 + (C + D*Y)*Y)
      RETURN
      END
      INCLUDE 'CALCUL.FOR'
*     *********************
      FUNCTION CODING(X,LP)
*     *********************
      LP=INT(ABS(X))/10
      CODING=(ABS(X)-10*LP)*SIGN(1.0,X)
      RETURN
      END
      INCLUDE 'CONPRE.FOR'
      INCLUDE 'COORD.FOR'
      INCLUDE 'CUBHAR.FOR'
*     *************************
      COMPLEX FUNCTION CWERF(Z)
*     *************************
*
*     Gautschi's algorithm for evaluating complex error function
*     for complex argument Z.
*
      COMPLEX VALUE,Z
      REAL LAMBDA
      LOGICAL B
      SAVE
*
*     TWOPI = 2/Sqrt(Pi)
*
      DATA TWOPI/1.12837 91671/
*
      XX=REAL(Z)
      YY=AIMAG(Z)
      X=ABS(XX)
      Y=ABS(YY)
      IF(Y .GE. 4.29 .OR. X .GE. 5.33) THEN
         H=0
         NC=0
         NU=8
         LAMBDA=0
         B=.TRUE.
      ELSE
         S=(1 - Y/4.29)*SQRT(1 - X*X/28.41)
         H=1.6*S
         H2=2*H
         NC=6+INT(23*S)
         NU=9+INT(21*S)
         LAMBDA=H2**NC
         B= LAMBDA .EQ. 0
      END IF
      R1=0
      R2=0
      S1=0
      S2=0
      N=NU+1
    3 N=N-1
      FN=N+1
      T1=Y+H+FN*R1
      T2=X-FN*R2
      C=0.5/(T1*T1+T2*T2)
      R1=C*T1
      R2=C*T2
      IF(H .LE. 0 .OR. N .GT. NC) GO TO 4
      T1=LAMBDA+S1
      S1=R1*T1-R2*S2
      S2=R2*T1+R1*S2
      LAMBDA=LAMBDA/H2
    4 IF(N .GT. 0) GO TO 3
      IF(B) THEN
         RS1=R1
         RS2=R2
      ELSE
         RS1=S1
         RS2=S2
      END IF
      RS1=TWOPI*RS1
      IF(Y .EQ. 0) RS1=EXP(-X*X)
      VALUE=CMPLX(RS1,TWOPI*RS2)
      IF(YY .GE. 0) THEN
         IF(XX .LT. 0) VALUE=CONJG(VALUE)
      ELSE
         VALUE=2*CEXP(-CMPLX(X,Y)**2)-VALUE
         IF(XX .GT. 0) VALUE=CONJG(VALUE)
      END IF
      CWERF=VALUE
      RETURN
      END
      INCLUDE 'DIRECT.FOR'
      INCLUDE 'DTSLIN.FOR'
*     *******************
      FUNCTION ERROR(A,B)
*     *******************
*
*     This function calculates the ESD of a function of 6 variables of
*     which the covariance matrix A and their derivatives B are given.
*
      REAL A(6,6),B(6)
      SUM=0
      DO 1 I=1,6
      DO 1 J=I,6
      X=2
      IF(I.EQ.J) X=1
    1 SUM=SUM+A(I,J)*B(I)*B(J)*X
      IF(SUM.LT.0) SUM=0
      ERROR=SQRT(SUM)
      RETURN
      END
*     ********************
      SUBROUTINE ESD(SM,V)
*     ********************
*
*     This subroutine calculates the E.S.D of the direct cell dimensions
*     from the inverse normal matrix SM and from the cell constants V.
*
      REAL SM(6,6),V(6),S(6),R(6),E(6)
      SAVE
      DATA DEGR/57.295779513/
*
      I=0
      DO 1 IA=1,3
      I=I+1
      IB=MOD(IA,3)+1
      IC=MOD(IA+1,3)+1
      ID=IA+3
      IE=IB+3
      IP=IC+3
      FNUM=4*V(IB)*V(IC)-V(ID)*V(ID)
      DEN=4*V(IA)*V(IB)*V(IC)-V(IA)*V(ID)*V(ID)-V(IB)*V(IE)*V(IE)-V(IC)
     1   *V(IP)*V(IP)+V(ID)*V(IE)*V(IP)
      S(I)=SQRT(FNUM/DEN)
      R(IA)=-S(I)**4
      DENL=1/(DEN*DEN)
      R(IB)=DENL*(4*V(IC)*DEN-(4*V(IA)*V(IC)-V(IE)*V(IE))*FNUM)
      R(IC)=DENL*(4*V(IB)*DEN-(4*V(IA)*V(IB)-V(IP)*V(IP))*FNUM)
      R(ID)=-DENL*(2*V(ID)*DEN-(2*V(IA)*V(ID)-V(IE)*V(IP))*FNUM)
      R(IE)=DENL*FNUM*(2*V(IB)*V(IE)-V(ID)*V(IP))
      R(IP)=DENL*FNUM*(2*V(IC)*V(IP)-V(ID)*V(IE))
      E(I)=ERROR(SM,R)
      E(I)=E(I)/(2*S(I))
      FNUM=V(IE)*V(IP)-2*V(IA)*V(ID)
      DEN=16*V(IA)*V(IA)*V(IB)*V(IC)+V(IE)*V(IE)*V(IP)*V(IP)-
     1    4*(V(IA)*V(IB)*V(IE)*V(IE)+V(IC)*V(IA)*V(IP)*V(IP))
      DENL=1/(DEN*DEN)
      FNUML=FNUM*FNUM
      S(I+3)=FNUML/DEN
      R(IA)=-4*FNUM*DENL*(V(ID)*DEN+FNUM*(8*V(IA)*V(IB)*V(IC)-
     1      V(IB)*V(IE)*V(IE)-V(IC)*V(IP)*V(IP)))
      R(IB)=-4*FNUML*DENL*(4*V(IA)*V(IA)*V(IC)-V(IA)*V(IE)*V(IE))
      R(IC)=-4*FNUML*DENL*(4*V(IA)*V(IA)*V(IB)-V(IA)*V(IP)*V(IP))
      R(ID)=-4*FNUM/DEN*V(IA)
      R(IE)=2*FNUM*DENL*(V(IP)*DEN-FNUM*(V(IE)*V(IP)*V(IP)-
     1      4*V(IA)*V(IB)*V(IE)))
      R(IP)=2*FNUM*DENL*(V(IE)*DEN-FNUM*(V(IE)*V(IE)*V(IP)-
     1      4*V(IC)*V(IA)*V(IP)))
      E(I+3)=ERROR(SM,R)
      IF(S(I+3).NE.0) THEN
         E(I+3)=E(I+3)/(2*SQRT(S(I+3)*(1-S(I+3))))*DEGR
         S(I+3)=ATAN(SQRT((1-S(I+3))/S(I+3)))*DEGR
      ELSE
         S(I+3)=90
         E(I+3)=0
      END IF
      IF(FNUM.LT.0) S(I+3)=180-S(I+3)
      IF(SM(ID,ID).EQ.0) E(ID)=0
    1 CONTINUE
      IF(ABS(S(1)-S(2)).LT.0.00008) E(1)=E(2)
      DO 2 I=1,6
      SM(I,I)=E(I)
    2 V(I)=S(I)
      RETURN
      END
*     ******************
      FUNCTION EXPERF(X)
*     ******************
*
*     This routine evaluates function EXP(X*X)*ERFC(X) with relative
*     error less than 1.3E-9 for negative X and less than 3.3E-11 for
*     positive x.
*
      T=1.0-7.5/(ABS(X)+3.75)
      Y=T*(+1.0260431255E-1+T*(-5.4802326329E-2+
     C  T*(+2.4143218599E-2+T*(-8.2206241220E-3+
     D  T*(+1.8029906156E-3+T*(-2.5522236047E-5+
     E  T*(-1.5254648703E-4+T*(+4.7864583725E-5+
     F  T*(+3.7322348606E-6+T*(-6.1287297159E-6+
     G  T*(+5.9065541351E-7+T*(+6.0603869357E-7+
     H  T*(-8.6338434635E-8+T*(-4.2166157960E-8))))))))))))))
      Y=   +1.4558972128E-1+T*(-2.7342193150E-1+
     A  T*(+2.2600806690E-1+T*(-1.6357189555E-1+Y)))
      EXPERF=Y
      IF(X.LT.0) EXPERF=2*EXP(X*X)-Y
      RETURN
      END
      INCLUDE 'EXPUT.FOR'
      INCLUDE 'GEOMET.FOR'
      INCLUDE 'GROUP.FOR'
      INCLUDE 'INPUT.FOR'
      INCLUDE 'ITER.FOR'
      INCLUDE 'MATRIX.FOR'
      INCLUDE 'PICK.FOR'
      INCLUDE 'PREPAR.FOR'
*     *****************
      SUBROUTINE PROFIL
*     *****************
*
*     This subroutine assigns the adjusted values to the profile
*     parameters and calculates direct cell dimensions from the cell
*     constants.
*
      INCLUDE 'BLANK.FOR'
      COMMON/PARAMS/AINSTR(10),ACRYST(22), ATOM(14,30),
     *              LINSTR(10),LCRYST(22),LATOM(14,30),
     *              PINSTR(10),PCRYST(22),PATOM(14,30)
      COMMON/SPACE/JUMP,A,B,C,COSALF,COSBET,COSGAM,SINALF,SINBET,SINGAM
      SAVE
*
      AINV=SQRT(PCRYST(1))
      BINV=SQRT(PCRYST(2))
      CINV=SQRT(PCRYST(3))
      AF=PCRYST(4)/(2*BINV*CINV)
      BF=PCRYST(5)/(2*AINV*CINV)
      CF=PCRYST(6)/(2*AINV*BINV)
      T1=PCRYST(1)*PCRYST(2)*PCRYST(3)
      T2=PCRYST(4)*PCRYST(5)*PCRYST(6)
      T3=PCRYST(1)*PCRYST(4)**2+
     *   PCRYST(2)*PCRYST(5)**2+
     *   PCRYST(3)*PCRYST(6)**2
      VL=SQRT(T1+(T2-T3)/4)
      T1=SQRT(1-AF**2)
      T2=SQRT(1-BF**2)
      T3=SQRT(1-CF**2)
*
*     Direct cell constants.
*
      A=BINV*CINV*T1/VL
      B=AINV*CINV*T2/VL
      C=AINV*BINV*T3/VL
      COSALF=(BF*CF-AF)/(T2*T3)
      COSBET=(AF*CF-BF)/(T1*T3)
      COSGAM=(AF*BF-CF)/(T1*T2)
      SINALF=SQRT(1-COSALF**2)
      SINBET=SQRT(1-COSBET**2)
      SINGAM=SQRT(1-COSGAM**2)
      RETURN
      END
*     ****************************************************************
      SUBROUTINE READER(UNIT,CHARS,KCHAR,REALS,KREAL,INTS,KINT,OPTION)
*     ****************************************************************
*
*     This subroutine works like a free format read statement with
*     following restrictions and extensions:
*
*     Restrictions:
*     -------------
*     A) Blanks may mean zero numbers but the rest of that line
*        must consist of zeros, blanks, commas or tabs.
*     B) If a line is blank or the end of it is blanks, all numbers
*        ment to be on that line or on part of it will become zeros.
*        Search is not continued on next line(s)!
*
*     Extensions:
*     -----------
*     A) Character texts may not be included between apostrophes
*     B) Real numbers may be given without a decimal point
*     C) Real numbers given in case where integers are needed,
*        are ROUNDED to nearest integer (e.g. -5.9 is -6)
*
      CHARACTER*(*) CHARS
      INTEGER KCHAR,KREAL,KINT,INTS(1),UNIT
      REAL REALS(1)
      LOGICAL OPTION
*
      CHARACTER*80 LINE,STORE,WORK,BLANK*1,COMMA*1,TAB*1,FMT*22
      INTEGER I,K,L,LENGTH,M
      REAL RESULT(20)
      SAVE
      DATA BLANK,COMMA,TAB/' ',',','	'/
*
      FMT='(F??.0)'
*
*                        A new line is read if so desired.
*
      IF(OPTION) THEN
         READ(UNIT,'(A)') LINE
         LENGTH=80
*
*                        Replace tabs with blanks
*
   10    I=INDEX(LINE,TAB)
         IF(I .GT. 0) THEN
            LINE(I:I)=BLANK
            GO TO 10
         END IF
*
*                        Replace commas with blanks
*
   20    I=INDEX(LINE,COMMA)
         IF(I .GT. 0) THEN
            LINE(I:I)=BLANK
            GO TO 20
         END IF
      STORE=LINE
      END IF
*
*                        Remove leading blanks
*
      I=1
   30 IF(LINE(I:I) .EQ. BLANK) THEN
         I=I+1
         IF(I .LT. 80) GO TO 30
      END IF
      WORK=LINE(I:LENGTH)
      LINE=WORK
      LENGTH=LENGTH-I+1
*
*                        Extract character part if existing
*
      IF(KCHAR.GT.0) THEN
         READ(LINE,'(A)') CHARS
         L=LEN(CHARS)
         WORK=LINE(L+1:80)
         LINE=WORK
         LENGTH=80-L
      END IF
*
*                        Search for all numbers wanted
*
      DO 50 M=1,KREAL+KINT
      RESULT(M)=0
      DO 40 I=1,LENGTH
      IF(LINE(I:I) .EQ. BLANK) GO TO 40
*
*                        A non-blank character was found.
*                        Search for next blank after it.
*
      WORK=LINE(I:LENGTH)
      LENGTH=LENGTH-I+1
      K=INDEX(WORK,BLANK)
      IF(K .LT. 2) GO TO 80
*
*                        Next blank was found on position K.
*                        Read the number before it.
*
      WRITE(FMT(3:4),'(I2)') K-1
      READ(LINE(I:I+K-2),FMT,ERR=90) RESULT(M)
      LINE=WORK(K:LENGTH)
      LENGTH=LENGTH-K+1
      GO TO 50
   40 CONTINUE
   50 CONTINUE
*
*                        All numbers are found. Copy them
*                        to proper arrays.
      DO 60 I=1,KREAL
   60 REALS(I)=RESULT(I)
      DO 70 I=1,KINT
   70 INTS(I)=NINT(RESULT(I+KREAL))
      RETURN
*
*                        Error: a line without blanks was found!
*
   80 WRITE(6,1000)
      STOP
*
*                        Error: Write the erroneous line out
*                               and mark the position where the
*                               error was encountered.
*
   90 WRITE(6,2000) ('1234567890',K=1,8),STORE
      FMT='(??X,''^'',/,1X,??(''=''))'
      WRITE(FMT(2:3),'(I2)') 96-LENGTH
      WRITE(FMT(15:16),'(I2)') 96-LENGTH
      WRITE(6,FMT)
      STOP
 1000 FORMAT(' No blanks on the line!')
 2000 FORMAT(16X,8A10,/,' Error in line: ',A)
      END
      INCLUDE 'SORT.FOR'
*     **************************
      FUNCTION STEP(RELAX,X,K,Z)
*     **************************
      INCLUDE 'BLANK.FOR'
      IF(K.NE.0) THEN
         STEP=V(K)*X*RELAX
         Z=SQRT(SMM(K,K))*ABS(X)
      ELSE
         STEP=0
         Z=0
      END IF
   10 RETURN
      END
      INCLUDE 'TEXTUR.FOR'
      INCLUDE 'TIMES.FOR'
      INCLUDE 'UPDATE.FOR'
      INCLUDE 'Y.FOR'
