C----PROGRAM LATCON                                                      XYLT
C----ADAPTED TO SYSTEM BY D. SCHWARZENBACH, UNI LAUSANNE, SWITZERLAND    XYLT
C    AUGUST 1975                                                         XYLT
C     PROGRAM XRY745                                                     XYHT
CUNIV SUBROUTINE XRY745                                                  XYHT
C                                                                        XYLT
C                                                                        XYHT
C---)A GENERAL PROGRAM FOR THE LS-REFINEMENT OF LATTICE CONSTANTS. THIS  XYLT
C    IS A ONE-PASS NON-ITERATIVE PROGRAM WHICH @ETERMINES THE LS-VALUES  XYLT
C    OF THE RECIPROCAL METRIC TENSOR AND ITS VARIANCE-COVARIANCE MATRIX. XYLT
C    CONSTRAINTS DUE TO SYMMETRY MAY BE INTRODUCED. INPUT IS H, K, L,    XYLT
C    2-THETA, SIGMA OF 2-THETA. FROM THIS, THE APPROPRIATE WEIGHT FOR    XYLT
C    1/D**2 (D = SPACING OF LATTICE PLANES) IS CALCULATED.               XYLT
C----WRITTEN BY D. SCHWARZENBACH ETH-ZURICH. THE MAIN PARTS OF THIS      XYLT
C    PROGRAM HAVE BEEN COPIED FROM A SIMILAR IBM-1620 PROGRAM BY         XYLT
C    D. SCHWARZENBACH (PITTSBURGH, PA. 1966).       CHRISTMAS 1971.      XYLT
C----COMMON PACKAGE FOR SYSTEM                                           XYLT
      COMMON /SYS/ IDCOMP(2),NTIN,NTOUT,NTAPEA,NTAPEB,NTAPEC,NTAPED,     XYLT
     *NTAPEE,NTAPEF,NTAPEG,NTAPEH,NTAPEI,NTAPEJ,NCDBUF,NTINM,NTOUTM,     XYLT
     *NTPCHM,IDPROG(2),ITOOB,NPAGE,LINMX,LINCT,NCRDS,ISHORT,             XYLT
     *TITLE(19),CARD(21),QTQ(100),ICYCLE,IQUIT,                          XYLT
     * JUNK,JUNK1,JUNK2,JUNK3,JUNK4,IOZY,IOZTM,ICDC(20)                  XYLT
      DIMENSION IQTQ(100),ICARD(21),NTIND(1)                             XYLT
C     MODEL 70 ADDITION
      INTEGER*2 IDENT(40),ENDIT
      CHARACTER*20 INPUTFILE,OUTFILE
C     ADDITION ENDS
      EQUIVALENCE (QTQ(1),IQTQ(1)),(NTAPEC,NTPCH),(NTCU,NTPCHM)          XYLT
      EQUIVALENCE (CARD(1),ICARD(1)),(NTIND(1),NTIN)                     XYLT
      COMMON/XYINFO/MAXCOR,MINCOR,MAXUTL,IDPREV(2)                       XYLT
C----COMMON PACKAGE FOR LATCON                                           XYLT
      COMMON/OMNIA/A(7,14),B(6,12),C(7),CC(3,6),D(9),                    XYLT
     *DERIV(6,6),ESD1(6),ESD2(9),G(5),IG(5),IS(5),I1,I2,V(7)             XYLT
C                                                                        XYLT
CD6   COMMON XYDATA(1)                                                   XYLT
C     MODEL 70 MOD  (1) TO (800) IN NEXT 2 LINES
      COMMON /XYDATA/ QX(800)                                              XYLT
      DIMENSION AA(800)                                                    XYLT
CD6   EQUIVALENCE (AA(1),XYDATA(5000))                                   XYLT
      EQUIVALENCE (AA(1),QX(1))                                          XYLT
C     LEVEL 2,     AA,QX                                                 XYLT
CUNIV                                                                    XYLT
C                                                                        XYLT
      DIMENSION IIS(45),IIG(45),JSYM(27),IFL1(2),IFL2(2),GII(45)         XYLT
C----DATA STATEMENTS (WASTEFUL, BUT HANDY)                               XYLT
      DATA IIS /0,0,0,0,0,5,6,0,0,0,4,6,0,0,0,4,5,0,0,0,4,5,6,0,0,       XYLT
     * 2,4,5,6,0,2,4,5,6,0,2,3,5,6,0,2,3,4,5,6 /                         XYLT
      DATA IIG / 25*0,1,0,0,0,0,1,1,0,0,0,1,1,4,4,0,1,1,0,0,0 /          XYLT
      DATA GII / 25*0.,1.,0.,0.,0.,0.,1.,.5,0.,0.,0.,1.,1.,1.,1.,0.,     XYLT
     * 1.,1.,0.,0.,0. /                                                  XYLT
      DATA JSYM /4HTRIC,4HLINI,4HC   ,4HMONO,4HCLIN,4HIC-C,4HMONO,       XYLT
     * 4HCLIN,4HIC-B,4HMONO,4HCLIN,4HIC-A,4HORTH,4HORHO,4HMBIC,4HTETR,   XYLT
     * 4HAGON,4HAL  ,4HHEXA,4HGONA,4HL   ,4HRHOM,4HBOHE,4HDRAL,          XYLT
     * 4HCUBI,4HC   ,4H     /                                            XYLT
      DATA IFL1,IFL2,IFMX / 2HTH,2HEN,4HETA ,4HD   ,2 /                  XYLT
      DATA CON1,CON2 / 0.8726647E-02,1.7453294E-02 /                     XYLT
C     MODEL 70 ADDITION NEXT LINE
      DATA ENDIT/2HEN/
