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 in        C
C   sin^2(theta) to (1/A vs sin^2(theta)), by minimizing weighted      C
C   differences between Fo^2 and Fc^2.                                 C
C   XABS2 calculates either 24, 18 or 12 coefficients from a least-    C
C   squares fit of (1/A vs sin^2(theta)) to either a cubic, quadratic  C
C   or linear equation in sin^2(theta) by minimization of Fo^2 and     C
C   Fc^2 differences.  The program reads a file in format              C
C   (3I4,2X,F10.2,2X,F10.2,2X,F8.2) named *.fcf containing             C
C   [h, k, l, Fc2, Fo2, sigma(Fo2)] created by LIST 4 in SHELXL-93,    C
C   and outputs a file in format (3I4,2F8.2) or (3I4,F8.0,F8.2) named  C
C   *.new which can be renamed *.hkl, to be read into SHELXTL-plus or  C
C   SHELXL-93 by the HKLF 4 instruction.                               C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   History of changes:                                                C
C                                                                      C
C   'Extinction' rejection based on F(000) ----------------- 93.04.18  C
C   Subroutine CHLSKY replaced IMSL canned routine --------- 93.07.14  C
C   Bug in EXT rejection criterion fixed ------------------- 94.02.17  C
C   Loops rearranged/optimized (5-fold speed increase) ----- 94.02.18  C
C   Various bugs fixed ------------------------------------- 94.08.18  C
C   Modified to read an unedited *.fcf file ---------------- 95.03.01  C
C   Modified to read CELL info from *.res in free format --- 95.03.01  C
C   Modified to allow choice of approximation used --------- 95.08.08  C
C   DO loops converted to Fortran 77 ----------------------- 95.08.08  C
C   Change to default values made -------------------------- 95.08.09  C
C   FORMAT of the *.fcf file changed to accept F12.2 ------- 96.01.05  C
C   Recently introduced bug associated with SCALE fixed ---- 96.01.05  C
C   Rejection of *.fcf file header bug fixed --------------- 96.10.25  C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   This program can be used free of charge by anyone PROVIDED they    C
C   make adequate reference to it.  The most suitable reference is:    C
C   Parkin, S., Moezzi, B. & Hope, H. (1995) J. Appl. Cryst. 28, 53-56 C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                                                        
      PROGRAM XABS2
      IMPLICIT NONE                                                     
                                                                        
      CHARACTER*80 NAME, LINE, FCFFIL, RESFIL, NEWFIL, ABSFIL, OUTPUT   
      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, PI        
                                                                        
      DOUBLE PRECISION AS, VOLUME, ST2, H2AS2, COSALP, COSBET, COSGAM,  
     &   W2, L2, WLFOFC, WLFCFC, MATRIX, VECTOR, NORMAL, RHSVEC,        
     &   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), NORMAL(24,24), RHSVEC(24), 
     &   MATRIX(24,24), VECTOR(24), ABSCOR(24), ABCOPT(6), DIAGNL(24)   
                                                                        
      PI = 3.14159265354                                         

C  Initialize all arrays with zeros

      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

100   FORMAT (4X, A76)                                                   
200   FORMAT (3I4, 2F12.2, F10.2)                            
300   FORMAT (A80)                                                      
400   FORMAT (A1)                                                       
500   FORMAT (3I4, 2F8.2)                                                
600   FORMAT (3I4, F8.0, F8.2)                                            
700   FORMAT (4F10.4)                                                   
800   FORMAT (F6.3, A5, F6.3)                                           
                                                                        
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 *.new contains corrected data (output).                        
C   File *.abs contains the coefficients, min/max correction            
C   factors, and R-values before and after correction.                  
                                                                        
      WRITE (*,'(1X,A$)') 'Input filename without extension: '          
      READ (*, 300) NAME                                                    
                                                                        
      RESFIL = NAME(1:INDEX(NAME,' ')-1)//'.res'                        
      FCFFIL = NAME(1:INDEX(NAME,' ')-1)//'.fcf'                        
      NEWFIL = NAME(1:INDEX(NAME,' ')-1)//'.new'                        
      ABSFIL = NAME(1:INDEX(NAME,' ')-1)//'.abs'                        
                                                                        
C   The next READ line just skips the TITL line in a SHELXL *.res       
C   file, and the one after it reads in wavelength and unit cell        
C   parameters.                                                         
                                                                        
      OPEN (11, FILE = RESFIL, STATUS = 'OLD', ERR = 9990)              
      READ (11, 300)                                                    
      READ (11, 100) LINE                                               
      CLOSE (11)                                                        
      OPEN (11, STATUS = 'SCRATCH')                                     
      WRITE (11,100) LINE                                               
      REWIND (11)                                                       
      READ (11,*) L, A(1), A(2), A(3), ALP, BET, GAM                    
      CLOSE (11)

      L2 = L * L                                                        
                                                                        
