      PROGRAM POWDER
C***  INPUT COORDINATES ON CARDS OR TAPE1
C***  INPUT REFLEXIONS ON CARDS OR TAPE2
C***  OUTPUT COORDINATES ON TAPE3
C***  OUTPUT H,K,L,FCALC,A,B ON TAPE 11
C***
C***  MAXIMUM DIMENSIONS
C***        9 ATOM TYPES
C***       60 ATOMS
C***       99 VARIABLES
C***
C***  2 DATA SETS
C***
      CHARACTER*8 TITLE(10),AN
      REAL*8 AM,PD,DIAG,V,PDI,SIGSQ
      INTEGER H
      COMMON NC,NV,IRT,IPR,ISYM,NF,NT,NA,NS,NP,IC,IWARN,I10,
     1 WVL(2),SC(2),BO(2),SINT,DIDBO(2),SINANG,NDT,KDT,
     2 PRCP(6),DR(9,2),DI(9,2),BI(9),DIDB(9),SF(9,32,2),
     3 TS(3,48),IS(2,3,48),KI(253),P(253),DA(253),
     4 AN(60),ISF(60),ITF(60),AI(60),DIDA(60),
     5 X(3,60),DIDX(3,60)
      COMMON /MATRIX/AM(4950)
      DIMENSION V(99),DV(99),DIAG(99),PD(99)
      DIMENSION H(3),PW(3,2),PLP(3),DIDS(2)
      LUC=60  
      LUR=60
      OPEN(3,FILE='COORDS.SCR',STATUS='UNKNOWN')
      OPEN(11,FILE='FOURIER.SCR',STATUS='UNKNOWN')
      OPEN(60,FILE=' ')
      OPEN(61,FILE=' ')
      OPEN(10,FILE='DIFFS10.SCR',STATUS='UNKNOWN')
      REWIND 3
      REWIND 10
      REWIND 11