C                                                                        XYLT
C                                                                        XYLT
C----THE PROGRAM MINIMIZES THE QUANTITY WEIGHT*(Q-OBS - Q-CALC)**2.      XYLT
C    Q-CALC IS NOT NECESSARILY THE -TRUE- Q-CALC, BUT IT INCLUDES A      XYLT
C    CORRECTION TERM  -  Q-CALC = SUM(H(I)*H(J)*TENSORELEMENT(I,J)) +    XYLT
C    DK*(4./(LAMBDA**2))*SIN(2-THETA)**N, N = 1 OR 2.                    XYLT
C   *N=1*  TO BE USED ONLY IF THE POSITION OF THE POWDER LINES HAS BEEN  XYLT
C    MEASURED ON ONE SIDE OF THE PRIMARY BEAM ONLY. THUS THERE IS AN     XYLT
C    UNCERTAINTY ABOUT THE LOCATION OF THE PRIMARY BEAM, I.E. 2-THETA-   XYLT
C    ZERO. DK EQUALS THETA-ZERO, I.E. (2-THETA)TRUE = (2-THETA)OBS -2*DK XYLT
C    SEVERAL LS-CYCLES MAY BE NECESSARY TO GET A PROPER 2-THETA-ZERO,    XYLT
C    BECAUSE THE EQUATIONS ARE NOT LINEAR FOR THIS PARAMETER.            XYLT
C    2-THETA-OBS WILL BE CORRECTED BEFORE EACH CYCLE, SO THAT THE FINAL  XYLT
C    2-THETA-OBS VALUES AS WELL AS THE CALCULATED QUANTITES ARE THOSE    XYLT
C    WITH THE BEST LOCATION OF THE PRIMARY BEAM.                         XYLT
C   *N=2*  TO BE USED ONLY IF THE POSITIONS OF CORRESPONDING POWDER      XYLT
C    LINES HAVE BEEN MEASURED, I.E. 2-THETA = 0.5*(2-THETA-MINUS +       XYLT
C    2-THETA-PLUS). ABSORPTION EFFECTS AND SPECIMEN ECCENTRICITY TEND TO XYLT
C    INTRODUCE AN ERROR INTO MEASURED D-VALUES WHICH VARIES WITH         XYLT
C    COS(THETA)**2, DELTA(D)=2*DK*D*COS(THETA)**2. THE SYSTEMATIC ERROR  XYLT
C    IN Q IS THEN DELTA(Q) = -(4./(LAMBDA**2))*DK*SIN(2-THETA)**2.       XYLT
C    SEE AZAROFF AND BUERGER, THE POWDER METHOD, MC GRAWW-HILL 1958      XYLT
C    PAGE 240. ONLY ONE CYCLE IS NECESSARY. THE COMPUTED Q AND 2-THETA   XYLT
C    VALUES ARE NOT THE -TRUE- ONES, BUT Q-CALC = Q-TRUE + DELTA(Q),     XYLT
C    AND 2-THETA-CALC IS DERIVED FROM Q-CALC.                            XYLT
C----THE PROGRAM COMPUTES THE RECIPROCAL AND DIRECT METRIC TENSORS, THE  XYLT
C    RECIPROCAL AND DIRECT LATTICE CONSTANTS, AND ALL THE CORRESPONDING  XYLT
C    VARIANCE-COVARIANCE MATRICES. IT MAY ALSO BE USED TO GET ALL METRIC XYLT
C    QUANTITIES WITH THEIR ERRORS (VAR.-COVAR.) FROM A GIVEN RECIPROCAL  XYLT
C    METRIC TENSOR AND ITS VAR.-COVAR. MATRIX.                           XYLT
C                                                                        XYLT
C*********************************************************************** XYLT
C*----INPUT CARDS                                                      * XYLT
C*----CALLING CARD                                                     * XYLT
C*    LATCON COMPID ICORR,ISYM,WLAM,SIGLAM,CORR                        * XYLT
C*    FORMAT(A2,A4,1X,A2,A4,2I5,3F10.5)                                * XYLT
C*    COLS. 14-18 (0)/(-1)/(N) FOR CORRECTION TERM TO Q NOT TO BE REFI-* XYLT
C*     NED)/(REFINE CORRECTION TERM PROPORTIONAL TO SIN(2-THETA)**2,   * XYLT
C*     ONE CYCLE)/(REFINE N CYCLES TO GET 2-THETA-ZERO. REFINEMENT WILL* XYLT
C*     STOP IF CHANGE IS SMALLER THAN 0.002 DEGREES)                   * XYLT
C*    COLS. 19-23 SYMMETRY INDICATOR                                   * XYLT
C*     1 = TRICLINIC                                                   * XYLT
C*     2 = MONOCLINIC FIRST SETTING (C)                                * XYLT
C*     3 = MONOCLINIC SECOND SETTING (B)                               * XYLT
C*     4 = MONOCLINIC THIRD SETTING (A)                                * XYLT
C*     5 = ORTHORHOMBIC                                                * XYLT
C*     6 = TETRAGONAL                                                  * XYLT
C*     7 = HEXAGONAL - TRIGONAL                                        * XYLT
C*     8 = RHOMBOHEDRAL                                                * XYLT
C*     9 = CUBIC                                                       * XYLT
C*    COLS. 24-33 WAVE LENGTH IN ANGSTROM                              * XYLT
C*    COLS. 34-43 E.S.D. OF WAVELENGTH. LEAVE NORMALLY BLANK, UNLESS   * XYLT
C*     YOU WANT TO CORRECT THE E.S.D.-S OF THE REFINED QUANTITIES TO   * XYLT
C*     BE CORRECTED FOR THE INCERTAINTY IN THE WAVELENGTH.             * XYLT
C*     DEVISED MAINLY FOR NEUTRON-PEOPLE                               * XYLT
C*    COLS. 44-53 INITIAL VALUE OF 2-THETA ZERO, I.E. 2-THETA(TRUE) =  * XYLT
C*     2-THETA(OBS) - CORR. CORR NEED NOT BE REFINED                   * XYLT
C*                                                                     * XYLT
C*----REFLECTION DATA CARDS                                            * XYLT
C*    THETA  H,K,L, 2-THETA, E.S.D. OF 2-THETA, WAVE LENGTH            * XYLT
C*    FORMAT(A2,A4,9X,3F4.0,3F10.6)                                    * XYLT
C*    COLS.  1-6 THETA                                                 * XYLT
C*    COLS. 16-19 INDEX H                                              * XYLT
C*    COLS. 20-23 INDEX K                                              * XYLT
C*    COLS. 24-27 INDEX L                                              * XYLT
C*    COLS. 28-37 MEASURED 2-THETA VALUE                               * XYLT
C*    COLS. 38-47 E.S.D. OF 2-THETA (DEFAULT = .1 DEGREES)             * XYLT
C*    COLS. 48-57 WAVELENGTH FOR THIS REFLECTION IF DIFFERENT FROM     * XYLT
C*     VALUE GIVEN ON THE CALLING CARD                                 * XYLT
C*                                                                     * XYLT
C*----END CARD                                                         * XYLT
C*                                                                     * XYLT
C*********************************************************************** XYLT
C*C----BINARY DATA FILES AND SCRATCH FILES ARE NOT USED                  XYLT
C                                                                        XYLT
C*********************************************************************** XYLT
C----IMPORTANT NOTICE .....                                              XYLT
C    THE CORRECTION TO Q IS PROPORTIONAL TO Q*COTAN(THETA) (ICORR = N)   XYLT
C    OR TO Q*COS(THETA)**2 (ICORR = -1). 2-THETA-ZERO IS THEREFORE TO BE XYLT
C    REFINED MAINLY WITH LOW ORDER REFLECTIONS, THE SIN(2-THETA)**2      XYLT
C    TERM MUST BE REFINED WITH AS MANY ORDERS AS POSSIBLE. THE CORRELA-  XYLT
C    TIONS TEND TO BE HIGH, BUT IF THE WRONG REFLECTIONS ARE USED,       XYLT
C    THEY ARE SIMPLY LUDICROUS...                                        XYLT
C*********************************************************************** XYLT
C                                                                        XYLT
C                                                                        XYLT
C                                                                        XYLT
C      MODEL 70 ADDITIONS TO SIMULATE MAIN?
      NCDBUF=1
      NTOUT=2
      WRITE(*,*) 'ENTER NAME OF INPUT FILE'
      READ(*,1001) INPUTFILE
      WRITE(*,*) 'ENTER NAME OF OUTPUT FILE'
      READ(*,1001) OUTFILE