C   The reflections used to calculate the coefficients can be           
C   controlled by altering the input parameters.  Sin(theta) min.       
C   and max. cutoffs allow low and high resolution reflections to       
C   be skipped.  The OMIT parameter allows reflections with poor        
C   counting statistics to be skipped.  The EXT option allows           
C   reflections at low angle with Fo^2 less than (EXT*Fc^2) to be       
C   skipped to exclude those reflections suffering from severe          
C   extinction.  The sin(theta) min. cutoff will do a similar job       
C   but is less selective.  Reflexions can also be skipped if           
C   F(obs) is greater than some fraction of F(000).                     
                                                                        
      WRITE (*,'(1X,A$)') 'Enter a rough value for F(000): '            
      READ (*,*) F03                                                       
      WRITE (*,'(1X,A$)') 'Enter a rough overall thermal parameter: '   
      READ (*,*) MEANU                                                     
      WRITE (*,*) 'Do you want to control the input parameters ?'       
      WRITE (*,'(1X,A$)') 'Only answer yes for good reason [N]: '       
      READ 400, CTRL                                                    
                                                                        
      IF ((CTRL.EQ.'Y').OR.(CTRL.EQ.'y')) THEN                          
                                                                        
         WRITE (*,'(1X,A$)') 'Enter cubic(C), quadratic(Q) or '         
         WRITE (*,'(1X,A$)') 'linear(L): '                              
1000     READ (*,400) 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                                                           
            WRITE (*,'(1X,A$)') 'Not acceptable, again please: '        
            GOTO 1000                                                   
         ENDIF                                                          
                                                                        
         WRITE (*,'(1X,A$)') 'Enter sin(theta) minimum cutoff: '        
1010     READ (*,*) STMIN                                                  
         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) maximum cutoff: '        
1020     READ (*,*) STMAX                                                  
         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): '                         
1030     READ (*,*) OMIT                                                   
         IF (OMIT.LT.0.0) THEN                                          
            WRITE (*,'(1X,A$)') 'Not acceptable, again please: '        
         GOTO 1030                                                      
         ENDIF                                                          
                                                                        
         WRITE (*,'(1X,A$)') 'Enter EXT, (reject if Fo2 < EXT(Fc2)): '  
1040     READ (*,*) EXT                                                    
         IF ((EXT.LT.0.0).OR.(EXT.GE.1.0)) THEN                         
            WRITE (*,'(1X,A$)') 'Not acceptable, again please: '        
         GOTO 1040                                                      
         ENDIF                                                          
                                                                        
         WRITE (*,'(1X,A$)') 'Enter sin(theta) max. for '               
         WRITE (*,'(1X,A$)') 'EXT criterion: '                          
1050     READ (*,*) EXTMAX                                                 
         IF ((EXTMAX.LT.0.0).OR.(EXTMAX.GT.STMAX)) THEN                 
            WRITE (*,'(1X,A$)') 'Not acceptable, again please: '        
         GOTO 1050                                                      
         ENDIF                                                          
                                                                        
         WRITE (*,'(1X,A$)') 'Enter maximum allowable Fobs/F(000): '    
1060     READ (*,*) FDF03                                                  
         IF (FDF03.GE.1.0) THEN                                         
            WRITE (*,'(1X,A$)') 'Not acceptable, again please: '        
         GOTO 1060                                                      
         ENDIF                                                          
         FDF032 = FDF03 ** 2                                            
                                                                        
         WRITE (*,*) ' '                                                
         WRITE (*,*) 'End of input'                                     
         WRITE (*,*) ' '                                                
                                                                        
C   These parameters are set to reasonable defaults.  The values were   
C   arrived at by analysis of many test cases, but They are not         
C   necessarily optimal for all situations.                             
                                                                        
      ELSE                                                              
         STMIN = 0.05                                                   
         STMAX = 1.0                                                    
         OMIT = 3.0                                                     
         EXT = 0.5                                                      
         EXTMAX = 0.35 * SQRT(L)                                        
         FDF032 = 0.25                                                  
         IF (L.LE.1.0) THEN                                             
            O1 = 12                                                     
            O2 = 2                                                      
            O3 = 1                                                      
         ELSE IF (L.GT.1.0) THEN                                        
            O1 = 24                                                     
            O2 = 4                                                      
            O3 = 3                                                      
         ENDIF                                                          
      ENDIF                                                             

      WRITE (*,*) 'Running...'                                          
                                                                        
C   Fortran wants angles in radians in order to calculate cosines.      
                                                                        
      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                           
                                                                        
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.                                          
                                                                        
            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 from each reflexion 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 reflexion, 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)                                                        
                                                                        
C   The normal equations are now sent to the subroutine CHLSKY to be    
C   solved.                                                             
                                                                        
      WRITE (*,*) 'Solving normal equations...'  

      CALL CHLSKY(NORMAL, DIAGNL, RHSVEC, ABSCOR, O1)                   
                                                                        
C   The 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                                                      
      SCALE = 1.0                                                       
      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                                                           
                                                                        
            DO 4030 X1 = 1, 6                                           
               ABCOPT(X1) = ((H2AS2(X1) * ABSCOR(X1)) +                 
     &         (H2AS2(X1) * ABSCOR(X1 + 6) / ST2))                      