C***  READ TITLE AND CONTROL CARDS
      READ(60,1000)(TITLE(I),I=1,10)
 1000 FORMAT(10A8)
      READ(60,1010)NC,NV,INC,INR,IW,ILP,IRT,IPR,IDP,NDT
 1010 FORMAT(10I3)
      IF(INC.NE.0)  LUC=1
      IF(INR.NE.0)  LUR=2
      NM=(NV*(NV+1))/2
      WRITE(61,1030)(TITLE(I),I=1,10)
 1030 FORMAT(1H1,10A8)
      IF(IRT-1)101,102,103
  101 WRITE(61,1141)  
      GOTO 110
 1141 FORMAT(/'0INPUT COORDINATES, AND ATOMS VARIED WILL BE PRINTED'//
     1)
  102 WRITE(61,1142)  
      GOTO 110
 1142 FORMAT(/'0ALL CRYSTAL DATA, AND ATOMS VARIED WILL BE PRINTED'//)
  103 WRITE(61,1143)  
      GOTO 110
 1143 FORMAT(/'0FULL PRINTOUT OF INPUT DATA; PARAMETERS EACH CYCLE'//)
  110 CONTINUE
      I10=0
      IF(IPR.GT.1)I10=1
      IF(IPR.GT.1)IPR=IPR-2
      IF(IPR.NE.0) WRITE(61,1144)
 1144 FORMAT(/'0REFLEXION DATA WILL BE PRINTED AFTER LAST CYCLE'//)
      IF(I10.NE.0)WRITE(61,1147)
 1147 FORMAT(/'0TAPE 10 FOR PROGRAMSUITE DIFFUSE WILL BE WRITTEN'//)
      DP=IDP/10.0
      IF(IDP.NE.0.AND.NC.GT.1)  WRITE(61,1145)DP
 1145 FORMAT('0PARAMETER SHIFTS MULTIPLIED BY'F4.2)
C***  READ STRUCTURE CARD
      READ(60,1150)ISYM,NF,NT,NA,NS
 1150 FORMAT(5I3)
      IF(IW.EQ.0.AND.ILP.EQ.0)  GO TO 160
      READ(60,1200)PW,PLP
 1200 FORMAT(8F9.4,1F8.4)
      IF(PW(1,1).LE.0.0) PW(1,1)=0.4
      IF(PW(2,1).LE.0.0) PW(2,1)=10.0
      WRITE(61,1210)PW,PLP
 1210 FORMAT(/'0SUBROUTINE PARAMETERS'/1H0,4X,'WEIGHTING SCHEME SET 1',
     +3F12.4/' ',25X,'2',3F12.4,/1H0,4X,'LP FUNCTION'3F12.4)
  160 CONTINUE
      IF(IRT.LT.1)  GO TO 170
      WRITE(61,1040)NC,NV
 1040 FORMAT(/1H0,I3,' REFINEMENT CYCLES'/1H0,I3,' PARAMETERS VARIED')
      IF(INC.NE.0)  GO TO 10
      WRITE(61,1050)
 1050 FORMAT('0COORDINATES WILL BE READ FROM CARDS')
      GO TO 20
   10 WRITE(61,1060)
 1060 FORMAT('0COORDINATES WILL BE READ FROM FILE ON LOGICAL UNIT 1')
   20 CONTINUE
      IF(INR.NE.0)  GO TO 30
      WRITE(61,1070)
 1070 FORMAT('0REFLEXION DATA WILL BE READ FROM CARDS AND COPIED ON TO L
     1OGICAL UNIT 2')
      GO TO 40
   30 WRITE(61,1080)
 1080 FORMAT('0REFLEXION DATA WILL BE READ FROM FILE ON LOGICAL UNIT 2')
   40 CONTINUE
      IF(NDT.GT.1) WRITE(61,1081)NDT
 1081 FORMAT('0',I4,' DATA SETS')
      IF(IW-1)50,60,70
   50 WRITE(61,1090)
 1090 FORMAT('0WEIGHTS TO BE READ ON INPUT')
      GO TO 80
   60 WRITE(61,1100)
 1100 FORMAT('0WEIGHTS CALCULATED BY SUBROUTINE')
      GO TO 80
   70 WRITE(61,1110)
 1110 FORMAT('0INPUT OBSERVATIONAL WEIGHTS MULTIPLIED BY SUBROUTINE')
   80 CONTINUE
      IF(LUR.EQ.60.AND.IW.NE.0) WRITE(61,1111)
 1111 FORMAT('  AND WRITTEN TO TAPE2')
      IF(ILP.NE.0)  GO TO 90
      WRITE(61,1120)
 1120 FORMAT('0LP FACTORS READ ON INPUT')
      GO TO 100
   90 WRITE(61,1130)
 1130 FORMAT('0LP FACTORS CALCULATED BY SUBROUTINE')
  100 CONTINUE
      IF(LUR.EQ.60.AND.ILP.NE.0) WRITE(61,1111)
      GO TO(120,130,140,145)ISYM
  120 WRITE(61,1160)
 1160 FORMAT('0STRUCTURE IS CENTROSYMMETRIC')
      GO TO 150
  130 WRITE(61,1170)
 1170 FORMAT('0STRUCTURE IS NON-CENTROSYMMETRIC')
      GO TO 150
  140 WRITE(61,1180)
 1180 FORMAT('0STRUCTURE IS CENTROSYMMETRIC WITH ANOMALOUS DISPERSION')
      GO TO 150
  145 WRITE(61,1185)
 1185 FORMAT('0STRUCTURE IS NON-CENTROSYMMETRIC WITH ANOMALOUS DISPERSIO
     1N')
  150 CONTINUE
      WRITE(61,1190)NF,NT,NA,NS
 1190 FORMAT(1H0,I3,' SCATTERING FACTORS'/1H0,I3,' TEMPERATURE FACTORS'
     1 /1H0,I3,' ATOMS IN ASYMMETRIC UNIT'/1H0,I3,' SYMMETRY CARDS')
  170 CONTINUE
      IF(INC.NE.0) OPEN(1,FILE=' ')
      IF(INR.EQ.0) OPEN(2,FILE='DATA.SCR',STATUS='NEW')
      IF(INR.NE.0) OPEN(2,FILE=' ')
      NC=NC+1
C***  READ AND PUT OUT STRUCTURE DATA
      CALL READFX
      IF(IWARN.EQ.1)  GO TO 9999
C***  READ AND PUT OUT TRIAL PARAMETERS
      CALL READPR(LUC)
C***  START REFINEMENT CYCLES
      NIC=NC
      DO 400  IIC=1,NIC
      IC=IIC
C***  SET INITIAL COUNTS
      KDT=1
      NO=0  
      RDEN=0.0
      SIGSQ=0.0
      IF(IC.LT.NC)  GO TO 190
      IF(IPR.NE.0) WRITE(61,1240)(TITLE(I),I=1,10),IC
 1240 FORMAT(1H1,10A8,/'0CALCULATED INTENSITIES BASED ON PARAMETERS BEFO
     1RE CYCLE'I4,//'0   H   K   L     M    LP FACT   SINTHETA    I(CALC
     2)'11X,'I(CALC)    I(OBS)    DELTA      WEIGHT      WO     IO'//)
  190 CONTINUE
C***  CLEAR MATRIX AND VECTOR ARRAYS EXCEPT ON FINAL CYCLE
      IF(IC.EQ.NC)  GO TO 220
      DO 200  I=1,NV
  200 V(I)=0.0
      DO 210  I=1,NM
  210 AM(I)=0.0
  220 CONTINUE
C***  ENTRY POINT FOR EACH NEW REFLEXION GROUP
      READ(LUR,1220,END=260)ISENT,(H(J),J=1,3),IO,YO,W,PL,M,NR,WO
 1220 FORMAT(1I1,4I3,1F7.2,2F9.5,2I3,F6.2)
      IF(WO.EQ.0.0.AND.IO.NE.0) WO=1.0
      IF(NR.EQ.0)  NR=1
      IF(ISENT.EQ.0)  GO TO 230
      IF(LUR.EQ.60)  WRITE(2,1220) ISENT
      KDT=KDT+1
      IF(KDT.LE.NDT) GO TO 220
      IF(KDT.GT.NDT)  GO TO 260
  230 ISENT=0
      SINT=0.0
      DIDBO(1)=0.0
      DIDBO(2)=0.0
      DO 233  I=1,NT
  233 DIDB(I)=0.0
      DO 235  I=1,NA
      DIDA(I)=0.0
      DO 235  J=1,3
      DIDX(J,I)=0.0
  235 CONTINUE
      IF(IO.NE.0)  NO=NO+1
      RHO=PRCP(1)*H(1)**2+PRCP(2)*H(2)**2+PRCP(3)*H(3)**2+
     1 PRCP(4)*H(2)*H(3)+PRCP(5)*H(1)*H(3)+PRCP(6)*H(1)*H(2)
      SINANG=WVL(KDT)*SQRT(RHO)
      IF(ILP.NE.0)  PL=PLFACT(SINANG,PLP,KDT)
      IF(IW.NE.0)  W=WEIGHT(IW,PW,YO,WO,SINANG,KDT)
      IF(LUR.EQ.60)  WRITE(2,1220)ISENT,(H(J),J=1,3),IO,YO,W,PL,M,NR,WO
      CALL CALCIS(H,M,PL,RHO,IO)
      IF(NR.EQ.1)  GO TO 250
      DO 240  ICR=2,NR
      READ(LUR,1230)(H(J),J=1,3),PL,M
 1230 FORMAT(1X,3I3,19X,F9.5,I3,3X)
      RHO=PRCP(1)*H(1)**2+PRCP(2)*H(2)**2+PRCP(3)*H(3)**2+
     1 PRCP(4)*H(2)*H(3)+PRCP(5)*H(1)*H(3)+PRCP(6)*H(1)*H(2)
      SINANG=WVL(KDT)*SQRT(RHO)
      IF(ILP.NE.0)  PL=PLFACT(SINANG,PLP,KDT)
      IF(LUR.EQ.60)  WRITE(2,1230)(H(J),J=1,3),PL,M
      CALL CALCIS(H,M,PL,RHO,IO)
  240 CONTINUE
  250 CONTINUE
      DIDS(KDT)=SINT/SC(KDT)
      DY=YO-SINT
      IF(IC.EQ.NC.AND.IPR.NE.0)  WRITE(61,1250)SINT,YO,DY,W,WO,IO
 1250 FORMAT(1H+,58X,3F10.2,1F12.6,1F10.4,1I5/116X)
      IF(IC.EQ.NC.AND.I10.NE.0) WRITE(10,1250)SINT,YO,DY,W,WO,IO
C***  ACCUMULATE SUMS FOR AGREEMENT FACTORS
      IF(IO.EQ.0)  GO TO 255
      RDEN=RDEN+W*YO**2
      SIGSQ=SIGSQ+W*DY**2
C***  COLLECT REQUIRED DERIVATIVES,EXCEPT ON FINAL CYCLE
      CALL COPTOP(DA,DIDS,DIDBO,DIDB,DIDA,DIDX,NA,NT)
      IF(IC.EQ.NC)  GO TO 255
      J=1
      DO 270  K=1,NP
      IF(KI(K).NE.1)  GO TO 270
      DV(J)=DA(K)
      J=J+1
  270 CONTINUE
C***  STORE MATRIX AND VECTOR
      JK=NM
      DO 300  J=1,NV
      IF(DV(J).NE.0.0)  GO TO 280
      JK=JK-NV+J-1
      GO TO 300
  280 CONTINUE
      DO 290  K=J,NV
      AM(JK)=AM(JK)+W*DV(J)*DV(K)
      JK=JK-1
  290 CONTINUE
      V(J)=V(J)+W*DY*DV(J)
  300 CONTINUE
  255 CONTINUE
C***  READ NEW REFLEXION GROUP
      GO TO 220
  260 CONTINUE
      IF(LUR.NE.60) GO TO 261
      ENDFILE 2
      IW=0
      ILP=0
  261 LUR=2
      REWIND 2
C***  PUT OUT AGREEMENT FACTORS FOR THIS CYCLE
      RTSIG=SQRT(SIGSQ/(NO-NV))
      RW=SQRT(SIGSQ/RDEN)*0.50
      WRITE(61,1260)IC,NO,SIGSQ,RTSIG,RW
 1260 FORMAT('1AGREEMENT FACTORS BEFORE CYCLE'I3,//1H0I3,' OBSERVATIONS
     1 READ'/19H0SUM(W*(IO-IC)**2)= E12.4,/
     234H0SQRT((SUM(W*(IO-IC)**2)/(NO-NV))=  F10.4,/
     3 '0WEIGHTED RESIDUAL FACTOR='F9.4)
      IF(IC.EQ.NC)  GO TO 400
C***  TEST FOR ZERO DIAGONAL ELEMENTS IN MATRIX
      ISING=0
      INM=NM  
      INV=NV
      DO 320  I=1,NV
      IF(AM(INM).NE.0.0)  GO TO 310
      ISING=1
      WRITE(61,1270)I
 1270 FORMAT('0MATRIX HAS ZERO DIAGONAL ELEMENT CORRESPONDING TO PARAMET
     1ER'I3,' VARIED')
  310 INM=INM-INV
  320 INV=INV-1
      IF(ISING.EQ.1)  GO TO 9999
      CALL SMI(AM,NV,ISING,PD,DIAG)
      IF(ISING.EQ.0)  GO TO 330
      WRITE(61,1280)
 1280 FORMAT('0SINGULAR MATRIX,JOB THEREFORE TERMINATED')
      GO TO 9999
  330 CONTINUE
C***  START VECTOR MULTIPLICATION
      DO 380  I=1,NV
      PDI=0.0
      IJ=NM-I+1
      IJD=NV-1
      DO 370  J=1,NV
      PDI=PDI+AM(IJ)*V(J)
      IF(J-I)340,350,360
  340 IJ=IJ-IJD
      IJD=IJD-1
      GO TO 370
  350 DIAG(I)=AM(IJ)
  360 IJ=IJ-1
  370 CONTINUE
      IF(IDP.NE.0)  PDI=PDI*DP
      PD(I)=PDI
      SIGSQ=SIGSQ-PDI*V(I)
  380 CONTINUE
C***  OUTPUT PARAMETERS AND CHANGES
      WRITE(61,1285)IC,SIGSQ
 1285 FORMAT(/'0ESTIMATED AGREEMENT FACTOR AFTER CYCLE'I3,/
     1 19H0SUM(W*(IO-IC)**2)= E12.4)
      IF(SIGSQ.LT.0.0) SIGSQ=0.0
      RTSIG=SQRT(SIGSQ/(NO-NV))
      J=1  
      K=0
      WRITE(61,1290)IC
 1290 FORMAT(/'0PARAMETERS AFTER CYCLE'I3/
     1/1H0,25X,'OLD'10X,'CHANGE'11X,'NEW'11X,'ERROR')
      IDEG=0  
      IMIN=NT+4
      KS=0
      DO 410  I=1,NP
      IF(I.EQ.1)  WRITE(61,1292)
 1292 FORMAT(/'0SCALE FACTORS')
      IF(I.EQ.3)  WRITE(61,1293)
 1293 FORMAT(/'0OVERALL TEMPERATURE FACTORS')
      IF(I.EQ.5)  WRITE(61,1295)
 1295 FORMAT(/'0INDIVIDUAL TEMPERATURE FACTOR(S)')
      IMOD=MOD(I-NT-5,4)
      IF(I.LE.IMIN.OR.IMOD.NE.0)  GO TO 385
      K=K+1
      KS=KI(I)+KI(I+1)+KI(I+2)+KI(I+3)
      IF(KS.EQ.0.AND.IRT.LT.2)  GO TO 410
      WRITE(61,1300)AN(K),ITF(K)
 1300 FORMAT(/1H0,A8,4X,'TF TYPE 'I3)
  385 CONTINUE
      IF(KI(I).EQ.0)  GO TO 390
      POLD=P(I)
      P(I)=P(I)+PD(J)
      IF(I.LE.IMIN)  GO TO 375
      TEST=0.25
      IF(IMOD.EQ.0)  TEST=2.0
      IF(ABS(PD(J)).GT.TEST)  IDEG=1
  375 ERR=SQRT(DIAG(J))*RTSIG
      WRITE(61,1310)POLD,PD(J),P(I),ERR
 1310 FORMAT(1H ,15X,4F15.7)
      IF(IMOD.EQ.0.AND.I.GT.NT)  WRITE(61,1315)
 1315 FORMAT(1X)
      J=J+1
      GO TO 410
  390 CONTINUE
      IF(KS.EQ.0.AND.IRT.LT.2)  GO TO 410
      WRITE(61,1320)P(I),P(I)
 1320 FORMAT(1H ,2(15X,F15.7))
      IF(IMOD.EQ.0.AND.I.GT.NT)  WRITE(61,1315)
  410 CONTINUE
      CALL COPEXP(P,SC,BO,BI,AI,X,NA,NT)
C***  CALL RESETX AND COPTOP
      CALL RESETX(AI,X)
      CALL COPTOP(P,SC,BO,BI,AI,X,NA,NT)
C***  TEST TEMPERATURE FACTORS
      DO 430  I=1,NT
      TF=BI(I)
      IF(TF.LT.0.0)  WRITE(61,1330)I
 1330 FORMAT('0TEMPERATURE FACTOR NUMBER'I4,' IS NOW INVALID')
  430 CONTINUE
      IF(IDEG.EQ.1)  GO TO 405
  400 CONTINUE
  405 CONTINUE
      IF(NC.EQ.1)  GO TO 9999
      IF(IDEG.EQ.1)  GO TO 439
      WRITE(3,1365)SC(1),SC(2)
 1365 FORMAT(E12.6)
      WRITE(3,1360)BO(1),BO(2)
      WRITE(3,1360)(BI(I),I=1,NT)
 1360 FORMAT(F9.6)
      DO 420  I=1,NA
  420 WRITE(3,1370)AN(I),ISF(I),ITF(I),AI(I),(X(J,I),J=1,3)
 1370 FORMAT(A8,2I3,4X,4F9.6)
      ENDFILE 3
C***  PRINT CORRELATION MATRIX
  439 WRITE(61,1340)
 1340 FORMAT(/'0CORRELATION MATRIX'///)
      DO 440  I=1,NV
  440 DIAG(I)=1.0/SQRT(DIAG(I))
      IJ=NM
      DO 470  I=1,NV
      DO 450  J=1,NV
  450 V(J)=0.0
      DO 460  J=I,NV
      V(J)=AM(IJ)*DIAG(I)*DIAG(J)
  460 IJ=IJ-1
      WRITE(61,1350)I,(V(J),J=1,NV)
 1350 FORMAT(1H0,I3,20F6.2,/(4X,20F6.2))
  470 CONTINUE
 9999 CONTINUE
      REWIND 10
      REWIND 11
      END