1001  FORMAT(A)
      MAXCOR=800
      OPEN(UNIT=NCDBUF,FILE=INPUTFILE,STATUS='UNKNOWN')
      OPEN(UNIT=NTOUT,FILE=OUTFILE,STATUS='UNKNOWN')
      DO 99888 I=1,4
      READ(NCDBUF,99887) IDENT
99887 FORMAT(40A2)
      WRITE(NTOUT,99887) IDENT
99888 CONTINUE
C     END OF ADDITIONS
      MARK = 1                                                           XYLT
      INCREF = 8                                                         XYLT
      MAXREF = MAXCOR / INCREF                                           XYLT
C                                                                        XYLT
C----DECODE THE CALLING CARD                                             XYLT
      READ(NCDBUF,10)ICORR,ISYM,WLAM,SIGLAM,CORR                         XYLT
   10 FORMAT(13X,2I5,3F10.5)                                             XYLT
      ICYC = 1                                                           XYLT
C                                                                        XYLT
C----IS THE SYMMETRY OK .....                                            XYLT
      IF (ISYM*(10-ISYM)) 1103,1103,1111                                 XYLT
 1103 WRITE(NTOUT,16)                                                    XYLT
C     PRINT       16                                                     XYLT
   16 FORMAT(44H0INVALID SYMMETRY INDICATOR, JOB TERMINATED.)            XYLT
      GO TO 9999                                                         XYLT
C                                                                        XYLT
 1111 JUNK1 = 3*ISYM-2                                                   XYLT
      WRITE (NTOUT,18) JSYM(JUNK1),JSYM(JUNK1+1),JSYM(JUNK1+2),          XYLT
     * WLAM,SIGLAM                                                       XYLT
   18 FORMAT(29H0LATTICE CONSTANT REFINEMENT./                           XYLT
     * 17H THE SYMMETRY IS ,3A4/18H THE WAVELENGTH IS,F10.5,             XYLT
     * 27H ANGSTROM WITH AN E.S.D. OF,F10.5)                             XYLT
      IF(ICORR) 1121,1123,1125                                           XYLT
 1121 WRITE(NTOUT,15)                                                    XYLT
   15 FORMAT(98H0CORRECTION TO Q-VALUE PROPORTIONAL TO (4/(LAMBDA**2))*S XYLT
     *IN(2-THETA)**2 TO BE REFINED. (ICORR = -1))                        XYLT
      GO TO 1123                                                         XYLT
 1125 WRITE(NTOUT,17) CORR                                               XYLT
C     PRINT       17, CORR                                               XYLT
   17 FORMAT(52H02-THETA-ZERO TO BE REFINED WITH A STARTING VALUE OF,    XYLT
     * F7.3)                                                             XYLT
      ICYC = ICORR                                                       XYLT
 1123 CONTINUE                                                           XYLT
C                                                                        XYLT
C----ENTER SYMMETRY CONSTRAINTS                                          XYLT
      ISYM = 5*ISYM-5                                                    XYLT
      DO 1113 I = 1,5                                                    XYLT
      ITEMP1=ISYM+I
      IS(I) = IIS(ITEMP1)                                                XYLT
      IG(I) = IIG(ITEMP1)                                                XYLT
 1113 G(I) = GII(ITEMP1)                                                 XYLT
      WLAMSQ = 4./(WLAM*WLAM)                                            XYLT
      MZ = 0                                                             XYLT
C                                                                        XYLT
C----READ REFLECTIONS, GET PROPER WEIGHTS, BUILD NORMAL EQUATIONS        XYLT
C     MODEL 70 ADDITION
1201  READ(NCDBUF,99920) IDENT(1),FH,FK,FL,THET2,SIGTH,VLAM
99920 FORMAT(A2,13X,3F4.0,3F10.6)
      KCD=1
      IF(IDENT(1).EQ.ENDIT) KCD=2
      GO TO(1211,1251,9999),KCD                                          XYLT
C                                                                        XYLT
 1211 CONTINUE                                                           XYLT
   20 FORMAT(15X,3F4.0,3F10.6)                                           XYLT
C                                                                        XYLT
      IF (THET2-CORR) 1201,1201,1217                                     XYLT
 1217 MZ = MZ + 1                                                        XYLT
      IF (MZ-MAXREF-1) 1221,1213,1201                                    XYLT
 1213 WRITE (NTOUT,22) MAXREF                                            XYLT
   22 FORMAT(15H0ONLY THE FIRST,I4,55H REFLECTIONS ARE USED. YOU MAY, HO XYLT
     *WEVER, CHANGE MAXCORE )                                            XYLT
      GO TO 1201                                                         XYLT
C----STORE REFLECTION                                                    XYLT
 1221 MM=(MZ-1)*INCREF + MARK - 1                                        XYLT
      AA(MM+1)=FH                                                        XYLT
      AA(MM+2)=FK                                                        XYLT
      AA(MM+3)=FL                                                        XYLT
      AA(MM+4)=THET2                                                     XYLT
      IF(SIGTH) 1225,1225,1227                                           XYLT
 1225 SIGTH = 0.1                                                        XYLT
 1227 AA(MM+5)=SIGTH                                                     XYLT
      AA(MM+6)=WLAMSQ                                                    XYLT
      IF(VLAM) 1201,1201,1205                                            XYLT
 1205 AA(MM+6)=4.0/(VLAM*VLAM)                                           XYLT
      GO TO 1201                                                         XYLT
