CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C                  XABS2 - (c) Sean Parkin, 92.12.18                   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  XABS2 is a program to calculate empirical absorption corrections.   C
C  The original XABS (H. Hope and B. Moezzi, Jul83, Feb84, May85)      C
C  obtained 12 coefficients from a linear least-squares fit between    C
C  sin^2(theta) and (1/A vs sin^2(theta)), by minimizing weighted      C
C  differences of Fo^2 and Fc^2.                                       C
C  XABS2 calculates either 24, 18, 12 or 6 coefficients from a least-  C
C  squares fit of (1/A vs sin^2(theta)) to either a cubic, quadratic   C
C  linear or 'flat' equation in sin^2(theta) by minimization of Fo^2   C
C  and Fc^2 differences.  The program reads a file in format           C
C  (3I4,F12.2,F12.2,F10.2) called filename.fcf containing              C
C  [h, k, l, Fc2, Fo2, sigma(Fo2)] created by LIST 4 in SHELXL-97,     C
C  and outputs a file in format (3I4,2F8.2) named *x.hkl to be read    C
C  by the HKLF 4 instruction in SHELXL-97.  It also writes a new       C
C  instructions file for SHELXL-97 named *x.ins                        C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This program can be used freely provided it is cited properly.      C
C                                                                      C
C  Parkin, S., Moezzi, B. & Hope, H. (1995) J. Appl. Cryst. 28, 53-56  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   How to compile this program:                                       C
C                                                                      C
C   Compilation for 64bit systems may need some additional compiler    C
C   flags.  On Mac OSX with gfortran, the following work well:         C
C                                                                      C
C   gfortran -O3 -o xabs2 xabs2.f                                      C
C   gfortran -m64 -O3 -o xabs2 xabs2.f                                 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   How to run this program:                                           C
C                                                                      C
C   You should have SHELXL-97 *.res and *.fcf files, with the latter   C
C   having been created with the 'LIST 4' command.  Also, in order to  C
C   sample as much of reciprocal space as possible, it is best to use  C
C   'MERG 0' for the refinement run prior to application of XABS2.     C
C                                                                      C
C   The program can be started from a command line in three different  C
C   ways depending on how much control you want, as follows:           C
C                                                                      C
C   xabs2                                                              C
C                                                                      C
C   or                                                                 C
C                                                                      C
C   xabs2 filename                                                     C
C                                                                      C
C   or                                                                 C
C                                                                      C
C   xabs2 filename X                                                   C
C                                                                      C
C   In the first case, it asks for a filename (without an extension).  C
C   Once it has the filename, it asks if you want to control some of   C
C   the input parameters.  If you answer 'N' for no it will do a       C
C   linear (12 parameters) or a cubic (24 parameters) fit, depending   C
C   on whether your data were collected with MoK(alpha) or CuK(alpha)  C
C   x-rays.  No other input is required.  If you choose 'Y' for yes,   C
C   it will give a choice of cubic, quadratic (18 parameters), linear  C
C   or so-called 'flat' (6 parameters) fits.  The 'flat' fit is maybe  C
C   not so useful for small-molecule work but it was easy to include   C
C   in the code, and may be useful for low-resolution data as is       C
C   common for biological macromolecules.                              C 
C                                                                      C
C   The other input parameters are sin(theta) min. and max. cutoffs,   C
C   and an OMIT parameter.  Since the *.fcf file written by SHELXL-97  C
C   has already had the EXTI extinction correction applied (i.e. if    C
C   it was used in SHELXL-97), it should be safe to include the full   C
C   resolution range of data (values 0.0 and 1.0 for sin(theta) min.   C
C   and max.  Be careful to ensure that reflections obscured by the    C
C   beamstop have been flagged for removal during the SHELXL-97 run.   C
C   The default for OMIT in XABS2 is 3.0, which ensures that the very  C
C   weakest data (those with Fo^2 < 3sigma(Fo^2) do not get included   C
C   when calculating the correction factors.  In practice it does not  C
C   seem to make a whole lot of difference.                            C
C                                                                      C
C   If you include a filename (i.e. second case above) on the command  C
C   line then it will use that filename and then continue on to ask    C
C   whether you want to control the input parameters, similar to the   C
C   first case above.                                                  C
C                                                                      C
C   The third case allows you to run either cubic, quadratic, linear   C
C   or flat fits fully automatically without further input.  In this   C
C   case, simply replace 'X' by either C,Q,L,F (or c,q,l,f) for cubic  C
C   quadratic, linear or flat.  Any other characters will be ignored   C
C   and will cause the program to ask for more input parameters, as    C
C   per the second case above.                                         C
C                                                                      C
C   The program writes a corrected dataset with 'x.hkl' appended to    C
C   the stem of the filename.  This corrected dataset has an extra     C
C   column of numbers containing the correction factor applied to      C
C   each reflection.  It also writes a new SHELXL-97 INS file with     C
C   'x.ins' appended to the filename stem.  This new file is almost    C
C   the same as the previous *.res file, except that any 'MERG 0'      C
C   statement is changed to 'REM MERG 0', and any EXTI command has     C
C   also been REM'd out.                                               C
C                                                                      C
C   If you find that there is no improvement after XABS2, that just    C
C   means that it could not model any slowly varying change in the     C
C   scale factor over reciprocal space.  This could happen if your     C
C   primary absorption correction (e.g. SADABS, SORTAV etc.) was ok    C
C   or that absorption was not a problem to begin with.  Generally     C
C   speaking, the best indication of a beneficial effect after XABS2   C
C   is a flatter difference map and a reduction in the pancake effect  C
C   that can plague ellipsoids for absorption affected data.  This     C
C   may also be coupled with slightly lower R-values.  If you see no   C
C   obvious improvement, there is no point in using the 'corrected'    C
C   dataset.                                                           C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                                                        
      PROGRAM XABS2
      IMPLICIT NONE                                                     
                                                                        
      CHARACTER*80 LINE, FCFFIL, RESFIL, NEWFIL, ABSFIL, INSFIL, 
     &   OUTPUT, KIND
      CHARACTER*1 CTRL, FIT                                             
      INTEGER H, HSUM, X1, Y1, Z1, X2, Y2, COUNT, REJECT, O1, O2, O3    
      REAL L, A, ALP, BET, GAM, FC2, FO2, SFO2, STMIN, STMAX, MAXFO2,   
     &   OMIT, EXT, EXTMAX, F03, MEANU, FDF03, FDF032, SCALE, OSF
      DOUBLE PRECISION AS, ST2, H2AS2,  
     &   W2, L2, WLFOFC, WLFCFC, MATRIX, VECTOR, NORMAL, RHSVEC, PI,
     &   ABSCOR, CFO2, CSFO2, ABCOFA, ABCOPT, ABCMIN, ABCMAX, SCFOFC,   
     &   SFOFC, SCFO, SFO, DIAGNL, SPDV1, SPDV2, COSALS, COSBES, COSGAS 
      DIMENSION H(3), A(3), AS(3), H2AS2(6), ABCOPT(6), RHSVEC(24), 
     &   VECTOR(24), ABSCOR(24), DIAGNL(24), NORMAL(24,24), 
     &   MATRIX(24,24)  
      PI = 3.14159265354

100   FORMAT (4X, A76)                                                   
200   FORMAT (3I4, 2F14.2, F13.2)                            
300   FORMAT (A80)                                                      
400   FORMAT (A1)                                                       
500   FORMAT (3I4, 2F8.2, F10.3)                                                
600   FORMAT (3I4, F8.0, F8.2, F10.3)                                            
700   FORMAT (A5, 4F10.4, A34)                                                   
710   FORMAT (A5, 3F10.4, A44)                                                   
720   FORMAT (A5, 2F10.4, A54)                                                   
730   FORMAT (A5, F10.4, A64)                                                   
800   FORMAT (F8.3, A5, F8.3, A15)                                           
850   FORMAT (F8.3, A28)                                           
                                                                        
C  Print some identifying information on the screen.

      CALL BANNER                                         

C  Initialize arrays with zeros.

      CALL INTLZE (H, A, AS, H2AS2, ABCOPT, RHSVEC, VECTOR, ABSCOR,
     &   DIAGNL, NORMAL, MATRIX)

C  Get the names of files to be used.

      CALL GETNAM (RESFIL, FCFFIL, NEWFIL, INSFIL, ABSFIL, KIND)
             
C  Get information needed for the run.

      CALL GETNFO (RESFIL, L, A, ALP, BET, GAM, L2, F03, MEANU, 
     &   O1, O2, O3, STMIN, STMAX, OMIT, EXT, EXTMAX, FDF032, OSF,
     &   FCFFIL, KIND) 
          
C  Get the reciprocal cell parameters, volume etc.

      CALL RSCELL (A, ALP, BET, GAM, AS, COSALS, COSBES, COSGAS)
            
C  Build the XABS2 normal equations.

      CALL BLDNML (FCFFIL, O1, O2, O3, STMIN, STMAX, OMIT, EXT, EXTMAX, 
     &   MEANU, FDF032, F03, AS, L2, COSALS, COSBES, COSGAS, NORMAL, 
     &   RHSVEC)
     
C   Solve the XABS2 normal equations.                                                             
                                                                        
      CALL CHLSKY(NORMAL, DIAGNL, RHSVEC, ABSCOR, O1)   
                                                                        
C   The 6, 12, 18 or 24 coefficients returned by CHLSKY are then used.               
                                                                        
      WRITE (*,*) '  Applying correction ...'            
      
      OPEN (12, FILE = FCFFIL, STATUS = 'OLD')                          
      OPEN (13, FILE = NEWFIL, STATUS = 'UNKNOWN')                      
                                                                        
2010  READ (12, 300) LINE
      IF (LINE.EQ.' ') GOTO 2010 
      READ (LINE, 200, ERR = 2010) H(1), H(2), H(3), FC2, FO2, SFO2

      HSUM = ABS(H(1)) + ABS(H(2)) + ABS(H(3))                          
      ABCMIN = 1.0                                                      
      ABCMAX = 1.0                                                      
c      SCALE = 1.0 
c Scale is now that from the SHELXL FVAR line, but it needs to be squared
c because the SHELXL scale factor is relative to F, not F^2.
      scale = 1.0/(osf*osf)

c 2010/06/08 Next line commented out so that the data does not get scaled 
c by an arbitrary amount of 100.  It now gets scaled using the SHELXL scale
c from the FVAR line of the RES file (actually the square on that value
c because SHELXL scales relative to F rather than F^2).
c      IF (MAXFO2.GT.999999) SCALE = 100.0                               
                                                                        
C   This loop applies corrections to each reflection in turn and        
C   creates a corrected data set in file *.new.  If the corrected Fo^2  
C   numbers are too large for the format required by SHELXTL-plus or    
C   SHELXL-93, they (and sigma(Fo^2)) are scaled by a factor of 100.    
                                                                        
4000  IF (HSUM.NE.0) THEN                                               
                                                                        
         H2AS2(1) = H(1) * H(1) * AS(1) * AS(1) / 4.0                   
         H2AS2(2) = H(2) * H(2) * AS(2) * AS(2) / 4.0                   
         H2AS2(3) = H(3) * H(3) * AS(3) * AS(3) / 4.0                   
         H2AS2(4) = H(1) * H(2) * AS(1) * AS(2) / 2.0                   
         H2AS2(5) = H(1) * H(3) * AS(1) * AS(3) / 2.0                   
         H2AS2(6) = H(2) * H(3) * AS(2) * AS(3) / 2.0                   
                                                                        
         ST2 = L2 * (H2AS2(1) + H2AS2(2) + H2AS2(3) + (H2AS2(4)         
     &      * COSGAS) + (H2AS2(5) * COSBES) + (H2AS2(6) * COSALS))      
                                                                        
C   These DO loops calculate the parts of the correction factor that    
C   correspond to the six components of the ellipsoidal approximation   
C   to the absorption surface.                                          
                                                                        
         IF (O1.EQ.24) THEN                                             
                                                                        
            DO 4010 X1 = 1, 6                                           
               ABCOPT(X1) = ((H2AS2(X1) * (ST2 ** 2) * ABSCOR(X1))      
     &            + (H2AS2(X1) * ST2 * ABSCOR(X1 + 6)) + (H2AS2(X1)     
     &            * ABSCOR(X1 + 12)) + (H2AS2(X1)                       
     &            * ABSCOR(X1 + 18) / ST2))                             
4010        CONTINUE                                                    
                                                                        
         ELSE IF (O1.EQ.18) THEN                                        

            DO 4020 X1 = 1, 6                                           
               ABCOPT(X1) = ((H2AS2(X1) * ST2 * ABSCOR(X1)) + (H2AS2(X1)
     &         * ABSCOR(X1 + 6)) + (H2AS2(X1) * ABSCOR(X1 + 12) / ST2)) 
4020        CONTINUE                                                    
                                                                        
         ELSE IF (O1.EQ.12) THEN                                                       
                                                                        
            DO 4030 X1 = 1, 6                                           
               ABCOPT(X1) = ((H2AS2(X1) * ABSCOR(X1)) +                 
     &         (H2AS2(X1) * ABSCOR(X1 + 6) / ST2))                      
4030        CONTINUE                                                    

         ELSE IF (O1.EQ.6) THEN
            DO 4040 X1 = 1, 6
               ABCOPT(X1) = H2AS2(X1) * ABSCOR(X1) / ST2
4040        CONTINUE
                                                                        
         ENDIF                                                          
                                                                                                                                                
c   The absorption correction factor is then the sum of these parts,    
c   modified by a scale factor if corrected F^2 get too big for the     
c   format of *.hkl files in SHELXL92/3/6/7.                                  
                                                                        
         ABCOFA = L2 * (ABCOPT(1) + ABCOPT(2) + ABCOPT(3)       
     &      + ABCOPT(4) + ABCOPT(5) + ABCOPT(6)) 

c   Apply the correction to both Fo2 and sigma(Fo2).   
c   Next line commented out because SCALE is now applied to the 
c   corrected values on the subsequent two lines.
c         ABCOFA = ABCOFA * SCALE
         CFO2 = FO2 / (ABCOFA * SCALE)                                     
         CSFO2 = SFO2 / (ABCOFA * SCALE)                                     
         FO2 = FO2 / SCALE                                              
         FC2 = FC2 / SCALE                                              

c 2010/06/10 Next two lines changed because SCALE is now applied 
c elsewhere, as it should have been from the beginning.
c         ABCMIN = DMIN1(ABCMIN,(ABCOFA / SCALE))                                  
c         ABCMAX = DMAX1(ABCMAX,(ABCOFA / SCALE))                                  
         ABCMIN = DMIN1(ABCMIN,ABCOFA)                                  
         ABCMAX = DMAX1(ABCMAX,ABCOFA)                                  
         SCFOFC = SCFOFC + ABS(SQRT(ABS(CFO2)) - SQRT(ABS(FC2)))        
         SFOFC = SFOFC + ABS(SQRT(ABS(FO2)) - SQRT(ABS(FC2)))           
         SCFO = SCFO + SQRT(ABS(CFO2))                                  
         SFO = SFO + SQRT(ABS(FO2))                                     
                                                                        
C   If the value of the corrected Fo2 is too large for format F8.2      
C   then it is written in format F8.0.                                  
                                                                        
         IF (CFO2.LT.99999) THEN                                        
            WRITE(13, 500) H(1), H(2), H(3), CFO2, CSFO2, ABCOFA                
         ELSE                                                           
            WRITE(13, 600) H(1), H(2), H(3), CFO2, CSFO2, ABCOFA            
         ENDIF                                                          
                                                                        
         READ (12, 200) H(1), H(2), H(3), FC2, FO2, SFO2                
         HSUM = ABS(H(1)) + ABS(H(2)) + ABS(H(3))                       
      GOTO 4000                                                         
      ENDIF                                                             
                                                                        
      WRITE (13,500) 0, 0, 0, 0.00, 0.00                                
      CLOSE (12)                                                        
      CLOSE (13)                                                        
                                                                        
      OPEN (14, FILE = ABSFIL, STATUS = 'UNKNOWN')                      
      WRITE (14,'(1X,A78)') '+-----------------------------------------'
     &   //'-----------------------------------+'
      WRITE (14,'(A2,76X,A1)') '|', '|'
      WRITE (14,'(1X,A78)') '|               XABS2 - A four-dimensional'
     &   //' scale factor program              |'
      WRITE (14,'(1X,A78)') '|        (approximate treatment of anisotr'
     &   //'opic absorption/extinction)        |'
      WRITE (14,'(A2,76X,A1)') '|', '|'
      write (14,'(1X,A78)') '|    Parkin, S., Moezzi, B. & Hope, H. (19'
     &   //'95). J. Appl. Cryst. 28, 53-56.    |'
      WRITE (14,'(A2,76X,A1)') '|', '|'
      WRITE (14,'(1X,A78)') '+========================================='
     &   //'===================================+'
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                 
      WRITE (14,'(1X,A78)') '|  The coefficients of the XABS2 tensor ar'
     &   //'e:                                 |'                             
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                
                                                                        
      DO 5000 X1 = 1,6                                                  
         IF (O1.EQ.24) THEN                                             
            WRITE (14,700) ' |   ', ABSCOR(X1), ABSCOR(X1+6),      
     &      ABSCOR(X1+12), ABSCOR(X1+18), '|'                                           
         ELSE IF (O1.EQ.18) THEN                                        
            WRITE (14,710) ' |   ', ABSCOR(X1), ABSCOR(X1+6), 
     &      ABSCOR(X1+12), '|'   
         ELSE IF (O1.EQ.12) THEN
            WRITE (14,720) ' |   ', ABSCOR(X1), ABSCOR(X1+6), '|'
         ELSE
            WRITE (14,730) ' |   ', ABSCOR(X1), '|'
         ENDIF                                                          
5000  CONTINUE                                                          
                                                                        
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                 
      WRITE (14, '(1X,A$)') '|  Relative min/max transmission factors: '
      WRITE (14, 800) ABCMIN, ' and ', ABCMAX, '           |'
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                 
      WRITE (14, '(1X,A$)') '|  Ratio of min/max transmission factors: '
      WRITE (14, 850) ABCMIN/ABCMAX, '                  |'
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                 
      WRITE (14, '(1X,A$)') '|  R-values before and after correction:  '
      WRITE (14, 800) SFOFC/SFO, ' and ', SCFOFC/SCFO, '      |'                   
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                 
      WRITE (14, '(1X,A78)') '|  A corrected data file *x.hkl has been '
     &   //'written.                            |' 
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                 

c 2010/06/10 Next lines not needed because data are now automatically 
c scaled using the square of the SHELXL scale factor from FVAR line.
c      IF (SCALE.NE.1.) THEN                                             
c         WRITE (14, *) 'The corrected data were scaled down by 100'     
c         WRITE (14, *) ' '                                              
c      ELSE                                                              
c         WRITE (14, *) ' '                                              
c         WRITE (14, *) ' '                                              
c      ENDIF                                                             

      CALL NEWINS (RESFIL, INSFIL)
                                                                              
      WRITE (14, '(1X,A78)') '|  A new instructions file *x.ins has bee'
     &   //'n written.                          |'
      WRITE (14,'(A2,76X,A1)') '|', '|'                                                 
      WRITE (14,'(1X,A78)') '+========================================='
     &   //'===================================+'

      CLOSE (14)         
      
      OPEN (14, FILE = ABSFIL, STATUS = 'OLD')                          
                                                                        
      DO 5010 COUNT = 1, 28
         READ (14, 300) OUTPUT                                          
         IF (COUNT.LE.8) GOTO 5010
         WRITE (*, *) OUTPUT(3:77)                                            
5010  CONTINUE                                                          
                                                                        
      CLOSE (14)                                                        
                                                                        
      STOP                                                              
                                                                        
9990  WRITE (*,'(1X,A$)') 'THERE IS NO SUCH FILE AS '//RESFIL
      WRITE (*,*)CHAR(32)
      STOP                                                              
9991  WRITE (*,'(1X,A$)') 'THERE IS NO SUCH FILE AS '//FCFFIL         
      WRITE (*,*)CHAR(32)
      STOP                                                              
                                                                        
      END                                                               
C=======================================================================                                                                        


C-----------------------------------------------------------------------

      SUBROUTINE BANNER

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This simple subroutine prints out some identifying information to   C
C  the screen.                                                         C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      WRITE (*,*)
      WRITE (*,'(1X,A78)') '+------------------------------------------'
     &   //'----------------------------------+'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      WRITE (*,'(1X,A78)') '|               XABS2 - A four-dimensional '
     &   //'scale factor program              |'
      WRITE (*,'(1X,A78)') '|        (approximate treatment of anisotro'
     &   //'pic absorption/extinction)        |'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      write (*,'(1X,A78)') '|    Parkin, S., Moezzi, B. & Hope, H. (199'
     &   //'5). J. Appl. Cryst. 28, 53-56.    |'
      WRITE (*,'(A2,76X,A1)') '|', '|'
      WRITE (*,'(1X,A78)') '+=========================================='
     &   //'==================================+'
      WRITE (*,*)

      RETURN
      END

C=======================================================================                                                                        


C-----------------------------------------------------------------------

      SUBROUTINE INTLZE (H, A, AS, H2AS2, ABCOPT, RHSVEC, VECTOR, 
     &   ABSCOR, DIAGNL, NORMAL, MATRIX)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This simple subroutine initializes a bunch of variables and arrays  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      INTEGER X1, H, Y1
      REAL A
      DOUBLE PRECISION AS, H2AS2, ABCOPT, RHSVEC, VECTOR, ABSCOR,
     &   DIAGNL, NORMAL, MATRIX
      DIMENSION H(3), A(3), AS(3), H2AS2(6), ABCOPT(6), RHSVEC(24),
     &   VECTOR(24), ABSCOR(24), DIAGNL(24), NORMAL(24,24), 
     &   MATRIX(24,24)

      DO 1 X1 = 1, 3
         H(X1) = 0
         A(X1) = 0.0
         AS(X1) = 0.0
1     CONTINUE

      DO 2 X1 = 1, 6
         H2AS2(X1) = 0.0
         ABCOPT(X1) = 0.0
2     CONTINUE

      DO 3 X1 = 1, 24
         RHSVEC(X1) = 0.0
         VECTOR(X1) = 0.0
         ABSCOR(X1) = 0.0
         DIAGNL(X1) = 0.0
         DO 4 Y1 = 1, 24
            NORMAL(X1,Y1) = 0.0
            MATRIX(X1,Y1) = 0.0
4        CONTINUE
3     CONTINUE

      RETURN
      END

C=======================================================================                                                                        


C-----------------------------------------------------------------------

      SUBROUTINE GETNAM (RESFIL, FCFFIL, NEWFIL, INSFIL, ABSFIL, KIND)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine gets the filename either from the command line or   C
C  by asking the user if there is no filename available on the input   C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

C   File *.res contains the wavelength and cell dimensions (input).     
C   File *.fcf contains h, k, l, Fc^2, Fo^2, sigma(Fo^2) (input).       
C   File *x.hkl contains corrected data (output).                        
C   File *x.abs contains the coefficients, min/max correction            
C   factors, and R-values before and after correction.  
c
c   2011/12/23: Added ability to start program with a filename and 
c   thereby bypass filename request.
                                                                        
      CHARACTER*80 NAME, RESFIL, FCFFIL, NEWFIL, INSFIL, ABSFIL, KIND
      
 100  FORMAT (A80)

      CALL GETARG(1, NAME)
      CALL GETARG(2, KIND)
      
      if ((kind.eq.'c').or.(kind.eq.'C')) then
         kind = 'C'
         goto 1
      endif
      if ((kind.eq.'q').or.(kind.eq.'Q')) then
         kind = 'Q'
         goto 1
      endif
      if ((kind.eq.'l').or.(kind.eq.'L')) then
         kind = 'L'
         goto 1
      endif
      if ((kind.eq.'f').or.(kind.eq.'F')) then
         kind = 'F'
         goto 1
      endif
      
      kind = 'N' 
      
1     continue
 
      IF (NAME.NE.'') THEN
         RESFIL = NAME(1:INDEX(NAME,' ')-1)//'.res'                        
         FCFFIL = NAME(1:INDEX(NAME,' ')-1)//'.fcf'                        
         NEWFIL = NAME(1:INDEX(NAME,' ')-1)//'x.hkl'                        
         ABSFIL = NAME(1:INDEX(NAME,' ')-1)//'x.abs'                        
         INSFIL = NAME(1:INDEX(NAME,' ')-1)//'x.ins'                        
      ELSE
         WRITE (*,'(1X,A$)') '  Input filename without extension: '          
         READ (*, 100) NAME                                                    
         RESFIL = NAME(1:INDEX(NAME,' ')-1)//'.res'                        
         FCFFIL = NAME(1:INDEX(NAME,' ')-1)//'.fcf'                        
         NEWFIL = NAME(1:INDEX(NAME,' ')-1)//'x.hkl'                        
         ABSFIL = NAME(1:INDEX(NAME,' ')-1)//'x.abs'                        
         INSFIL = NAME(1:INDEX(NAME,' ')-1)//'x.ins'                        
      ENDIF

      RETURN
      END
                                                                        
C=======================================================================                                                                        

 
C-----------------------------------------------------------------------

      SUBROUTINE GETNFO (RESFIL, L, A, ALP, BET, GAM, L2, F03, MEANU, 
     &   O1, O2, O3, STMIN, STMAX, OMIT, EXT, EXTMAX, FDF032, OSF, 
     &   FCFFIL, KIND)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine gets information needed for the run.  Some of it    C
C  is extracted from the *.res and *.fcf files, or is input by user    C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 RESFIL, FCFFIL, LINE, WORD, KIND
      CHARACTER*1 CTRL, FIT
      
      INTEGER O1, O2, O3
      REAL L, A, ALP, BET, GAM, F03, MEANU, STMIN, STMAX, OMIT, EXT,
     &   EXTMAX, FDF03, FDF032, OSF
      DOUBLE PRECISION L2
      DIMENSION A(3)

100   FORMAT (A80)
200   FORMAT (4X, A76)
250   FORMAT (20X, A60)
300   FORMAT (A1)

C  The next lines skip the TITL line in a SHELXL *.res file, and      
C  then read in wavelength and cell parameters.  
C  S.Parkin 161109: Fix for SHELXL2016/6.                                                       
                                                                        
      OPEN (11, FILE = RESFIL, STATUS = 'OLD', ERR = 9990)              

1     READ (11, 100) LINE
      IF (LINE(1:4).NE.'CELL') goto 1

      line = line(5:80)                                           
      CLOSE (11)                                                        
      OPEN (11, STATUS = 'SCRATCH')                                     
      WRITE (11,200) LINE                                               
      REWIND (11)                                                       
      READ (11,*) L, A(1), A(2), A(3), ALP, BET, GAM     
      
      CLOSE (11)
      OPEN (11, FILE = RESFIL, STATUS = 'OLD', ERR = 9990)              

C Get the overall scale factor from the FVAR line in the RES file.

800   READ (11, 100) LINE
      IF (LINE(1:4).NE.'FVAR') GOTO 800
      READ (LINE, 200) WORD
      READ (WORD, *) OSF
      REWIND (11)
      CLOSE (11)

      L2 = L * L                                                        
                                                                        
C The reflections used to calculate the coefficients can be           
C controlled by altering the input parameters: 
C Sin(theta) min. and max. cutoffs allow low and high resolution 
C reflections to be skipped.  
C OMIT allows reflections with poor counting statistics to be 
C skipped.  
C EXT allows reflections at low angle with Fo^2 less than (EXT*Fc^2) 
C to be skipped to exclude those reflections suffering from severe          
C extinction.  The sin(theta) min. cutoff will do a similar job but      
C is less selective.  
C Reflections can also be skipped if F(obs) is greater than some 
C fraction of F(000).  
C Note: The SHELXL97 *.FCF file has Fo^2 values that have already 
C had an isotropic extinction correction applied.  This means that 
C including them here will force XABS2 to calculate an approximate 
C treatment of any anisotropy in the extinction correction.  This
C is especially true if a reasonable absorption correction (such 
C as SADABS, SORTAV, similar or better) has already been performed.
C This treatment for anisotropic extinction is an approximation,
C but it seems to do a passable job.
                                                                        
C 111223: Removed next two lines because it is less effort to just 
C automatically get the value of F(000) from the FCF file.
C      WRITE (*,'(1X,A$)') 'Enter a rough value for F(000): '            
C      READ (*,*) F03     
      
      OPEN (11, FILE = FCFFIL, STATUS = 'OLD', ERR = 9991)      
900   READ (11,100) LINE
      IF (LINE(1:20).NE.'_exptl_crystal_F_000') GOTO 900
      READ (LINE, 250) WORD
      READ (WORD, *) F03
      REWIND (11)
      CLOSE (11)
      
C 111223: Moved request for a mean isotropic U value to after the 
C request for control over the input parameters.  Default set at 0.01, 
C which should mean it does very little if anything.
C      WRITE (*,'(1X,A$)') 'Enter a rough overall thermal parameter: '   
C      READ (*,*) MEANU   

      MEANU = 0.005
            
      if (kind.ne.'N') then
         fit = kind(1:1)
         IF (FIT.EQ.'C'.OR.FIT.EQ.'c') THEN                         
            O1 = 24                                                     
            O2 = 4                                                      
            O3 = 3                                                      
         ELSE IF (FIT.EQ.'Q'.OR.FIT.EQ.'q') THEN                    
            O1 = 18                                                     
            O2 = 3                                                      
            O3 = 2                                                      
         ELSE IF (FIT.EQ.'L'.OR.FIT.EQ.'l') THEN                    
            O1 = 12                                                     
            O2 = 2                                                      
            O3 = 1
         ELSE IF (FIT.EQ.'F'.OR.FIT.EQ.'f') THEN
            O1 = 6
            O2 = 1
            O3 = 0
         ENDIF                                                          
         STMIN = 0.0                                                   
         STMAX = 1.0                                                    
         OMIT = 3.0                                                     
         EXT = 0.5                                                      
         EXTMAX = 0.35 * SQRT(L)                                        
         FDF032 = 0.25                                                  
         goto 2
      endif
                     
      WRITE (*,'(1X,A$)') '  Do you want to control the input'
      WRITE (*,'(1X,A$)') 'parameters ? [N]: '   
      
c      WRITE (*,'(1X,A$)') 'Only answer yes for good reason [N]: '       
      READ (*,300), CTRL                                                    
      write (*,*)
                                                                        
      IF (CTRL.EQ.'Y'.OR.CTRL.EQ.'y') THEN  
      
c 120116: After lots of tests it seems better to not bother with MEANU as
c as a user definable parameter at all.  Probably this is because LIST 4  
c format *.fcf file is already corrected for extinction.
c         WRITE (*,'(1X,A$)') 'Enter a rough overall thermal parameter: '   
c         READ (*,*) MEANU   
                                                                        
         WRITE (*,'(1X,A$)') '     Enter cubic(C), quadratic(Q) or '         
         WRITE (*,'(1X,A$)') 'linear(L) or flat(F) '                              
1000     READ (*,300) FIT

         IF (FIT.EQ.'C'.OR.FIT.EQ.'c') THEN                         
            O1 = 24                                                     
            O2 = 4                                                      
            O3 = 3                                                      
         ELSE IF (FIT.EQ.'Q'.OR.FIT.EQ.'q') THEN                    
            O1 = 18                                                     
            O2 = 3                                                      
            O3 = 2                                                      
         ELSE IF (FIT.EQ.'L'.OR.FIT.EQ.'l') THEN                    
            O1 = 12                                                     
            O2 = 2                                                      
            O3 = 1
         ELSE IF (FIT.EQ.'F'.OR.FIT.EQ.'f') THEN
            O1 = 6
            O2 = 1
            O3 = 0
         ELSE                                                           
            WRITE (*,'(1X,A$)') '  Not acceptable, again please: '        
            GOTO 1000                                                   
         ENDIF                                                          

         stmin = 0.0
         stmax = 1.0
         omit = 3.0
                                                                                 
         WRITE (*,'(1X,A$)') '     Enter sin(theta) min. cutoff [0.0]: '        
1010     READ (*,'(F8.0)', err=1011) STMIN 
         goto 1012
1011     stmin = 0.0
         write (*,*) '  Nonsense input, using 0.0 for sin(theta)min. '
1012     continue
         IF (STMIN.LT.0.0.OR.STMIN.GE.1.0) THEN                     
            WRITE (*,'(1X,A$)') '  Not acceptable, again please: '        
         GOTO 1010                                                      
         ENDIF                                                          

         WRITE (*,'(1X,A$)') '     Enter sin(theta) max. cutoff [1.0]: '        
1020     READ (*,'(F8.0)', err=1021) STMAX  
         goto 1022
1021     stmax = 1.0
         write (*,*) '  Nonsense input, using 1.0 for sin(theta)max. '
1022     continue
         if (stmax.eq.0.0) stmax = 1.0
         IF (STMAX.LE.STMIN.OR.STMAX.GT.1.0) THEN                   
            WRITE (*,'(1X,A$)') '  Not acceptable, again please: '        
         GOTO 1020                                                      
         ENDIF   
                                                                        
         WRITE (*,'(1X,A$)') '     Enter OMIT, (reject if Fo2 <'            
         WRITE (*,'(1X,A$)') 'OMIT(sigmaFo2) [3.0]: '                         
1030     READ (*,'(F8.0)', err=1031) OMIT  
         goto 1032
1031     omit = 3.0
         write (*,*) '  Nonsense input, using 3.0 for OMIT. '
1032     continue
         if (omit.eq.0.0) omit = 3.0
         IF (OMIT.LT.0.0) THEN                                          
            WRITE (*,'(1X,A$)') '  Not acceptable, again please: '        
         GOTO 1030                                                      
         ENDIF  
                                                                         
c 120116: Since *.fcf is already corrected for extinction, this parameter
c does very little if anything.  Best to just set it to some value that 
c will do no damage.
c          WRITE (*,'(1X,A$)') 'Enter EXT, (reject if Fo2 < EXT(Fc2)): '  
c 1040     READ (*,*) EXT                                                    
c          IF (EXT.LT.0.0.OR.EXT.GE.1.0) THEN                         
c             WRITE (*,'(1X,A$)') '  Not acceptable, again please: '        
c          GOTO 1040                                                      
c          ENDIF   

         EXT = 0.5
                                                                        
c 120116: Since *.fcf is already corrected for extinction, this parameter
c does very little if anything.  Best to just set it to some value that 
c will do no damage.
c          WRITE (*,'(1X,A$)') 'Enter sin(theta) max. for '               
c          WRITE (*,'(1X,A$)') 'EXT criterion: '                          
c 1050     READ (*,*) EXTMAX                                                 
c          IF (EXTMAX.LT.0.0.OR.EXTMAX.GT.STMAX) THEN                 
c             WRITE (*,'(1X,A$)') '  Not acceptable, again please: '        
c          GOTO 1050                                                      
c          ENDIF  

         EXTMAX = 0.35 * SQRT(L)
                                                                        
c 120116: Since *.fcf is already corrected for extinction, this parameter
c does very little if anything.  Best to just set it to some value that 
c will do no damage.
c          WRITE (*,'(1X,A$)') 'Enter maximum allowable Fobs/F(000): '    
c 1060     READ (*,*) FDF03                                                  
c          IF (FDF03.GE.1.0) THEN                                         
c             WRITE (*,'(1X,A$)') '  Not acceptable, again please: '        
c          GOTO 1060                                                      
c          ENDIF                                                          
c          FDF032 = FDF03 ** 2    

         FDF032 = 0.25
                                                                        
         WRITE (*,*) ' '                                                
         WRITE (*,*) '  End of input'                                     
         WRITE (*,*) ' '                                                
                                                                        
C   These parameters are set to reasonable defaults.  The values were   
C   arrived at by looking at a bunch of test cases, but they may not be        
C   optimal for all situations.                             
                                                                        
      ELSE                                                              
         STMIN = 0.0                                                   
         STMAX = 1.0                                                    
         OMIT = 3.0                                                     
         EXT = 0.5                                                      
         EXTMAX = 0.35 * SQRT(L)                                        
         FDF032 = 0.25                                                  
         IF (L.LE.1.0) THEN
            fit = 'l'
            O1 = 12                                                     
            O2 = 2                                                      
            O3 = 1                                                      
         ELSE IF (L.GT.1.0) THEN 
            fit = 'c'
            O1 = 24                                                     
            O2 = 4                                                      
            O3 = 3                                                      
         ENDIF                                                          
      ENDIF                                                             

2     continue

      if ((fit.eq.'c').or.(fit.eq.'C')) 
     &   kind = ' cubic (24 parameters)    '
      if ((fit.eq.'q').or.(fit.eq.'Q')) 
     &   kind = ' quadratic (18 parameters)'
      if ((fit.eq.'l').or.(fit.eq.'L')) 
     &   kind = ' linear (12 parameters)   '
      if ((fit.eq.'f').or.(fit.eq.'F')) 
     &   kind = ' flat (6 parameters)      '
      
      write (*,*) '  Running with the following parameters: '
      write (*,*)
      write (*,'(A23, A26)') '      Type of fit    =  ', kind
      write (*,'(A23, F5.2)') '      Sin(theta)min. = ', stmin
      write (*,'(A23, F5.2)') '      Sin(theta)max. = ', stmax
      write (*,'(A23, F5.2)') '      OMIT parameter = ', omit
      write (*,*)

      RETURN

9990  WRITE (*,'(1X,A$)') 'THERE IS NO SUCH FILE AS ', RESFIL           
      STOP                                                              
9991  WRITE (*,'(1X,A$)') 'THERE IS NO SUCH FILE AS ', FCFFIL           
      STOP                                                              

      END

C=======================================================================                                                                        


C-----------------------------------------------------------------------

      SUBROUTINE RSCELL (A, ALP, BET, GAM, AS, cosals, cosbes, cosgas)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine calculates the reciprocal cell parameters from the  C
C  real-space unit cell parameters                                     C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      REAL A, ALP, BET, GAM
      DOUBLE PRECISION AS, COSALP, COSBET, COSGAM, COSALS, COSBES,
     &   COSGAS, VOLUME, PI
      DIMENSION A(3), AS(3)

C   Fortran wants angles in radians in order to calculate cosines.      

      PI = 3.14159265354
      ALP = ALP * PI / 180.0                                            
      BET = BET * PI / 180.0                                            
      GAM = GAM * PI / 180.0                                            
      COSALP = COS(ALP)                                                 
      COSBET = COS(BET)                                                 
      COSGAM = COS(GAM)                                                 
      COSALS = ((COSBET * COSGAM) - COSALP) / (SIN(BET) * SIN(GAM))     
      COSBES = ((COSALP * COSGAM) - COSBET) / (SIN(ALP) * SIN(GAM))     
      COSGAS = ((COSALP * COSBET) - COSGAM) / (SIN(ALP) * SIN(BET))     
                                                                        
C   Calculate cell volume, a*, b* and c*.                               
                                                                        
      VOLUME = A(1) * A(2) * A(3) * SQRT(1.0 - (COSALP*COSALP)          
     &   - (COSBET * COSBET) - (COSGAM * COSGAM)                        
     &   + (2.0 * COSALP * COSBET * COSGAM))                            
      AS(1) = A(2) * A(3) * SIN(ALP) / VOLUME                           
      AS(2) = A(1) * A(3) * SIN(BET) / VOLUME                           
      AS(3) = A(1) * A(2) * SIN(GAM) / VOLUME  
            
      RETURN
      END

C=======================================================================                                                                        


C-----------------------------------------------------------------------

      SUBROUTINE BLDNML (FCFFIL, O1, O2, O3, STMIN, STMAX, OMIT, EXT, 
     &   EXTMAX, MEANU, FDF032, F03, AS, L2, COSALS, COSBES, COSGAS, 
     &   NORMAL, RHSVEC)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  This subroutine builds the normal equations                         C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      
      CHARACTER*80 FCFFIL, LINE
      INTEGER H, HSUM, REJECT, X1, Y1, Z1, X2, Y2, O1, O2, O3 
      REAL FC2, FO2, SFO2, MAXFO2, OMIT, STMIN, STMAX, MEANU, FDF032, 
     &   FO3, EXT, EXTMAX 
      DOUBLE PRECISION AS, H2AS2, ST2, L2, COSALS, COSBES, COSGAS, W2, 
     &   WLFOFC, WLFCFC, MATRIX, VECTOR, SPDV1, SPDV2, RHSVEC, NORMAL
      DIMENSION H(3), AS(3), H2AS2(6), RHSVEC(24), VECTOR(24), 
     &   NORMAL(24,24), MATRIX(24,24)
      
200   FORMAT (3I4, 2F14.2, F13.2)                             
300   FORMAT (A80)                                                      

      WRITE (*,*) '  Building the normal equations ...'                                          

C   Read in the first line of data from the *.fcf file.  HSUM is just   
C   a counter that looks for the end of the data file.  All the 'CIF'   
C   junk at the top of the file is skipped automatically.               
                                                                        
      OPEN (12, FILE = FCFFIL, STATUS = 'OLD', ERR = 9991)              

2000  READ (12, 300) LINE
      IF (LINE.EQ.' ') GOTO 2000 
      READ (LINE, 200, ERR = 2000) H(1), H(2), H(3), FC2, FO2, SFO2

      HSUM = ABS(H(1)) + ABS(H(2)) + ABS(H(3))   
            
      REJECT = 0                                                        
      MAXFO2 = FO2                                                      
                                                                        
C   This loop reads in each reflection and constructs the normal        
C   equations.  If the reflection fails to meet the input criteria,     
C   it is skipped.                                                      
                                                                        
3000  IF (HSUM.NE.0) THEN                                               
                                                                        
         H2AS2(1) = H(1) * H(1) * AS(1) * AS(1) / 4.0                   
         H2AS2(2) = H(2) * H(2) * AS(2) * AS(2) / 4.0                   
         H2AS2(3) = H(3) * H(3) * AS(3) * AS(3) / 4.0                   
         H2AS2(4) = H(1) * H(2) * AS(1) * AS(2) / 2.0                   
         H2AS2(5) = H(1) * H(3) * AS(1) * AS(3) / 2.0                   
         H2AS2(6) = H(2) * H(3) * AS(2) * AS(3) / 2.0                   
                                                                        
         ST2 = L2 * (H2AS2(1) + H2AS2(2) + H2AS2(3) + (H2AS2(4)         
     &      * COSGAS) + (H2AS2(5) * COSBES) + (H2AS2(6) * COSALS))      
                                                                        
C   The following series of IF's reject the current reflexion if it     
C   fails to meet the inclusion criteria.                               
                                                                        
         IF (FO2.LT.OMIT * SFO2) REJECT = 1                             
         IF (ST2.LT.STMIN ** 2) REJECT = 1                            
         IF (ST2.GT.STMAX ** 2) REJECT = 1                            
         IF (FO2 / EXP(-78.96 * 2 * MEANU * ST2 / L2)                 
     &      .GT.FDF032 * (F03 ** 2)) REJECT = 1                       
         IF (ST2.LE.EXTMAX ** 2.AND.FO2.LT.EXT * FC2) REJECT = 1
                                                                        
         IF (REJECT.NE.1) THEN                                          
                                                                        
C   W2 is the weight assigned to each reflexion when its contribution   
C   to the normal equations is being calculated.                        
                                                                        
            W2 = 1. / (SFO2 * SFO2)                                     
            WLFOFC = W2 * L2 * FO2 * FC2                                
            WLFCFC = W2 * L2 * FC2 * FC2                                
                                                                        
C   The next two groups of nested DO loops calculate the contribution   
C   to the normal equations (MATRIX) and right hand side (VECTOR) of    
C   each individual reflexion.  SPDV1 and SPDV2 are dummy variables 
C   that enable shifting stuff out of the inner loop.  This resulted  
C   in a roughly five fold speed up.
                                                                        
            DO 3010 X1 = 1, 6                                           
               SPDV1 = WLFCFC * H2AS2(X1)                               
               DO 3020 Z1 = 1, O2                                       
                  X2 = X1 + ((Z1 - 1) * 6)                              
                  SPDV2 = SPDV1 * ST2 ** (O3 - Z1)                      
                  DO 3030 Y1 = 1, O1                                    
                     MATRIX(Y1,X2) = SPDV2                              
3030              CONTINUE                                              
3020           CONTINUE                                                 
3010        CONTINUE                                                    
                                                                        
            DO 3040 Y1 = 1, 6                                           
               DO 3050 Z1 = 1, O2                                       
                  Y2 = Y1 + ((Z1 - 1) * 6)                              
                  SPDV2 = H2AS2(Y1) * ST2 ** (O3 - Z1)                  
                  VECTOR(Y2) = WLFOFC * SPDV2                           
                  SPDV1 = L2 * SPDV2                                    
                  DO 3060 X1 = 1, O1                                    
                     MATRIX(Y2,X1) = MATRIX(Y2,X1) * SPDV1              
3060              CONTINUE                                              
3050           CONTINUE                                                 
3040        CONTINUE                                                    
                                                                        
C   The next group of nested loops builds up the normal equations       
C   by adding MATRIX and VECTOR for each reflection to the complete     
C   normal matrix (NORMAL) and right hand side vector (RHSVEC).         
                                                                        
            DO 3070 Y1 = 1, O1                                          
               RHSVEC(Y1) = RHSVEC(Y1) + VECTOR(Y1)                     
               DO 3080 X1 = 1, O1                                       
                  NORMAL(Y1,X1) = NORMAL(Y1,X1) + MATRIX(Y1,X1)         
3080           CONTINUE                                                 
3070        CONTINUE                                                    
                                                                        
         ENDIF                                                          
                                                                        
C   Read in the next reflection, change the value of MAXFO2 if           
C   necessary, and then complete the loop.                              
                                                                        
         READ (12, 200) H(1), H(2), H(3), FC2, FO2, SFO2
         HSUM = ABS(H(1)) + ABS(H(2)) + ABS(H(3))                       
         REJECT = 0                                                     
         MAXFO2 = AMAX1(FO2,MAXFO2)                                     
                                                            
      GOTO 3000                                                         
      ENDIF                                                             
      CLOSE (12)   

      RETURN

      STOP

9991  WRITE (*,'(1X,A$)') 'THERE IS NO SUCH FILE AS ', FCFFIL           
      STOP                                                              
      
      END
      
C=======================================================================                                                                        


C-----------------------------------------------------------------------

      SUBROUTINE CHLSKY(NORMAL, DIAGNL, RHSVEC, ABSCOR, O1)             

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   CHLSKY solves the normal equations.  Adapted from a routine in     C  
C   Numerical Recipes in Fortran by Press, Teukolsky, Vetterling and   C
C   Flannery.                                                          C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                                                                                                                                
      INTEGER I, J, K, O1                                               
      DOUBLE PRECISION NORMAL, DIAGNL, RHSVEC, ABSCOR, SUM              
      DIMENSION NORMAL(24,24), DIAGNL(24), RHSVEC(24), ABSCOR(24)       
                                                                        
      WRITE (*,*) '  Solving normal equations ...'  

      DO 1003 I = 1, O1                                                 
         DO 1002 J = I, O1                                              
            SUM = NORMAL(I,J)                                           
            DO 1001 K = I-1, 1, -1                                      
               SUM = SUM - NORMAL(I,K) * NORMAL(J,K)                    
 1001       CONTINUE                                                    
            IF (I.EQ.J) THEN  
c               IF (SUM.LE.0.0) PAUSE 'CHOLESKY FAILED'
c               write (*,*) i,j,sum
               DIAGNL(I) = SQRT(SUM)                                    
            ELSE                                                        
               NORMAL(J,I) = SUM / DIAGNL(I)                            
            ENDIF                                                       
 1002    CONTINUE                                                       
 1003 CONTINUE                                                          
                                                                        
      DO 1005 I = 1, O1                                                 
         SUM = RHSVEC(I)                                                
         DO 1004 K = I-1, 1, -1                                         
            SUM = SUM - NORMAL(I,K) * ABSCOR(K)                         
 1004    CONTINUE                                                       
         ABSCOR(I) = SUM / DIAGNL(I)                                    
 1005 CONTINUE                                                          
                                                                        
      DO 1007 I = O1, 1, -1                                             
         SUM = ABSCOR(I)                                                
         DO 1006 K = I+1, O1                                            
            SUM = SUM - NORMAL(K,I) * ABSCOR(K)                         
 1006    CONTINUE                                                       
         ABSCOR(I) = SUM / DIAGNL(I)                                    
 1007 CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               

C=======================================================================   


C-----------------------------------------------------------------------

      SUBROUTINE NEWINS (RESFIL, INSFIL)             

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  Writes a new instructions file *x.ins from the old *.res file, but  C
C  with the 'MERG 0' instruction (if used) changed to 'REM MERG 0'     C
C  and any 'EXTI' line also REM'ed out.                                C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CHARACTER*80 RESFIL, INSFIL, LINE
      
100   FORMAT (A80)

      OPEN (11, FILE = RESFIL, STATUS = 'OLD', ERR = 9990)              
      OPEN (12, FILE = INSFIL, STATUS = 'UNKNOWN')              

1     READ (11, 100) LINE
            
      IF (LINE(1:6).EQ.'MERG 0') THEN
         LINE = 'REM '//LINE
      ENDIF
      
      if (line(1:4).eq.'EXTI') then
         line = 'REM '//line
      endif

      if (line(1:4).eq.'HKLF') then
         line = line(1:4)//' 4'
      endif
      
      WRITE (12, 100) LINE
      
      IF (LINE(1:4).EQ.'END ') THEN
         GOTO 2
      ENDIF
      
      GOTO 1
      
2     CLOSE (12)
      CLOSE (11)
      
      RETURN

9990  WRITE (*,'(1X,A$)') 'THERE IS NO SUCH FILE AS ', RESFIL           
      STOP                                                              

      RETURN                                                            
      END                                                               

C=======================================================================   


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C  History of changes:                                                 C
C                                                                      C
C  Extinction rejection based on F(000) ------------------ 1993.04.18  C
C  Subroutine CHLSKY replaced IMSL canned routine -------- 1993.07.14  C
C  Bug in EXT rejection criterion fixed ------------------ 1994.02.17  C
C  Loops rearranged/optimized (5-fold speed increase) ---- 1994.02.18  C
C  Minor bugs fixed -------------------------------------- 1994.08.18  C
C  Modified to read unedited *.fcf file ------------------ 1995.03.01  C
C  Modified to read CELL from *.res in free format ------- 1995.03.01  C
C  Modified to allow choice of approximation used -------- 1995.08.08  C
C  DO loops standardized --------------------------------- 1995.08.08  C
C  Small change to default values ------------------------ 1995.08.09  C
C  FORMAT of the *.fcf file changed to accept F12.2 ------ 1996.01.05  C
C  Bug associated with SCALE fixed ----------------------- 1996.01.05  C
C  Rejection of *.fcf file header bug fixed -------------- 1996.10.25  C
C  Changed to use the SHELXL overall SCALE from FVAR ----- 2010.06.10  C
C  Ability to start program with "xabs2 filename" -------- 2011.12.23  C
C  Get F(000) from *.fcf file ---------------------------- 2011.12.23  C
C  By default do not ask for mean isotropic U value ------ 2011.12.23  C
C  Normal equation build moved to subroutine BLDNML ------ 2011.12.27  C
C  Removed some not so useful user inputs ---------------- 2012.01.16  C
C  Added NEWINS subroutine to write new *x.ins file ------ 2012.01.16  C
C  Some output formatting changes ------------------------ 2012.01.16  C
C  Added ability to set kind of fit from command line ---- 2012.01.17  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
