      SUBROUTINE RFINE
C
C             STRUCTURE FACTOR CALCULATION AND LEAST SQUARES REFINEMENT
C             MAXIMUM OF QQATM ATOMS, 10 SCALE FACTORS,
C              20 SCATTERING FACTORS, QQPRM LEAST-SQUARES PARAMETERS
C              AND 24 SYMMETRY OPERATIONS
C
C             JULY  1966    BY L. W. FINGER
C             MODIFIED NOV 7,1966 FOR BLOCKED STRUCTURE FACTOR RECORDS
C             MODIFIED JAN 30,1967 FOR EXTINCTION CORRECTION
C             SPECIAL PARAMETER REFINEMENT ADDED FEB. 1972
C             ROBUST/RESISTANT FEATURES ADDED MAY 1977
C             TWIN COMPONENTS, R-FACTORS BY LEVELS, BECKER-COPPENS
C                EXTINCTION AND ANALYTICAL ABSORPTION ADDED MAY 1987
C
      INCLUDE 'RFINE.PRM'
      INCLUDE 'A.INC'
      INCLUDE 'C.INC'
      INCLUDE 'CNSTRS.INC'
      INCLUDE 'D.INC'
      INCLUDE 'E.INC'
      INCLUDE 'F.INC'
      INCLUDE 'FBUF.INC'
      INCLUDE 'G.INC'
      INCLUDE 'LSMAT.INC'
      INCLUDE 'EXTN.INC'
      INCLUDE 'TWIN.INC'
      INCLUDE 'XPAR.INC'
      INCLUDE 'ABSORB.INC'
C
      DIMENSION KREJ(4),NEW(8),R(36),DIN(3),DOUT(3)
      DIMENSION EGGY(6),BBETA(6,QQATM),SIGB(6,QQATM)
      DIMENSION SIGO(QQATM),SIGXYZ(3,QQATM),TFB(6)
C
      EQUIVALENCE (R,DADX)
      EQUIVALENCE (HIB,DADO),(KB,DADO(21)),(JB,DADO(41))
C
      DATA KREJ/2H  ,2H* ,2H  ,2H**/