C                                                                        XYLT
C----COME HERE AFTER ALL INPUT                                           XYLT
C----COMPUTE Q-VALUES                                                    XYLT
 1251 WRITE(NTOUT,21) MZ                                                 XYLT
   21 FORMAT(1H0,I4,17H REFLECTIONS READ)                                XYLT
      MZ = MIN0(MZ,MAXREF)                                               XYLT
      MAXUTL = MZ * INCREF                                               XYLT
      MARK1 = MAXUTL + MARK - 1                                          XYLT
 1253 CONTINUE                                                           XYLT
      IF(ICYC) 1257,1257,1259                                            XYLT
 1259 DO 1115 I = 1,7                                                    XYLT
      DO 1117 J = 1,14                                                   XYLT
 1117 A(I,J) = 0.0                                                       XYLT
 1115 V(I) = 0.0                                                         XYLT
 1257 CONTINUE                                                           XYLT
      DO 1255 K=MARK,MARK1,INCREF                                        XYLT
      FH=AA(K)                                                           XYLT
      FK=AA(K+1)                                                         XYLT
      FL=AA(K+2)                                                         XYLT
      THET2=AA(K+3)-CORR                                                 XYLT
      SIGTH=AA(K+4)                                                      XYLT
      WLAMSQ=AA(K+5)                                                     XYLT
      Q = WLAMSQ*(SIN(CON1*THET2)**2)                                    XYLT
      SINT = SIN(CON2*THET2)                                             XYLT
C----SIGLAM -MUST- NOT BE INCLUDED AT THIS POINT                         XYLT
      SIGSQ = (WLAMSQ*SINT*CON1*SIGTH)**2                                XYLT
      IF(SIGSQ) 1223,1223,1229                                           XYLT
 1223 SIGSQ = 1.0E-04*Q*Q                                                XYLT
 1229 AA(K+6)=Q                                                          XYLT
      WEIGHT = 1./SIGSQ                                                  XYLT
      AA(K+7)=WEIGHT                                                     XYLT
      IF(ICYC) 1255,1255,1247                                            XYLT
 1247 CONTINUE                                                           XYLT
C----BUILD MATRIX                                                        XYLT
      C(1) = FH*FH                                                       XYLT
      C(2) = FK*FK                                                       XYLT
      C(3) = FL*FL                                                       XYLT
      C(4) = 2.*FH*FK                                                    XYLT
      C(5) = 2.*FH*FL                                                    XYLT
      C(6) = 2.*FK*FL                                                    XYLT
      C(7) = 0.                                                          XYLT
      IF(ICORR) 1241,1243,1245                                           XYLT
 1241 C(7) = WLAMSQ*SINT*SINT                                            XYLT
      GO TO 1243                                                         XYLT
 1245 C(7) = WLAMSQ*SINT                                                 XYLT
 1243 DO 1231 I = 1,7                                                    XYLT
      DO 1233 J = 1,7                                                    XYLT
 1233 A(I,J) = A(I,J) + WEIGHT*C(I)*C(J)                                 XYLT
 1231 V(I) = V(I) + WEIGHT*C(I)*Q                                        XYLT
 1255 CONTINUE                                                           XYLT
      IF(ICYC) 1711,1711,1301                                            XYLT
C                                                                        XYLT
C----INTRODUCE SYMMETRY CONSTRAINTS INTO MATRIX                          XYLT
 1301 DO 1303 I = 1,7                                                    XYLT
 1303 A(I,I+7) = 1.0                                                     XYLT
      NR = 7                                                             XYLT
      DO 1313 I = 1,5                                                    XYLT
      K = IS(I)                                                          XYLT
      J = IG(I)                                                          XYLT
      IF(K) 1314,1313,1314                                               XYLT
 1314 IF(J) 1315,1316,1315                                               XYLT
 1315 DO 1317 L = 1,7                                                    XYLT
 1317 A(L,J) = A(L,J) + G(I)*A(L,K)                                      XYLT
      DO 1318 L = 1,7                                                    XYLT
 1318 A(J,L) = A(J,L) + G(I)*A(K,L)                                      XYLT
      V(J) = V(J) + G(I)*V(K)                                            XYLT
 1316 DO 1319 L = 1,7                                                    XYLT
      A(L,K) = 0.0                                                       XYLT
 1319 A(K,L) = 0.0                                                       XYLT
      V(K) = 0.0                                                         XYLT
      A(K,K+7) = 0.0                                                     XYLT
      NR = NR -1                                                         XYLT
 1313 CONTINUE                                                           XYLT
      IF(ICORR) 1323,1325,1323                                           XYLT
 1325 A(7,7) = 1.0                                                       XYLT
      NR = NR - 1                                                        XYLT
 1323 CONTINUE                                                           XYLT
C----ARE THERE ENOUGH OBSERVATIONS....                                   XYLT
      IF(MZ-NR) 1321,1401,1401                                           XYLT
 1321 WRITE (NTOUT,23) MZ                                                XYLT
   23 FORMAT(49H0NOT ENOUGH OBSERVATIONS TO SOLVE THIS PROBLEM - ,       XYLT
     * I2,18H REFLECTIONS READ.)                                         XYLT
      GO TO 9999                                                         XYLT
C----INVERT A(I,J), TAKING INTO ACCOUNT ZERO ROWS AND COLUMS OF DEPENDEN XYLT
C    T VARIABLES.                                                        XYLT
 1401 CONTINUE                                                           XYLT
      DO 1420 N = 1,7                                                    XYLT
      DO 1426 K = 1,5                                                    XYLT
      IF(IS(K)-N) 1426,1420,1426                                         XYLT
 1426 CONTINUE                                                           XYLT
      IF(A(N,N)) 1429,1431,1429                                          XYLT
 1431 WRITE (NTOUT,24)                                                   XYLT
   24 FORMAT(49H0LS-MATRIX IS SINGULAR. CHECK YOUR 2-THETA CARDS.)       XYLT
      GO TO 9999                                                         XYLT
 1429 BB = 1./A(N,N)                                                     XYLT
      DO 1421 J = 1,14                                                   XYLT
 1421 A(N,J) = A(N,J)*BB                                                 XYLT
      DO 1422 I = 1,6                                                    XYLT
      M = N+I                                                            XYLT
      IF(M-7) 1424,1424,1423                                             XYLT
 1423 M = M - 7                                                          XYLT
 1424 DO 1427 K = 1,5                                                    XYLT
      IF(IS(K)-M) 1427,1422,1427                                         XYLT
 1427 CONTINUE                                                           XYLT
      BB = A(M,N)                                                        XYLT
      DO 1425 J = 1,14                                                   XYLT
 1425 A(M,J) = A(M,J) - A(N,J)*BB                                        XYLT
 1422 CONTINUE                                                           XYLT
 1420 CONTINUE                                                           XYLT