4030        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 SHELXL93.                                  
                                                                        
         ABCOFA = L2 * (ABCOPT(1) + ABCOPT(2) + ABCOPT(3)       
     &      + ABCOPT(4) + ABCOPT(5) + ABCOPT(6))                        
         ABCMIN = DMIN1(ABCMIN,ABCOFA)                                  
         ABCMAX = DMAX1(ABCMAX,ABCOFA)                                  
                                                                        
C   Apply the correction to both Fo2 and sigma(Fo2).                    
                                                                        
         CFO2 = FO2 / (ABCOFA * SCALE)                                           
         CSFO2 = SFO2 / (ABCOFA * SCALE)                                       
         FO2 = FO2 / SCALE                                              
         FC2 = FC2 / SCALE                                              
         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                
         ELSE                                                           
            WRITE(13, 600) H(1), H(2), H(3), CFO2, CSFO2                
         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, *) ' '                                                 
      WRITE (14, *) 'The coefficients are:'                             
      WRITE (14, *) ' '                                                 
                                                                        
      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,700) ABSCOR(X1), ABSCOR(X1+6), ABSCOR(X1+12)      
         ELSE                                                           
            WRITE (14,700) ABSCOR(X1), ABSCOR(X1+6)                     
         ENDIF                                                          
5000  CONTINUE                                                          
                                                                        
      WRITE (14, *) ' '                                                 
      WRITE (14, '(1X,A$)') 'Relative min/max transmission factors were'
      WRITE (14, 800) ABCMIN, ' and ', ABCMAX
      print*, ABCMIN, ABCMAX                           
      WRITE (14, *) ' '                                                 
      WRITE (14, '(1X,A$)') 'R-values before and after correction are'  
      WRITE (14, 800) SFOFC/SFO, ' and ', SCFOFC/SCFO                   
      WRITE (14, *) ' '                                                 
      WRITE (14, *) 'The corrected data file has the suffix .new'       
      WRITE (14, *) ' '                                                 
      IF (SCALE.NE.1.) THEN                                             
         WRITE (14, *) 'The corrected data were scaled down by 100'     
         WRITE (14, *) ' '                                              
      ELSE                                                              
         WRITE (14, *) ' '                                              
         WRITE (14, *) ' '                                              
      ENDIF                                                             
      CLOSE (14)                                                        
                                                                        
      OPEN (14, FILE = ABSFIL, STATUS = 'OLD')                          
                                                                        
      DO 5010 COUNT = 1,18                                              
         READ (14, 300) OUTPUT                                          
         WRITE (*, *) OUTPUT                                            
5010  CONTINUE                                                          
                                                                        
      CLOSE (14)                                                        
                                                                        
      STOP                                                              
                                                                        
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                                                               
                                                                        
                                                                        
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                                                      C
C   Subroutine CHLSKY solves the normal equations by Cholesky          C
C   decomposition followed by backsubstitution.                        C
C                                                                      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
                                                                        
      SUBROUTINE CHLSKY(NORMAL, DIAGNL, RHSVEC, ABSCOR, O1)             
                                                                        
      INTEGER I, J, K, O1                                               
                                                                        
      DOUBLE PRECISION NORMAL, DIAGNL, RHSVEC, ABSCOR, SUM              
                                                                        
      DIMENSION NORMAL(24,24), DIAGNL(24), RHSVEC(24), ABSCOR(24)       
                                                                        
      DO 1000 I = 1, O1                                                 
         DO 1010 J = I, O1                                              
            SUM = NORMAL(I,J)                                           
            DO 1020 K = I-1, 1, -1                                      
               SUM = SUM - NORMAL(I,K) * NORMAL(J,K)                    
1020        CONTINUE                                                    
            IF (I.EQ.J) THEN                                            
               DIAGNL(I) = SQRT(SUM)                                    
            ELSE                                                        
               NORMAL(J,I) = SUM / DIAGNL(I)                            
            ENDIF                                                       
1010     CONTINUE                                                       
1000  CONTINUE                                                          
                                                                        
      DO 1030 I = 1, O1                                                 
         SUM = RHSVEC(I)                                                
         DO 1040 K = I-1, 1, -1                                         
            SUM = SUM - NORMAL(I,K) * ABSCOR(K)                         
1040     CONTINUE                                                       
         ABSCOR(I) = SUM / DIAGNL(I)                                    
1030  CONTINUE                                                          
                                                                        
      DO 1050 I = O1, 1, -1                                             
         SUM = ABSCOR(I)                                                
         DO 1060 K = I+1, O1                                            
            SUM = SUM - NORMAL(K,I) * ABSCOR(K)                         
1060     CONTINUE                                                       
         ABSCOR(I) = SUM / DIAGNL(I)                                    
1050  CONTINUE                                                          
                                                                        
      RETURN                                                            
      END                                                               