C
   10 FORMAT(18A4)
   12 FORMAT(15I1,2F5.0,3I1,3I2,4I1,F5.0,30I1)
   14 FORMAT(26H0 CONTROL PARAMETERS ARE  /15I2,2F7.2,3I2,3I3,4I2,F7.2,
     $30I2)
   16 FORMAT(9F8.0,2A4)
   18 FORMAT(1H1,18A4)
   24 FORMAT(2A3,3F6.4/2I2,4X,9F7.5,I2)
   26 FORMAT(80I1)
   28 FORMAT(1H1,18A4,5X,17HSTRUCTURE FACTORS,10X,4HPAGE,I3/1H0,7X,1HH,
     1 4X,1HK,4X,1HL,6X,6HF(OBS),5X,7HF(CALC),5X,7HA(CALC),5X,7HB(CALC),
     2 5X,7HDELTA F,3X,'DELTA/SIGMA  EXTIN     ABS    BIWT'/)
   30 FORMAT(4X,3I5,F12.3,1X,A1,F10.3,3F12.3,F14.4,3F8.4,2X,2A2)
   32 FORMAT(6E12.6)
   34 FORMAT( /10X,8HSUM FCAL,20X,F15.2/ 10X,31HSTANDARD DEV OF UNIT WEI
     1GHT OBS,F12.2)
   36 FORMAT(46H0INPUT PARAMETERS READ FROM FILE WITH LABEL **,3A4,2H**)
   38 FORMAT(I4)
   40 FORMAT(6H IREJ=,I2,8H DELMAX=,F7.2)
   42 FORMAT(37H0FORMAT FOR STRUCTURE FACTOR INPUT IS/1X,18A4)
   44 FORMAT(49H0SIGMA WILL BE COMPUTED BY SIGMA=(SQRT(SIGMA**2+(,F6.4,
     17H*F)**2))
   48 FORMAT('ATOM',3X,2A3,3F8.4/'BIJ',4X,2A3,6F8.4,1X,'4')
   50 FORMAT(10F8.4)
   52 FORMAT(9F8.6,F6.4,I2)
   54 FORMAT(10F8.6)
   56 FORMAT(' WRONG VALUE OF NFLAG,ACTUAL VALUE IS',I3)
   58 FORMAT(3I1,F8.5,7E10.4)
   60 FORMAT(I2/(6E12.5))
C             INITIALIZE
      PI=ACOS(-1.0)
      TEMPZZ=8.*PI*PI
      NPAR=0
      FISGN=1.0
      MAXPG=53
      IFIRST=-1
C     IF(IWROTE.NE.0) FOURIER FILE ALREADY EXISTS.
      IWROTE=0
      REWIND ISC1
      REWIND ISC2
C             STARTING POINT FOR CYCLE
C
  300 READ(IN,10)TITLE
      WRITE(IOUT,18)TITLE
      READ(IN,12)MFSEL,MSF,MCORR,MBODAN,MELVIB,MREF,IFILE,NEW,(BUFF(I),
     1I=1,2),(IBUFF(I),I=3,12),BUFF(13),(IBUFF(I),I=15,32),IEF,ILVG
     2,ITWIN,JABS
      IF(IEF.NE.0)IWROTE=1
      WRITE(IOUT,14)MFSEL,MSF,MCORR,MBODAN,MELVIB,MREF,IFILE,NEW,(BUFF(
     1I),I=1,2),(IBUFF(I),I=3,12),BUFF(13),(IBUFF(I),I=15,32),IEF,ILVG,
     2ITWIN,JABS
      IFIRST=IFIRST+1
      NBUFF=0
      NFBR=0
C             CHECK IF FIRST CYCLE OF RUN
      IF(IFIRST)340,310,340
  310                        DO 320 I=1,8
      IF(NEW(I))320,330,320
  320 CONTINUE
C             ALL PARAMETERS NEW, DO NOT READ SAVE UNIT
      GO TO 340
C             READ SAVE UNIT TO START OF REFLECTION LIST
  330 READ(ISC1)
     1KCYL,LABEL,RECELL,G,VOL,NSYM,NATOM,NSCAT,JCELL,MRANK,NPARA,NSF,DMI
     2N,DMAX,DELMAX,ICENT,(((RMAT(I,J,K),I=1,3),J=1,3),K=1,NSYM),((TRANS
     3(I,J),I=1,3),J=1,NSYM),((CELLP(I,J),I=1,3),J=1,JCELL),(IPARA(I),I=
     41,NPARA),((FCURVE(I,J),I=1,9),J=1,NSCAT),(SCALE(I),I=1,NSF),IREJ,F
     5OMAX,FOMIN,RHOMN,RHOMX,(JSCRF(I),I=1,NSF),NRSF
     9,NPAR,(PAR(I),(PARNAM(J,I),J=1,2),IPAR(I),IRPAR(I),I=1,NPAR),FISGN
     O,IFSGN,EGGY,KEXTRF
      READ(ISC1 )
     1((ANOM(J,I),J=1,2),I=1,NSCAT),IEXT,EXTIN,((TAG(J,I),J=1,2),ISOT(I)
     2,(ISCAT(J,I),J=1,2),OCCA(I),OCCUP(I),SITE(I),(XYZ(J,I),J=1,3),
     3(BETA(J,I),J=1,6),(CR(J,I),J=1,10),(DR(J,I),J=1,15),I=1,NATOM),
     4 LDEGN,(KDEGN(I),I=1,LDEGN),NT,TFPAR,XWIN,WAVE,UBMAT,(TFRAC(I)
     5,((TWINTM(J,K,I),J=1,3),K=1,3),JTWNMT(I),I=1,NT),BDF,WVLNTH
     6,HELPE,HELPF,IZACH,NIEXTA,NABSRN,KABSRF,ABSCF,JTFLG
      NDEGN=LDEGN/2
      WRITE(IOUT,36)LABEL
      DELTMP=DELMAX
C             CHECK FOR NEW PARAMETERS
  340 CALL LOADNW(NEW,EGGY,IFIRST,NSF,NPARA,ITWIN,JABS)
      IF(NDEGN.NE.0.AND.(MFSEL.NE.0.OR.IFIRST.EQ.0))CALL DERIVS(NEW(7))
      IF(NDEGN.EQ.0)MDEGN=0
      IF(MFSEL.NE.1)GO TO 380
C             IF LEAST-SQUARES, CLEAR MATRICES
      IC=MDEGN
      IF(IWROTE.NE.0.AND.IREJ.GE.2)CALL MEDWT(DELMAX)
                             DO 360 IA=1,MRANK
      B(IA)=0.0
                             DO 360 IB=IA,MRANK
      IC=IC+1
  360 A(IC)=0.0
C             INITIALIZE FOR CYCLE
  380 IF(MFSEL.EQ.2)GO TO 940
      NEXTRA=13
      NSC9=NSCAT+NEXTRA
      LBUFF=1022/NSC9
      LBUFF=NSC9*LBUFF+1
      NPRINT=0
      SUMFC=0.0
      RMINUS=0.0
      FLOSS=0.
      NFLOSS=0
      EM=0.
      SSUM=0.
      CALL RCALC(SSUM,0)
      IF(MSF.LE.1)GO TO 390
C***********************************************************************
C        THIS SECTION SHOULD PREPARE TO WRITE A STRUCTURE FACTOR FILE
C          FOR A FOURIER PROGRAM
C             PREPARE BINARY STRUCTURE FACTOR FILE
      REWIND ISC4
      WRITE(ISC4)NSYM,VOL,(((RMAT(I,J,K),I=1,3),J=1,3),K=1,NSYM),((TRANS
     1(J,I),J=1,3),I=1,NSYM)
C***********************************************************************
  390 DELTMP=0.0
      IF(IBUFF(12).LT.5)GO TO 400
C             CHANGE REJECTION MODE
         IREJ=IBUFF(12)-5
         IF(BUFF(13).GT.0.)DELTMP=BUFF(13)
  400 IF(DELTMP.GT.0.0)DELMAX=DELTMP
      IF(IBUFF(12).GE.5.OR.((IREJ.EQ.2.OR.IREJ.EQ.3).AND.IWROTE.EQ.0))
     1WRITE(IOUT,40)IREJ,DELMAX
  410 IF(NEW(8))450,420,450
C             NEW STRUCTURE FACTORS NOT TO BE READ
  420 IF(IFIRST)440,430,440
C             FIRST CYCLE, STRUCTURE FACTORS FROM SAVE FILE
  430 ASSIGN 570 TO INTYPE
      GO TO 570
C             NOT FIRST CYCLE, STRUCTURE FACTORS FROM INTERMEDIATE UNIT
  440 ASSIGN 710 TO INTYPE
      GO TO 710
C             NEW STRUCTURE FACTORS
  450 IREJ=IBUFF(12)
      IF(IBUFF(11))470,460,470
C             NEW STRUCTURE FACTORS FROM STANDARD INPUT UNIT
  460 INN=IN
      GO TO 480
C             NEW STRUCTURE FACTORS FROM UNIT 4 (ISC3)
  470 INN=ISC3
  480 DELMAX=BUFF(13)
      IWT=IBUFF(10)+1
      FOMAX=0.0
      FOMIN=1.0E+6
      RHOMN=1.0E+6
      RHOMX=0.0
      READ(IN,10)FRMT
      WRITE(IOUT,42)FRMT
C  IF IWT.EQ.6 AN 'IGNORANCE FACTOR' WILL BE READ IN FOR COMPUTING
C    WEIGHTS
      IF(IWT.EQ.6)READ(IN,16)WTFAC
      IF(IWT.EQ.6)WRITE(IOUT,44)WTFAC
C  NEW(8) SPECIFIES THE INPUT FORMAT AS FOLLOWS. 
C           1=OUTPUT OF ABSORB, GEOLAB VERSION
C           2=XRAY67 HKL CARD SEQUENCE
C           3=BETA TO BE COMPUTED
C           4=HKL CARD SEQUENCE, XRAY BETA TO BE COMPUTED.
      IF(NEW(8).LT.3)GO TO 490
      BETFAC=(2.*WVLNTH**2/VOL**2)*1.0E+8
      IF(NEW(8).EQ.4)BETFAC=BETFAC*.0795
  490 ASSIGN 500 TO INTYPE
      IBUFF(1)=0
      NBUFF=1
  500 IEF=NEW(8)
      CALL INPUT(IEF)
      IF(LBUFF-NBUFF)510,510,520
  510 WRITE(ISC2)(IBUFF(I),I=1,LBUFF)
      NBUFF=1
  520 IBUFF(NBUFF+1)=IH(1)
      IBUFF(NBUFF+2)=IH(2)
      IBUFF(NBUFF+3)=IH(3)
      BUFF(NBUFF+4)=FOBS
C  HERE WT = SIGMA F
      BUFF(NBUFF+5)=WT
      BUFF(NBUFF+6)=RHO
      IBUFF(NBUFF+7)=ISCL
      BUFF(NBUFF+8)=EXBETA
      IBUFF(NBUFF+9)=MREJ
      DO 530 I=10,NEXTRA
  530 BUFF(NBUFF+I)=EXANG(I-9)
                             DO 540 I=1,NSCAT
  540 BUFF(I+NBUFF+NEXTRA)=SCAT(I)
      NBUFF=NBUFF+NSC9
      IF(IEF)550,790,550
  550 IBUFF(1)=1
      IBUFF(NBUFF-NSC9+1)=3000
      WRITE(ISC2)(IBUFF(I),I=1,LBUFF)
  560 END FILE ISC2
      IF(IFILE.LT.2)GO TO 900
C     READ SAVED SIGMA. FIRST SKIP OVER CONSTRAINTS AND VAR-COV MATRIX.
      READ(ISC1)IBUFF(1)
      READ(ISC1)IBUFF(1)
      READ(ISC1)((BBETA(J,I),J=1,6),SIGO(I),(SIGXYZ(J,I),J=1,3),(SIGB(J,
     1I),J=1,6),I=1,NATOM)
      GO TO 900
C             READ STRUCTURE FACTORS FROM SAVE UNIT
  570 IF(NBUFF)650,580,650
  580 READ(ISC1)(IBUFF(I),I=1,LBUFF)
      IF(NEW(4))590,640,590
C             NEW SCATTERING FACTORS
  590 ITEMP=7-NSC9
  600 ITEMP=ITEMP+NSC9
      IF(ITEMP-LBUFF)610,640,640
  610 RHO=BUFF(ITEMP)
                             DO 630 I=1,NSCAT
      TEMP=FCURVE(9,I)+ANOM(1,I)
                             DO 620 J=1,4
  620 TEMP=TEMP+FCURVE(2*J-1,I)*EXP(-FCURVE(2*J,I)*RHO)
  630 BUFF(I+ITEMP+6)=TEMP
      GO TO 600
  640 WRITE(ISC2)(IBUFF(I),I=1,LBUFF)
      NBUFF=1
  650 IH(1)=IBUFF(NBUFF+1)
      IH(2)=IBUFF(NBUFF+2)
      IH(3)=IBUFF(NBUFF+3)
      FOBS=BUFF(NBUFF+4)
C  HERE WT = SIGMA F
      WT=BUFF(NBUFF+5)
      RHO=BUFF(NBUFF+6)
      ISCL=IBUFF(NBUFF+7)
      EXBETA=BUFF(NBUFF+8)
      MREJ=IBUFF(NBUFF+9)
      DO 660 I=10,NEXTRA
  660 EXANG(I-9)=BUFF(NBUFF+I)
                             DO 670 I=1,NSCAT
  670 SCAT(I)=BUFF(I+NBUFF+NEXTRA)
      NBUFF=NBUFF+NSC9
      IF(LBUFF-NBUFF)680,680,690
  680 NBUFF=0
  690                        DO 700 I=1,3
  700 H(I)=IH(I)
      IF(IH(1)-3000)790,560,790
C             READ FROM INTERMEDIATE FILE
  710 IF(NBUFF)730,720,730
  720 READ(ISC2)(IBUFF(I),I=1,LBUFF)
      NBUFF=1
  730 IH(1)=IBUFF(NBUFF+1)
      IH(2)=IBUFF(NBUFF+2)
      IH(3)=IBUFF(NBUFF+3)
      FOBS=BUFF(NBUFF+4)
      WT=BUFF(NBUFF+5)
      RHO=BUFF(NBUFF+6)
      ISCL=IBUFF(NBUFF+7)
      EXBETA=BUFF(NBUFF+8)
      MREJ=IBUFF(NBUFF+9)
      DO 740 I=10,NEXTRA
  740 EXANG(I-9)=BUFF(NBUFF+I)
                             DO 750 I=1,NSCAT
  750 SCAT(I)=BUFF(I+NBUFF+NEXTRA)
      NBUFF=NBUFF+NSC9
      IF(LBUFF-NBUFF)760,760,770
  760 NBUFF=0
  770                        DO 780 I=1,3
  780 H(I)=IH(I)
      IF(IH(1)-3000)790,900,790
C             COMPUTE STRUCTURE FACTOR
  790 CALL SFAC
      FPAREN=(ACALC**2+BCALC**2)
      IF(NT.EQ.0.AND.FPAREN.LT.1.0E-4)GO TO 895
      FCALC=FPAREN*TFPAR
      FCALCA=FCALC
      IF(NT.GT.0)THEN
        DO 800 I=1,NT
          FCALC=FCALC+TFRAC(I)*(ATWN(I)**2+BTWN(I)**2)
  800	CONTINUE
      END IF
      FCALC=ABS(FCALC)
      IF(IBDF(1).GE.4.OR.NABSRN.GT.0) THEN
C CALCULATE DIRECTION COSINES OF INCIDENT AND DIFFRACTED BEAM
        SCHI = SIN(EXANG(3))
        CCHI = COS(EXANG(3))
        STPO = SIN(EXANG(1)+EXANG(2))
        CTPO = COS(EXANG(1)+EXANG(2))
        STMO = SIN(EXANG(1)-EXANG(2))
        CTMO = COS(EXANG(1)-EXANG(2))
        SPHI = SIN(EXANG(4))
        CPHI = COS(EXANG(4))
C
C---- VECTOR IN DIRECTION OF INCIDENT BEAM (IN THE PHI SYSTEM)
C
        DIN(1) = -STPO*CCHI*CPHI-CTPO*SPHI
        DIN(2) = -STPO*CCHI*SPHI+CTPO*CPHI
        DIN(3) = -STPO*SCHI
C
C---- VECTOR IN DIRECTION OF REFLECTED BEAM (PHI SYSTEM)
C
        DOUT(1) = STMO*CCHI*CPHI-CTMO*SPHI
        DOUT(2) = STMO*CCHI*SPHI+CTMO*CPHI
        DOUT(3) = STMO*SCHI
      END IF
C        CALCULATE EXTINCTION FACTOR
      IF(IZACH.NE.0)THEN
        EXX=EXBETA*EXTIN(1)*FCALC*1.0E-4
        EXY=1.0/SQRT(1.0+EXX)
        DRDE(1)=1.0E-4
      ELSE
        CALL EXTICA(NONPOS,INEUT,DIN,DOUT)
      END IF
      CALL ABSRP(DIN,DOUT)
      FCALC=SQRT(FCALC*EXY*ABSFCT)
      FCALCA=SQRT(FCALCA*EXY*ABSFCT)
      IF(IFSGN.NE.0)FMINUS=SQRT(EXY*(AMINUS**2+BMINUS**2)*ABSFCT)
      SUMFC=SUMFC+FCALC
      NREJ=0
      FOBS=FOBS/SCALE(ISCL)
      SIGF=WT/SCALE(ISCL)
      IF(IREJ.EQ.0)GO TO 830
      GO TO (810,810,830,820),IREJ
C             REJECT IF LESS THAN MIN OBS
  810 NREJ=MREJ
      GO TO 830
C             USE SPECIAL ROUTINE REJECT
  820 CALL REJECT
C             CALCULATE CONTRIBUTION TO R
C* * * RCALC ALSO CONVERTS 'WT' TO A WEIGHT  (SIGMA F TO HERE)
  830 CALL RCALC(SCALE(ISCL),1)
      IF(NREJ.NE.0.OR.IFSGN.EQ.0)GO TO 840
      TEMP=FOBS-FMINUS
      IF(MREF.NE.0)TEMP=FOBS**2-FMINUS**2
      RMINUS=RMINUS+WT*TEMP**2
  840 IF(MSF.EQ.0)GO TO 890
      IF(MSF.EQ.3)GO TO 880
C             WRITE PRINTED OUTPUT FILE
      NPAGE=NPRINT/MAXPG+1
      IF(MAXPG*(NPAGE-1)-NPRINT)860,850,860
C             PAGE OUTPUT
  850 WRITE(IOUT,28)TITLE,NPAGE
C             PRINT STRUCTURE FACTOR
  860 TEMP=FOBS-FCALC
      TEMP1=DELTA/SIGF
      WRITE(IOUT,30)IH,FOBS,KREJ(MREJ+1),FCALC,ACALC,BCALC,TEMP,TEMP1,
     1 EXY,ABSFCT,BIW,KREJ(NREJ+1),KREJ(MREJ+3)
      NPRINT=NPRINT+1
  870 IF(MSF.LT.2)GO TO 890
C***********************************************************************
C        THIS SECTION SHOULD ADD THE CURRENT REFLECTION TO THE FOURIER
C          BINARY STRUCTURE FACTOR FILE
  880 FTEMP=FCALC
      IF(FTEMP.EQ.0.0)FTEMP=100.0
      IFBUF(NFBR+1)=IH(1)
      IFBUF(NFBR+2)=IH(2)
      IFBUF(NFBR+3)=IH(3)
      FBUF(NFBR+4)=FOBS*FCALCA/FTEMP
      FBUF(NFBR+5)=FCALCA
      FBUF(NFBR+6)=ACALC
      FBUF(NFBR+7)=BCALC
      IFBUF(NFBR+8)=NREJ
      FBUF(NFBR+9)=SIGF
      FBUF(NFBR+10)=EXY
      NFBR=NFBR+10
      IF(NFBR.LT.500)GO TO 890
      WRITE(ISC4)FBUF
      NFBR=0
C***********************************************************************
C  ADD TO LEAST SQUARES MATRIX IF REFINED AND REFLECTION NOT
C     REJECTED.
  890 IF(MFSEL.NE.0.AND.NREJ.EQ.0)CALL MATRIX(ILVG)
C             GO TO GET NEW STRUCTURE FACTOR
  895 GO TO INTYPE,(710,570,500)
C             END OF LIST FOUND, REWIND FILES
  900 REWIND ISC1
      REWIND ISC2
      IF(ILVG.NE.0)REWIND 50
      IF(MSF.LT.2)GO TO 910
C***********************************************************************
C THIS SECTION ENDS THE FOURIER FILE
      IFBUF(NFBR+1)=99
      WRITE(ISC4)FBUF
      IWROTE=1
C***********************************************************************
C             COMPUTE AND PRINT OUT R FACTORS
  910 CALL RCALC(SCALE(1),2)
      TEMP=SQRT(SSUM/FLOAT(NFLOSS-MRANK))
      IF(IWROTE.EQ.0)DELTMP=7.0337*TEMP
      IF(IREJ.NE.2.AND.IREJ.NE.3)GO TO 920
      EM=EM/NFLOSS
      TEMP=TEMP*(1.+(FLOAT(MRANK)*(1.-EM))/(FLOAT(NFLOSS)*EM))/EM
  920 WRITE(IOUT,34)SUMFC,TEMP
      IF(MFSEL.EQ.0)GO TO 940
      CALL INVERT(MCORR,SIGB,BBETA,SIGXYZ,SIGO,EGGY,NSF,TEMP)
C             CHECK TO SEE IF AUXILLARY FUNCTIONS TO BE DONE
  940 IF(MBODAN.NE.0)CALL BODAN1(MBODAN)
      IF(MELVIB.NE.0)CALL ELVIB1(MELVIB,MBODAN)
      IF((MBODAN.EQ.0.AND.MELVIB.EQ.0).OR.NDEGN.EQ.0)GO TO 990
      REWIND ISC5
      READ(ISC5)LDEGN,(KDEGN(I),I=1,LDEGN)
  990 READ(IN,26)IEF,NFLAG,ITEM1,ITEM2
      IF(ITEM1.NE.0.OR.ITEM2.NE.0)ISC1=10*ITEM1+ITEM2
C        CHECK FOR LEGAL RANGE OF 'NFLAG'
      IF(NFLAG.LT.0.OR.NFLAG.GT.3)GO TO 1050
C        NO CARDS IF NFLAG=0
      IF(NFLAG.EQ.0)GO TO 1060
C        PUNCH XRAY67 CARDS IF NFLAG=1 OR 3
      IF(NFLAG.EQ.2)GO TO 1020
      DO 1010 I=1,NATOM
      DO 1000 J=1,6
 1000 TFB(J)=BETA(J,I)*EGGY(J)
 1010 WRITE(IPUN,48)(TAG(J,I),J=1,2),(XYZ(J,I),J=1,3),(TAG(J,I),J=1,2)
     $,TFB
 1020 IF(NFLAG.EQ.1)GO TO 1060
C        PUNCH 'RFINE' CARDS IF NFLAG=2 OR 3
C        PUNCH SCALE FACTOR AND EXTINCTION CARDS
      WRITE(IPUN,50)(SCALE(I),I=1,NSF)
      IF(IZACH.NE.0)THEN
        WRITE(IPUN,58)IBDF(1),IBDF(2),INEUT,WVLNTH,EXTIN
      ELSE
        NV=1
        IF(IBDF(1).EQ.1)IS=3
        IF(IBDF(1).EQ.2)IS=4
        IF(IBDF(1).EQ.3)THEN
          IS=3
          NV=2
        END IF
        IF(IBDF(1).EQ.4.OR.IBDF(1).EQ.6)THEN
          IS=4
          NV=6
        END IF
        IF(IBDF(1).EQ.5)THEN
          IS=10
          IV=6
        END IF
        IF(IBDF(1).NE.6)THEN
          TEMP=0.0
        ELSE
          TEMP=BDF(4)
        END IF
        WRITE(IPUN,58)IBDF(1),IBDF(2),INEUT,WVLNTH,(BDF(I+IS-1),I=1,NV)
     1      ,TEMP
        IF(INEUT.EQ.2)WRITE(IPUN,50)DIHED,THETM
      END IF
      IF(NABSRN.GE.1)THEN
        N=NABSRN*(NABSRN+2)
        WRITE(IPUN,60)NABSRN,(ABSCF(I),I=1,N)
      END IF
      DO 1030 I=1,NATOM
      K=0
      IF(ISOT(I).GT.1)K=1
      IFBUF(I)=(-ISOT(I)**3+6*ISOT(I)**2-2*ISOT(I))/3
      WRITE(IPUN,24)(TAG(J,I),J=1,2),SITE(I),OCCUP(I),OCCA(I),(ISCAT(J,I
     $),J=1,2),(XYZ(J,I),J=1,3),(BETA(J,I),J=1,6),K
      IF(K.EQ.0) GO TO 1030
C        PUNCH 3RD CUMULANT
      K=ISOT(I)-2
      WRITE(IPUN,52)(CR(J,I),J=1,10),K
      IF(K.EQ.0)GO TO 1030
C        PUNCH 4TH CUMULANT
      WRITE(IPUN,54)(DR(J,I),J=1,15)
 1030 CONTINUE
      IF(NPAR.EQ.0)GO TO 1040
C        PUNCH THE SPECIAL PARAMETERS
      WRITE(IPUN,38)NPAR
      CALL PUNPAR
C        PUNCH THE EXTINCTION AND TEMPERATURE FACTOR TYPES
 1040 WRITE(IPUN,26)IEXT,(IFBUF(I),I=1,NATOM)
      GO TO 1060
 1050 WRITE(IOUT,56)NFLAG
 1060 CONTINUE
      IF(MSF-2)1080,1070,1070
 1070 WRITE(ISC4)NATOM,((TAG(J,I),J=1,2),ISOT(I),(ISCAT(J,I),J=1,2),
     $OCCA(I),OCCUP(I),SITE(I),(XYZ(J,I),J=1,3),(BBETA(J,I),J=1,6),
     $I=1,NATOM),(SIGO(I),(SIGXYZ(J,I),J=1,3),(SIGB(J,I),J=1,6),I=1,NATO
     $M)
      END FILE ISC4
 1080 CONTINUE
      IF(IEF)300,1090,300
C             END OF CALCULATION
 1090 REWIND ISC1
      IF(MOD(IFILE,2))1140,1100,1140
C             WRITE SAVE UNIT
 1100 LDEGN=2*NDEGN
      WRITE(ISC1)
     1KCYL,LABEL,RECELL,G,VOL,NSYM,NATOM,NSCAT,JCELL,MRANK,NPARA,NSF,DMI
     2N,DMAX,DELMAX,ICENT,(((RMAT(I,J,K),I=1,3),J=1,3),K=1,NSYM),((TRANS
     3(I,J),I=1,3),J=1,NSYM),((CELLP(I,J),I=1,3),J=1,JCELL),(IPARA(I),I=
     41,NPARA),((FCURVE(I,J),I=1,9),J=1,NSCAT),(SCALE(I),I=1,NSF),IREJ,F
     5OMAX,FOMIN,RHOMN,RHOMX,(JSCRF(I),I=1,NSF),NRSF
     9,NPAR,(PAR(I),(PARNAM(J,I),J=1,2),IPAR(I),IRPAR(I),I=1,NPAR),FISGN
     O,IFSGN,EGGY,KEXTRF
      WRITE(ISC1)
     1((ANOM(J,I),J=1,2),I=1,NSCAT),IEXT,EXTIN,((TAG(J,I),J=1,2),ISOT(I)
     2,(ISCAT(J,I),J=1,2),OCCA(I),OCCUP(I),SITE(I),(XYZ(J,I),J=1,3),
     3(BETA(J,I),J=1,6),(CR(J,I),J=1,10),(DR(J,I),J=1,15),I=1,NATOM),
     4 LDEGN,(KDEGN(I),I=1,LDEGN),NT,TFPAR,XWIN,WAVE,UBMAT,(TFRAC(I)
     5,((TWINTM(J,K,I),J=1,3),K=1,3),JTWNMT(I),I=1,NT),BDF,WVLNTH
     6,HELPE,HELPF,IZACH,NIEXTA,NABSRN,KABSRF,ABSCF,JTFLG
 1110 READ(ISC2)(IBUFF(I),I=1,LBUFF)
      WRITE(ISC1)(IBUFF(I),I=1,LBUFF)
      IF(IBUFF(1))1120,1110,1120
 1120 REWIND ISC2
      IF(IFIRST.EQ.0.AND.MFSEL.EQ.0)GO TO 1130
      IK=MDEGN+1
      IJ=MDEGN+(MRANK*(MRANK+1))/2
      WRITE(ISC1)MDEGN,(A(J),J=1,MDEGN)
      WRITE(ISC1)(A(J),J=IK,IJ)
      WRITE(ISC1)((BBETA(J,I),J=1,6),SIGO(I),(SIGXYZ(J,I),J=1,3),
     1(SIGB(J,I),J=1,6),I=1,NATOM)
 1130 END FILE ISC1
      REWIND ISC1
C        END OF THIS ROUTINE IF NEW PARAMETERS ARE TO BE SAVED
      RETURN
C        END OF RUN FOR NEW PARAMETERS NOT SAVED
 1140 CALL EXIT
      END