C                                                                        XYLT
C----INVERSE MATRIX NOW IN A(I,J),I=1,7,J=8,14                           XYLT
C                                                                        XYLT
C----CALCULATE NOW RECIPROCAL METRIC TENSOR                              XYLT
      DO 1530 I = 1,7                                                    XYLT
      C(I) = 0.                                                          XYLT
      DO 1530 J = 1,7                                                    XYLT
 1530 C(I) = C(I) + A(I,J+7)*V(J)                                        XYLT
C----COMPLETION OF SOLUTION VECTOR                                       XYLT
      DO 1531 I = 1,5                                                    XYLT
      K = IS(I)                                                          XYLT
      J = IG(I)                                                          XYLT
      IF(K) 1532,1531,1532                                               XYLT
 1532 IF(J) 1533,1531,1533                                               XYLT
 1533 C(K) = G(I)*C(J)                                                   XYLT
      DO 1534 L = 1,7                                                    XYLT
      A(K,L+7) = G(I)*A(J,L+7)                                           XYLT
 1534 A(L,K+7) = G(I)*A(L,J+7)                                           XYLT
      A(K,K+7) = A(J,J+7)*G(I)*G(I)                                      XYLT
 1531 CONTINUE                                                           XYLT
C----FIRST PART OF 1620-PROGRAM FINISHED. IN C WE HAVE THE RECIPROCAL    XYLT
C    METRIC TENSOR                                                       XYLT
      ESD = A(7,14)                                                      XYLT
      DO 1631 I = 1,6                                                    XYLT
      V(I) = 0.0                                                         XYLT
      AJUNK = ESD*A(I,I+7)                                               XYLT
      IF(AJUNK) 1631,1631,1633                                           XYLT
 1633 V(I) = A(7,I+7)/FSQRT(AJUNK)                                        XYLT
 1631 CONTINUE                                                           XYLT
C----IS THERE AN APPRECIABLE CHANGE IN 2-THETA-ZERO......                XYLT
C    THIS IS A SHORTENED ITERATIVE PROCEDURE TO GET THIS ONE PARAMETER.  XYLT
      IF(ICORR) 1711,1711,1621                                           XYLT
 1621 DK = C(7)/CON1                                                     XYLT
      CORR = CORR + DK                                                   XYLT
      WRITE (NTOUT,25) DK,CORR,(V(I),I=1,6)                              XYLT
   25 FORMAT(47H02-THETA-ZERO HAS BEEN CHANGED BY THE AMOUNT OF,  F7.3,  XYLT
     *3H TO, F7.3, 8H DEGREES/   71H THE CORRELATIONS WITH THE ELEMENTS  XYLT
     *OF THE RECIPROCAL METRIC TENSOR ARE, 6F7.3 )                       XYLT
      ICYC = ICYC - 1                                                    XYLT
      IF(ABS(DK) - 0.002) 1625,1625,1253                                 XYLT
 1625 ICYC = 0                                                           XYLT
      GO TO 1253                                                         XYLT
C----CALCULATE Q-CALC, 2-THETA-CALC AND WEIGHTED RESIDUALS               XYLT
 1711 CONTINUE                                                           XYLT
      R = 0.0                                                            XYLT
      DK = C(7)                                                          XYLT
      WRITE (NTOUT,26) CORR                                              XYLT
   26 FORMAT(/82H0MEASURED AND COMPUTED 2-THETA ANGLES AND Q-VALUES. Q I XYLT
     *S (2*SIN(THETA)/LAMBDA)**2./40H 2THETA-TRUE = 2THETA-OBS - CORR, C XYLT
     *ORR =,F7.3/   107H FOR ICORR = -1, 2THETA-CALC AND Q-CALC ARE NOT  XYLT
     *THE -TRUE- VALUES BUT DERIVED FROM Q-TRUE + CORRECTION TERM)       XYLT
      WRITE (NTOUT,32)                                                   XYLT
   32 FORMAT(1H ,3X,1HH,4X,1HK,4X,1HL,3X,10H2THETA-OBS,2X,11H2THETA-TRUE XYLT
     *,2X,11H2THETA-CALC,1X,12HDELTA-2THETA,1X,12HE.S.D.2THETA,          XYLT
     * 1X,12HDELTA/E.S.D.,8X,5HQ-OBS,7X,6HQ-CALC  )                      XYLT
C                                                                        XYLT
      DO 2001 I = 1,3                                                    XYLT
 2001 CC(I,I) = C(I)                                                     XYLT
      CC(1,2) = C(4)                                                     XYLT
      CC(2,1) = C(4)                                                     XYLT
      CC(1,3) = C(5)                                                     XYLT
      CC(3,1) = C(5)                                                     XYLT
      CC(2,3) = C(6)                                                     XYLT
      CC(3,2) = C(6)                                                     XYLT
      DO 2199 I=MARK,MARK1,INCREF                                        XYLT
      QCALC = 0.0                                                        XYLT
      IH=AA(I)                                                           XYLT
      IK=AA(I+1)                                                         XYLT
      IL=AA(I+2)                                                         XYLT
      WLAMSQ=AA(I+5)                                                     XYLT
      THET2=AA(I+3)-CORR                                                 XYLT
      DO 2101 J = 1,3                                                    XYLT
      DO 2101 K = 1,3                                                    XYLT
      ITEMP1=I+J-1
      ITEMP2=I+K-1
 2101 QCALC = QCALC + CC(J,K)*AA(ITEMP1)*AA(ITEMP2)                      XYLT
      IF(ICORR) 2103,2105,2105                                           XYLT
 2103 QCALC = QCALC + WLAMSQ*DK*(SIN(CON2*THET2)**2)                     XYLT
 2105 CONTINUE                                                           XYLT
      SINT = FSQRT(QCALC/WLAMSQ)                                          XYLT
      THETC = ASIN(SINT)/CON1                                            XYLT
      DELTA = THET2 - THETC                                              XYLT
      DESD=DELTA/AA(I+4)                                                 XYLT
      R=R+AA(I+7)*((AA(I+6)-QCALC)**2)                                   XYLT
C----WRITE THIS REFLECTION                                               XYLT
      GO TO 2107
C      IF(NUC005(1)) 2107,2107,2109                                       XYLT
C????????
C2109 WRITE (NTOUT,32)                                                   XYLT
 2107 CONTINUE                                                           XYLT
      WRITE(NTOUT,28)IH,IK,IL,AA(I+3),THET2,THETC,DELTA,AA(I+4),         XYLT
     * DESD, AA(I+6), QCALC                                              XYLT
   28 FORMAT(3I5,6(2X,F11.3),2F13.8)                                     XYLT
C                                                                        XYLT
 2199 CONTINUE                                                           XYLT
C----FINALLY GET THE VARIANCE-COVARIANCE MATRIX                          XYLT
      ZM = MZ - NR                                                       XYLT
      IF(ZM) 2221,2221,2223                                              XYLT
 2223 R = R/ZM                                                           XYLT
 2221 NR = ZM                                                            XYLT
      DO 2201 I = 1,7                                                    XYLT
      DO 2201 J = 1,7                                                    XYLT
 2201 A(I,J) = R*A(I,J+7)                                                XYLT
      R = FSQRT(R)                                                        XYLT
      WRITE (NTOUT,30) R,NR                                              XYLT
   30 FORMAT(13H0ERROR OF FIT,11X,F10.4/                                 XYLT
     * 29H NUMBER OF DEGREES OF FREEDOM, I5   )                          XYLT
      ESD = FSQRT(A(7,7))                                                 XYLT
      IF(ICORR) 2301,3001,2303                                           XYLT
 2303 ESD = ESD/CON1                                                     XYLT
      WRITE (NTOUT,27) ESD                                               XYLT
C     PRINT        27, ESD                                               XYLT
   27 FORMAT(30H0THE E.S.D. OF 2-THETA-ZERO IS,F10.5)                    XYLT
      GO TO 3001                                                         XYLT
 2301 WRITE (NTOUT,29) DK,ESD,(V(I),I=1,6)                               XYLT
C     PRINT        29, DK,ESD,(V(I),I=1,6)                               XYLT
   29 FORMAT(43H0THE CORRECTION COEFFICIENT (ICORR = -1) IS,E14.6,       XYLT
     * 12H WITH E.S.D.,F10.5/ 34H AND WITH CORRELATION COEFFICIENTS,     XYLT
     * 6F7.3)                                                            XYLT
C                                                                        XYLT
C----NOW WE HAVE RECIPROCAL METRIC TENSOR AND VAR-COVAR MATRIX.          XYLT
C----AN UNCERTAINTY IN THE WAVELENGTH -SIGLAM- INTRODUCES A RELATIVE     XYLT
C    ERROR IN Q, NAMELY 2*(SIGLAM/WLAM). THE SAME RELATIVE ERROR IS INT- XYLT
C    RODUCED IN THE COMPONENTS OF THE METRIC TENSOR, I.E. THE VARIANCES  XYLT
C    HAVE TO BE CORRECTED BY AN ADDITIVE TERM (2*(SIGLAM/WLAM)*COMPONENT XYLT
C    ))**2. THE COVARIANCES ARE CORRECTED SUCH THAT THE CORRELATIONS     XYLT
C    REMAIN UNCHANGED.                                                   XYLT
      SIGLQ = 0.0                                                        XYLT
 3001 IF(SIGLAM) 3101,3101,3015                                          XYLT
 3015 SIGLQ = (2.*SIGLAM)/WLAM                                           XYLT
      R = (SIGLAM/WLAM)*100.                                             XYLT
      WRITE (NTOUT,33) R                                                 XYLT
   33 FORMAT(20H0A RELATIVE ERROR OF,F6.3,70H PERCENT IN THE WAVELENGTH  XYLT
     *IS USED TO CORRECT THE STANDARD DEVIATIONS.)                       XYLT
      DO 3017 I = 1,6                                                    XYLT
      V(I) = 0.                                                          XYLT
      IF(A(I,I)) 3017,3017,3021                                          XYLT
 3021 V(I) = FSQRT(A(I,I))                                                XYLT
      A(I,I) = A(I,I) + (SIGLQ*C(I))**2                                  XYLT
      V(I) = FSQRT(A(I,I))/V(I)                                           XYLT
 3017 CONTINUE                                                           XYLT
      DO 3019 I = 1,5                                                    XYLT
      N = I + 1                                                          XYLT
      DO 3019 J = N,6                                                    XYLT
      A(I,J) = A(I,J)*V(I)*V(J)                                          XYLT
 3019 A(J,I) = A(I,J)                                                    XYLT
C                                                                        XYLT
C----CALCULATE RECIPROCAL VECTORS AND THEIR VAR-COVAR MATRIX             XYLT
 3101 IF(LINCT.LT.12) LINCT = 0                                          XYLT
      WRITE (NTOUT,36)                                                   XYLT
   36 FORMAT(1H ,120(1H-)/21H0RECIPROCAL ELEMENTS.)                      XYLT
      I1 = -2                                                            XYLT
      I2 = -1                                                            XYLT
C----CALL METRIC                                                         XYLT
      CALL XRY746                                                        XYLT
C                                                                        XYLT
C----INVERT RECIPROCAL METRIC TENSOR, GET DIRECT LATTICE                 XYLT
      DO 4011 N = 1,3                                                    XYLT
      DO 4016 M = 4,6                                                    XYLT
 4016 CC(N,M) = 0.                                                       XYLT
 4011 CC(N,N+3) = 1.                                                     XYLT
      DO 4015 N = 1,3                                                    XYLT
      IF(CC(N,N)) 4009,4007,4009                                         XYLT
 4007 WRITE (NTOUT,38)                                                   XYLT
C     PRINT        38,                                                   XYLT
   38 FORMAT(64H0RECIPROCAL METRIC TENSOR IS A SINGULAR MATRIX. CHECK YO XYLT
     *UR DATA.)                                                          XYLT
      GO TO 9999                                                         XYLT
 4009 BB = 1./CC(N,N)                                                    XYLT
      DO 4012 J = 1,6                                                    XYLT
 4012 CC(N,J) = CC(N,J)*BB                                               XYLT
      DO 4015 I = 1,2                                                    XYLT
      M = N + I                                                          XYLT
      IF(M-3) 4014,4014,4013                                             XYLT
 4013 M = M - 3                                                          XYLT
 4014 BB = CC(M,N)                                                       XYLT
      DO 4015 J = 1,6                                                    XYLT
 4015 CC(M,J) = CC(M,J) - CC(N,J)*BB                                     XYLT
C                                                                        XYLT
      DO 4021 I = 1,3                                                    XYLT
 4021 C(I) = CC(I,I+3)                                                   XYLT
      C(4) = CC(1,5)                                                     XYLT
      C(5) = CC(1,6)                                                     XYLT
      C(6) = CC(2,6)                                                     XYLT
C                                                                        XYLT
C----COMPUTE THE DERIVATIVES OF THE METRIC TENSOR ELEMENTS WRT RECIPRO-  XYLT
C    CAL TENSOR ELEMENTS.                                                XYLT
      DO 5020 I = 1,3                                                    XYLT
      DO 5021 K = 1,3                                                    XYLT
 5021 DERIV(I,K) = -CC(I,K+3)**2                                         XYLT
      DO 5020 J = 1,2                                                    XYLT
      M = J+1                                                            XYLT
      DO 5020 N = M,3                                                    XYLT
      K = J+N+1                                                          XYLT
      DERIV(I,K) = -CC(I,J+3)*CC(I,N+3)                                  XYLT
 5020 DERIV(K,I) = 2.0*DERIV(I,K)                                        XYLT
      DERIV(4,4) = -C(4)**2-C(1)*C(2)                                    XYLT
      DERIV(5,5) = -C(5)**2 - C(1)*C(3)                                  XYLT
      DERIV(6,6) = -C(6)**2 - C(2)*C(3)                                  XYLT
      DERIV(4,5) = -C(4)*C(5) - C(1)*C(6)                                XYLT
      DERIV(4,6) = -C(4)*C(6) - C(2)*C(5)                                XYLT
      DERIV(5,6) = -C(5)*C(6) - C(3)*C(4)                                XYLT
      DO 5026 I = 4,6                                                    XYLT
      DO 5026 K = I,6                                                    XYLT
 5026 DERIV(K,I) = DERIV(I,K)                                            XYLT
C----DERIV(I,K) = DERIVATIVE OF METRIC TENSOR ELEMENT (K) WRT RECIPROCAL XYLT
C    METRIC TENSOR ELEMENT (I), WHERE THE ORDER OF THE ELEMENTS IS 11,   XYLT
C    22, 33, 12, 13, 23.                                                 XYLT
C                                                                        XYLT
C----COMPUTE VAR-COVAR MATRIX OF DIRECT METRIC TENSOR                    XYLT
      DO 5130 I = 1,6                                                    XYLT
      DO 5130 K = I,6                                                    XYLT
      J = K+6                                                            XYLT
      A(I,J) = 0.0                                                       XYLT
      DO 5131 M = 1,6                                                    XYLT
      DO 5131 N = 1,6                                                    XYLT
 5131 A(I,J) = A(I,J) + (DERIV(M,I)*DERIV(N,K)*A(M,N))                   XYLT
 5130 CONTINUE                                                           XYLT
      DO 5133 I = 1,6                                                    XYLT
      DO 5133 K = I,6                                                    XYLT
      A(I,K) = A(I,K+6)                                                  XYLT
 5133 A(K,I) = A(I,K)                                                    XYLT
C                                                                        XYLT
C----READY AGAIN TO COMPUTE THE VECTORS AND THEIR VAR-COVAR MATRICES     XYLT
      WRITE (NTOUT,40)                                                   XYLT
   40 FORMAT(1H0,120(1H-)/17H0DIRECT ELEMENTS.)                          XYLT
      I1 = 2                                                             XYLT
      I2 = 1                                                             XYLT
C----CALL METRIC                                                         XYLT
      CALL XRY746                                                        XYLT
      WRITE (NTOUT,42)                                                   XYLT
   42 FORMAT(1H0,120(1H-)/11H0NORMAL END)                                XYLT
C                                                                        XYLT
C----MADE IT ....                                                        XYLT
      GO TO 9991                                                         XYLT
 9999 IQUIT = 1                                                          XYLT
 9991 CONTINUE                                                           XYLT
      END                                                                XYLT
      SUBROUTINE XRY746                                                  LTME
C----SUBROUTINE METRIC                                                   LTME
C                                                                        LTME
C----THIS SUBROUTINE COMPUTES THE VECTORS FROM A METRIC TENSOR, THEIR    LTME
C    E.S.D.-S, THE VARIANCE-COVARIANCE MATRICES AND THE CORRELATION      LTME
C    MATRICES.                                                           LTME
C                                                                        LTME
C----COMMON PACKAGE FOR SYSTEM                                           LTME
      COMMON /SYS/ IDCOMP(2),NTIN,NTOUT,NTAPEA,NTAPEB,NTAPEC,NTAPED,     LTME
     $NTAPEE,NTAPEF,NTAPEG,NTAPEH,NTAPEI,NTAPEJ,NCDBUF,NTINM,NTOUTM,     LTME
     $NTPCHM,IDPROG(2),ITOOB,NPAGE,LINMX,LINCT,NCRDS,ISHORT,             LTME
     $TITLE(19),CARD(21),QTQ(100),ICYCLE,IQUIT,                          LTME
     $ JUNK,JUNK1,JUNK2,JUNK3,JUNK4,IOZY,IOZTM,ICDC(20)                  LTME
      DIMENSION IQTQ(100),ICARD(21),NTIND(1)                             LTME
      EQUIVALENCE (QTQ(1),IQTQ(1)),(NTAPEC,NTPCH),(NTCU,NTPCHM)          LTME
      EQUIVALENCE (CARD(1),ICARD(1)),(NTIND(1),NTIN)                     LTME
C----COMMON STATEMENTS                                                   LTME
      COMMON/OMNIA/A(7,14),B(6,12),C(7),CC(3,6),D(9),                    LTME
     *DERIV(6,6),ESD1(6),ESD2(9),G(5),IG(5),IS(5),I1,I2,V(7)             LTME
C                                                                        LTME
      DATA CON /57.295780/                                               LTME
C                                                                        LTME
C                                                                        LTME
C----COMPUTE E.S.D.-S AND CORRELATION MATRIX OF METRIC TENSOR            LTME
C    A = VAR-COVAR MATRIX AND CORRELEATION MATRIX OF METRIC TENSOR       LTME
C    B = VAR-COVAR MATRIX AND CORRELATION MATRIX OF VECTORS AND ANGLES   LTME
C    DERIV = DERIVATIVES                                                 LTME
C    C = METRIC TENSOR IN ORDER 11, 22, 33, 12, 13 , 23                  LTME
C    D = CORRESPONDING LENGTHS AND ANGLES                                LTME
      DO 97 I = 1,6                                                      LTME
      DO 99 J = 7,12                                                     LTME
      A(I,J) = 0.0                                                       LTME
   99 B(I,J) = 0.0                                                       LTME
      A(I,I+6) = 1.0                                                     LTME
   97 B(I,I+6) = 1.0                                                     LTME
      DO 101 I = 1,6                                                     LTME
  101 ESD1(I) = FSQRT(A(I,I))                                             LTME
      DO 103 I = 1,6                                                     LTME
      DO 103 J = 1,6                                                     LTME
      ESD = ESD1(I)*ESD1(J)                                              LTME
      IF(ESD) 103,103,105                                                LTME
  105 A(I,J+6) = A(I,J)/ESD                                              LTME
  103 CONTINUE                                                           LTME
C                                                                        LTME
C----COMPUTE VECTORS AND THEIR DERIVATIVES                               LTME
      DO 111 I = 1,6                                                     LTME
      DO 111 J = 1,6                                                     LTME
  111 DERIV(I,J) = 0.0                                                   LTME
      DO 113 I = 1,3                                                     LTME
      D(I) = FSQRT(C(I))                                                  LTME
  113 DERIV(I,I) = 0.5/D(I)                                              LTME
      DO 115 I = 1,2                                                     LTME
      N = I + 1                                                          LTME
      DO 115 J = N,3                                                     LTME
      M = I + J + 1                                                      LTME
      D(M) = C(M)/(D(I)*D(J))                                            LTME
      DERIV(I,M) = -0.5*D(M)/C(I)                                        LTME
      DERIV(J,M) = -0.5*D(M)/C(J)                                        LTME
  115 DERIV(M,M) = 1./(D(I)*D(J))                                        LTME
C----DERIV(I,J) CONTAINS THE DERIVATIVE OF D(J) WRT C(I)                 LTME
C                                                                        LTME
C----CALCULATE VAR-COVAR MATRIX INTO B                                   LTME
      DO 121 I = 1,6                                                     LTME
      DO 121 K = I,6                                                     LTME
      B(I,K) = 0.0                                                       LTME
      DO 123 M = 1,6                                                     LTME
      DO 123 N = 1,6                                                     LTME
  123 B(I,K) = B(I,K) + (DERIV(M,I)*DERIV(N,K)*A(M,N))                   LTME
  121 B(K,I) = B(I,K)                                                    LTME
C----TREAT ROUNDING OFF ERRORS IN VAR-COVAR OF COS OF ANGLES.            LTME
      DO 125 I = 1,6                                                     LTME
      DO 125 J = 4,6                                                     LTME
      IF(ABS(B(I,J))-1.0E-15) 127,127,125                                LTME
  127 B(I,J) = 0.0                                                       LTME
      B(J,I) = 0.0                                                       LTME
  125 CONTINUE                                                           LTME
C                                                                        LTME
C----COMPUTE E.S.D.-S AND CORRELEATION MATRIX                            LTME
      DO 131 I = 1,6                                                     LTME
  131 ESD2(I) = FSQRT(B(I,I))                                             LTME
      DO 133 I = 1,6                                                     LTME
      DO 133 J = 1,6                                                     LTME
      ESD = ESD2(I)*ESD2(J)                                              LTME
      IF(ESD) 133,133,135                                                LTME
  135 B(I,J+6) = B(I,J)/ESD                                              LTME
  133 CONTINUE                                                           LTME
C                                                                        LTME
C----COMPUTE ANGLES AND THEIR VARIANCES                                  LTME
      DO 141 I = 4,6                                                     LTME
      AJUNK = ACOS(D(I))                                                 LTME
      D(I+3) = CON*AJUNK                                                 LTME
      AJUNK = SIN(AJUNK)                                                 LTME
  141 ESD2(I+3) = CON*ESD2(I)/AJUNK                                      LTME
C                                                                        LTME
C----WRITE THE RESULTS                                                   LTME
      WRITE (NTOUT,10) I1,(C(I),I=1,6),(ESD1(I),I=1,6)                   LTME
   10 FORMAT(91H0METRIC TENSOR ELEMENTS AND THEIR E.S.D-S LISTED IN THE  LTME
     *ORDER 11, 22, 33, 12, 13, 23 IN A**,I2/ 6(1X,E13.6) / 6(1X,E13.6)/ LTME
     * 45H0VARIANCE-COVARIANCE AND CORRELATION MATRICES)                 LTME
      DO 201 I = 1,6                                                     LTME
      WRITE (NTOUT,12) (A(I,J),J=1,12)                                   LTME
   12 FORMAT(6(1X,E13.6),2X,6F7.4)                                       LTME
  201 CONTINUE                                                           LTME
      DO 205 I = 1,12                                                    LTME
      AJUNK = B(6,I)                                                     LTME
      B(6,I) = B(4,I)                                                    LTME
      B(4,I) = AJUNK                                                     LTME
  205 CONTINUE                                                           LTME
      WRITE (NTOUT,14) I2,D(1),D(2),D(3),D(6),D(5),D(4),D(9),D(8),D(7),  LTME
     * ESD2(1),ESD2(2),ESD2(3),ESD2(6),ESD2(5),ESD2(4),ESD2(9),ESD2(8),  LTME
     * ESD2(7)                                                           LTME
   14 FORMAT(15H0VECTORS IN A**,I2,69H  -- A, B, C, COSA, COSB, COSC, AL LTME
     *PHA, BETA, GAMMA, AND THEIR E.S.D-S/9F14.6/9F14.6/                 LTME
     * 45H0VARIANCE-COVARIANCE AND CORRELATION MATRICES)                 LTME
      DO 203 I = 1,6                                                     LTME
      WRITE (NTOUT,12) B(I,1),B(I,2),B(I,3),B(I,6),B(I,5),B(I,4),        LTME
     * B(I,7),B(I,8),B(I,9),B(I,12),B(I,11),B(I,10)                      LTME
  203 CONTINUE                                                           LTME
C                                                                        LTME
C----MADE IT...                                                          LTME
      RETURN                                                             LTME
      END                                                                LTME
      FUNCTION FSQRT(FNUM)
      IF(FNUM.GT.0.0) GO TO 2
      FSQRT=0.0
      IF(FNUM.EQ.0.0) GO TO 3
      WRITE(*,100)
100   FORMAT(' Square Root of -ve Number  Answer set to 0.0')
      GO TO 3
2     FSQRT=SQRT(FNUM)
3     RETURN
      END
