/.JANA4S LOGON
/SYSFILE SYSLST=C.JANA4S
/ER *
/EXEC $FOR1
COMOPT MSGLEVEL=E
COMOPT PAD
COMOPT OPT=2
COMOPT LIST=(SOURCE,XREF)
      PROGRAM JANA4S
      COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF,
     1              CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ
      COMMON/UNITS/ W,M40,M50
      COMMON/REFL/  YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708),
     1              DC(9708),DYP,IQ
      COMMON/IND/   H(24),K(24),L(24),T(24),HQ(24),KQ(24),LQ(24),HK2(24)
     1             ,HL2(24),KL2(24),FSMM(24),EPS(24),EPSM(24),MM,PHFP
      COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15)
      COMMON/CM4/   ATOM(100),AI(100),ISF(100),ITF(100),X(3,100),
     1              BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100)
     2             ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100),
     3              AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100),
     4              KFX(100),KFB(100)
      COMMON/SYMM/  IS(2,4,24),TS(3,24),TAU(24),FSM(24)
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      COMMON/VAZBY/ NVX,NX(10),NVB,NB(10),RMX(3,3,10),RMB(6,6,10),
     1              NVAI,NAI(2,10),NAIXB(10),SUMAI(10)
      COMMON/STALY/ PI2
      COMMON/MODUL/ QUI(3),QUIR(3)
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON/TEMPAR/SQRP(3),RHOM
      COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708),
     1       AM(75000)
      EQUIVALENCE (NPSS(10),NPSM)
      LOGICAL CSYM,NULOVA
      INTEGER W
      REAL K,L,KQ,LQ,KL2,KP,LP
      CHARACTER*1 ISPACE,ISTAR,INR,NSPEC,KSPEC
      CHARACTER*8 MOLNAME,ATOM
      CHARACTER*2 SYMBOL
      DATA ISPACE,ISTAR,INR/' ','*','#'/
      DIMENSION HH(3),HHP(3),RNUM(6),RZNUM(6),RDEN(6),RZDEN(6),WRNUM(6),
     1          WRZNUM(6),WRDEN(6),WRZDEN(6)
      CALL PRELIM
      EFPIP(7)=1.5*RR
      DO 9000IU=1,NCYKL+1
      CALL NEWPG
      IF(IU.EQ.NCYKL+1) GO TO 400
      WRITE(W,107) IU,TLUM
      GO TO 450
  400 WRITE(W,110)
  450 CALL CONTR(1,0)
      CALL CONTRM(1,0)
      IF(NMOL.GT.0) CALL SETAMP
      DO 500I=1,NPSM
  500 PS(I)=0.
      DO 700I=1,NAMM
  700 AM(I)=0.
      DO 800I=1,6
      RNUM(I)=0.
      RZNUM(I)=0.
      RDEN(I)=0.
      RZDEN(I)=0.
      WRNUM(I)=0.
      WRZNUM(I)=0.
      WRDEN(I)=0.
  800 WRZDEN(I)=0.
      NO=0
      NZ=0
      NVYH=0
 1000 READ(91,333) ISENT,IH1,IK1,IL1,IM1,YO,SIGYO,IQ,SINTHL,NXX
  333 FORMAT(I1,4I4,2F9.4,I4,F9.6,I4)
      IF(ISENT.NE.0) GO TO 8000
      NO=NO+1
      MM=IM1
      MMABS=IABS(MM)
      PHFP=-MMABS*(MMABS-1)
      IF(KIM.LE.-2) GO TO 1100
      IF(KIM.EQ.-1.AND.MM.NE.0) GO TO 1100
      IF(MMABS.EQ.KIM) GO TO 1100
      GO TO 1000
 1100 IF(NXX.GT.NC) GO TO 1000
      IF(SINTHL.LT.SNLMN.OR.SINTHL.GT.SNLMX) GO TO 1000
      MMABS=MMABS+2
      NZ=NZ+1
      HH(1)=IH1
      HH(2)=IK1
      HH(3)=IL1
      FMM=MM
      DO 1110I=1,3
 1110 HH(I)=HH(I)+FMM*(QUI(I)+QUIR(I))
      POM=1./SQRT(HH(1)**2+HH(2)**2+HH(3)**2)
      DO 2111I=1,3
      EFPIP(I)=POM*HH(I)
 2111 EFPIP(I+3)=EFPIP(I)
      SIGYO=SQRT(SIGYO**2+BLKOEF*YO**2)
      RHO=SINTHL**2*RHOM
      NULOVA=NXX.GT.0
      NSPEC=ISPACE
      IF(NULOVA) NSPEC=ISTAR
C
      DO 2500I=1,NS
      DO 2200J=1,3
 2200 HHP(J)=0.
      DO 2400J=1,3
      DO 2300M=1,2
      ISP=IS(M,J,I)
      ISPA=IABS(ISP)
      IF(ISP) 2230,2300,2250
 2230 HHP(ISPA)=HHP(ISPA)-HH(J)
      GO TO 2300
 2250 HHP(ISPA)=HHP(ISPA)+HH(J)
 2300 CONTINUE
 2400 CONTINUE
C
      EPS(I)=1.
      IF(IS(1,4,I).LT.0) EPS(I)=-1.
      IF(IS(1,4,I).EQ.0.AND.IS(2,4,I).LT.0) EPS(I)=-1.
      FSMM(I)=FSM(I)*FMM
      EPSM(I)=EPS(I)*FMM
      HP=HHP(1)
      KP=HHP(2)
      LP=HHP(3)
      H(I)=HP*PI2
      K(I)=KP*PI2
      L(I)=LP*PI2
      T(I)=(TS(1,I)*HH(1)+TS(2,I)*HH(2)+TS(3,I)*HH(3))*PI2
      HP=HP*SQRP(1)
      KP=KP*SQRP(2)
      LP=LP*SQRP(3)
      HQ(I)=HP**2
      KQ(I)=KP**2
      LQ(I)=LP**2
      HK2(I)=2.*HP*KP
      HL2(I)=2.*HP*LP
 2500 KL2(I)=2.*KP*LP
      PT=SINTHL/.05+1.
      IPT=PT
      FRACPT=PT-FLOAT(IPT)
      DO 2800I=1,NF
      FX(I)=FF(IPT,I)+(FF(IPT+1,I)-FF(IPT,I))*FRACPT+FF(29,I)
 2800 FY(I)=FF(31,I)/FX(I)
      CALL CALC
      IF(IU.GT.NCYKL) THEN
      WRITE(80) ISENT,IH1,IK1,IL1,IM1,YO,YC,A,B,SINTHL,NXX
      GO TO 2990
      ENDIF
      IF(NVAI.EQ.0) GO TO 2819
      DO 2818I=1,NVAI
      J1=NAI(1,I)*92-73
      J2=NAI(2,I)*92-73
      IF(NAIXB(I).LE.0) GO TO 2810
      DER(J1)=DER(J1)-DER(J2)
      DER(J2)=0.
      GO TO 2815
 2810 DO 2811J=10,18
      DER(J1+J)=DER(J1+J)-DER(J2+J)
 2811 DER(J2+J)=0.
 2815 IF(IABS(NAIXB(I)).EQ.3) GO TO 2818
      JK=9
      IF(IABS(NAIXB(I)).EQ.2) JK=3
      DO 2817J=1,JK
      DER(J1+J)=DER(J1+J)+DER(J2+J)
 2817 DER(J2+J)=0.
      IF(JK.NE.9) GO TO 2818
      DO 12817J=19,91
      DER(J1+J)=DER(J1+J)+DER(J2+J)
12817 DER(J2+J)=0.
 2818 CONTINUE
 2819 IF(NVX.LE.0.AND.NVB.LE.0) GO TO 2900
      IF(NVX.LE.0) GO TO 2850
      DO 2820I=1,NVX
      M=NX(I)*92-72
 2820 CALL ZMDER(DER(M),RMX(1,1,I),3)
 2850 IF(NVB.LE.0) GO TO 2900
      DO 2860I=1,NVB
      M=NB(I)*92-69
 2860 CALL ZMDER(DER(M),RMB(1,1,I),6)
 2900 IF(NMOL.EQ.0) GO TO 2990
      DO 2910I=NPMP,NPMK
 2910 DER(I)=0.
      J1=-73
      DO 2960I=1,NA
      J1=J1+92
      IF(KMODX(I).GE.0) GO TO 2960
      N=-KMODX(I)
      M=KMODXM(N)
      IF(M.LE.0) GO TO 2960
      J2=NPMP+49*(N-1)
      JX=J1+1
      JM=J1+19
      DO 2920JJ=1,M
      CALL CULTM(DER(JM),DURDX(1,JJ,N),DER(JX),1,6,3)
 2920 JM=JM+6
      JM=J2-1
      JX=J1+18
      JJ=6*M
      DO 2930J=1,JJ
 2930 DER(JM+J)=DER(JM+J)+DER(JX+J)
      JJ=J1+19
      JD=J2-J1+5
      DO 2950JX=JJ,JJ+6*(M-1),6
      ISP=JX+JD
      CALL CULTM(DER(JX),DURDR(1,I),DER(ISP),1,3,3)
 2950 CALL CULTM(DER(JX+3),DURDR(1,I),DER(ISP+3),1,3,3)
      JM=J2+48
      DER(JM)=DER(JM)+DER(J1+91)
 2960 CONTINUE
 2990 CALL DSETA(DER(19),NA)
      IF(NMOL.GT.0) CALL DSETM(DER(NPMP),NMOL)    
C
      GO TO (3000,3100,3200),IWQ
 3000 WT=1./SIGYO**2
      GO TO 3300
 3100 WT=1.
      SIGYO=1.
      GO TO 3300
 3200 WT=1./(YOMIN+YO+YO**2*YOMAX)
C
 3300 DY=YO-YC
      DYP=DY
      YOP=YO**2
      IF(IFSQ.NE.2) GO TO 3350
      DYP=YOP-YC**2
      WT=WT/(4.*YOP)
      YOP=4.*YOP**2
      IF(IU.GT.NCYKL) GO TO 3350
      DO 3310I=1,NPMK
 3310 DER(I)=DER(I)*2.*YC
 3350 WDY=SQRT(WT)*DYP
      WDYQ=WDY**2
      WYOQ=WT*YOP
      ADY=ABS(DY)
      AYO=ABS(YO)
      KSPEC=ISPACE
      IF(ABS(WDY).LE.VYH) GO TO 4500
      KSPEC=INR
      NVYH=NVYH+1
C
 4500 WRZNUM(1)=WRZNUM(1)+WDYQ
      RZNUM(1)=RZNUM(1)+ADY
      RZDEN(1)=RZDEN(1)+AYO
      WRZDEN(1)=WRZDEN(1)+WYOQ
      IF(NULOVA) GO TO 4700
      RNUM(1)=RNUM(1)+ADY
      RDEN(1)=RDEN(1)+AYO
      WRNUM(1)=WRNUM(1)+WDYQ
      WRDEN(1)=WRDEN(1)+WYOQ
 4700 WRZNUM(MMABS)=WRZNUM(MMABS)+WDYQ
      RZNUM(MMABS)=RZNUM(MMABS)+ADY
      RZDEN(MMABS)=RZDEN(MMABS)+AYO
      WRZDEN(MMABS)=WRZDEN(MMABS)+WYOQ
      IF(NULOVA) GO TO 4750
      RNUM(MMABS)=RNUM(MMABS)+ADY
      RDEN(MMABS)=RDEN(MMABS)+AYO
      WRNUM(MMABS)=WRNUM(MMABS)+WDYQ
      WRDEN(MMABS)=WRDEN(MMABS)+WYOQ
C
 4750 CALL CONTR(2,MM)
      CALL CONTRM(2,MMABS-1)
      IF(IU.LE.NCYKL) CALL SUMA
      IF(IABS(NOWR).EQ.2.AND.IU.NE.1.AND.IU.NE.NCYKL+1) GO TO 1000
      IF(NOWR) 4800,1000,4900
 4800 IF(KSPEC.EQ.ISPACE) GO TO 1000
      IF(MOD(NVYH,50).NE.1) GO TO 4950
      GO TO 4940
 4900 IF(MOD(NZ,50).NE.1) GO TO 4950
 4940 CALL NEWPG
      WRITE(W,100)
 4950 WRITE(W,101) IH1,IK1,IL1,IM1,YO,YC,A,B,DY,SIGYO,WDY,NO,NSPEC,KSPEC
     1            ,SINTHL,IQ,EXTKOR
C
      GO TO 1000
C
 8000 Q=WRZNUM(1)
      IF(IU.GT.NCYKL) WRITE(80) ISENT,IH1,IK1,IL1,IM1,YO,YC,A,B,SINTHL,
     1                NXX
      CALL NEWPG
      WRITE(W,108)
      DO 8500I=1,6
      IF(RZDEN(I).EQ.0.) GO TO 8500
      IF(I.GE.2) WRITE(W,109) I-2,I-2
      R=RNUM(I)/RDEN(I)
      RZ=RZNUM(I)/RZDEN(I)
      WRNUM(I)=SQRT(WRNUM(I))
      WRDEN(I)=SQRT(WRDEN(I))
      WRZNUM(I)=SQRT(WRZNUM(I))
      WRZDEN(I)=SQRT(WRZDEN(I))
      WR=WRNUM(I)/WRDEN(I)
      WRZ=WRZNUM(I)/WRZDEN(I)
      WRITE(W,102)
      WRITE(W,103) RZNUM(I),RZDEN(I),RZ
      WRITE(W,104) RNUM(I),RDEN(I),R
      WRITE(W,105) WRZNUM(I),WRZDEN(I),WRZ
      WRITE(W,106) WRNUM(I),WRDEN(I),WR
 8500 CONTINUE
      CALL CONTR(3,0)
      CALL CONTRM(3,0)
      IF(IU.GT.NCYKL) GO TO 9000
C
      CALL INVER
      CALL INTER
C
 9000 REWIND 91
      STOP
C
  100 FORMAT('0',3X,'H',3X,'K',3X,'L',3X,'M',5X,'YO',9X,'YC',10X,'A',10X
     1,'B',9X,'DY',8X,'SIGYO',7X,'WDY',4X,'NREF',4X,'SINTHL',4X,'IQ'
     2,6X,'EXT'//)
  101 FORMAT(1X,4I4,7F11.4,I5,1X,2A1,F9.5,I5,F9.5)
  102 FORMAT(1H0,57X,'NUMERATOR    DENOMINATOR       R')
  103 FORMAT(' R FACTOR INCLUDING UNOBSERVED REFLECTIONS',10X,2F15.3,
     1       F10.4)
  104 FORMAT(' R FACTOR  OMITTING UNOBSERVED REFLECTIONS',10X,2F15.3,
     1       F10.4)
  105 FORMAT(' WEIGHTED R FACTOR INCLUDING UNOBSERVED REFLECTIONS ',
     1       2F15.3,F10.4)
  115 FORMAT(' WEIGHTED R FACTOR INCLUDING UNOBSERVED REFLECTIONS ',
     1       2F15.3,2F10.4)
  106 FORMAT(' WEIGHTED R FACTOR  OMITTING UNOBSERVED REFLECTIONS ',
     1       2F15.3,F10.4)
  107 FORMAT(1X,//////////////////////'0',35X,53('*')/36X,'* REFINEMENT 
     1CYCLE #',I2,' WITH DUMPING FACTOR : ',F5.2,'  *'/36X,53('*'))
  108 FORMAT('0TABLE OF R FACTORS FOR ALL REFLECTIONS')
  109 FORMAT('0TABLE OF R FACTORS FOR REFLECTIONS WITH SATELLITE INDEX :
     1 ',    I1,' OR -',I1)
  110 FORMAT(1X,//////////////////////'0',51X,21('*')/52X,'* FINAL CALCU
     1LATION *'/1X,51X,21('*'))
      END
      BLOCK DATA
      COMMON/UNITS/ W,M40,M50
      COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10
      COMMON/TITLE/ TITLE(20),NAME(20),PAGE
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      COMMON/STALY/ PI2
      INTEGER W,TITLE,PAGE
      CHARACTER*12 L1,L2,L5,L7,L8,L9,L10
      DATA L1/'OCCUPATION  '/,L2/'X           ',    'Y           ',
     1        'Z           '/,L5/'U ISO       '/,
     2     L7/'RHO ISO     '/,L8/'G ISO       '/,L9/'PHASON      '/,
     3     L10/'A0          '/,
     4     W,M40,M50,PAGE/6,40,50,0/,TORAD/0.017453293/
      DATA PI2/6.2831853/
      END
      SUBROUTINE CONTR(ICONT,MM)
      DIMENSION SSQWDS(8),SSQWDF(9),SINMEZ(8,3),FMEZ(8,3),
     1          NLSP(8,3),NLSN(8,3),NLS(8),NLFP(9,3),NLFN(9,3),NLF(9),
     2          RNUMSP(8,3),RNUMSN(8,3),RNUMS(8),RDENS(8,3),RFAKS(8),
     3          RNUMFP(9,3),RNUMFN(9,3),RNUMF(9),RDENF(9,3),RFAKF(9),
     4          SPQWDS(8,3),SNQWDS(8,3),SPQWDF(9,3),SNQWDF(9,3)
      COMMON/UNITS/ W,M40,M50
      COMMON/REFL/  YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708),
     1              DC(9708),DYP,IQ
      DIMENSION SLIM(3),FLIM(3)
      LOGICAL NULOVA,CTENO
      INTEGER W
      CHARACTER*21 IDENT(3)
      DATA IDENT/'ALL REFLECTIONS      ','MAIN REFLECTIONS    ',
     1           'SATELLITE REFLECTIONS'/
      DATA CTENO/.FALSE./
      GO TO (5,15,50),ICONT
    5 DO 12J=1,3
      DO 10I=1,9
      SPQWDF(I,J)=0.0
      SNQWDF(I,J)=0.0
      RNUMFP(I,J)=0.0
      RNUMFN(I,J)=0.0
      RDENF(I,J)=0.0
      NLFP(I,J)=0
      NLFN(I,J)=0
      IF(I.EQ.9) GO TO 10
      SPQWDS(I,J)=0.0
      SNQWDS(I,J)=0.0
      RNUMSP(I,J)=0.0
      RNUMSN(I,J)=0.0
      RDENS(I,J)=0.0
      NLSP(I,J)=0
      NLSN(I,J)=0
   10 CONTINUE
   12 CONTINUE
      IF(CTENO) RETURN
      READ(M50,75)((SINMEZ(I,J),I=1,8),(FMEZ(I,J),I=1,8),J=1,3)
   75 FORMAT (8F9.5/8F9.1)
      DO 13J=1,3
      SLIM(J)=SINMEZ(8,J)
   13 FLIM(J)=FMEZ(8,J)
      CTENO=.TRUE.
      RETURN
   15 DO 30I=1,3
      IF(MM.EQ.0.AND.I.EQ.3) GO TO 30
      IF(MM.NE.0.AND.I.EQ.2) GO TO 30
      IF(YO.GT.FLIM(I).OR.SINTHL.GT.SLIM(I)) RETURN
      LF=9
      IF(NULOVA) GO TO 21
      LS=0
   17 LS=LS+1
      IF(SINTHL.GT.SINMEZ(LS,I)) GO TO 17
      LF=0
   20 LF=LF+1
      IF(YO.GT.FMEZ(LF,I)) GO TO 20
   21 IF(DY) 22,24,24
   22 RNUMFN(LF,I)=RNUMFN(LF,I)+DY
      NLFN(LF,I)=NLFN(LF,I)+1
      SNQWDF(LF,I)=SNQWDF(LF,I)+WDYQ
      IF(NULOVA) GO TO 25
      RNUMSN(LS,I)=RNUMSN(LS,I)+DY
      NLSN(LS,I)=NLSN(LS,I)+1
      SNQWDS(LS,I)=SNQWDS(LS,I)+WDYQ
      GO TO 25
   24 RNUMFP(LF,I)=RNUMFP(LF,I)+DY
      NLFP(LF,I)=NLFP(LF,I)+1
      SPQWDF(LF,I)=SPQWDF(LF,I)+WDYQ
      IF(NULOVA) GO TO 25
      RNUMSP(LS,I)=RNUMSP(LS,I)+DY
      NLSP(LS,I)=NLSP(LS,I)+1
      SPQWDS(LS,I)=SPQWDS(LS,I)+WDYQ
   25 RDENF(LF,I)=RDENF(LF,I)+YO
      IF(NULOVA) RETURN
      RDENS(LS,I)=RDENS(LS,I)+YO
   30 CONTINUE
      RETURN
   50 DO 70J=1,3
      SNUMP=0.0
      SNUMN=0.0
      FNUMP=0.0
      FNUMN=0.0
      SDEN=0.0
      FDEN=0.0
      LSP=0
      LSN=0
      LFP=0
      LFN=0
      DO 60I=1,9
      RNUMF(I)=RNUMFP(I,J)-RNUMFN(I,J)
      IF(RDENF(I,J).NE.0.0) RFAKF(I)=100.0*RNUMF(I)/RDENF(I,J)
      NLF(I)=NLFP(I,J)+NLFN(I,J)
      SSQWDF(I)=SPQWDF(I,J)+SNQWDF(I,J)
      IF(NLF(I).NE.0)  SSQWDF(I)=SSQWDF(I)/FLOAT(NLF(I))
      IF(I.EQ.9) GO TO 60
      RNUMS(I)=RNUMSP(I,J)-RNUMSN(I,J)
      IF(RDENS(I,J).NE.0.0) RFAKS(I)=100.0*RNUMS(I)/RDENS(I,J)
      SNUMP=SNUMP+RNUMSP(I,J)
      SNUMN=SNUMN+RNUMSN(I,J)
      FNUMP=FNUMP+RNUMFP(I,J)
      FNUMN=FNUMN+RNUMFN(I,J)
      SDEN=SDEN+RDENS(I,J)
      FDEN=FDEN+RDENF(I,J)
      LSP=LSP+NLSP(I,J)
      LSN=LSN+NLSN(I,J)
      LFP=LFP+NLFP(I,J)
      LFN=LFN+NLFN(I,J)
      NLS(I)=NLSP(I,J)+NLSN(I,J)
      SSQWDS(I)=SPQWDS(I,J)+SNQWDS(I,J)
      IF(NLS(I).NE.0) SSQWDS(I)=SSQWDS(I)/FLOAT(NLS(I))
   60 CONTINUE
      LS=LSP+LSN
      LF=LFP+LFN
      IF(LF.LE.0) GO TO 70
      CALL NEWPG
      WRITE(W,54) IDENT(J)
      WRITE(W,55)
      WRITE(W,56)(SINMEZ(I,J),I=1,8),(NLSP(I,J),I=1,8),(NLSN(I,J),I=1,8)
     1          ,NLS,SSQWDS,(SPQWDS(I,J),I=1,8),(SNQWDS(I,J),I=1,8),
     2           (RNUMSP(I,J),I=1,8),(RNUMSN(I,J),I=1,8),RNUMS,
     3           (RDENS(I,J),I=1,8),RFAKS
      WRITE(W,57)
      WRITE(W,59)(FMEZ(I,J),I=1,8),(NLFP(I,J),I=1,9),(NLFN(I,J),I=1,9)
     1          ,NLF,SSQWDF,(SPQWDF(I,J),I=1,9),(SNQWDF(I,J),I=1,9),
     2           (RNUMFP(I,J),I=1,9),(RNUMFN(I,J),I=1,9),RNUMF,
     3           (RDENF(I,J),I=1,9),RFAKF
      SNUM=SNUMP-SNUMN
      FNUM=FNUMP-FNUMN
      RFS=100.0*SNUM/SDEN
      RFF=100.0*FNUM/FDEN
      WRITE(W,58) LSP,LFP,LSN,LFN,LS,LF,SNUMP,FNUMP,SNUMN,FNUMN,
     1                   SNUM,FNUM,SDEN,FDEN,RFS,RFF
   70 CONTINUE
   54 FORMAT('0STATISTICS AS A FUNCTION OF SIN(TH)/LAMBDA AND STRUCTURE 
     1FACTORS FOR ',A21)
   55 FORMAT(///1X,'SIN(TH)/LAMBDA')
   56 FORMAT(1H+,17X,'LIMITS',5X,8F11.5//18X,'NUMBER +',3X,8I11/25X,'-',
     13X,8I11/18X,'TOGETHER',3X,8I11//18X,'AVERAGE WDY',8F11.4/18X,
     2'SSQWD +',4X,8F11.4/24X,'-',4X,8F11.4///18X,'NUMERATOR +',8F11.1/
     328X,'-',8F11.1//18X,'TOGETHER',3X,8F11.1/18X,'DENOMINATOR',8F11.1
     4//18X,'R FACTOR',3X,8F11.2)
   57 FORMAT(///' STRUCT. FACTORS')
   59 FORMAT(1H+,17X,'LIMITS',5X,8F11.1,6X,'UNOBS'//18X,'NUMBER +',3X,
     19I11/25X,'-',3X,9I11/18X,'TOGETHER',3X,9I11//18X,'AVERAGE WDY',
     29F11.4/18X,'SSQWD +',4X,9F11.4/24X,'-',4X,9F11.4///18X,'NUMERATOR 
     3+',9F11.1/28X,'-',9F11.1//18X,'TOGETHER',3X,9F11.1/18X,'DENOMINATO
     4R',9F11.1//18X,'R FACTOR',3X,9F11.2)
   58 FORMAT(///32X,16HFINAL CHECK     ,15X,14HSIN(TH)/LAMBDA,5X,'STRUCT
     1URE FACTORS'//47X,'NUMBER +',9X,I9,13X,I9/54X,1H-,9X,I9,13X,I9/47X
     2,'TOGETHER',9X,I9,13X,I9//47X,'NUMERATOR +',5X,F11.1,9X,F13.1/57X,
     3'-',5X,F11.1,9X,F13.1//47X,'TOGETHER',8X,F11.1,9X,F13.1/
     447X,'DENOMINATOR',5X,F11.1,9X,F13.1///47X,'R-FACTOR',8X,F11.2,9X,
     5F13.2)
      RETURN
      END
      SUBROUTINE CONTRM(ICONT,M)
      DIMENSION SSQWD(5),SNQWD(5),SPQWD(5),NLP(5),NLN(5),NL(5),RNUMP(5),
     1          RNUMN(5),RNUM(5),RDEN(5),RFAK(5)
      COMMON/UNITS/ W,M40,M50
      COMMON/REFL/  YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708),
     1              DC(9708),DYP,IQ
      INTEGER W
      LOGICAL NULOVA
      GO TO (1000,2000,3000),ICONT
 1000 DO 1500I=1,5
      SPQWD(I)=0.
      SNQWD(I)=0.
      RNUMP(I)=0.
      RNUMN(I)=0.
      RDEN(I)=0.
      NLP(I)=0
      NLN(I)=0.
 1500 CONTINUE
      RETURN
 2000 IF(NULOVA) RETURN
      IF(DY) 2100,2200,2200
 2100 RNUMN(M)=RNUMN(M)+DY
      NLN(M)=NLN(M)+1
      SNQWD(M)=SNQWD(M)+WDYQ
      GO TO 2300
 2200 RNUMP(M)=RNUMP(M)+DY
      NLP(M)=NLP(M)+1
      SPQWD(M)=SPQWD(M)+WDYQ
 2300 RDEN(M)=RDEN(M)+YO
      RETURN
 3000 DO 3100I=1,5
      RNUM(I)=RNUMP(I)-RNUMN(I)
      IF(RDEN(I).NE.0.) RFAK(I)=100.*RNUM(I)/RDEN(I)
      NL(I)=NLN(I)+NLP(I)
      SSQWD(I)=SNQWD(I)+SPQWD(I)
      IF(NL(I).NE.0) SSQWD(I)=SSQWD(I)/FLOAT(NL(I))
 3100 CONTINUE
      CALL NEWPG
      WRITE(W,100)(I,-I,I=1,4)
      WRITE(W,101) NLP,NLN,NL,SSQWD,SPQWD,SNQWD,RNUMP,RNUMN,RNUM,RDEN,
     1             RFAK
      RETURN
  100 FORMAT('0STATISTICS AS A FUNCTION OF SATELLITE INDEX'//'0SATELLITE
     1 INDEX  :',22X,'0',10X,4(I2,'  OR  ',I2,6X)//)
  101 FORMAT(19X,'NUMBER +',5I16/26X,'-',5I16/19X,'TOGETHER',5I16//19X,
     1       'AVERAGE WDY',5F16.4/19X,'SSQWD +',4X,5F16.4/25X,'-',4X,
     2       5F16.4///19X,'NUMERATOR +',5F16.1/29X,'-',5F16.1//19X,
     3       'TOGETHER   ',5F16.1/19X,'DENOMINATOR',5F16.1//19X,'R FACTO
     4R'    ,3X,5F16.2)
      END
      SUBROUTINE PRELIM
      COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF,
     1              CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ
      COMMON/UNITS/ W,M40,M50
      COMMON/CM4/   ATOM(100),AI(100),ISF(100),ITF(100),X(3,100),
     1              BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100)
     2             ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100),
     3              AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100),
     4              KFX(100),KFB(100)
      COMMON/CSM4/  SAI(100),SX(3,100),SB(6,100),SUX(3,4,100),
     1              SUY(3,4,100),SPHF(100),SA0(100),SAX(4,100),
     2              SAY(4,100),SBX(6,4,100),SBY(6,4,100)
      COMMON/TITLE/ TITLE(20),NAME(20),PAGE
      COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15)
      COMMON/SYMM/  IS(2,4,24),TS(3,24),TAU(24),FSM(24)
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      COMMON/VAZBY/ NVX,NX(10),NVB,NB(10),RMX(3,3,10),RMB(6,6,10),
     1              NVAI,NAI(2,10),NAIXB(10),SUMAI(10)
      COMMON/MODUL/ QUI(3),QUIR(3)
      COMMON/STALY/ PI2
      COMMON/METRICS/ G(9),GI(9),TR(9),VOL
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON/COVC/  ICOV
      COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10
      COMMON/TEMPAR/SQRP(3),RHOM
      COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708),
     1       AM(75000)
      DIMENSION RMT(3,10),URCP(6),RCP(6),PRCP(6)
      LOGICAL CSYM
      INTEGER W,TITLE,PAGE
      CHARACTER*8 ATOM,AT,ATO,MOLNAME
      CHARACTER*1 IC(7),ICENT,IQS(11),IQUIR,SMB(2),PLUS,MINUS,MEZERA,
     1            SMBT(5)
      CHARACTER*12 L1,L2,L5,L7,L8,L9,L10
      EQUIVALENCE (A,RCP(1)),(B,RCP(2)),(C,RCP(3)),(CSA,RCP(4)),
     1            (CSB,RCP(5)),(CSG,RCP(6)) 
      DATA IC/'P','A','B','C','I','R','F'/,IQS/'P','A','B','C','L','M',
     1        'N','U','V','W','R'/PLUS,MINUS,MEZERA,SMBT/'+','-',' ',
     2        '1','S','T','Q','H'/,TPISQ,EPISQ/19.7392088,78.95683521/
C
      READ(M50,100) NAME
      READ(M50,129) A,B,C,ALFA,BETT,GAMA,FLAME
      CSA=COS(TORAD*ALFA)
      CSB=COS(TORAD*BETT)
      CSG=COS(TORAD*GAMA)
      SNB=SQRT(1.-CSB**2)
      SNG=SQRT(1.-CSG**2)
      POM=(CSA-CSB*CSG)/SNG
      TR(1)=A
      TR(2)=0.
      TR(3)=0.
      TR(4)=B*CSG
      TR(5)=B*SNG
      TR(6)=0.
      TR(7)=C*CSB
      TR(8)=C*POM
      TR(9)=C*SNB*SQRT(1.-(POM/SNB)**2)
      VOL=A*B*C*SQRT(1.-CSA**2-CSB**2-CSG**2+2.*CSA*CSB*CSG)
      VOLUME=VOL**2*12.593E-4/FLAME**3
      G(1)=A**2
      G(2)=A*B*CSG
      G(3)=A*C*CSB
      G(4)=G(2)
      G(5)=B**2
      G(6)=B*C*CSA
      G(7)=G(3)
      G(8)=G(6)
      G(9)=C**2
      CALL SMI3(G,GI,CENTR)
      CALL RECIP(RCP,RCP)
      DO 250I=1,3
  250 PRCP(I)=.25*RCP(I)**2
      PRCP(4)=.25*RCP(1)*RCP(2)*RCP(6)
      PRCP(5)=.25*RCP(1)*RCP(3)*RCP(5)
      PRCP(6)=.25*RCP(2)*RCP(3)*RCP(4)
      DO 300I=1,3
  300 URCP(I)=TPISQ*RCP(I)**2
      URCP(4)=TPISQ*RCP(1)*RCP(2)
      URCP(5)=TPISQ*RCP(1)*RCP(3)
      URCP(6)=TPISQ*RCP(2)*RCP(3)
      RHOM=EPISQ
      DO 301I=1,3
  301 SQRP(I)=SQRT(URCP(I))
C
      READ(M50,113) QUI
      READ(M50,101) NS,NI,ICENT,IQUIR
      CSYM=NI.EQ.1
      DO 1000NC=1,7
      IF(ICENT.EQ.IC(NC)) GO TO 1100
 1000 CONTINUE
      NC=1
 1100 GO TO (1200,1210,1210,1210,1210,1220,1230),NC
 1200 CENTR=1.
      GO TO 1300
 1210 CENTR=2.
      GO TO 1300
 1220 CENTR=3.
      GO TO 1300
 1230 CENTR=4.
 1300 CENTR=CENTR*FLOAT(3-NI)
      DO 1310I=1,11
      IF(IQUIR.EQ.IQS(I)) GO TO 1320
 1310 CONTINUE
      I=1
 1320 IF(I.GE.8) GO TO 1330
      DO 1321J=1,3
 1321 QUIR(J)=0.
      IF(I.EQ.1) GO TO 1350
      I=I-1
      P=.5
      IF(I.LT.4) GO TO 1322
      I=I-3
      P=1.
 1322 QUIR(I)=P
      GO TO 1350
 1330 IF(I.EQ.11) GO TO 1340
      DO 1335J=1,3
 1335 QUIR(J)=.5
      I=I-7
      QUIR(I)=0.
      GO TO 1350
 1340 QUIR(1)=.333333333
      QUIR(2)=.333333333
      QUIR(3)=0.
C
 1350 DO 1360I=1,NS
      READ(M50,102)(TS(J,I),(IS(K,J,I),K=1,2),J=1,3),SMB,
     1                      (IS(K,4,I),K=1,2) 
      P=1.
      J=2
      IF(SMB(1).EQ.MEZERA.OR.SMB(1).EQ.PLUS) GO TO 1352
      IF(SMB(1).EQ.MINUS) GO TO 1351
      J=1
      GO TO 1352
 1351 P=-1.
 1352 DO 1353K=1,5
      IF(SMB(J).EQ.SMBT(K)) GO TO 1354
 1353 CONTINUE
      K=1
 1354 IF(K.EQ.1) P=0.
      IF(K.EQ.5) K=6
      TAU(I)=P/FLOAT(K)
 1360 FSM(I)=PI2*(TAU(I)-QUI(1)*TS(1,I)-QUI(2)*TS(2,I)-QUI(3)*TS(3,I))
C
      READ(M50,100)
      READ(M50,100)
      READ(M50,100)
      READ(M50,100)
      READ(M50,101) NF
      READ(M50,103)((FF(J,I),J=1,32),I=1,NF)
      READ(M50,103)(AW(I),I=1,NF)
      CALL CONTR(1,0)
      CALL NAJDI(M50)
      READ(M50,100) TITLE
      CALL NEWPG
      WRITE(W,108) ICENT
      IF(CSYM) GO TO 1400
      WRITE(W,105)
      GO TO 1450
 1400 WRITE(W,106)
 1450 WRITE(W,107)
      WRITE(W,109)((TS(J,I),(IS(K,J,I),K=1,2),J=1,3),TAU(I),
     1                      (IS(K,4,I),K=1,2),I=1,NS)
      WRITE(W,156) QUI
      WRITE(W,157) QUIR,IQUIR
      CALL NEWPG
      WRITE(W,110)
      WRITE(W,111)(I,(FF(J,I),J=1,32),I=1,NF)
      READ(M50,112) IWQ,NC,NOWR,ICOV,LITE,IFSQ
      READ(M50,112) IEXT,ITYPEX,IDISTR,IMOSAIC,IAB
      READ(M50,112) KIM,KIW
      IF(KIW.NE.1) GO TO 1470
      DO 1460I=1,NF
 1460 AW(I)=1.
 1470 IF(IEXT.GT.1) ITYPEX=IEXT-1
      IF(IEXT.LE.0) IAB=0
      IDISTR=MIN0(IDISTR,2)
      IDISTR=MAX0(IDISTR,1)
      IMOSAIC=MIN0(IMOSAIC,1)
      IMOSAIC=MAX0(IMOSAIC,0)
      READ(M50,113) YOMIN,YOMAX,VYH,BLKOEF,RR,UHMON,SNLMN,SNLMX
      IF(SNLMX.LE.0.) SNLMX=10.
      IF(IAB.NE.0) READ(M50,130) AA,BA,DADMI
      CALL NEWPG
      IF(IFSQ.NE.2) WRITE(W,180)
      IF(IFSQ.EQ.2) WRITE(W,181)
      IWQ=IWQ+1
      GO TO (1500,1510,1520),IWQ
 1500 WRITE(W,114) BLKOEF
      BLKOEF=(BLKOEF*.01)**2
      GO TO 2000
 1510 WRITE(W,115)
      GO TO 2000
 1520 WRITE(W,116) YOMIN,YOMAX
      YOMIN=2.*YOMIN
      YOMAX=2./YOMAX
 2000 WRITE(W,117) VYH
      WRITE(W,118) NC
      IF(KIM.EQ.0)  WRITE(W,158)
      IF(KIM.EQ.-1) WRITE(W,159)
      IF(KIM.GT.0)  WRITE(W,160) KIM,-KIM
      WRITE(W,150) SNLMN,SNLMX
      IF(NOWR) 12000,22000,32000
12000 WRITE(W,119)
      GO TO 32000
22000 WRITE(W,173)
      GO TO 52000
32000 I=ABS(NOWR)
52000 IF(I.EQ.1) WRITE(W,174)
      IF(I.EQ.2) WRITE(W,175)  
      IF(ICOV.NE.0) WRITE(W,171)
      READ(M40,112) NA,NMOL
      IF(NMOL.LT.0) NMOL=0
      IF(NMOL.EQ.0) GO TO 2004
      IF(KIW-1) 2001,2002,2003
 2001 WRITE(W,165)
      GO TO 2004
 2002 WRITE(W,166)
      WRITE(W,167)
      GO TO 2004
 2003 WRITE(W,166)
      WRITE(W,168)
 2004 READ(M40,131) SC,(KI(M),M=1,6)
      READ(M40,131) TO
      READ(M40,131)(EC(M),M=1,6),(KI(M),M=7,12),(EC(M),
     1                          M=7,12),(KI(M),M=13,18)
      L=18
      DO 2010I=1,NA
      K=L+1
      L=L+10
      READ(M40,121) ATOM(I),ISF(I),KFB(I),KFX(I),ITF(I),KMODX(I),
     1              KMODS(I),AI(I),(X(J,I),J=1,3),(BETA(J,I),J=1,6),
     2              (KI(M),M=K,L)
      KMODB(I)=IABS(KMODX(I)/10)
      M=MOD(IABS(KMODX(I)),10)
      IF(KMODX(I).GE.0) THEN
      KMODX(I)=M
      ELSE
      KMODX(I)=-M
      ENDIF
      QCNTM(I)=-PI2*((QUI(1)+QUIR(1))*X(1,I)+(QUI(2)+QUIR(2))*X(2,I)+
     1               (QUI(3)+QUIR(3))*X(3,I))
      K=L+1
      L=L+1
      A0(I)=1.
      KI(L)=0
      IF(KMODS(I).GT.0) READ(M40,151) A0(I),KI(L)
      DO 2005J=1,4
      K=L+1
      L=L+2
      IF(J.GT.KMODS(I)) THEN
      AX(J,I)=0.
      AY(J,I)=0.
      KI(K)=0
      KI(L)=0
      ELSE
      READ(M40,182) AX(J,I),AY(J,I),KI(K),KI(L)
      ENDIF
 2005 CONTINUE
      DO 2009J=1,4
      K=L+1
      L=L+6
      IF(J.GT.KMODX(I)) THEN
      DO 2006M=1,3
      UX(M,J,I)=0.
 2006 UY(M,J,I)=0.
      DO 2007M=K,L
 2007 KI(M)=0
      ELSE
      READ(M40,131)(UX(M,J,I),M=1,3),(UY(M,J,I),M=1,3),(KI(M),M=K,L)
      ENDIF
 2009 CONTINUE
      DO 2008J=1,4
      K=L+1
      L=L+12
      IF(J.GT.KMODB(I)) THEN
      DO 12008M=1,6
      BX(M,J,I)=0.
12008 BY(M,J,I)=0.
      DO 22008M=K,L
22008 KI(M)=0
      ELSE
      READ(M40,131)(BX(M,J,I),M=1,6),(KI(M),M=K,K+5),
     1             (BY(M,J,I),M=1,6),(KI(M),M=K+6,L)
      ENDIF
 2008 CONTINUE
      K=L+1
      L=L+1
      PHF(I)=0.
      KI(L)=0
      IF(KMODX(I)+KMODS(I)+KMODB(I).GT.0) READ(M40,151) PHF(I),KI(L)
 2010 CONTINUE
      IF(NMOL.EQ.0) GO TO 2017
      DO 2016I=1,NMOL
      READ(M40,163) MOLNAME(I),NPG(I),KMODXM(I),(XM(J,I),J=1,3)
      NPG(I)=MAX0(1,NPG(I))
      KPG(1,I)=1
      IF(NPG(I).GT.1) READ(M40,112)(KPG(J,I),J=2,NPG(I))
      DO 2014J=1,4
      K=L+1
      L=L+6
      IF(J.LE.KMODXM(I)) GO TO 2013
      DO 2011M=1,3
      UTX(M,J,I)=0.
 2011 UTY(M,J,I)=0.
      DO 2012M=K,L
 2012 KI(M)=0
      GO TO 2014
 2013 READ(M40,131)(UTX(M,J,I),M=1,3),(UTY(M,J,I),M=1,3),(KI(M),M=K,L)
 2014 CONTINUE
      DO 12014J=1,4
      K=L+1
      L=L+6
      IF(J.LE.KMODXM(I)) GO TO 12013
      DO 12011M=1,3
      URX(M,J,I)=0.
12011 URY(M,J,I)=0.
      DO 12012M=K,L
12012 KI(M)=0
      GO TO 12014
12013 READ(M40,131)(URX(M,J,I),M=1,3),(URY(M,J,I),M=1,3),(KI(M),M=K,L)
12014 CONTINUE
      K=L+1
      L=L+1
      PHFM(I)=0.
      KI(L)=0
      IF(KMODXM(I).GT.0) READ(M40,151) PHFM(I),KI(L)
      CALL CENTER(I)
      QCNTMM=-PI2*((QUI(1)+QUIR(1))*XM(1,I)+(QUI(2)+QUIR(2))*XM(2,I)+
     1             (QUI(3)+QUIR(3))*XM(3,I))
      DO 62015M=1,NA
      IF(KMODX(M).NE.-I) GO TO 62015
      QCNTM(M)=QCNTMM
62015 CONTINUE
 2016 CONTINUE
 2017 DO 2019I=1,NA
      READ(M40,120,END=2020,ERR=2020) SAI(I),(SX(J,I),J=1,3),
     1                                (SB(J,I),J=1,6)
      M=KMODS(I)
      IF(M.LE.0) GO TO 22017
      READ(M40,182,END=2020,ERR=2020) SA0(I)
      DO 12017J=1,M
12017 READ(M40,182,END=2020,ERR=2020) SAX(J,I),SAY(J,I)
22017 DO 2018J=1,4
      IF(J.GT.KMODX(I)) GO TO 32017
 2018 READ(M40,131,END=2020,ERR=2020)(SUX(K,J,I),K=1,3),
     1                               (SUY(K,J,I),K=1,3)
32017 DO 12018J=1,4
      IF(J.GT.KMODB(I)) GO TO 22018
      READ(M40,131,END=2020,ERR=2020)(SBX(K,J,I),K=1,6)
12018 READ(M40,131,END=2020,ERR=2020)(SBY(K,J,I),K=1,6)
22018 IF(KMODS(I)+KMODX(I)+KMODB(I).GT.0) READ(M40,131,END=2020,
     1                       ERR=2020) SPHF(I)
 2019 CONTINUE
      IF(NMOL.EQ.0) GO TO 2050
      DO 32019I=1,NMOL
      READ(M40,120,END=2020,ERR=2020)
      K=KMODXM(I)
      IF(K.LE.0) GO TO 32019
      DO 12019J=1,K
12019 READ(M40,131,END=2020,ERR=2020)(SUTX(M,J,I),M=1,3),
     1                               (SUTY(M,J,I),M=1,3)
      DO 22019J=1,K
22019 READ(M40,131,END=2020,ERR=2020)(SURX(M,J,I),M=1,3),
     1                               (SURY(M,J,I),M=1,3)
      READ(M40,131,END=2020,ERR=2020) SPHFM(I)
32019 CONTINUE
      GO TO 2050
 2020 DO 2030I=1,NA
      SAI(I)=0.
      DO 2021J=1,3
 2021 SX(J,I)=0.
      DO 2022J=1,6
 2022 SB(J,I)=0.
      DO 2024K=1,3
      DO 2023J=1,4
      SUX(K,J,I)=0.
 2023 SUY(K,J,I)=0.
 2024 CONTINUE
      DO 2026K=1,6
      DO 2025J=1,4
      SBX(K,J,I)=0.
 2025 SBY(K,J,I)=0.
 2026 CONTINUE
      SPHF(I)=0.
 2030 CONTINUE
      IF(NMOL.EQ.0) GO TO 2050
      DO 2040I=1,NMOL
      DO 2032K=1,3
      DO 2031J=1,4
      SUTX(K,J,I)=0.
      SUTY(K,J,I)=0.
      SURX(K,J,I)=0.
 2031 SURY(K,J,I)=0.
 2032 CONTINUE
      SPHFM(I)=0.
 2040 CONTINUE
 2050 IF(TO.EQ.0.) GO TO 2200
C
      DO 2100I=1,NA
      IF(ITF(I).EQ.1) BETA(1,I)=BETA(1,I)+TO
 2100 CONTINUE
 2200 F=0.
C
      DO 2300I=1,NA
      J=ISF(I)
 2300 F=F+FF(1,J)*AI(I)
      F=F*FLOAT(NS)*CENTR
C
      WRITE(W,122) F
      CALL NEWPG
      WRITE(W,123)
      IF(LITE.LT.1) GO TO 2350
      WRITE(W,176)
      IF(LITE.EQ.1) WRITE(W,177)
      GO TO 2360
 2350 WRITE(W,178)
      IF(LITE.LE.-1) WRITE(W,179)
 2360 WRITE(W,162) SC,(KI(I),I=1,6)
      J=1
      L=18
C
      DO 2450I=1,NA
      J=J+1
      IF(MOD(J,50).NE.1) GO TO 2400
      CALL NEWPG
      WRITE(W,123)
 2400 N=IABS(KMODX(I))+10*KMODB(I)
      IF(KMODX(I).LT.0) N=-N
      WRITE(W,125) ATOM(I),ISF(I),KFB(I),KFX(I),ITF(I),N,KMODS(I),AI(I),
     1             (X(M,I),M=1,3)
      J=J+1
      IF(MOD(J,50).NE.1) GO TO 2410
      CALL NEWPG
      WRITE(W,123)
 2410 K=L+1
      L=L+10
      WRITE(W,161) (BETA(M,I),M=1,6),(KI(M),M=K,L)
      IF(LITE.EQ.0.OR.LITE.GE.2) GO TO 2418
      IF(LITE.LT.0) GO TO 2415
      IF(ITF(I).NE.2) GO TO 2412
      DO 2411M=1,6
      DO 12410J=1,KMODB(I)
      BX(M,J,I)=BX(M,J,I)/URCP(M)
12410 BY(M,J,I)=BY(M,J,I)/URCP(M)
 2411 BETA(M,I)=BETA(M,I)/URCP(M)
      GO TO 2418
 2412 BETA(1,I)=BETA(1,I)/EPISQ
      GO TO 2418
 2415 IF(ITF(I).NE.2) GO TO 2417
      DO 2416M=1,6
      DO 12415J=1,KMODB(I)
      BX(M,J,I)=BX(M,J,I)*URCP(M)
12415 BY(M,J,I)=BY(M,J,I)*URCP(M)
 2416 BETA(M,I)=BETA(M,I)*URCP(M)
      GO TO 2418
 2417 BETA(1,I)=BETA(1,I)*EPISQ
 2418 K=L+1
      L=L+1
      IF(KMODS(I).GT.0) WRITE(W,124) A0(I),KI(L)
      DO 2419J=1,4
      K=L+1
      L=L+2
      IF(J.GT.KMODS(I)) GO TO 2419
      WRITE(W,183) AX(J,I),AY(J,I),KI(K),KI(L)
 2419 CONTINUE
      DO 2430N=1,4
      K=L+1
      L=L+6
      IF(N.GT.KMODX(I)) GO TO 2430
      J=J+1
      IF(MOD(J,50).EQ.1) THEN
      CALL NEWPG
      WRITE(W,123)
      ENDIF
      WRITE(W,162)(UX(M,N,I),M=1,3),(UY(M,N,I),M=1,3),(KI(M),M=K,L)
 2430 CONTINUE
      DO 2432N=1,4
      K=L+1
      L=L+12
      IF(N.GT.KMODB(I)) GO TO 2432
      J=J+1
      IF(MOD(J,50).EQ.1) THEN
      CALL NEWPG
      WRITE(W,123)
      ENDIF
      WRITE(W,162)(BX(M,N,I),M=1,6),(KI(M),M=K,K+5)
      J=J+1
      IF(MOD(J,50).EQ.1) THEN
      CALL NEWPG
      WRITE(W,123)
      ENDIF
      WRITE(W,162)(BY(M,N,I),M=1,6),(KI(M),M=K+6,L)
 2432 CONTINUE
      K=L+1
      L=L+1
      IF(KMODS(I)+KMODX(I)+KMODB(I).LE.0) GO TO 2450
      J=J+1
      IF(MOD(J,50).EQ.1) THEN
      CALL NEWPG
      WRITE(W,123)
      ENDIF
      WRITE(W,124) PHF(I),KI(K)
 2450 CONTINUE
C
      IF(NMOL.EQ.0) GO TO 2500
      DO 2499I=1,NMOL
      J=J+1
      IF(MOD(J,50).NE.1) GO TO 2451
      CALL NEWPG
      WRITE(W,123)
 2451 WRITE(W,164) MOLNAME(I),NPG(I),KMODXM(I),(XM(J,I),J=1,3)
      DO 2460N=1,4
      K=L+1
      L=L+6
      IF(N.GT.KMODXM(I)) GO TO 2460
      J=J+1
      IF(MOD(J,50).NE.1) GO TO 2455
      CALL NEWPG
      WRITE(W,123)
 2455 WRITE(W,162)(UTX(M,N,I),M=1,3),(UTY(M,N,I),M=1,3),(KI(M),M=K,L)
 2460 CONTINUE
      DO 2470N=1,4
      K=L+1
      L=L+6
      IF(N.GT.KMODXM(I)) GO TO 2470
      J=J+1
      IF(MOD(J,50).NE.1) GO TO 2465
      CALL NEWPG
      WRITE(W,123)
 2465 WRITE(W,162)(URX(M,N,I),M=1,3),(URY(M,N,I),M=1,3),(KI(M),M=K,L)
 2470 CONTINUE
      K=L+1
      L=L+1
      IF(KMODXM(I).LE.0) GO TO 2499
      J=J+1
      IF(MOD(J,50).NE.1) GO TO 2480
      CALL NEWPG
      WRITE(W,123)
 2480 WRITE(W,124) PHFM(I),KI(K)
 2499 CONTINUE
C
 2500 IF(LITE.LT.0.OR.LITE.GT.2) LITE=2
      IF(LITE.EQ.1) LITE=0
      IF(LITE.EQ.0) GO TO 22500
      L5='B'//L5(2:12)
      RHOM=1.
      DO 12500M=1,3
12500 SQRP(M)=1.
C
22500 DO 72500I=1,NA
      IF(ITF(I).NE.3) GO TO 62500
      IF(LITE.EQ.2) GO TO 42500
      DO 32500J=2,3
32500 BETA(J,I)=BETA(1,I)
      BETA(4,I)=BETA(1,I)*RCP(6)
      BETA(5,I)=BETA(1,I)*RCP(5)
      BETA(6,I)=BETA(1,I)*RCP(4)
      ITF(I)=2
      GO TO 62500
42500 POM=BETA(1,I)
      DO 52500J=1,6
52500 BETA(J,I)=POM*PRCP(J)
      ITF(I)=2
62500 IF(ITF(I).NE.1.AND.ITF(I).NE.2) ITF(I)=1
72500 CONTINUE
C
      NVX=0
      NVB=0
      NVAI=0
 2501 READ(M50,144) AT,K,ATO
      DO 2502M=1,NA
      IF(ATOM(M).EQ.AT) GO TO 2503
 2502 CONTINUE
      GO TO 2550
 2503 DO 2504I=1,NA
      IF(ATOM(I).EQ.ATO) GO TO 2505
 2504 CONTINUE
      GO TO 2512
 2505 NVAI=NVAI+1
      IF(M.LT.I) GO TO 2506
      WRITE(W,170) AT,ATO
      NVAI=NVAI-1
      GO TO 2501
 2506 NAI(1,NVAI)=M
      NAI(2,NVAI)=I
      NAIXB(NVAI)=K
      SUMAI(NVAI)=AI(M)+AI(I)
      IF(K.LE.0) SUMAI(NVAI)=A0(M)+A0(I)
      GO TO 2501
 2512 IF(K.EQ.2) GO TO 2520
      NVX=NVX+1
      NX(NVX)=M
      READ(M50,145)((RMX(I,J,NVX),J=1,3),I)
      DO 2515I=1,3
      POM=0.
      DO 2514J=1,3
 2514 POM=POM+RMX(I,J,NVX)*X(J,M)
 2515 RMT(I,NVX)=X(I,M)-POM
      GO TO 2501
 2520 NVB=NVB+1
      NB(NVB)=M
      READ(M50,146)((RMB(I,J,NVB),J=1,6),I=1,6)
      GO TO 2501
 2550 IF(NVX.LE.0.AND.NVB.LE.0.AND.NVAI.LE.0) GO TO 2585
      CALL NEWPG
      WRITE(W,147)
      IF(NVX.LE.0) GO TO 2570
      DO 2560M=1,NVX
      L=NX(M)
 2560 WRITE(W,148) ATOM(L),((RMX(I,J,M),J=1,3),RMT(I,M),I=1,3)
 2570 IF(NVB.LE.0) GO TO 2580
      DO 2571M=1,NVB
      L=NB(M)
 2571 WRITE(W,149) ATOM(L),((RMB(I,J,M),J=1,6),I=1,6)
 2580 IF(NVAI.LE.0) GO TO 2585
      DO 2584M=1,NVAI
      J1=NAI(1,M)
      J2=NAI(2,M)
      IF(NAIXB(M).LT.0) GO TO 12583
      WRITE(W,152) ATOM(J1),ATOM(J2),SUMAI(M)
12580 GO TO (2581,2582,2583),IABS(NAIXB(M))
 2581 WRITE(W,153)
      GO TO 2584
 2582 WRITE(W,154)
      GO TO 2584
 2583 WRITE(W,155)
      GO TO 2584
12583 WRITE(W,184) ATOM(J1),ATOM(J2)
      GO TO 12580
 2584 CONTINUE
C
 2585 IF(NMOL.LE.0) GO TO 2599
      J=0
      DO 2588I=1,NMOL
      IF(NPG(I).LE.1) GO TO 2588
      IF(J.NE.0) GO TO 2587
      CALL NEWPG
      WRITE(W,169)
 2587 WRITE(W,172) MOLNAME(I),(KPG(J,I),J=1,NPG(I))
 2588 CONTINUE
C
 2599 IF(IEXT.LE.0) GO TO 2900
      CALL NEWPG
      WRITE(W,132) RR,UHMON
      C2THMQ=COS(2.*UHMON*TORAD)**2
      IF(IEXT.EQ.1) WRITE(W,133)
      IF(IEXT.NE.1) WRITE(W,134)
      GO TO (2600,2610,2620),ITYPEX
 2600 WRITE(W,135)
      GO TO 2630
 2610 WRITE(W,136)
      GO TO 2630
 2620 WRITE(W,137)
 2630 IF(ITYPEX.EQ.2) GO TO 2640
      IF(IDISTR.EQ.1) WRITE(W,138)
      IF(IDISTR.EQ.2) WRITE(W,139)
 2640 IF(IEXT.EQ.1.OR.IEXT.EQ.3) GO TO 2700
      IF(IMOSAIC.EQ.0) WRITE(W,140)
      IF(IMOSAIC.EQ.1) WRITE(W,141)
 2700 IF(IAB.EQ.0) GO TO 2800
      WRITE(W,142) AA,BA,DADMI
 2800 WRITE(W,143) (EC(I),I=1,6),(KI(I),I=7,12),(EC(I),I=7,12),
     1                      (KI(I),I=13,18)
C
 2900 NP=NA*92+18
      NPMP=NP+1
      NPMK=NP+NMOL*49
      J=0
      NAMM=0
      NPSS(1)=0
      DO 5000I=1,9
      NPSI=0
      DO 3000M=1,NPMK
      IF(KI(M).EQ.I) NPSI=NPSI+1
 3000 CONTINUE
      IF(NPSI.LE.200) GO TO 4100
      WRITE(W,126) I,NPSI
      STOP
 4100 NPS(I)=NPSI
      J=J+NPSI
      NPSS(I+1)=J
      NAM(I)=NPSI*(NPSI+1)/2
      NAMM=NAMM+NAM(I)
 5000 CONTINUE
      IF(NAMM.LE.75000) GO TO 7000
      WRITE(W,127)
      STOP
C
 7000 CALL NEWPG
      WRITE(W,128) (I,I=1,9),NPS,J,NAM,NAMM
      READ(M50,104) TLUM,NCYKL
      IF(TLUM.LE.0.0.OR.TLUM.GT.1.0) TLUM=1.0
      READ(M50,113) DIF
      RETURN
  100 FORMAT(20A4)
  101 FORMAT(2I5,4X,A1,4X,A1)
  102 FORMAT(3(F11.6,2I2),3X,2A1,1X,2I2)
  103 FORMAT(8F9.3)
  104 FORMAT(F5.2,I5)
  105 FORMAT('0NON-CENTROSYMMETRIC SUPER-SPACE GROUP')
  106 FORMAT('0CENTROSYMMETRIC SUPER-SPACE GROUP')
  107 FORMAT('0SYMMETRY CARDS'//)
  108 FORMAT('0CENTERING OF CELL : ',A1)
  109 FORMAT(4(F20.6,2I2))
  110 FORMAT('0ATOMIC SCATTERING TABLES :')
  111 FORMAT('0',I5,5X,8F9.3/11X,8F9.3
     1                      /11X,8F9.3
     2                      /11X,8F9.3)
  112 FORMAT(11I5)
  113 FORMAT(8F10.5)
  114 FORMAT('0WEIGHT 1/SIGMA**2 COEFFICIENT OF UNSTABILITY IS ',F8.2)
  115 FORMAT('0UNIT WEIGHT')
  116 FORMAT('0WEIGHT AFTER CRUIKSHANCK WITH YOMIN =',F7.2,' YOMAX =',
     1        F9.2)
  117 FORMAT('0REFLECTIONS WITH ABS(YO-YC)>',F5.2,'*SIGYO WILL BE SYMBOL
     1IZED IN OUTPUT BY #')
  118 FORMAT('0REFLECTIONS WITH NXX>',I2,' WILL BE OMITTED')
  119 FORMAT('0ONLY REFLECTIONS SYMBOLISED BY # (SEE ABOVE) WILL BE PRIN
     1TED')
  120 FORMAT(18X,4F9.6/6F9.6)
  121 FORMAT(A8,I3,3I1,I3,I1,4F9.6/6F9.6,6X,10I1)
  122 FORMAT('0F000 =',F8.0)
  123 FORMAT('0INPUT PARAMETERS AND THEIR KEYS FOR REFINEMENT'//)
  124 FORMAT(1X,F9.6,48X,'=>',I1,'<=')
  125 FORMAT(1X,A8,I3,3I1,I3,I1,4F9.6)
  126 FORMAT('0ERROR - NUMBER PARAMETERS TO BE REFINED IN ',I2,'.BLOCK I
     1S :',I5,' THAT IS GREATER THAN 200')
  127 FORMAT('0NUMBER OF ELEMENTS OF MATRIX IS',I10,' THAT IS GREATER TH
     1AN 75000')
  128 FORMAT('0BLOCK :',9I8,'  TOGETHER'/'0NPS   :',9I8,I10/' NAM   :',
     19I8,I10)
  129 FORMAT(3F10.3,3F10.2,F10.5)
  130 FORMAT(10F7.4)
  131 FORMAT(6F9.6,6X,6I1)
  132 FORMAT('0INFORMATION ON EXTINCTION CORRECTION :'/39('=')//'0RADIUS
     1 OF SHERICAL SAMPLE :',F9.4,' CM'/'0MONOCHROMATOR ANGLE :',F8.2,
     2' DEGREES')
  133 FORMAT('0ISOTROPIC CORRECTION')
  134 FORMAT('0ANISOTROPIC CORRECTION')
  135 FORMAT('0EXTINCTION TYP I')
  136 FORMAT('0EXTINCTION TYP II')
  137 FORMAT('0GENERAL TYP OF EXTINCTION')
  138 FORMAT('0GAUSSIAN DISTRIBUTION')
  139 FORMAT('0LORENTZIAN DISTRIBUTION')
  140 FORMAT('0ANISOTROPY OF MOSAIC SPREAD ACCORDING COPPENS AND HAMILTO
     1N')
  141 FORMAT('0ANISOTROPY OF MOSAIC SPREAD ACCORDING NELMES AND THORNLEY
     1')
  142 FORMAT('0MODIFICATION FOR ABSORPTION'/'0A(SIN(TH)) STEP 0.1 IN SIN
     1(TH)     : ',10F8.2                  /' B(SIN(TH)) STEP 0.1 IN SIN
     2(TH)     : ',10F8.2                  /' DADMI(TH)  STEP 5 DEGREES
     3 IN TH   : ',10F8.4/45X,9F8.4)
  143 FORMAT('0LIST OF PARAMETERS OF EXTINCTION AND THEIR KEYS OF REFINE
     1MENT' /'0RHO    : ',6F9.6,5X,'=>',6I1,'<='/
     2       '0G(I,J) : ',6F9.6,5X,'=>',6I1,'<=')
  144 FORMAT(A8,I2,2X,A8)
  145 FORMAT(3F5.0)
  146 FORMAT(6F5.0)
  147 FORMAT('0RESTRINCTIONS :'/16('='))
  148 FORMAT('0RESTRINCTION ON COORDINATES OF ATOM :',A8/
     1       '0X = ',F5.1,'*X +',F5.1,'*Y +',F5.1,'*Z +',F10.5/
     2       ' Y = ',F5.1,'*X +',F5.1,'*Y +',F5.1,'*Z +',F10.5/
     3       ' Z = ',F5.1,'*X +',F5.1,'*Y +',F5.1,'*Z +',F10.5)
  149 FORMAT('0RESTRINCTION ON TEMPERATURE PARAMETERS OF ATOM :',A8/
     1       '0B(1,1) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3
     2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/
     1       ' B(2,2) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3
     2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/
     1       ' B(3,3) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3
     2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/
     1       ' B(1,2) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3
     2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/
     1       ' B(1,3) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3
     2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)'/
     1       ' B(2,3) = ',F5.1,'*B(1,1) +',F5.1,'*B(2,2) +',F5.1,'*B(3,3
     2) +',F5.1,'*B(1,2) +',F5.1,'*B(1,3) +',F5.1,'*B(2,3)')
  150 FORMAT('0SIN(TH)/LAMBDA LIMITS FOR ACCEPTANCE OF REFLECTION',
     1       2F10.5)
  151 FORMAT(F9.6,51X,I1)
  152 FORMAT('0OCCUPATIONAL RESTRINCTION OF ATOMS : ',A8,' AND : ',A8,
     1        ' SUMMA OF OCCUPATIONS WILL BE : ',F9.6)
  153 FORMAT(38X,'COORDINATES ,TEMPERATURE AND POSITIONAL MODULATION PAR
     1AMETERS OF THESE ATOMS WILL BE IDENTICAL')
  154 FORMAT(38X,'COORDINATE OF THESE ATOMS WILL BE IDENTICAL')
  155 FORMAT(38X,'WITHOUT RESTRINCTION ON COORDINATES, TEMPERATURE AND M
     1ODULATION PARAMETERS OF THESE ATOMS')
  156 FORMAT('0IRRATIONAL PART OF MODULATION VECTOR Q : ',3F9.5)
  157 FORMAT('0  RATIONAL PART OF MODULATION VECTOR Q : ',3F9.5,' DERIVE
     1D FROM SYMBOL : ',A1)
  158 FORMAT('0ONLY MAIN REFLECTIONS WILL BE ACCEPTED INTO CALCULATION')
  159 FORMAT('0ONLY SATELLITE REFLECTIONS WILL BE ACCEPTED INTO CALCULAT
     1ION')
  160 FORMAT('0ONLY REFLECTIONS WITH SATELLITE INDEX',I2,' OR',I2,'WILL
     1ACCEPTED INTO CALCULATION')
  161 FORMAT(1X,6F9.6,3X,'=>',10I1,'<=')
  162 FORMAT(1X,6F9.6,3X,'=>', 6I1,'<=')
  163 FORMAT(A8,3X,2I3,10X,3F9.6)
  164 FORMAT(1X,A8,3X,2I3,10X,3F9.6)
  165 FORMAT('0MOLECULAR CENTERS WILL NOT BE CHANGED')
  166 FORMAT('0MOLECULAR CENTERS WILL BE CHANGED TO BE')
  167 FORMAT('+',40X,'A GEOMETRICAL CENTER OF ATOMS')
  168 FORMAT('+',40X,'A CENTER OF MASS OF ATOMS')
  169 FORMAT('0INFORMATION ON A CRYSTALLOGRAPHIC POINT SYMMETRY OF RIGID
     1 GROUPS (MOLECULES)'/'0ELEMENTS OF POINT SYMMETRY ARE GIVEN BY NUM
     2BER OF SYMMETRY CARD'/' NEGATIVE NUMBER MEANS THE ELEMENT HAS TO B
     3E COMBINED WITH INVERSION CENTER'//)
  170 FORMAT('OERROR ATOM : ',A8,' IS PLACED BEHIND ATOM : ',A8,' RESTRI
     1NTION WILL NOT BE USED')
  171 FORMAT('0ELEMENTS OF COVARIATION MATRIX LARGER THAN 0.9 WILL BE PR
     1INTED')
  172 FORMAT(' GROUP : ',A8,' HAS ELEMENTS :',24I4)
  173 FORMAT('0PRINT OF REFLECTIONS SUPRESSED')
  174 FORMAT('0PRINT OF REFLECTIONS IN EVERY CYCLE OF REFINEMENT')
  175 FORMAT('0PRINT OF REFLECTIONS BEFORE FIRST AND AFTER LAST CYCLE OF
     1 REFINEMENT')
  176 FORMAT('0INPUT TEMPERATURE PARAMETERS ARE BETA')
  177 FORMAT('                              AND WILL BE CHANGED TO BE U'
     1) 
  178 FORMAT('0INPUT TEMPERATURE PARAMETERS ARE U')
  179 FORMAT('                              AND WILL BE CHANGED TO BE BE
     1TA')
  180 FORMAT('0REFINEMENT BASED ON F')
  181 FORMAT('0REFINEMENT BASED ON F**2')
  182 FORMAT(2F9.6,42X,2I1)
  183 FORMAT(1X,2F9.6,39X,'=>',2I1,'<=')
  184 FORMAT('0OCCUPATIONAL MODULATION OF ATOM : ',A8,' AND ATOM : ',A8,
     1       ' HAS TO BE OPPOSITE')
      END
      SUBROUTINE NAJDI(M)
      CHARACTER*8 NZ1,NZ2
      DATA NZ1/'JANA    '/
 1000 READ(M,100) NZ2
  100 FORMAT(A8)
      IF(NZ2.NE.NZ1) GO TO 1000
      RETURN
      END
      SUBROUTINE NEWPG
      COMMON/UNITS/ W,M40,M50
      COMMON/TITLE/ TITLE(20),NAME(20),PAGE
      INTEGER TITLE,PAGE,W
      PAGE=PAGE+1
      WRITE(W,100) PAGE,NAME,TITLE
  100 FORMAT('1REFINEMENT OF MODULATED STRUCTURE',75X,'PAGE =',I3/
     1       ' STRUCTURE : ',20A4/13X,20A4/)
      RETURN
      END
      SUBROUTINE CALC
      COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF,
     1              CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ
      COMMON/REFL/  YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708),
     1              DC(9708),DYP,IQ
      COMMON/CM4/   ATOM(100),AI(100),ISF(100),ITF(100),X(3,100),
     1              BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100)
     2             ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100),
     3              AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100),
     4              KFX(100),KFB(100)
      COMMON/IND/   H(24),K(24),L(24),T(24),HQ(24),KQ(24),LQ(24),HK2(24)
     1             ,HL2(24),KL2(24),FSMM(24),EPS(24),EPSM(24),MM,PHFP
      COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15)
      COMMON/SYMM/  IS(2,4,24),TS(3,24),TAU(24),FSM(24)
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON/STALY/ PI2
      LOGICAL CSYM,NULOVA,IZO,HARM,EPSL,MAIN
      DIMENSION DADUXX(4),DADUXY(4),DADUXZ(4),DADUYX(4),DADUYY(4),
     1          DADUYZ(4),DBDUXX(4),DBDUXY(4),DBDUXZ(4),DBDUYX(4),
     2          DBDUYY(4),DBDUYZ(4),UXX(4),UXY(4),UXZ(4),UYX(4),UYY(4),
     3          UYZ(4),BESP(9,4),DBFDU1(9,4),DBFDU2(9,4),DCHIDU1(4),
     4          DCHIDU2(4),DSDU1(4),DSDU2(4),CHI(4),DADU1(4),DADU2(4),
     5          DBDU1(4),DBDU2(4),DSDCHI(4),SB(5),DSBDAX(4),DSBDAY(4),
     6          KSI(4),DKSIDAX(4),DKSIDAY(4),DADAX(4),DADAY(4),DBDAX(4)
     7         ,DBDAY(4)
      DIMENSION DADBX1(4),DADBX2(4),DADBX3(4),DADBX4(4),DADBX5(4),
     1          DBDBX1(4),DBDBX2(4),DBDBX3(4),DBDBX4(4),DBDBX5(4),
     2          DADBY1(4),DADBY2(4),DADBY3(4),DADBY4(4),DADBY5(4),
     3          DBDBY1(4),DBDBY2(4),DBDBY3(4),DBDBY4(4),DBDBY5(4),
     4          DADBX6(4),DBDBX6(4),DADBY6(4),DBDBY6(4),BX1(4),BX2(4),
     5          BX3(4),BX4(4),BX5(4),BX6(4),BY1(4),BY2(4),BY3(4),BY4(4),
     6          BY5(4),BY6(4),BESB(5,4),DBFDB1(5,4),DBFDB2(5,4),
     7          DETADB1(4),DETADB2(4),DSDB1(4),DSDB2(4),ETA(4),DADB1(4),
     8          DADB2(4),DBDB1(4),DBDB2(4),DSDETA(4),DSDB3(4),DSDU3(4),
     9          DBFDB3(5,4),DBFDU3(5,4)
      CHARACTER*8 ATOM,MOLNAME
      REAL K,L,KQ,LQ,KL2,KJ,LJ,KQJ,LQJ,KL2J,KSI
C
      A=0.
      B=0.
      AANOM=0.
      BANOM=0.
      M=19
      DO 500I=1,6
  500 DER(I)=0.
      EXTKOR=1.
      MAIN=MM.EQ.0
      DO 6000I=1,NA
      AII=AI(I)
      QCNTMI=QCNTM(I)
      KMODXI=KMODX(I)
      KFXI=KFX(I)
      KMODBI=KMODB(I)
      KFBI=KFB(I)
      KMODSI=KMODS(I)
      IF(KMODXI.LT.0) KMODXI=KMODXM(-KMODXI)
      HARM=KMODXI.LE.0.AND.KMODSI.LE.0.AND.KMODBI.LE.0
      IF(AII.NE.0..AND.(.NOT.HARM.OR.MAIN)) GO TO 910
      DO 900J=M,M+91
      DER(J)=0.
  900 DC(J)=0.
      GO TO 6000
  910 XI=X(1,I)
      YI=X(2,I)
      ZI=X(3,I)
      J=ISF(I)
      FXI=FX(J)
      FYI=FY(J)
      IZO=ITF(I).EQ.1
      B1=BETA(1,I)
      EXPIJ=1.
      IF(IZO) GO TO 1001
      B2=BETA(2,I)
      B3=BETA(3,I)
      B4=BETA(4,I)
      B5=BETA(5,I)
      B6=BETA(6,I)
 1001 IF(HARM) GO TO 1003
      SB(1)=A0(I)
      DO 8J=1,4
      IF(J.LE.KMODSI) GO TO 6
    5 SB(J+1)=0.
      KSI(J)=0.
      DSBDAX(J)=0.
      DSBDAY(J)=0.
      GO TO 8
    6 UQR=AX(J,I)**2+AY(J,I)**2
      IF(UQR.LE.0) GO TO 5
      U=SQRT(UQR)
      UQR=1./UQR
      POM=U*.5
      SB(J+1)=POM
      DSBDAX(J)=AX(J,I)/U*.5
      DSBDAY(J)=AY(J,I)/U*.5
      KSI(J)=ATAN2(-AX(J,I),AY(J,I))
      DKSIDAX(J)=-AY(J,I)*UQR
      DKSIDAY(J)= AX(J,I)*UQR
    8 CONTINUE
      M1P=0
      M1K=0
      M2P=0
      M2K=0
      M3P=0
      M3K=0
      M4P=0
      M4K=0
      M5P=0
      M5K=0
      M6P=0
      M6K=0
      M7P=0
      M7K=0
      M8P=0
      M8K=0
      M9P=-KMODSI
      M9K= KMODSI
      IF(KMODXI.EQ.0) GO TO 10
      M1P=-4
      M1K= 4
      IF(KMODXI.EQ.1) GO TO 10
      M2P=-4
      M2K= 4
      IF(KMODXI.EQ.2) GO TO 10
      M3P=-4
      M3K= 4
      IF(KMODXI.EQ.3) GO TO 10
      M4P=-4
      M4K= 4
   10 IF(KMODBI.EQ.0) GO TO 15
      M5P=-4
      M5K= 4
      IF(KMODBI.EQ.1) GO TO 15
      M6P=-4
      M6K= 4
      IF(KMODBI.EQ.2) GO TO 15
      M7P=-4
      M7K= 4
      IF(KMODBI.EQ.3) GO TO 15
      M8P=-4
      M8K= 4
   15 LP=0
      LK=4
      IF(KMODXI+KMODSI+KMODBI.NE.1) GO TO 20
      LP=IABS(MM)
      LK=IABS(MM)
   20 DO 1002J=1,KMODXI
      UXX(J)=UX(1,J,I)
      UXY(J)=UX(2,J,I)
      UXZ(J)=UX(3,J,I)
      UYX(J)=UY(1,J,I)
      UYY(J)=UY(2,J,I)
 1002 UYZ(J)=UY(3,J,I)
      DO 11002J=1,KMODBI
      BX1(J)=BX(1,J,I)
      BX2(J)=BX(2,J,I)
      BX3(J)=BX(3,J,I)
      BX4(J)=BX(4,J,I)
      BX5(J)=BX(5,J,I)
      BX6(J)=BX(6,J,I)
      BY1(J)=BY(1,J,I)
      BY2(J)=BY(2,J,I)
      BY3(J)=BY(3,J,I)
      BY4(J)=BY(4,J,I)
      BY5(J)=BY(5,J,I)
11002 BY6(J)=BY(6,J,I)
      PHFI=PHF(I)
 1003 DADFII=0.
      IF(CSYM) GO TO 1005
      DBDFII=0.
 1005 DADXI=0.
      DADYI=0.
      DADZI=0.
      IF(IZO) GO TO 1010
      DADB1I=0.
      DADB2I=0.
      DADB3I=0.
      DADB4I=0.
      DADB5I=0.
      DADB6I=0.
 1010 IF(CSYM) GO TO 1015
      DBDXI=0.
      DBDYI=0.
      DBDZI=0.
      IF(IZO) GO TO 1015
      DBDB1I=0.
      DBDB2I=0.
      DBDB3I=0.
      DBDB4I=0.
      DBDB5I=0.
      DBDB6I=0.
 1015 DADA0=0.
      DBDA0=0.
      IF(KMODSI.LE.0) GO TO 1020
      DO 1017J=1,KMODSI
      DADAX(J)=0.
      DADAY(J)=0.
      IF(CSYM) GO TO 1017
      DBDAX(J)=0.
      DBDAY(J)=0.
 1017 CONTINUE
 1020 IF(KMODXI.LE.0) GO TO 1022
      DO 1021J=1,KMODXI
      DADUXX(J)=0.
      DADUXY(J)=0.
      DADUXZ(J)=0.
      DADUYX(J)=0.
      DADUYY(J)=0.
      DADUYZ(J)=0.
      IF(CSYM) GO TO 1021
      DBDUXX(J)=0.
      DBDUXY(J)=0.
      DBDUXZ(J)=0.
      DBDUYX(J)=0.
      DBDUYY(J)=0.
 1021 DBDUYZ(J)=0.
 1022 DO 1023J=1,KMODBI
      DADBX1(J)=0.
      DADBX2(J)=0.
      DADBX3(J)=0.
      DADBX4(J)=0.
      DADBX5(J)=0.
      DADBX6(J)=0.
      DADBY1(J)=0.
      DADBY2(J)=0.
      DADBY3(J)=0.
      DADBY4(J)=0.
      DADBY5(J)=0.
      DADBY6(J)=0.
      IF(CSYM) GO TO 1023
      DBDBX1(J)=0.
      DBDBX2(J)=0.
      DBDBX3(J)=0.
      DBDBX4(J)=0.
      DBDBX5(J)=0.
      DBDBX6(J)=0.
      DBDBY1(J)=0.
      DBDBY2(J)=0.
      DBDBY3(J)=0.
      DBDBY4(J)=0.
      DBDBY5(J)=0.
 1023 DBDBY6(J)=0.
 1025 DADPHF=0.
      DBDPHF=0.
      DO 5000J=1,NS
      HJ  =H(J)
      KJ  =K(J)
      LJ  =L(J)
      EPSJ=EPS(J)
      EPSL=EPSJ.LE.0.
      IF(IZO) GO TO 2001
      HQJ =HQ(J)
      KQJ =KQ(J)
      LQJ =LQ(J)
      HK2J=HK2(J)
      HL2J=HL2(J)
      KL2J=KL2(J)
 2001 FI   =    T(J)+HJ*XI+KJ*YI+LJ*ZI
      IF(HARM) GO TO 2070
      DO 2020KK=1,KMODXI
      U1=UXX(KK)*HJ+UXY(KK)*KJ+UXZ(KK)*LJ
      IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN
      U2=UYX(KK)*HJ+UYY(KK)*KJ+UYZ(KK)*LJ
      UQR=U1**2+U2**2
      IF(UQR.LE.0.) GO TO 2006
      U=SQRT(UQR)
      RU=1./U
      UQR=1./UQR
      CHI(KK)=ATAN2(U2,U1)
      DO 2002LL=LP,LK
      POM=BESJ(U,LL,DIF,0)
      BESP(LL+5,KK)=POM
      IF(LL.EQ.0) GO TO 2002
      IF(MOD(LL,2).EQ.1) POM=-POM
      BESP(-LL+5,KK)=POM
 2002 CONTINUE
      POMB=BESJ(U,LK+1,DIF,0)*RU
      DO 2005LL=LK,LP,-1
      POMA=BESP(LL+5,KK)*RU
      POM=-POMB
      IF(LL.EQ.0) GO TO 2004
      POM=POM+FLOAT(LL)*RU*POMA
      POMU1=POM*U1
      POMU2=POM*U2
      DBFDU1(LL+5,KK)=POMU1
      DBFDU2(LL+5,KK)=POMU2
      IF(MOD(LL,2).EQ.0) GO TO 2003
      POMU1=-POMU1
      POMU2=-POMU2
 2003 DBFDU1(-LL+5,KK)=POMU1
      DBFDU2(-LL+5,KK)=POMU2
      GO TO 2005
 2004 DBFDU1(5,KK)=POM*U1
      DBFDU2(5,KK)=POM*U2
 2005 POMB=POMA
C
      DCHIDU1(KK)=-UQR*U2
      DCHIDU2(KK)= UQR*U1
      GO TO 2020
 2006 DO 2010LL=LP,LK
      IF(LL.EQ.0) THEN
      BESP(LL+5,KK)=1.
      DBFDU1(LL+5,KK)=0.
      DBFDU2(LL+5,KK)=0.
      ELSE
      IF(LL.EQ.1) THEN
      BESP( LL+5,KK)=0.
      BESP(-LL+5,KK)=0.
      DBFDU1( LL+5,KK)= .5
      DBFDU1(-LL+5,KK)=-.5
      DBFDU2( LL+5,KK)= .5
      DBFDU2(-LL+5,KK)=-.5
      ELSE
      BESP( LL+5,KK)=0.
      BESP(-LL+5,KK)=0.
      DBFDU1( LL+5,KK)=0.
      DBFDU1(-LL+5,KK)=0.
      DBFDU2( LL+5,KK)=0.
      DBFDU2(-LL+5,KK)=0.
      ENDIF
      ENDIF
 2010 CONTINUE
      CHI(KK)=0.
      DCHIDU1(KK)=0.
      DCHIDU2(KK)=0.
      ELSE
      CHI(KK)=PI2*UYX(KK)
      DO 2015LL=LP,LK
      CALL XMOD(BESP(LL+5,KK),DBFDU1(LL+5,KK),DBFDU2(LL+5,KK),
     1          DBFDU3(LL+5,KK),U1,UYY(KK),UYZ(KK),LL,KFXI)
      CALL XMOD(BESP(-LL+5,KK),DBFDU1(-LL+5,KK),DBFDU2(-LL+5,KK),
     1          DBFDU3(-LL+5,KK),U1,UYY(KK),UYZ(KK),-LL,KFXI)
 2015 CONTINUE
      DCHIDU1(KK)=PI2
      ENDIF
 2020 CONTINUE
      DO 12020KK=1,KMODBI
      BB1=-HQJ *BX1(KK)-KQJ *BX2(KK)-LQJ*BX3(KK)-HK2J*BX4(KK)
     1    -HL2J*BX5(KK)-KL2J*BX6(KK)
      IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN
      BB2=-HQJ *BY1(KK)-KQJ *BY2(KK)-LQJ*BY3(KK)-HK2J*BY4(KK)
     1    -HL2J*BY5(KK)-KL2J*BY6(KK)
      UQR=BB1**2+BB2**2
      IF(UQR.LE.0.) GO TO 12006
      U=SQRT(UQR)
      RU=1./U
      UQR=1./UQR
      ETA(KK)=ATAN2(BB2,BB1)
      DO 12002LL=LP,LK
12002 BESB(LL+1,KK)=BESJ(U,LL,DIF,1)
      POMB=BESJ(U,LK+1,DIF,1)*RU
      DO 12005LL=LK,LP,-1
      POMA=BESB(LL+1,KK)*RU
      POM=POMB
      IF(LL.NE.0) POM=POM+FLOAT(LL)*RU*POMA
      DBFDB1(LL+1,KK)=POM*BB1
      DBFDB2(LL+1,KK)=POM*BB2
12005 POMB=POMA
      DETADB1(KK)=-UQR*BB2
      DETADB2(KK)= UQR*BB1
      GO TO 12020
12006 DO 12010LL=LP,LK
      IF(LL.EQ.0) THEN
      BESB(1,KK)=1.
      DBFDB1(1,KK)=0.
      DBFDB2(1,KK)=0.
      ELSE
      IF(LL.EQ.1) THEN
      BESB(2,KK)=0.
      DBFDB1(2,KK)=.5
      DBFDB2(2,KK)=.5
      ELSE
      BESB(LL+1,KK)=0.
      DBFDB1(LL+1,KK)=0.
      DBFDB2(LL+1,KK)=0.
      ENDIF
      ENDIF
12010 CONTINUE
      ELSE
      ETA(KK)=PI2*BY1(KK)
      DO 22002LL=LP,LK
      CALL BMOD(BESB(LL+1,KK),DBFDB1(LL+1,KK),DBFDB2(LL+1,KK),
     1          DBFDB3(LL+1,KK),BB1,BY2(KK),BY3(KK),LL,KFBI)
22002 CONTINUE
      DETADB1(KK)=PI2
      ENDIF
12020 CONTINUE
C
      FIP=FI+FSMM(J)+EPSM(J)*QCNTMI
      COSIJ=0.
      SINIJ=0.
      DO 2021KK=1,KMODXI
      DADU1(KK)=0.
      DADU2(KK)=0.
      IF(CSYM) GO TO 2021
      DBDU1(KK)=0.
      DBDU2(KK)=0.
 2021 CONTINUE
      DADU3=0.
      DADU4=0.
      DBDU3=0.
      DBDU4=0.
      DO 12021KK=1,KMODBI
      DADB1(KK)=0.
      DADB2(KK)=0.
      IF(CSYM) GO TO 12021
      DBDB1(KK)=0.
      DBDB2(KK)=0.
12021 CONTINUE
      DADB3=0.
      DBDB3=0.
      DADB4=0.
      DBDB4=0.
      DO 2069M9=M9P,M9K
      DO 2068M8=M8P,M8K
      DO 2067M7=M7P,M7K
      DO 2066M6=M6P,M6K
      DO 2065M5=M5P,M5K
      DO 2060M4=M4P,M4K
      DO 2050M3=M3P,M3K
      DO 2040M2=M2P,M2K
      M1=MM-2*M2-3*M3-4*M4-M5-2*M6-3*M7-4*M8-M9
      IF(KFXI.NE.0) THEN
      IF(KMODXI.EQ.2) M1=M1+M2
      IF(KMODXI.EQ.3) M1=M1+2*M3
      IF(KMODXI.EQ.4) M1=M1+3*M4
      ENDIF
      IF(KFBI.NE.0) THEN
      IF(KMODBI.EQ.2) M1=M1+M6
      IF(KMODBI.EQ.3) M1=M1+2*M7
      IF(KMODBI.EQ.4) M1=M1+3*M8
      ENDIF
      IF(M1.LT.M1P.OR.M1.GT.M1K) GO TO 2040
      MS=M1+M2+M3+M4
      IF(KFXI.NE.0) THEN
      IF(KMODXI.EQ.1) MS=MS-M1
      IF(KMODXI.EQ.2) MS=MS-M2
      IF(KMODXI.EQ.3) MS=MS-M3
      IF(KMODXI.EQ.4) MS=MS-M4
      ENDIF
      CHIS=0.
      S=1.
      IF(MOD(IABS(MS),2).EQ.1.AND..NOT.EPSL) S=-S
      IF(KK.EQ.KMODXI.AND.KFXI.NE.0.AND.EPSL) S=-S
C
      DO 2030KK=1,KMODXI
      GO TO (2022,2023,2024,2025),KK
 2022 FMP=-M1
      MP=M1+5
      GO TO 2026
 2023 FMP=-M2
      MP=M2+5
      GO TO 2026
 2024 FMP=-M3
      MP=M3+5
      GO TO 2026
 2025 FMP=-M4
      MP=M4+5
 2026 IF(KK.EQ.KMODXI.AND.KFXI.NE.0) FMP=-FMP
      IF(EPSL) FMP=-FMP
      DSDU1(KK)=S*DBFDU1(MP,KK)
      DSDU2(KK)=S*DBFDU2(MP,KK)
      POM=BESP(MP,KK)
      S=S*POM
      IF(ABS(S).LT.DIF) GO TO 2040
      IF(KK.LE.1) GO TO 2028
      DO 2027LL=1,KK-1
      DSDU1(LL)=DSDU1(LL)*POM
 2027 DSDU2(LL)=DSDU2(LL)*POM
 2028 DSDCHI(KK)=FMP
      CHIS=CHIS+FMP*CHI(KK)
 2030 CONTINUE
      DO 12030KK=1,KMODBI
      GO TO (12022,12023,12024,12025),KK
12022 FMP=-M5
      MP=IABS(M5)+1
      GO TO 12026
12023 FMP=-M6
      MP=IABS(M6)+1
      GO TO 12026
12024 FMP=-M7
      MP=IABS(M7)+1
      GO TO 12026
12025 FMP=-M8
      MP=IABS(M8)+1
12026 IF(KK.EQ.KMODBI.AND.KFBI.NE.0) FMP=-FMP
      IF(EPSL) FMP=-FMP
      DSDB1(KK)=S*DBFDB1(MP,KK)
      DSDB2(KK)=S*DBFDB2(MP,KK)
      IF(KK.EQ.KMODBI.AND.KFBI.NE.0) DSDB3(KK)=S*DBFDB3(MP,KK)
      POM=BESB(MP,KK)
      S=S*POM
      IF(ABS(S).LT.DIF) GO TO 2040
      DO 12027LL=1,KK-1
      DSDB1(LL)=DSDB1(LL)*POM
      IF(KK.EQ.KMODBI.AND.KFBI.NE.0) DSDB3(LL)=DSDB3(LL)*POM
12027 DSDB2(LL)=DSDB2(LL)*POM
      DO 22028LL=1,KMODXI
      DSDU1(LL)=DSDU1(LL)*POM
22028 DSDU2(LL)=DSDU2(LL)*POM
12028 DSDETA(KK)=FMP
      CHIS=CHIS+FMP*ETA(KK)
      IF(KK.NE.KMODBI.OR.KFBI.EQ.0) CHIS=CHIS-.25*FMP*PI2
12030 CONTINUE
      IF(KMODSI.LE.0) GO TO 2034
      MP=IABS(M9)
      IF(MP.EQ.0) GO TO 2031
      DSDAX=S*DSBDAX(MP)
      DSDAY=S*DSBDAY(MP)
      FMP=-SIGN(1.,FLOAT(M9))
      IF(EPSL) FMP=-FMP
      CHIS=CHIS+KSI(MP)*FMP
      DSDKSI=FMP
      GO TO 12031
 2031 DSDA0=S
12031 POM=SB(MP+1)
      S=S*POM
      DO 2032KK=1,KMODXI
      DSDU1(KK)=DSDU1(KK)*POM
 2032 DSDU2(KK)=DSDU2(KK)*POM
      DO 2033KK=1,KMODBI
      DSDB1(KK)=DSDB1(KK)*POM
      IF(KK.EQ.KMODBI.AND.KFBI.NE.0) DSDB3(KK)=DSDB3(KK)*POM
 2033 DSDB2(KK)=DSDB2(KK)*POM
C
 2034 POM=FIP+CHIS
      COSP=COS(POM)
      SINP=SIN(POM)
      COSPS=COSP*S
      SINPS=SINP*S
      COSIJ=COSIJ+COSPS
      SINIJ=SINIJ+SINPS
      DO 2035KK=1,KMODXI
      POMA=-DSDCHI(KK)*SINPS
      POMB= DSDCHI(KK)*COSPS
      IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN
      DADU1(KK)=DADU1(KK)+POMA*DCHIDU1(KK)+COSP*DSDU1(KK)
      DADU2(KK)=DADU2(KK)+POMA*DCHIDU2(KK)+COSP*DSDU2(KK)
      ELSE
      DADU1(KK)=DADU1(KK)+COSP*DSDU1(KK)
      DADU2(KK)=DADU2(KK)+COSP*DSDU2(KK)
      DADU3    =DADU3    +COSP*DSDU3(KK)
      DADU4    =DADU4    +POMA*DCHIDU1(KK)
      ENDIF
      IF(CSYM) GO TO 2035
      IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN
      DBDU1(KK)=DBDU1(KK)+POMB*DCHIDU1(KK)+SINP*DSDU1(KK)
      DBDU2(KK)=DBDU2(KK)+POMB*DCHIDU2(KK)+SINP*DSDU2(KK)
      ELSE
      DBDU1(KK)=DBDU1(KK)+SINP*DSDU1(KK)
      DBDU2(KK)=DBDU2(KK)+SINP*DSDU2(KK)
      DBDU3    =DBDU3    +SINP*DSDU3(KK)
      DBDU4    =DBDU4    +POMB*DCHIDU1(KK)
      ENDIF
 2035 CONTINUE
      DO 2036KK=1,KMODBI
      POMA=-DSDETA(KK)*SINPS
      POMB= DSDETA(KK)*COSPS
      IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN
      DADB1(KK)=DADB1(KK)+POMA*DETADB1(KK)+COSP*DSDB1(KK)
      DADB2(KK)=DADB2(KK)+POMA*DETADB2(KK)+COSP*DSDB2(KK)
      ELSE
      DADB1(KK)=DADB1(KK)+COSP*DSDB1(KK)
      DADB2(KK)=DADB2(KK)+COSP*DSDB2(KK)
      DADB3    =DADB3    +COSP*DSDB3(KK)
      DADB4    =DADB4    +POMA*DETADB1(KK)
      ENDIF
      IF(CSYM) GO TO 2036
      IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN
      DBDB1(KK)=DBDB1(KK)+POMB*DETADB1(KK)+SINP*DSDB1(KK)
      DBDB2(KK)=DBDB2(KK)+POMB*DETADB2(KK)+SINP*DSDB2(KK)
      ELSE
      DBDB1(KK)=DBDB1(KK)+SINP*DSDB1(KK)
      DBDB2(KK)=DBDB2(KK)+SINP*DSDB2(KK)
      DBDB3    =DBDB3    +SINP*DSDB3(KK)
      DBDB4    =DBDB4    +POMB*DETADB1(KK)
      ENDIF
 2036 CONTINUE
      IF(KMODSI.LE.0) GO TO 2040
      IF(MP.EQ.0) GO TO 2038
      POMA=-DSDKSI*SINPS
      POMB= DSDKSI*COSPS
      DADAX(MP)=DADAX(MP)+POMA*DKSIDAX(MP)+COSP*DSDAX
      DADAY(MP)=DADAY(MP)+POMA*DKSIDAY(MP)+COSP*DSDAY
      IF(CSYM) GO TO 2040
      DBDAX(MP)=DBDAX(MP)+POMB*DKSIDAX(MP)+SINP*DSDAX
      DBDAY(MP)=DBDAY(MP)+POMB*DKSIDAY(MP)+SINP*DSDAY
      GO TO 2040
 2038 DADA0=DADA0+DSDA0*COSP
      IF(CSYM) DBDA0=DBDA0+DSDA0*SINP
 2040 CONTINUE
 2050 CONTINUE
 2060 CONTINUE
 2065 CONTINUE
 2066 CONTINUE
 2067 CONTINUE
 2068 CONTINUE
 2069 CONTINUE
C
      GO TO 2075
 2070 SINIJ=SIN(FI)
      COSIJ=COS(FI)
 2075 IF(IZO) GO TO 2076
      EXPIJ=EXP(-HQJ*B1-KQJ*B2-LQJ*B3-HK2J*B4-HL2J*B5-KL2J*B6)
      COSIJ=COSIJ*EXPIJ
      SINIJ=SINIJ*EXPIJ
C
 2076 IF(HARM) GO TO 2080
      DO 2077KK=1,KMODXI
      DADU1(KK)=EXPIJ*DADU1(KK)
      DADU2(KK)=EXPIJ*DADU2(KK)
      IF(KK.EQ.KMODXI.AND.KFXI.NE.0) THEN
      DADU3=EXPIJ*DADU3
      DADU4=EXPIJ*DADU4
      ENDIF
      IF(CSYM) GO TO 2077
      DBDU1(KK)=EXPIJ*DBDU1(KK)
      DBDU2(KK)=EXPIJ*DBDU2(KK)
      IF(KK.EQ.KMODXI.AND.KFXI.NE.0) THEN
      DBDU3=EXPIJ*DBDU3
      DBDU4=EXPIJ*DBDU4
      ENDIF
 2077 CONTINUE
      DO 12077KK=1,KMODBI
      DADB1(KK)=EXPIJ*DADB1(KK)
      DADB2(KK)=EXPIJ*DADB2(KK)
      IF(KK.EQ.KMODBI.AND.KFBI.NE.0) THEN
      DADB3=EXPIJ*DADB3
      DADB4=EXPIJ*DADB4
      ENDIF
      IF(CSYM) GO TO 12077
      DBDB1(KK)=EXPIJ*DBDB1(KK)
      DBDB2(KK)=EXPIJ*DBDB2(KK)
      IF(KK.EQ.KMODBI.AND.KFBI.NE.0) THEN
      DBDB3=EXPIJ*DBDB3
      DBDB4=EXPIJ*DBDB4
      ENDIF
12077 CONTINUE
      DO 2078KK=1,KMODSI
      DADAX(KK)=EXPIJ*DADAX(KK)
      DADAY(KK)=EXPIJ*DADAY(KK)
      IF(CSYM) GO TO 2078
      DBDAX(KK)=EXPIJ*DBDAX(KK)
      DBDAY(KK)=EXPIJ*DBDAY(KK)
 2078 CONTINUE
 2080 DADFII=DADFII+COSIJ
      IF(CSYM) GO TO 2090
      DBDFII=DBDFII+SINIJ
 2090 DADXI=DADXI-HJ*SINIJ
      DADYI=DADYI-KJ*SINIJ
      DADZI=DADZI-LJ*SINIJ
      IF(IZO) GO TO 2100
      DADB1I=DADB1I-HQJ *COSIJ
      DADB2I=DADB2I-KQJ *COSIJ
      DADB3I=DADB3I-LQJ*COSIJ
      DADB4I=DADB4I-HK2J*COSIJ
      DADB5I=DADB5I-HL2J*COSIJ
      DADB6I=DADB6I-KL2J*COSIJ
 2100 IF(CSYM) GO TO 2200
      DBDXI=DBDXI+HJ*COSIJ
      DBDYI=DBDYI+KJ*COSIJ
      DBDZI=DBDZI+LJ*COSIJ
      IF(IZO) GO TO 2200
      DBDB1I=DBDB1I-HQJ*SINIJ
      DBDB2I=DBDB2I-KQJ*SINIJ
      DBDB3I=DBDB3I-LQJ *SINIJ
      DBDB4I=DBDB4I-HK2J*SINIJ
      DBDB5I=DBDB5I-HL2J*SINIJ
      DBDB6I=DBDB6I-KL2J*SINIJ
C
 2200 IF(HARM) GO TO 5000
      DO 3000KK=1,KMODXI
      DADUXX(KK)=DADUXX(KK)+HJ*DADU1(KK)
      DADUXY(KK)=DADUXY(KK)+KJ*DADU1(KK)
      DADUXZ(KK)=DADUXZ(KK)+LJ*DADU1(KK)
      IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN
      DADUYX(KK)=DADUYX(KK)+HJ*DADU2(KK)
      DADUYY(KK)=DADUYY(KK)+KJ*DADU2(KK)
      DADUYZ(KK)=DADUYZ(KK)+LJ*DADU2(KK)
      ELSE
      DADUYX(KK)=DADUYX(KK)+DADU4
      DADUYY(KK)=DADUYY(KK)+DADU2(KK)
      DADUYZ(KK)=DADUYZ(KK)+DADU3
      ENDIF
      IF(CSYM) GO TO 3000
      DBDUXX(KK)=DBDUXX(KK)+HJ*DBDU1(KK)
      DBDUXY(KK)=DBDUXY(KK)+KJ*DBDU1(KK)
      DBDUXZ(KK)=DBDUXZ(KK)+LJ*DBDU1(KK)
      IF(KK.NE.KMODXI.OR.KFXI.EQ.0) THEN
      DBDUYX(KK)=DBDUYX(KK)+HJ*DBDU2(KK)
      DBDUYY(KK)=DBDUYY(KK)+KJ*DBDU2(KK)
      DBDUYZ(KK)=DBDUYZ(KK)+LJ*DBDU2(KK)
      ELSE
      DBDUYX(KK)=DBDUYX(KK)+DBDU4
      DBDUYY(KK)=DBDUYY(KK)+DBDU2(KK)
      DBDUYZ(KK)=DBDUYZ(KK)+DBDU3
      ENDIF
 3000 CONTINUE
      DO 4000KK=1,KMODBI
      DADBX1(KK)=DADBX1(KK)-HQJ*DADB1(KK)
      DADBX2(KK)=DADBX2(KK)-KQJ*DADB1(KK)
      DADBX3(KK)=DADBX3(KK)-LQJ*DADB1(KK)
      DADBX4(KK)=DADBX4(KK)-HK2J*DADB1(KK)
      DADBX5(KK)=DADBX5(KK)-HL2J*DADB1(KK)
      DADBX6(KK)=DADBX6(KK)-KL2J*DADB1(KK)
      IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN
      DADBY1(KK)=DADBY1(KK)-HQJ*DADB2(KK)
      DADBY2(KK)=DADBY2(KK)-KQJ*DADB2(KK)
      DADBY3(KK)=DADBY3(KK)-LQJ*DADB2(KK)
      DADBY4(KK)=DADBY4(KK)-HK2J*DADB2(KK)
      DADBY5(KK)=DADBY5(KK)-HL2J*DADB2(KK)
      DADBY6(KK)=DADBY6(KK)-KL2J*DADB2(KK)
      ELSE
      DADBY1(KK)=DADBY1(KK)+DADB4
      DADBY2(KK)=DADBY2(KK)+DADB2(KK)
      DADBY3(KK)=DADBY3(KK)+DADB3
      ENDIF
      IF(CSYM) GO TO 4000
      DBDBX1(KK)=DBDBX1(KK)-HQJ*DBDB1(KK)
      DBDBX2(KK)=DBDBX2(KK)-KQJ*DBDB1(KK)
      DBDBX3(KK)=DBDBX3(KK)-LQJ*DBDB1(KK)
      DBDBX4(KK)=DBDBX4(KK)-HK2J*DBDB1(KK)
      DBDBX5(KK)=DBDBX5(KK)-HL2J*DBDB1(KK)
      DBDBX6(KK)=DBDBX6(KK)-KL2J*DBDB1(KK)
      IF(KK.NE.KMODBI.OR.KFBI.EQ.0) THEN
      DBDBY1(KK)=DBDBY1(KK)-HQJ*DBDB2(KK)
      DBDBY2(KK)=DBDBY2(KK)-KQJ*DBDB2(KK)
      DBDBY3(KK)=DBDBY3(KK)-LQJ*DBDB2(KK)
      DBDBY4(KK)=DBDBY4(KK)-HK2J*DBDB2(KK)
      DBDBY5(KK)=DBDBY5(KK)-HL2J*DBDB2(KK)
      DBDBY6(KK)=DBDBY6(KK)-KL2J*DBDB2(KK)
      ELSE
      DBDBY1(KK)=DBDBY1(KK)+DBDB4
      DBDBY2(KK)=DBDBY2(KK)+DBDB2(KK)
      DBDBY3(KK)=DBDBY3(KK)+DBDB3
      ENDIF
 4000 CONTINUE
 5000 CONTINUE
      PH=FXI
      IF(PHFP.EQ.0..OR.PHF(I).EQ.0..OR.HARM) GO TO 5100
      PH=FXI*EXP(PHFP*PHF(I))
 5100 IF(IZO) GO TO 5200
      FTA=AII*PH
      FTACOS=FTA*DADFII
      DER(M  )=PH*DADFII
      DER(M+4)=DADB1I*FTA
      DER(M+5)=DADB2I*FTA
      DER(M+6)=DADB3I*FTA
      DER(M+7)=DADB4I*FTA
      DER(M+8)=DADB5I*FTA
      DER(M+9)=DADB6I*FTA
      GO TO 5300
 5200 TFI=EXP(-RHO*B1)
      FTA=PH*AII*TFI
      FTACOS=FTA*DADFII
      DADFII=TFI*DADFII
      DER(M  )=PH*DADFII
      DER(M+4)=-RHO*FTACOS
 5300 A=A+FTACOS
      AANOM=AANOM+FYI*FTACOS
      DER(M+1)=DADXI*FTA
      DER(M+2)=DADYI*FTA
      DER(M+3)=DADZI*FTA
      IF(HARM) GO TO 5400
      IF(KMODSI.LE.0) GO TO 5311
      DER(M+10)=DADA0*FTA
      LL=M+9
      DO 5310KK=1,KMODSI
      LL=LL+2
      DER(LL)  =DADAX(KK)*FTA
 5310 DER(LL+1)=DADAY(KK)*FTA
 5311 LL=M+13
      DO 5350KK=1,KMODXI
      LL=LL+6
      DER(LL  )=DADUXX(KK)*FTA
      DER(LL+1)=DADUXY(KK)*FTA
      DER(LL+2)=DADUXZ(KK)*FTA
      DER(LL+3)=DADUYX(KK)*FTA
      DER(LL+4)=DADUYY(KK)*FTA
 5350 DER(LL+5)=DADUYZ(KK)*FTA
      LL=M+31
      DO 5360KK=1,KMODBI
      LL=LL+12
      DER(LL   )=DADBX1(KK)*FTA
      DER(LL+1 )=DADBX2(KK)*FTA
      DER(LL+2 )=DADBX3(KK)*FTA
      DER(LL+3 )=DADBX4(KK)*FTA
      DER(LL+4 )=DADBX5(KK)*FTA
      DER(LL+5 )=DADBX6(KK)*FTA
      DER(LL+6 )=DADBY1(KK)*FTA
      DER(LL+7 )=DADBY2(KK)*FTA
      DER(LL+8 )=DADBY3(KK)*FTA
      DER(LL+9 )=DADBY4(KK)*FTA
      DER(LL+10)=DADBY5(KK)*FTA
 5360 DER(LL+11)=DADBY6(KK)*FTA
      DER(M+91)=PHFP*FTACOS
 5400 IF(CSYM) GO TO 6000
      FTASIN=FTA*DBDFII
      B=B+FTASIN
      BANOM=BANOM+FYI*FTASIN
      DC(M+1)=DBDXI*FTA
      DC(M+2)=DBDYI*FTA
      DC(M+3)=DBDZI*FTA
      IF(IZO) GO TO 5600
      DC(M   )=PH*DBDFII
      DC(M+4 )=DBDB1I*FTA
      DC(M+5 )=DBDB2I*FTA
      DC(M+6 )=DBDB3I*FTA
      DC(M+7 )=DBDB4I*FTA
      DC(M+8 )=DBDB5I*FTA
      DC(M+9 )=DBDB6I*FTA
      GO TO 5700
 5600 DBDFII=TFI*DBDFII
      DC(M   )=PH*DBDFII
      DC(M+4 )=-RHO*FTASIN
 5700 IF(HARM) GO TO 6000
      DC(M+10)=DBDA0*FTA
      LL=M+9
      DO 5710KK=1,KMODSI
      LL=LL+2
      DC(LL)  =DBDAX(KK)*FTA
 5710 DC(LL+1)=DBDAY(KK)*FTA
      LL=M+13
      DO 5750KK=1,KMODXI
      LL=LL+6
      DC(LL  )=DBDUXX(KK)*FTA
      DC(LL+1)=DBDUXY(KK)*FTA
      DC(LL+2)=DBDUXZ(KK)*FTA
      DC(LL+3)=DBDUYX(KK)*FTA
      DC(LL+4)=DBDUYY(KK)*FTA
 5750 DC(LL+5)=DBDUYZ(KK)*FTA
      LL=M+31
      DO 5460KK=1,KMODBI
      LL=LL+12
      DC(LL   )=DBDBX1(KK)*FTA
      DC(LL+1 )=DBDBX2(KK)*FTA
      DC(LL+2 )=DBDBX3(KK)*FTA
      DC(LL+3 )=DBDBX4(KK)*FTA
      DC(LL+4 )=DBDBX5(KK)*FTA
      DC(LL+5 )=DBDBX6(KK)*FTA
      DC(LL+6 )=DBDBY1(KK)*FTA
      DC(LL+7 )=DBDBY2(KK)*FTA
      DC(LL+8 )=DBDBY3(KK)*FTA
      DC(LL+9 )=DBDBY4(KK)*FTA
      DC(LL+10)=DBDBY5(KK)*FTA
 5460 DC(LL+11)=DBDBY6(KK)*FTA
      DC(M+91)=PHFP*FTASIN
 6000 M=M+92
      A=A-BANOM
      B=B+AANOM
      A=A*CENTR
      B=B*CENTR
      SQRTAB=SQRT(A**2+B**2)
      DER(IQ)=SQRTAB
      YC=SC(IQ)*SQRTAB
      IF(SQRTAB.EQ.0.) GO TO 6100
      SQRTAB=SC(IQ)/SQRTAB
 6100 C1=A*SQRTAB*CENTR
      C2=B*SQRTAB*CENTR
      M=18
      IF(CSYM) GO TO 9500
      DO 9200I=1,NA
      KMODXI=KMODX(I)
      KMODBI=KMODB(I)
      KMODSI=KMODS(I)
      IF(KMODXI.LT.0) KMODXI=KMODXM(-KMODXI)
      J=ISF(I)
      FYI=FY(J)
      CA=C1+FYI*C2
      CB=C2-FYI*C1
      J=ITF(I)*5
      DO 9100N=1,J
      M=M+1
 9100 DER(M)=CA*DER(M)+CB*DC(M)
      IF(J.NE.10) M=M+5
      IF(KMODSI.GT.0) GO TO 9120
      M=M+9
      GO TO 9130
 9120 M=M+1
      DER(M)=CA*DER(M)+CB*DC(M)
      KK=KMODSI*2
      DO 9125N=1,KK
      M=M+1
 9125 DER(M)=CA*DER(M)+CB*DC(M)
      M=M+8-KK
 9130 IF(KMODXI.GT.0) GO TO 9140
      M=M+24
      GO TO 9160
 9140 KK=KMODXI*6
      DO 9150N=1,KK
      M=M+1
 9150 DER(M)=CA*DER(M)+CB*DC(M)
      M=M+24-KK
 9160 IF(KMODBI.GT.0) GO TO 9170
      M=M+49
      GO TO 9190
 9170 KK=KMODBI*12
      DO 9180N=1,KK
      M=M+1
 9180 DER(M)=CA*DER(M)+CB*DC(M)
      M=M+49-KK
 9190 DER(M)=CA*DER(M)+CB*DC(M)
 9200 CONTINUE
      GO TO 9900
 9500 DO 9300I=1,NA
      KMODXI=KMODX(I)
      KMODBI=KMODB(I)
      KMODSI=KMODS(I)
      IF(KMODXI.LT.0) KMODXI=KMODXM(-KMODXI)
      J=ISF(I)
      CA=C1+FY(J)*C2
      J=ITF(I)*5
      DO 9400N=1,J
      M=M+1
 9400 DER(M)=CA*DER(M)
      IF(J.NE.10) M=M+5
      IF(KMODSI.GT.0) GO TO 9420
      M=M+9
      GO TO 9430
 9420 M=M+1
      DER(M)=CA*DER(M)
      KK=KMODSI*2
      DO 9425N=1,KK
      M=M+1
 9425 DER(M)=CA*DER(M)
      M=M+8-KK
 9430 IF(KMODXI.GT.0) GO TO 9440
      M=M+24
      GO TO 9460
 9440 KK=KMODXI*6
      DO 9450N=1,KK
      M=M+1
 9450 DER(M)=CA*DER(M)
      M=M+24-KK
 9460 IF(KMODBI.GT.0) GO TO 9470
      M=M+49
      GO TO 9490
 9470 KK=KMODBI*12
      DO 9480N=1,KK
      M=M+1
 9480 DER(M)=CA*DER(M)
      M=M+49-KK
 9490 DER(M)=CA*DER(M)
 9300 CONTINUE
C
 9900 IF(IEXT.LE.0) GO TO 9999
      CALL PHILIP(YC,EXTINC,EXTINK,DER(7),SINTHL,SC(IQ))
      EXTKOR=SQRT(EXTINC)
      YC=YC*EXTKOR
      DER(IQ)=DER(IQ)*EXTKOR
      FI=EXTKOR*EXTINK/EXTINC
      DO 9910I=19,NP
 9910 DER(I)=DER(I)*FI
 9999 RETURN
      END
      SUBROUTINE PHILIP(F,EXTINC,EXTINK,DFDC,SINTHL,SC)
      DIMENSION DFDC(12)
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      YY(A,B,C,FX)=1./(1.+C*FX+A*FX**2/(1.+B*FX))
      ZZ(A,B,C,FX,FY)=-0.5*FY*(C+A*FX*(2.+B*FX)/(1.+B*FX)**2)
      DO 30I=1,12
   30 DFDC(I)=0.
      IF(F.LE.0.0) RETURN
      FSQR=(F/SC)**2
      TBAR=EFPIP(7)/VOLUME
      GAMMA=1.5E-4*FLAME/VOLUME
      HROOT=EC(1)
      IF(ITYPEX.EQ.1) HROOT=EC(7)
      IF(IEXT.EQ.1) GO TO 4
      N=0
      M=6
      IF(IEXT.NE.3) GO TO 3
      N=3
      M=0
    3 HROOT=EFPIP(N+1)**2*EC(M+1)+EFPIP(N+2)**2*EC(M+2)+EFPIP(N+3)**2
     1      *EC(M+3)+2.*(EFPIP(N+1)*EFPIP(N+2)*EC(M+4)+EFPIP(N+1)*
     2      EFPIP(N+3)*EC(M+5)+EFPIP(N+2)*EFPIP(N+3)*EC(M+6))
      HROOT=SQRT(HROOT)
      IF(IEXT.EQ.3.OR.IMOSAIC.EQ.1) HROOT=1./HROOT
    4 SINT=SINTHL*FLAME
      TH=ASIN(SINT)/TORAD
      COST=SQRT(1.-SINT**2)
      IF(IAB.EQ.0) GO TO 41
      I=TH*.2
      XI=TH*.2-FLOAT(I)
      I=I+1
      CP=.6666667*(DADMI(I)+(DADMI(I+1)-DADMI(I))*XI)
      TBAR=TBAR*CP
      GAMMA=GAMMA*CP
   41 SIN2T=2.*SINT*COST
      COS2T=COST**2-SINT**2
      AP=0.20+0.45*COS2T
      BP=0.22-0.12*(0.5-COS2T)**2
      CP=2.
      IF(IAB.EQ.0) GO TO 42
      I=SINT*10.
      XI=SINT*10.-FLOAT(I)
      I=I+1
      AG=AA(I)+(AA(I+1)-AA(I))*XI
      BG=BA(I)+(BA(I+1)-BA(I))*XI
      AL=AG
      BL=BG
      GO TO 43
   42 AG=0.58+0.48*COS2T+0.24*COS2T**2
      BG=0.020-0.025*COS2T
      AL=0.025+0.285*COS2T
      BL=0.15-0.2*(0.75-COS2T)**2
      IF(COS2T.LT.0.) BL=-0.45*COS2T
   43 IF(ITYPEX.NE.2) GO TO 130
      AS=AP
      BS=BP
      CS=2.
      GO TO 140
  130 AS=AG
      BS=BG
      CS=2.12
      IF(IDISTR.EQ.1) GO TO 140
      AS=AL
      BS=BL
      CS=2.
  140 IF(ITYPEX.NE.2) GO TO 100
C  HERE FOR TYPE 2
      RHO=HROOT
      BETA=RHO
      DBRHO=1.
      GO TO 120
  100 IF(ITYPEX.NE.1)GO TO 110
C  HERE FOR TYPE 1
      G=HROOT
      BETA=G/SIN2T
      DBG=1./SIN2T
      RHO=0.
      DBRHO=0.
      GO TO 120
C  HERE FOR GENERAL CASE
  110 RHO=EC(1)
      G=EC(7)
      IF(IEXT.EQ.1) GO TO 150
      G=HROOT
  150 BETA=RHO/SQRT(1.+(RHO*SIN2T/G)**2)
      DBRHO=(BETA/RHO)**3
      DBG=(BETA/G)**3*SIN2T**2
      IF(IDISTR.EQ.1)GO TO 120
      BETA=RHO/(1.+RHO*SIN2T/G)
      DBRHO=(BETA/RHO)**2
      DBG=(BETA/G)**2*SIN2T
  120 CONTINUE
      X2=TBAR*BETA*FSQR
      X1=0.
C  INTRODUCE PRIMARY EXTINCTION HERE FOR GENERAL CASE
      IF(ITYPEX.EQ.3) X1=GAMMA*FSQR*RHO**2
      XP=X1
      YP=YY(AP,BP,CP,XP)
      ZP=ZZ(AP,BP,CP,XP,YP)
      YP=SQRT(YP)
      IF(ITYPEX.NE.3) ZP=0.
      XS=X2*YP
      YS=YY(AS,BS,CS,XS)
      ZS=ZZ(AS,BS,CS,XS,YS)
      YS=SQRT(YS)
C  XRAY CASE TREATED HERE
      POL=COS2T**2
      FACT=1./(C2THMQ+POL)
      XPP=X1*POL
      YPP=YY(AP,BP,CP,XPP)
      ZPP=ZZ(AP,BP,CP,XPP,YPP)
      IF(ITYPEX.NE.3) ZPP=0.
      YPP=SQRT(YPP)
      XPS=YPP*X2*POL
      YPS=YY(AS,BS,CS,XPS)
      ZPS=ZZ(AS,BS,CS,XPS,YPS)
      YPS=SQRT(YPS)
      YPYS=C2THMQ*YP*YS
      YPPYPS=POL*YPP*YPS
      EXTINC=FACT*(YPYS+YPPYPS)
      OM1=FACT*(YPYS*(1.+XP*ZP+XS*ZS*(1.+XP*ZP))+YPPYPS*
     1                            (1.+XPP*ZPP+XPS*ZPS*(1.+XPP*ZPP)))
      OM2=FACT*(YPYS*ZS*YP+POL*YPPYPS*ZPS*YPP)
      OM3=FACT*(YPYS*(2.*GAMMA*RHO*ZP+ZS*YP*(TBAR*DBRHO +2.*X2*GAMMA*
     1    RHO*ZP))+POL*YPPYPS*(2.*GAMMA*RHO*ZPP+ZPS*YPP*(TBAR*DBRHO
     2    +2.*POL*X2*GAMMA*RHO*ZPP)))
      EXTINK=OM1
      TUMS=0.5*(SC*FSQR)**2/(F*SQRT(EXTINC))
      TUMS1=TUMS*OM3
      TUMS2=TUMS*TBAR*OM2*DBG
      DFDC(1)=TUMS1
      IF(IEXT.NE.1) GO TO 7
      DFDC(7)=TUMS2
      GO TO 70
    7 DIVID=0.5/HROOT
      IF(IEXT.EQ.3.OR.IMOSAIC.EQ.1) DIVID=-0.5*HROOT**3
      IF(ITYPEX.EQ.2) TUMS2=TUMS1
      DFDC(M+1)=TUMS2*EFPIP(N+1)**2*DIVID
      DFDC(M+2)=TUMS2*EFPIP(N+2)**2*DIVID
      DFDC(M+3)=TUMS2*EFPIP(N+3)**2*DIVID
      DFDC(M+4)=2.*TUMS2*EFPIP(N+1)*EFPIP(N+2)*DIVID
      DFDC(M+5)=2.*TUMS2*EFPIP(N+1)*EFPIP(N+3)*DIVID
      DFDC(M+6)=2.*TUMS2*EFPIP(N+2)*EFPIP(N+3)*DIVID
   70 RETURN
      END
      SUBROUTINE SUMA
      COMMON/REFL/  YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708),
     1              DC(9708),DYP,IQ
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708),
     1       AM(75000)
      LOGICAL NULOVA
      DIMENSION IPSA(9)
      DO 500IP=1,9
  500 IPSA(IP)=NPSS(IP)
      DO 700IP=1,NPMK
      N=KI(IP)
      IF(N.EQ.0) GO TO 700
      IPS=IPSA(N)+1
      DC(IPS)=DER(IP)
      IPSA(N)=IPS
  700 CONTINUE
      IM=0
      IPK=0
      DO 4000N=1,9
      NPSN=NPS(N)
      IF(NPSN.LE.0) GO TO 4000
      IPZ=IPK+1
      IPK=IPZ+NPSN-1
      DO 3000IP=IPZ,IPK
      DERI=DC(IP)*WT
      IF(DERI.NE.0.0) GO TO 1000
      IM=IM+IP-IPZ+1
      GO TO 3000
 1000 DO 2000JP=IPZ,IP
      IM=IM+1
      DERJ=DC(JP)
      IF(DERJ.EQ.0.) GO TO 2000
      AM(IM)=AM(IM)+DERI*DERJ
 2000 CONTINUE
      PS(IP)=PS(IP)+DERI*DYP
 3000 CONTINUE
 4000 CONTINUE
      RETURN
      END
      SUBROUTINE INVER
C
      COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF,
     1              CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ
      COMMON/UNITS/ W,M40,M50
      COMMON/REFL/  YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708),
     1              DC(9708),DYP,IQ
      COMMON/CM4/   ATOM(100),AI(100),ISF(100),ITF(100),X(3,100),
     1              BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100)
     2             ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100),
     3              AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100),
     4              KFX(100),KFB(100)
      COMMON/CSM4/  SAI(100),SX(3,100),SB(6,100),SUX(3,4,100),
     1              SUY(3,4,100),SPHF(100),SA0(100),SAX(4,100),
     2              SAY(4,100),SBX(6,4,100),SBY(6,4,100)
      COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      COMMON/VAZBY/ NVX,NX(10),NVB,NB(10),RMX(3,3,10),RMB(6,6,10),
     1              NVAI,NAI(2,10),NAIXB(10),SUMAI(10)
      COMMON/MODUL/ QUI(3),QUIR(3)
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON/STALY/ PI2
      COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708),
     1       AM(75000)
      EQUIVALENCE (NPSS(10),NPSM)
      DIMENSION ISING(9),ECP(12),SCP(6)
      LOGICAL CSYM
      INTEGER W
      CHARACTER*1 CU,CV,CW
      CHARACTER*8 ATOM,AT,MOLNAME
      CHARACTER*12 L1,L2,L5,L7,L8,L9,L10,LP
      DATA AT/'        '/,CU,CV,CW/'U','V','W'/
      IPS=1
      IAM=1
      Q=Q/FLOAT(NZ-NPSM-1)
      DO 3000I=1,9
      N=NPS(I)
      IF(N.LE.0) GO TO 3000
C
      CALL SMI(AM(IAM),N,ISING(I),I)
      CALL SKOVEJ(AM(IAM),Q,I,KI)
      IF(ISING(I).EQ.0) GO TO 1500
C
      KK=IPS+N-1
      DO 1200K=IPS,KK
 1200 DC(K)=0.
      WRITE(W,101) I
      GO TO 2800
 1500 CALL NASOB(AM(IAM),PS(IPS),DC(IPS),N)
 2800 IAM=IAM+NAM(I)
      IPS=IPS+N
 3000 CONTINUE
      Q=SQRT(Q)
      DO 3100I=1,NPMK
      PS(I)=0.
 3100 DER(I)=0.
      IAM=0
      IPS=0
      DO 3500I=1,9
      N=NPS(I)
      IF(N.EQ.0) GO TO 3500
      IF(ISING(I).EQ.0) GO TO 3200
      IAM=IAM+NAM(I)
      IPS=IPS+N
      GO TO 3500
 3200 K=0
C
      DO 3300IP=1,NPMK
      IF(KI(IP).NE.I) GO TO 3300
      K=K+1
      IAM=IAM+K
      IPS=IPS+1
      PS(IP)=DC(IPS)
      DER(IP)=SQRT(AM(IAM))*Q
 3300 CONTINUE
 3500 CONTINUE
      IF(NVX.LE.0) GO TO 3525
      DO 3521I=1,NVX
      K=NX(I)*92-72
 3521 CALL ZMZPET(PS(K),DER(K),RMX(1,1,I),3)
 3525 IF(NVB.LE.0) GO TO 3529
      DO 3526I=1,NVB
      K=NB(I)*92-69
 3526 CALL ZMZPET(PS(K),DER(K),RMB(1,1,I),6)
 3529 ISPASA=0
      TLUMP=TLUM
      DO 3530I=1,6
 3530 SCP(I)=SC(I)
      DO 3501I=1,12
 3501 ECP(I)=EC(I)
 3503 CALL ZMENA(L1,POM,POM,1,IT,AT,IP)
      IT=1
      CALL PIS(IT)
      DO 3504I=1,6
      ENCODE(12,100,LP) I
 3504 CALL ZMENA(LP,SC(I),POM,2,IT,AT,IP)
      IF(IEXT.LE.0) GO TO 3690
      IF(ITYPEX.EQ.1) GO TO 3600
      IF(IEXT.EQ.3) GO TO 3550
      CALL ZMENA(L7,EC(1),POM,2,IT,AT,IP)
      IF(EC(1).LT.0.) GO TO 3680
      IF(ITYPEX.EQ.3) GO TO 3600
      GO TO 3690
 3550 DO 3560J=1,6
      CALL INDEX(J,K,L)
      ENCODE(12,110,LP) K,L
 3560 CALL ZMENA(LP,EC(J),POM,2,IT,AT,IP)
      IF(EC(1).LE.0.) GO TO 3680
      IF(EC(1)*EC(2)-EC(4)**2.LE.0.) GO TO 3680
      IF(EC(1)*EC(2)*EC(3)+2.*EC(4)*EC(5)*EC(6)-EC(1)*EC(6)**2
     1  -EC(2)*EC(5)**2-EC(3)*EC(4)**2.LE.0.) GO TO 3680
      GO TO 3690
 3600 IP=12
      IF(IEXT.GT.1) GO TO 3650
      CALL ZMENA(L8,EC(7),POM,2,IT,AT,IP)
      IF(EC(7).LE.0.) GO TO 3680
      GO TO 3690
 3650 DO 3660J=1,6
      CALL INDEX(J,K,L)
      ENCODE(12,111,LP) K,L
 3660 CALL ZMENA(LP,EC(J+6),POM,2,IT,AT,IP)
      IF(EC(7).LE.0.) GO TO 3680
      IF(EC(7)*EC(8)-EC(10)**2.LE.0.) GO TO 3680
      IF(EC(7)*EC(8)*EC(9)+2.*EC(10)*EC(11)*EC(12)-EC(7)*EC(12)**2
     1  -EC(8)*EC(11)**2-EC(9)*EC(10)**2.LE.0.) GO TO 3680
      GO TO 3690
 3680 ISPASA=ISPASA+1
      DO 13680I=1,6
13680 SC(I)=SCP(I)
      DO 3681I=1,12
 3681 EC(I)=ECP(I)
      TLUM=TLUM*.5
      IF(ISPASA.GT.5) STOP
      WRITE(W,109) TLUMP,TLUM
      GO TO 3503
C
C .... CALCULATION OF CHANGES OF ATOMIC PARAMETERS
C
 3690 IP=18
      DO 4000I=1,NA
      DO 3691K=1,NVAI
      IF(I.EQ.NAI(2,K)) GO TO 3692
 3691 CONTINUE
      GO TO 3699
 3692 N=NAI(1,K)
      IF(NAIXB(K).LT.0) GO TO 13692
      AI(I)=SUMAI(K)-AI(N)
      GO TO 33692
13692 A0(I)=SUMAI(K)-A0(N)
      DO 23692KK=1,KMODS(I)
      AX(KK,I)=-AX(KK,N)
23692 AY(KK,I)=-AY(KK,N)
33692 IF(IABS(NAIXB(K)).EQ.3) GO TO 3699
      DO 3693KK=1,3
 3693 X(KK,I)=X(KK,N)
      IF(IABS(NAIXB(K)).EQ.2) GO TO 3699
      DO 3694KK=1,6
 3694 BETA(KK,I)=BETA(KK,N)
      DO 3696KK=1,KMODX(I)
      DO 3695L=1,3
      UX(L,KK,I)=UX(L,KK,N)
 3695 UY(L,KK,I)=UY(L,KK,N)
 3696 CONTINUE
      DO 3698KK=1,KMODB(I)
      DO 3697L=1,6
      BX(L,KK,I)=BX(L,KK,N)
 3697 BY(L,KK,I)=BY(L,KK,N)
 3698 CONTINUE
      PHF(I)=PHF(N)
C
 3699 CALL PIS(IT)
      WRITE(W,102) ATOM(I)
      CALL ZMENA(L1,AI(I),SAI(I),2,IT,ATOM(I),IP)
      CALL ZMENA(L2(1),X(1,I),SX(1,I),2,IT,ATOM(I),IP)
      CALL ZMENA(L2(2),X(2,I),SX(2,I),2,IT,ATOM(I),IP)
      CALL ZMENA(L2(3),X(3,I),SX(3,I),2,IT,ATOM(I),IP)
      IF(ITF(I).NE.1) GO TO 3700
      CALL ZMENA(L5,BETA(1,I),SB(1,I),2,IT,ATOM(I),IP)
      IP=IP+5
      GO TO 3800
 3700 DO 3710J=1,6
      CALL INDEX(J,K,L)
      IF(LITE.EQ.0) GO TO 3705
      ENCODE(12,112,LP) K,L
      GO TO 3710
 3705 ENCODE(12,117,LP) K,L     
 3710 CALL ZMENA(LP,BETA(J,I),SB(J,I),2,IT,ATOM(I),IP)
C
 3800 IF(KMODS(I).LE.0) GO TO 3808
      CALL ZMENA(L10,A0(I),SA0(I),2,IT,ATOM(I),IP)
      DO 3802J=1,4
      IF(J.GT.KMODS(I)) GO TO 3801
      ENCODE(12,118,LP) J
      CALL ZMENA(LP,AX(J,I),SAX(J,I),2,IT,ATOM(I),IP)
      ENCODE(12,119,LP) J
      CALL ZMENA(LP,AY(J,I),SAY(J,I),2,IT,ATOM(I),IP)
      GO TO 3802
 3801 IP=IP+2
 3802 CONTINUE
      GO TO 3809
 3808 IP=IP+9
 3809 IF(KMODX(I).LE.0.AND.KMODB(I).LE.0) GO TO 3995
      DO 3900J=1,4
      IF(J.GT.KMODX(I)) GO TO 3850
      DO 3810K=1,3
      ENCODE(12,113,LP) CU,K,J
 3810 CALL ZMENA(LP,UX(K,J,I),SUX(K,J,I),2,IT,ATOM(I),IP)
      DO 3820K=1,3
      ENCODE(12,114,LP) CU,K,J
 3820 CALL ZMENA(LP,UY(K,J,I),SUY(K,J,I),2,IT,ATOM(I),IP)
      GO TO 3900
 3850 IP=IP+6
 3900 CONTINUE
      DO 3990J=1,4
      IF(J.GT.KMODB(I)) GO TO 3950
      DO 3910K=1,6
      ENCODE(12,113,LP) 'T',K,J
 3910 CALL ZMENA(LP,BX(K,J,I),SBX(K,J,I),2,IT,ATOM(I),IP)
      DO 3920K=1,6
      ENCODE(12,114,LP) 'T',K,J
 3920 CALL ZMENA(LP,BY(K,J,I),SBY(K,J,I),2,IT,ATOM(I),IP)
      GO TO 3990
 3950 IP=IP+12
 3990 CONTINUE
      CALL ZMENA(L9,PHF(I),SPHF(I),2,IT,ATOM(I),IP)
      GO TO 4000
 3995 IP=IP+73
 4000 CONTINUE
      CALL RESETX(X)
      CALL RESETB(BETA)
      CALL RESETA(A0,AX,AY)
      CALL RESETD(UX,UY)
      CALL RESETC(BX,BY)
      IF(NMOL.EQ.0) GO TO 4200
      DO 4100I=1,NMOL
      CALL PIS(IT)
      WRITE(W,102) MOLNAME(I)
      IF(KMODXM(I).LE.0) GO TO 4090
      DO 4040J=1,4
      IF(J.GT.KMODXM(I)) GO TO 4030
      DO 4010K=1,3
      ENCODE(12,113,LP) CV,K,J
 4010 CALL ZMENA(LP,UTX(K,J,I),SUTX(K,J,I),2,IT,MOLNAME(I),IP)
      DO 4020K=1,3
      ENCODE(12,114,LP) CV,K,J
 4020 CALL ZMENA(LP,UTY(K,J,I),SUTY(K,J,I),2,IT,MOLNAME(I),IP)
      GO TO 4040
 4030 IP=IP+6
 4040 CONTINUE
      DO 4080J=1,4
      IF(J.GT.KMODXM(I)) GO TO 4070
      DO 4050K=1,3
      ENCODE(12,113,LP) CW,K,J
 4050 CALL ZMENA(LP,URX(K,J,I),SURX(K,J,I),2,IT,MOLNAME(I),IP)
      DO 4060K=1,3
      ENCODE(12,114,LP) CW,K,J
 4060 CALL ZMENA(LP,URY(K,J,I),SURY(K,J,I),2,IT,MOLNAME(I),IP)
      GO TO 4080
 4070 IP=IP+6
 4080 CONTINUE
      CALL ZMENA(L9,PHFM(I),SPHFM(I),2,IT,MOLNAME(I),IP)
      GO TO 4100
 4090 IP=IP+49
 4100 CONTINUE
      CALL RESETM(UTX,UTY,URX,URY)
 4200 CALL ZMENA(L1,POM,POM,3,IT,AT,IP)
      DO 4500I=1,NA
      IF(BETA(1,I).GT.0.) GO TO 4510
      WRITE(W,107) ATOM(I)
      GO TO 4500
 4510 IF(ITF(I).EQ.1) GO TO 4500
      IF(BETA(1,I)*BETA(2,I)-BETA(4,I)**2.GT.0.) GO TO 4520
      WRITE(W,107) ATOM(I)
      GO TO 4500
 4520 IF(BETA(1,I)*BETA(2,I)*BETA(3,I)+2.*BETA(4,I)*BETA(5,I)*BETA(6,I)
     1  -BETA(3,I)*BETA(4,I)**2-BETA(1,I)*BETA(6,I)**2-BETA(2,I)*
     2   BETA(5,I)**2.GT.0.) GO TO 4500
      WRITE(W,107) ATOM(I)
 4500 CONTINUE
      REWIND M40
      POM=0.0
      WRITE(M40,103) NA,NMOL
      WRITE(M40,108) SC,(KI(I),I=1,6)
      WRITE(M40,104) POM
      WRITE(M40,108)(EC(I),I=1,6),(KI(I),I=7,12),
     1              (EC(I),I=7,12),(KI(I),I=13,18)
      K=18
      DO 5000I=1,NA
      KK=K+1
      K=K+10
      J=IABS(KMODX(I))+10*KMODB(I)
      IF(KMODX(I).LT.0) J=-J
      WRITE(M40,105)  ATOM(I),ISF(I),KFB(I),KFX(I),ITF(I),J,KMODS(I),
     1                AI(I),(X(N,I),N=1,3),(BETA(N,I),N=1,6),
     2                (KI(N),N=KK,K)
      KK=K+1
      K=K+1
      IF(KMODS(I).GT.0) WRITE(M40,104) A0(I),KI(K)
      DO 4800J=1,4
      KK=K+1
      K=K+2
      IF(J.GT.KMODS(I)) GO TO 4800
      WRITE(M40,120) AX(J,I),AY(J,I),KI(KK),KI(K)
 4800 CONTINUE
      DO 4900J=1,4
      KK=K+1
      K=K+6
      IF(J.GT.KMODX(I)) GO TO 4900
      WRITE(M40,108)(UX(N,J,I),N=1,3),(UY(N,J,I),N=1,3),(KI(N),N=KK,K)
 4900 CONTINUE
      DO 4950J=1,4
      KK=K+1
      K=K+12
      IF(J.GT.KMODB(I)) GO TO 4950
      WRITE(M40,108)(BX(N,J,I),N=1,6),(KI(N),N=KK,KK+5),
     1              (BY(N,J,I),N=1,6),(KI(N),N=KK+6,K)
 4950 CONTINUE
      KK=K+1
      K=K+1
      IF(KMODX(I)+KMODS(I)+KMODB(I).GT.0) WRITE(M40,104) PHF(I),KI(KK)
 5000 CONTINUE
      IF(NMOL.EQ.0) GO TO 5200
      DO 5100I=1,NMOL
      CALL CENTER(I)
      QCNTMM=-PI2*((QUI(1)+QUIR(1))*XM(1,I)+(QUI(2)+QUIR(2))*XM(2,I)+
     1             (QUI(3)+QUIR(3))*XM(3,I))
      DO 5003N=1,NA
      IF(KMODX(N).NE.-I) GO TO 5003
      QCNTM(N)=QCNTMM
 5003 CONTINUE
 5005 WRITE(M40,116) MOLNAME(I),NPG(I),KMODXM(I),(XM(N,I),N=1,3)
      IF(NPG(I).GT.1) WRITE(M40,103)(KPG(J,I),J=2,NPG(I))
      DO 5010J=1,4
      KK=K+1
      K=K+6
      IF(J.GT.KMODXM(I)) GO TO 5010
      WRITE(M40,108)(UTX(N,J,I),N=1,3),(UTY(N,J,I),N=1,3),(KI(N),N=KK,K)
 5010 CONTINUE
      DO 5020J=1,4
      KK=K+1
      K=K+6
      IF(J.GT.KMODXM(I)) GO TO 5020
      WRITE(M40,108)(URX(N,J,I),N=1,3),(URY(N,J,I),N=1,3),(KI(N),N=KK,K)
 5020 CONTINUE
      KK=K+1
      K=K+1
      IF(KMODXM(I).GT.0) WRITE(M40,104) PHFM(I),KI(KK)
 5100 CONTINUE
 5200 DO 6000I=1,NA
      WRITE(M40,106) ATOM(I),SAI(I),(SX(N,I),N=1,3),(SB(N,I),N=1,6)
      K=KMODS(I)
      IF(K.LE.0) GO TO 5400
      WRITE(M40,104) SA0(I)
      DO 5300J=1,K
 5300 WRITE(M40,120) SAX(J,I),SAY(J,I)
 5400 DO 5500J=1,KMODX(I)
 5500 WRITE(M40,108)(SUX(N,J,I),N=1,3),(SUY(N,J,I),N=1,3)
      DO 5600J=1,KMODB(I)
      WRITE(M40,108)(SBX(N,J,I),N=1,6)
 5600 WRITE(M40,108)(SBY(N,J,I),N=1,6)
      IF(KMODX(I)+KMODS(I)+KMODB(I).GT.0) WRITE(M40,108) SPHF(I)
 6000 CONTINUE
      IF(NMOL.EQ.0) GO TO 7000
      DO 6100I=1,NMOL
      K=KMODXM(I)
      WRITE(M40,106) MOLNAME(I)
      IF(K.LE.0) GO TO 6100
      DO 6010J=1,K
 6010 WRITE(M40,108)(SUTX(N,J,I),N=1,3),(SUTY(N,J,I),N=1,3)
      DO 6020J=1,K
 6020 WRITE(M40,108)(SURX(N,J,I),N=1,3),(SURY(N,J,I),N=1,3)
      WRITE(M40,108) SPHFM(I)
 6100 CONTINUE
 7000 IF(ISPASA.NE.0) TLUM=TLUMP
  100 FORMAT('SCALE-',I1)
  101 FORMAT(' MATRIX',I2,'.BLOKU IS SINGULAR')
  102 FORMAT('0',A8)
  103 FORMAT(12I5)
  104 FORMAT(F9.6,51X,I1)
  105 FORMAT(A8,I3,3I1,I3,I1,4F9.6/6F9.6,6X,10I1)
  106 FORMAT(A8,10X,4F9.6/6F9.6)
  107 FORMAT('TTEPMERATURE FACTOR ATOMU : ',A8,' IS NOT POSITIVE DEFINIT
     1E')
  108 FORMAT(6F9.6,6X,6I1)
  109 FORMAT('0*** BEWARE ||| ***'/'0DUMPING FACTOR WAS CHANGED FROM'
     1,F10.4,' TO',F10.4/' IN ORDER TO SAVE POSITIVE DEFINITY OF EXTINCT
     2IONAL PARAMETERS')
  110 FORMAT('RHO(',I1,',',I1,')    ')
  111 FORMAT('G(',I1,',',I1,')      ')
  112 FORMAT('BETA(',I1,',',I1,')   ')
  113 FORMAT(A1,I1,'*SIN(',I1,'*X4)')
  114 FORMAT(A1,I1,'*COS(',I1,'*X4)')
  116 FORMAT(A8,3X,2I3,10X,3F9.6)
  117 FORMAT('U(',I1,',',I1,')      ')
  118 FORMAT('A*SIN(',I1,'*X4) ')
  119 FORMAT('A*COS(',I1,'*X4) ')
  120 FORMAT(2F9.6,42X,2I1)
      RETURN
      END
      SUBROUTINE PIS(IT)
      COMMON/UNITS/ W,M40,M50
      INTEGER W
      IF(MOD(IT,50).NE.1) GO TO 9999
      CALL NEWPG
      WRITE(W,100)
  100 FORMAT('0PRINT OF CHANGES OF PARAMETERS  :'/'0PARAMETER',9X,'  OLD
     1'    10X,'CHANGE',7X,'NEW',10X,'E.S.D.',1X,'CHANGE/E.S.D. BLOCK'/)
 9999 IT=IT+1
      RETURN
      END
      SUBROUTINE ZMENA(L,X,SX,K,IT,AT,IP)
      COMMON/UNITS/ W,M40,M50
      COMMON/REFL/  YO,YC,A,B,DY,WDYQ,NULOVA,SINTHL,RHO,WT,DER(9708),
     1              DC(9708),DYP,IQ
      COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(9708),PS(9708),
     1       AM(75000)
      EQUIVALENCE (NPSS(10),NPSM)
      LOGICAL NULOVA
      INTEGER W
      CHARACTER*8 AT,ATM
      CHARACTER*12 L,LM
      GO TO(1000,2000,6000),K
 1000 IP=0
      RXS=0.
      RXM=0.
      GO TO 9999
 2000 IP=IP+1
      CALL PIS(IT)
      DX=PS(IP)
      IF(KI(IP).EQ.0.AND.DX.EQ.0.) GO TO 5000
      SX=DER(IP)
      RX=ABS(DX/SX)
      IF(RX.LE.RXM) GO TO 4000
      RXM=RX
      ATM=AT
      LM=L
 4000 RXS=RXS+RX
      XNEW=X+TLUM*DX
      WRITE(W,100)  L,X,DX,XNEW,SX,RX,KI(IP)
      X=XNEW
      GO TO 9999
 5000 WRITE(W,101)  L,X,X
      GO TO 9999
 6000 CALL NEWPG
      RXS=RXS/FLOAT(NPSM)
      WRITE(W,102) RXS,RXM,ATM,LM
  100 FORMAT(1X,A12,1X,4F13.6,F10.2,I7)
  101 FORMAT(1X,A12,1X, F13.6,13X,F13.6)
  102 FORMAT('0AVERAGE CHANGE/E.S.D. IS',F7.2/'0MAXIMAL CHANGE/E.S.D. IS
     1',F7.2,' FOR ',A8,1X,A12)
 9999 RETURN
      END
      SUBROUTINE SMI(A,N,ISING,IBLOCK)
      DIMENSION A(1),H(200),D(200)
      COMMON/UNITS/ W,M40,M50
      COMMON/COVC/  ICOV
      INTEGER W
      M=(N*(N+1))/2
      ISING=0
      IF(N.GT.1) GO TO 100
      A(1)=1./A(1)
      RETURN
  100 IJ=0
      DO 500I=1,N
      IJ=IJ+I
      P=A(IJ)
      IF(P.LE.0.) GO TO 1010
  500 D(I)=1./SQRT(P)
      CALL ZNORM(A,D,N)
      K=N+1
 1000 K=K-1
      IF(K.LE.0) GO TO 9000
      P=A(1)
      IF(P.GT.0.)GO TO 1100
 1010 ISING=1
      RETURN
 1100 P=1./P
      II=1
      DO 3000I=2,N
      M=II
      II=II+I
      Q=A(M+1)
      HI=Q*P
      IF(I.LE.K) HI=-HI
      H(I)=HI
      M2=M+2
      DO 2000IJ=M2,II
      IJI=IJ-I
      IJM=IJ-M
 2000 A(IJI)=A(IJ)+Q*H(IJM)
 3000 CONTINUE
      M=M-1
      A(II)=P
      DO 4000I=2,N
      IJM=M+I
 4000 A(IJM)=H(I)
      GO TO 1000
 9000 CALL ZNORM(A,D,N)
      IF(ICOV.EQ.0) GO TO 9999
      IJ=0
      DO 9001I=1,N
      IJ=IJ+I
      P=A(IJ)
 9001 D(I)=1./SQRT(P)
      CALL ZNORM(A,D,N)
      IC=0
      K=0
      DO 9020J=1,N
      DO 9010I=1,J
      K=K+1
      IF(I.EQ.J) GO TO 9010
      IF(ABS(A(K)).LE..9) GO TO 9010
      IC=IC+1
      IF(MOD(IC,50).NE.1) GO TO 9005
      CALL NEWPG
      WRITE(W,101) IBLOCK
  101 FORMAT('0ELEMENTS OF COVARIATION MATRIX (BLOCK #',I1,
     1       ') LARGER THEN 0.9')
 9005 CALL COVA(I,J,A(K),IBLOCK)
 9010 CONTINUE
 9020 CONTINUE
      DO 9030I=1,N
 9030 D(I)=1./D(I)
      CALL ZNORM(A,D,N)
 9999 RETURN
      END
      SUBROUTINE ZNORM(A,D,N)
      DIMENSION A(1),D(1)
      K=0
      DO 2000J=1,N
      DJ=D(J)
      DO 1000I=1,J
      K=K+1
 1000 A(K)=A(K)*DJ*D(I)
 2000 CONTINUE
      RETURN
      END
      SUBROUTINE NASOB(A,P,X,N)
      DIMENSION A(1),P(1),X(1)
      IJP=0
      DO 2000I=1,N
      XI=0.0
      IJ=IJP
      ID=1
      DO 1000J=1,N
      IJ=IJ+ID
      XI=XI+A(IJ)*P(J)
      IF(J.GE.I) ID=J
 1000 CONTINUE
      X(I)=XI
 2000 IJP=IJP+I
      RETURN
      END
      SUBROUTINE ZMDER(DER,RM,N)
      DIMENSION DER(1),RM(1),DR(6)
      K=0
      DO 2000I=1,N
      DR(I)=0.
      DO 1000J=1,N
      K=K+1
 1000 DR(I)=DR(I)+RM(K)*DER(J)
 2000 CONTINUE
      DO 3000I=1,N
 3000 DER(I)=DR(I)
      RETURN
      END
      SUBROUTINE ZMZPET(PS,DER,RM,N)
      DIMENSION PS(1),DER(1),RM(1),P(6),DR(6)
      DO 2000J=1,N
      K=J
      P(J)=0.
      DR(J)=0.
      DO 1000I=1,N
      P(J)=P(J)+RM(K)*PS(I)
      DR(J)=DR(J)+(RM(K)*DER(I))**2
 1000 K=K+N
 2000 CONTINUE
      DO 3000I=1,N
      PS(I)=P(I)
 3000 DER(I)=SQRT(DR(I))
      RETURN
      END
      FUNCTION BESJ(X,N,D,K)
      M=IABS(N)
      Z=X*.5
      FN=1.
      IF(M.LE.1) GO TO 2000
      DO 1000I=2,M
 1000 FN=FN*FLOAT(I)
 2000 A=Z**M/FN
      Z=Z**2
      S=A
      DO 3000I=1,100
      A=A*Z/FLOAT(I*(M+I))
      IF(K.NE.1) A=-A
      S=S+A
      IF(ABS(A).LT.D) GO TO 4000
 3000 CONTINUE
 4000 IF(MOD(M,2).NE.0.AND.N.LT.0.AND.K.NE.1) S=-S
      BESJ=S
      RETURN
      END
      SUBROUTINE SMI3(A,AI,D)
      DIMENSION A(3,3),AI(3,3)
      DO 1I=1,3
      K=1+MOD(I+1,3)
      M=1+MOD(I+3,3)
      DO 1J=1,3
      L=1+MOD(J+1,3)
      N=1+MOD(J+3,3)
    1 AI(I,J)=A(L,K)*A(N,M)-A(L,M)*A(N,K)
      D=0.0
      DO 2I=1,3
    2 D=D+A(I,1)*AI(1,I)
      DO 3I=1,3
      DO 3J=1,3
    3 AI(I,J)=AI(I,J)/D
      RETURN
      END
      SUBROUTINE SETAMP
      COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF,
     1              CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ
      COMMON/CM4/   ATOM(100),AI(100),ISF(100),ITF(100),X(3,100),
     1              BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100)
     2             ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100),
     3              AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100),
     4              KFX(100),KFB(100)
      COMMON/METRICS/ G(9),GI(9),TR(9),VOL
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      CHARACTER*8 ATOM,MOLNAME
      DIMENSION A(3)
      DO 2000I=1,NMOL
      K=KMODXM(I)
      IF(K.LE.0) GO TO 2000
      DO 1000J=1,K
      A(1)=0.
      A(2)= URX(3,J,I)*VOL
      A(3)=-URX(2,J,I)*VOL
      CALL MULTM(GI,A,DURDX(1,J,I),3,3,1)
      A(2)= URY(3,J,I)*VOL
      A(3)=-URY(2,J,I)*VOL
      CALL MULTM(GI,A,DURDX(4,J,I),3,3,1)
      A(1)=-URX(3,J,I)*VOL
      A(2)=0.
      A(3)= URX(1,J,I)*VOL
      CALL MULTM(GI,A,DURDX(7,J,I),3,3,1)
      A(1)=-URY(3,J,I)*VOL
      A(3)= URY(1,J,I)*VOL
      CALL MULTM(GI,A,DURDX(10,J,I),3,3,1)
      A(1)= URX(2,J,I)*VOL
      A(2)=-URX(1,J,I)*VOL
      A(3)=0.
      CALL MULTM(GI,A,DURDX(13,J,I),3,3,1)
      A(1)= URY(2,J,I)*VOL
      A(2)=-URY(1,J,I)*VOL
 1000 CALL MULTM(GI,A,DURDX(16,J,I),3,3,1)
 2000 CONTINUE
      DO 9000I=1,NA
      IF(KMODX(I).GE.0) GO TO 9000
      N=-KMODX(I)
      K=KMODXM(N)
      IF(K.LE.0) GO TO 9000
      PHF(I)=PHFM(N)
      X1=X(1,I)-XM(1,N)
      X2=X(2,I)-XM(2,N)
      X3=X(3,I)-XM(3,N)
      DO 8000J=1,K
      DO 3000L=1,3
      UX(L,J,I)=UTX(L,J,N)
 3000 UY(L,J,I)=UTY(L,J,N)
      A(1)=-X2*URX(3,J,N)+X3*URX(2,J,N)
      A(2)=-X3*URX(1,J,N)+X1*URX(3,J,N)
      A(3)=-X1*URX(2,J,N)+X2*URX(1,J,N)
      DO 4000L=1,3
 4000 A(L)=A(L)*VOL
      CALL CULTM(GI,A,UX(1,J,I),3,3,1)
      A(1)=-X2*URY(3,J,N)+X3*URY(2,J,N)
      A(2)=-X3*URY(1,J,N)+X1*URY(3,J,N)
      A(3)=-X1*URY(2,J,N)+X2*URY(1,J,N)
      DO 5000L=1,3
 5000 A(L)=A(L)*VOL
 8000 CALL CULTM(GI,A,UY(1,J,I),3,3,1)
      A(1)=0.
      A(2)=-X3*VOL
      A(3)= X2*VOL
      CALL MULTM(GI,A,DURDR(1,I),3,3,1)
      A(1)= X3*VOL
      A(2)=0.
      A(3)=-X1*VOL
      CALL MULTM(GI,A,DURDR(4,I),3,3,1)
      A(1)=-X2*VOL
      A(2)= X1*VOL
      A(3)=0.
      CALL MULTM(GI,A,DURDR(7,I),3,3,1)
 9000 CONTINUE
      RETURN
      END
      SUBROUTINE CULTM(A,B,C,N1,N2,N3)
      DIMENSION A(1),B(1),C(1)
      DO 3000I=1,N1
      JKP=1
      IK=I
      DO 2000K=1,N3
      IJ=I
      JK=JKP
      P=C(IK)
      DO 1000J=1,N2
      P=P+A(IJ)*B(JK)
      IJ=IJ+N1
 1000 JK=JK+1
      C(IK)=P
      IK=IK+N1
 2000 JKP=JKP+N2
 3000 CONTINUE
      RETURN
      END
      SUBROUTINE MULTM(A,B,C,N1,N2,N3)
      DIMENSION A(1),B(1),C(1)
      DO 3000I=1,N1
      JKP=1
      IK=I
      DO 2000K=1,N3
      IJ=I
      JK=JKP
      P=0.
      DO 1000J=1,N2
      P=P+A(IJ)*B(JK)
      IJ=IJ+N1
 1000 JK=JK+1
      C(IK)=P
      IK=IK+N1
 2000 JKP=JKP+N2
 3000 CONTINUE
      RETURN
      END
      SUBROUTINE INDEX(J,K,L)
      K=J
      L=J
      IF(J.LE.3) RETURN
      K=J/3
      L=(J+4)/3
      RETURN
      END
      SUBROUTINE COVA(IA,IB,COV,IBLOCK)
      COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF,
     1              CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ
      COMMON/UNITS/ W,M40,M50
      COMMON/CM4/   ATOM(100),AI(100),ISF(100),ITF(100),X(3,100),
     1              BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100)
     2             ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100),
     3              AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100),
     4              KFX(100),KFB(100)
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      COMMON/LABEL/ L1,L2(3),L5,L7,L8,L9,L10
      COMMON NPS(9),NAM(9),NPSS(10),NAMM,NP,NZ,TLUM,Q,KI(6883),PS(6883),
     1       AM(75000)
      INTEGER W
      LOGICAL CSYM
      CHARACTER*1 CU,CV,CW,CP
      CHARACTER*4 CQ,S4
      CHARACTER*8 MOLNAME,ATOM,S8,ATO
      CHARACTER*12 L1,L2,L5,L7,L8,L9,L10,LI
      CHARACTER*85 IVEN
      DIMENSION IP(2)
      DATA CU,CV,CW/'U','V','W'/,S8/'        '/,S4/'    '/
      ENCODE(85,100,IVEN) IA,IB,COV
  100 FORMAT(1X,'(',I3,',',I3,')',F7.3,' COVARIATION BETWEEN',48X)
      J=38
      K=0
      DO 1000I=1,NPMK
      IF(KI(I).NE.IBLOCK) GO TO 1000
      K=K+1
      IF(K.EQ.IA) IP(1)=I
      IF(K-IB) 1000,500,1500
  500 IP(2)=I
 1000 CONTINUE
 1500 DO 3000I=1,2
      ATO=S8
      CQ=S4
      K=IP(I)
      IF(K.GT.6) GO TO 2100
      ENCODE(12,102,LI) K
      GO TO 2900
 2100 IF(K.GT.18) GO TO 2200
      K=K-6
      IF(K.NE.1.AND.K.NE.7) GO TO 2140
      IF(IEXT.NE.1) GO TO 2140
      IF(K.EQ.7) GO TO 2120
      LI=L7
      GO TO 2900
 2120 LI=L8
      GO TO 2900
 2140 IF(K.GT.6) GO TO 2150
      CALL INDEX(K,II,JJ)
      ENCODE(12,110,LI) II,JJ
      GO TO 2900
 2150 CALL INDEX(K-6,II,JJ)
      ENCODE(12,111,LI) II,JJ
      GO TO 2900
 2200 IF(K.GE.NPMP) GO TO 2700
      II=(K-19)/92+1
      ATO=ATOM(II)
      K=MOD(K-19,92)+1
      IF(K.GT.1) GO TO 2300
      LI=L1
      GO TO 2900
 2300 IF(K.GT.4) GO TO 2400
      LI=L2(K-1)
      GO TO 2900
 2400 IF(K.GT.10) GO TO 2500
      CALL INDEX(K-4,II,JJ)
      IF(LITE.EQ.1) GO TO 2450
      ENCODE(12,112,LI) II,JJ
      GO TO 2900
 2450 ENCODE(12,117,LI) II,JJ
      GO TO 2900
 2500 IF(K.GT.19) GO TO 2530
      IF(K.GT.11) GO TO 2510
      LI=L10
      GO TO 2900
 2510 II=(K-12)/2+1
      JJ=MOD(K-12,2)+1
      IF(JJ.GT.1) GO TO 2520
      ENCODE(12,115,LI) II
      GO TO 2900
 2520 ENCODE(12,116,LI) II
      GO TO 2900
 2530 IF(K.GT.43) GO TO 2590
      II=(K-20)/6+1
      JJ=MOD(K-20,6)+1
      IF(JJ.GT.3) GO TO 2580
      ENCODE(12,113,LI) CU,JJ,II
      GO TO 2900
 2580 ENCODE(12,114,LI) CU,JJ-3,II
      GO TO 2900
 2590 IF(K.GT.91) GO TO 2600
      II=(K-44)/12+1
      JJ=MOD(K-44,12)+1
      IF(JJ.LE.6) THEN
      ENCODE(12,113,LI) 'T',JJ,II
      ELSE
      ENCODE(12,113,LI) 'T',JJ-6,II
      ENDIF
      GO TO 2900
 2600 LI=L9
      GO TO 2900
 2700 II=(K-NPMP)/49+1
      ATO=MOLNAME(II)
      K=MOD(K-NPMP,49)+1
      IF(K.GT.48) GO TO 2800
      CP=CV
      IF(K.LE.24) GO TO 2720
      K=K-24
      CP=CW
 2720 II=(K-1)/6+1
      JJ=MOD(K-1,6)+1
      IF(JJ.GT.3) GO TO 2730
      ENCODE(12,113,LI) CP,JJ,II
      GO TO 2900
 2730 ENCODE(12,114,LI) CP,JJ-3,II
      GO TO 2900
 2800 LI=L9
 2900 IF(I.EQ.1) GO TO 2910
      IVEN=IVEN(1:J)//' AND '
      J=J+5
 2910 K=IDEL(LI)
      IF(K.NE.0) IVEN=IVEN(1:J)//LI(1:K)
      J=J+K+1
      K=IDEL(ATO)
      IF(K.NE.0) IVEN=IVEN(1:J)//ATO(1:K)
      J=J+K+1
      J=IDEL(IVEN)
 3000 CONTINUE
      WRITE(W,101) IVEN
      RETURN
  101 FORMAT(A85)
  102 FORMAT('SCALE-',I1)
  110 FORMAT('RHO(',I1,',',I1,')    ')
  111 FORMAT('G(',I1,',',I1,')      ')
  112 FORMAT('BETA(',I1,',',I1,')   ')
  113 FORMAT(A1,I1,'*SIN(',I1,'*X4)')
  114 FORMAT(A1,I1,'*COS(',I1,'*X4)')
  115 FORMAT('A*SIN(',I1,'*X4) ')
  116 FORMAT('A*COS(',I1,'*X4) ')
  117 FORMAT('U(',I1,',',I1,')      ')
      END
      FUNCTION IDEL(L)
      CHARACTER*(*) L
      IL=LEN(L)
  500 IF(L(IL:IL).NE.' ') GO TO 1000
      IL=IL-1
      IF(IL.GT.0) GO TO 500
 1000 IDEL=IL
      RETURN
      END
      SUBROUTINE CENTER(I)
      COMMON/KLICE/ NA,NS,NF,CSYM,VYH,NC,NOWR,IWQ,YOMIN,YOMAX,BLKOEF,
     1              CENTR,NCYKL,SNLMN,SNLMX,KIM,DIF,KIW,LITE,IFSQ
      COMMON/ATOMY/ FF(32,15),FX(15),FY(15),AW(15)
      COMMON/CM4/   ATOM(100),AI(100),ISF(100),ITF(100),X(3,100),
     1              BETA(6,100),SC(6),KMODX(100),UX(3,4,100),UY(3,4,100)
     2             ,PHF(100),QCNTM(100),KMODS(100),A0(100),AX(4,100),
     3              AY(4,100),KMODB(100),BX(6,4,100),BY(6,4,100),
     4              KFX(100),KFB(100)
      COMMON/SYMM/  IS(2,4,24),TS(3,24),TAU(24),FSM(24)
      COMMON/METRICS/ G(9),GI(9),TR(9),VOL
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON/UNITS/ W,M40,M50
      CHARACTER*8 ATOM,MOLNAME
      INTEGER W
      LOGICAL CSYM
      DIMENSION XP(3),A(3),TI(6),CP(3),TRM(9)
      DO 910J=1,6
  910 TI(J)=0.
      IF(KIW.LE.0) GO TO 3500
      FN=0.
      DO 900J=1,3
  900 CP(J)=0.
      DO 3000J=1,NPG(I)
      DO 2000K=1,NA
      IF(KMODX(K).NE.-I) GO TO 2000
      AIW=AI(K)*AW(ISF(K))
      IF(J.GT.1) GO TO 1000
      DO 950L=1,3
  950 XP(L)=X(L,K)
      GO TO 1100
 1000 M=KPG(J,I)
      MM=IABS(M)
      DO 1010L=1,3
      XP(L)=0.
 1010 A(L)=X(L,K)-XM(L,I)
      DO 1013L=1,3
      DO 1012LL=1,2
      IND=IS(LL,L,MM)
      IF(IND.EQ.0) GO TO 1012
      Z=SIGN(1.,FLOAT(IND))
      IND=IABS(IND)
      XP(L)=XP(L)+Z*A(IND)
 1012 CONTINUE
 1013 CONTINUE
      IF(M.GE.0) GO TO 1015
      DO 1014L=1,3
 1014 XP(L)=-XP(L)
 1015 DO 1020L=1,3
 1020 XP(L)=XP(L)+XM(L,I)
 1100 FN=FN+AIW
      DO 1500L=1,3
 1500 CP(L)=CP(L)+AIW*XP(L)
 2000 CONTINUE
 3000 CONTINUE
      DO 3100J=1,3
 3100 XM(J,I)=CP(J)/FN
 3500 DO 6000J=1,NPG(I)
      DO 5000K=1,NA
      IF(KMODX(K).NE.-I) GO TO 5000
      AIW=AI(K)*AW(ISF(K))
      IF(J.GT.1) GO TO 4000
      DO 3550L=1,3
 3550 XP(L)=X(L,K)-XM(L,I)
      GO TO 4100
 4000 M=KPG(J,I)
      MM=IABS(M)
      DO 4010L=1,3
      XP(L)=0.
 4010 A(L)=X(L,K)-XM(L,I)
      DO 4013L=1,3
      DO 4012LL=1,2
      IND=IS(LL,L,MM)
      IF(IND.EQ.0) GO TO 4012
      Z=SIGN(1.,FLOAT(IND))
      IND=IABS(IND)
      XP(L)=XP(L)+Z*A(IND)
 4012 CONTINUE
 4013 CONTINUE
      IF(M.GE.0) GO TO 4100
      DO 4014L=1,3
 4014 XP(L)=-XP(L)
 4100 CALL MULTM(TR,XP,A,3,3,1)
      TI(1)=TI(1)+AIW*(A(2)**2+A(3)**2)
      TI(2)=TI(2)-AIW*A(1)*A(2)
      TI(3)=TI(3)+AIW*(A(1)**2+A(3)**2)
      TI(4)=TI(4)-AIW*A(1)*A(3)
      TI(5)=TI(5)-AIW*A(2)*A(3)
      TI(6)=TI(6)+AIW*(A(1)**2+A(2)**2)
 5000 CONTINUE
 6000 CONTINUE
      CALL QL(TI,TRM,A,3)
      CALL TRMAT(TRM,TRMA(1,I),3)
      RETURN
      END
      SUBROUTINE INTER
      COMMON/UNITS/ W,M40,M50
      COMMON/EXTIN/ IEXT,ITYPEX,IDISTR,IMOSAIC,IAB,FLAME,VOLUME,AA(10),
     1              BA(10),DADMI(19),TORAD,RR,EFPIP(7),EXTKOR,C2THMQ,
     2              EC(12)
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      COMMON/METRICS/ G(9),GI(9),TR(9),VOL
      DIMENSION X(3),Y(3),Z(3),SV(3),SM(9),S1(3),S2(3),TRMAP(9),SX(3),
     1          SY(3),SZ(3)
      INTEGER W
      CHARACTER*8 MOLNAME
      IF(NMOL.LE.0) GO TO 9999
      CALL NEWPG
      WRITE(W,100)
  100 FORMAT('0INTERPRETATION OF DISPLACEMENT VECTORS'/39('=')//)
      WRITE(W,106)
  106 FORMAT('0COMPONENTS OF THE AMPLITUDE OF THE DISPLACEMENT VECTORS A
     1RE EXPRESSED AS :'/'0X(I) .... I=1,2,3   THE FRACTIONAL COORDINATE
     2 SYSTEM (NORMALISED TO BE THE MAXIMUM COORDINATE EQUAL 1) '/'0Y(I)
     3 .... I=1,2,3   THE CARTESIAN COORDINATE CRYSTAL SYSTEM (NORMALISE
     4 TO BE LENGTH EQUAL 1)'/'0Z(I) .... I=1,2,3   THE CARTESIAN COORDI
     5NATE TENSOR INERTIA SYSTEM (NORMALISE TO BE LENGTH EQUAL 1)')
      WRITE(W,107)
  107 FORMAT(/' MAGNITUDE OF DISPLACEMENT VECTORS IS DONE IN ANGSTROMS F
     1OR TRANSLATION AND IN DEGREES FOR ROTATIONAL DISPLACEMENT'/)
      DO 9000I=1,NMOL
      WRITE(W,101) MOLNAME(I)
  101 FORMAT(//' MOLECULE : ',A8/20('=')/)
      WRITE(W,102)
  102 FORMAT(/' TRANSLATION DISPLACEMENT VECTORS'/)
      WRITE(W,108)
  108 FORMAT('  WAVE    MAGNITUDE',4X,'X(1)',4X,'X(2)',4X,'X(3)',
     14X,'Y(1)',4X,'Y(2)',4X,'Y(3)',4X,'Z(1)',4X,'Z(2)',4X,'Z(3)'/
     2  1X,91('-'))
      CALL MULTM(TRMA(1,I),TR,TRMAP,3,3,3)
      DO 1300K=1,KMODXM(I)
      CALL NORM(UTX(1,K,I),X,DX)
      IF(DX.EQ.0.) GO TO 1100
      DO 500L=1,3
  500 SX(L)=SUTX(L,K,I)/DX
      CALL MULTM(TR,UTX(1,K,I),Y,3,3,1)
      D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2)
      DO 1000L=1,3
 1000 Y(L)=Y(L)/D
      CALL DODEJ(TR,D,Y,SV,SM)
      IP=1
      DO 1001L=1,3
      CALL NASOB(UTXM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SY(L),1,3,1)
 1001 IP=IP+3
      CALL NASOB(UTXM(1,K,I),SV,S1,3)
      CALL MULTM(SV,S1,SD,1,3,1)
      CALL MULTM(TRMAP,UTX(1,K,I),Z,3,3,1)
      DO 1010L=1,3
 1010 Z(L)=Z(L)/D
      CALL DODEJ(TRMAP,D,Z,SV,SM)
      IP=1
      DO 1011L=1,3
      CALL NASOB(UTXM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SZ(L),1,3,1)
 1011 IP=IP+3
      DO 1020L=1,3
      SY(L)=SQRT(SY(L))
 1020 SZ(L)=SQRT(SZ(L))
      SD=SQRT(SD)
      WRITE(W,103) K,D,X,Y,Z,SD,SX,SY,SZ
  103 FORMAT(' SIN(',I1,'X4)',F10.5,9F8.3/9X,F10.5,9F8.3)
 1100 CALL NORM(UTY(1,K,I),X,DX)
      IF(DX.EQ.0.) GO TO 1300
      DO 1110L=1,3
 1110 SX(L)=SUTY(L,K,I)/DX
      CALL MULTM(TR,UTY(1,K,I),Y,3,3,1)
      D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2)
      DO 1120L=1,3
 1120 Y(L)=Y(L)/D
      CALL DODEJ(TR,D,Y,SV,SM)
      IP=1
      DO 1021L=1,3
      CALL NASOB(UTYM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SY(L),1,3,1)
 1021 IP=IP+3
      CALL NASOB(UTYM(1,K,I),SV,S1,3)
      CALL MULTM(SV,S1,SD,1,3,1)
      CALL MULTM(TRMAP,UTY(1,K,I),Z,3,3,1)
      DO 1130L=1,3
 1130 Z(L)=Z(L)/D
      CALL DODEJ(TRMAP,D,Z,SV,SM)
      IP=1
      DO 1131L=1,3
      CALL NASOB(UTYM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SZ(L),1,3,1)
 1131 IP=IP+3
      DO 1140L=1,3
      SY(L)=SQRT(SY(L))
 1140 SZ(L)=SQRT(SZ(L))
      SD=SQRT(SD)
      WRITE(W,110) K,D,X,Y,Z,SD,SX,SY,SZ
  110 FORMAT(' COS(',I1,'X4)',F10.5,9F8.3/9X,F10.5,9F8.3)
 1300 CONTINUE
      WRITE(W,104)
  104 FORMAT(/' ROTATIONAL DISPLACEMENT VECTORS'/)
      WRITE(W,108)
      DO 2300K=1,KMODXM(I)
      CALL NORM(URX(1,K,I),X,DX)
      IF(DX.EQ.0.) GO TO 2100
      DO 1500L=1,3
 1500 SX(L)=SURX(L,K,I)/DX
      CALL MULTM(TR,URX(1,K,I),Y,3,3,1)
      D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2)
      DO 2000L=1,3
 2000 Y(L)=Y(L)/D
      CALL DODEJ(TR,D,Y,SV,SM)
      IP=1
      DO 2001L=1,3
      CALL NASOB(URXM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SY(L),1,3,1)
 2001 IP=IP+3
      CALL NASOB(URXM(1,K,I),SV,S1,3)
      CALL MULTM(SV,S1,SD,1,3,1)
      CALL MULTM(TRMAP,URX(1,K,I),Z,3,3,1)
      DO 2010L=1,3
 2010 Z(L)=Z(L)/D
      CALL DODEJ(TRMAP,D,Z,SV,SM)
      IP=1
      DO 2011L=1,3
      CALL NASOB(URXM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SZ(L),1,3,1)
 2011 IP=IP+3
      DO 2020L=1,3
      SY(L)=SQRT(SY(L))
 2020 SZ(L)=SQRT(SZ(L))
      SD=SQRT(SD)/TORAD
      D=D/TORAD
      WRITE(W,109) K,D,X,Y,Z,SD,SX,SY,SZ
  109 FORMAT(' SIN(',I1,'X4)',F10.2,9F8.3/9X,F10.2,9F8.3)
 2100 CALL NORM(URY(1,K,I),X,DX)
      IF(DX.EQ.0.) GO TO 2300
      DO 2110L=1,3
 2110 SX(L)=SURY(L,K,I)/DX
      CALL MULTM(TR,URY(1,K,I),Y,3,3,1)
      D=SQRT(Y(1)**2+Y(2)**2+Y(3)**2)
      DO 2120L=1,3
 2120 Y(L)=Y(L)/D
      CALL DODEJ(TR,D,Y,SV,SM)
      IP=1
      DO 2121L=1,3
      CALL NASOB(URYM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SY(L),1,3,1)
 2121 IP=IP+3
      CALL NASOB(URYM(1,K,I),SV,S1,3)
      CALL MULTM(SV,S1,SD,1,3,1)
      CALL MULTM(TRMAP,URY(1,K,I),Z,3,3,1)
      DO 2130L=1,3
 2130 Z(L)=Z(L)/D
      CALL DODEJ(TRMAP,D,Z,SV,SM)
      IP=1
      DO 2131L=1,3
      CALL NASOB(URYM(1,K,I),SM(IP),S1,3)
      CALL MULTM(SM(IP),S1,SZ(L),1,3,1)
 2131 IP=IP+3
      DO 2140L=1,3
      SY(L)=SQRT(SY(L))
 2140 SZ(L)=SQRT(SZ(L))
      SD=SQRT(SD)/TORAD
      D=D/TORAD
      WRITE(W,105) K,D,X,Y,Z,SD,SX,SY,SZ
  105 FORMAT(' COS(',I1,'X4)',F10.2,9F8.3/9X,F10.2,9F8.3)
 2300 CONTINUE
 9000 CONTINUE
 9999 RETURN
      END
      SUBROUTINE NORM(X,Y,D)
      DIMENSION X(3),Y(3)
      D=ABS(X(1))
      DO 1000 I=2,3
      IF(D.LT.ABS(X(I))) D=ABS(X(I))
 1000 CONTINUE
      IF(D.EQ.0.) RETURN
      DO 2000 I=1,3
 2000 Y(I)=X(I)/D
      RETURN
      END
      SUBROUTINE DODEJ(T,D,Z,SV,SM)
      DIMENSION T(9),Z(3),SV(3),SM(9),SMP(9)
      DO 1000I=1,9
 1000 SMP(I)=T(I)
      CALL MULTM(Z,T,SV,1,3,3)
      DO 1500I=1,3
 1500 SV(I)=-SV(I)
      CALL CULTM(SV,Z,SMP,3,1,3)
      DO 2000I=1,9
 2000 SMP(I)=SMP(I)/D
      CALL TRMAT(SMP,SM,3)
      RETURN
      END
      SUBROUTINE SKOVEJ(AM,Q,IBLOCK,KI)
      COMMON/MOLEC/ NMOL,NPMP,NPMK,MOLNAME(10),XM(3,10),KMODXM(10),
     1              UTX(3,4,10),UTY(3,4,10),URX(3,4,10),URY(3,4,10),
     2              PHFM(10),SUTX(3,4,10),SUTY(3,4,10),SURX(3,4,10),
     3              SURY(3,4,10),SPHFM(10),DURDX(18,4,10),DURDR(9,100),
     4              NPG(10),KPG(24,10),TRMA(9,10),UTXM(6,4,10),
     5              UTYM(6,4,10),URXM(6,4,10),URYM(6,4,10)
      DIMENSION AM(1),KI(1)
      CHARACTER*8 MOLNAME
      NK=0
      DO 500I=1,NPMP-1
      IF(KI(I).EQ.IBLOCK) NK=NK+1
  500 CONTINUE
      IP=NPMP
      DO 5000I=1,NMOL
      DO 3000K=1,4
      IF(KI(IP  ).NE.IBLOCK.AND.KI(IP+1).NE.IBLOCK
     1                     .AND.KI(IP+2).NE.IBLOCK) GO TO 1000
      CALL KTERE(KI,IBLOCK,NK,IP,I1,I2,I3)
      CALL ZAPLN(AM,I1,I2,I3,Q,UTXM(1,K,I))
 1000 IF(KI(IP+3).NE.IBLOCK.AND.KI(IP+4).NE.IBLOCK
     1                     .AND.KI(IP+5).NE.IBLOCK) GO TO 2000
      CALL KTERE(KI,IBLOCK,NK,IP+3,I1,I2,I3)
      CALL ZAPLN(AM,I1,I2,I3,Q,UTYM(1,K,I))
 2000 IF(KI(IP+24).NE.IBLOCK.AND.KI(IP+25).NE.IBLOCK
     1                      .AND.KI(IP+26).NE.IBLOCK) GO TO 2500
      CALL KTERE(KI,IBLOCK,NK,IP+24,I1,I2,I3)
      CALL ZAPLN(AM,I1,I2,I3,Q,URXM(1,K,I))
 2500 IF(KI(IP+27).NE.IBLOCK.AND.KI(IP+28).NE.IBLOCK
     1                      .AND.KI(IP+29).NE.IBLOCK) GO TO 3000
      CALL KTERE(KI,IBLOCK,NK,IP+27,I1,I2,I3)
      CALL ZAPLN(AM,I1,I2,I3,Q,URYM(1,K,I))
 3000 IP=IP+6
 5000 IP=IP+25
      RETURN
      END
      SUBROUTINE KTERE(KI,IBLOCK,NK,IP,I1,I2,I3)
      DIMENSION KI(1)
      I1=0
      IF(KI(IP).NE.IBLOCK) GO TO 1000
      NK=NK+1
      I1=NK
 1000 I2=0
      IF(KI(IP+1).NE.IBLOCK) GO TO 2000
      NK=NK+1
      I2=NK
 2000 I3=0
      IF(KI(IP+2).NE.IBLOCK) GO TO 3000
      NK=NK+1
      I3=NK
 3000 RETURN
      END
      SUBROUTINE ZAPLN(A,I1,I2,I3,Q,U)
      DIMENSION A(1),U(1)
      IMAT(II,JJ)=II+(JJ*(JJ-1))/2
      IF(I1.EQ.0) GO TO 2000
      U(1)=A(IMAT(I1,I1))*Q
      IF(I2.EQ.0) GO TO 1000
      U(2)=A(IMAT(I1,I2))*Q
 1000 IF(I3.EQ.0) GO TO 2000
      U(4)=A(IMAT(I1,I3))*Q
 2000 IF(I2.EQ.0) GO TO 3000
      U(3)=A(IMAT(I2,I2))*Q
      IF(I3.EQ.0) GO TO 9999
      U(5)=A(IMAT(I2,I3))*Q
 3000 IF(I3.EQ.0) GO TO 9999
      U(6)=A(IMAT(I3,I3))*Q
 9999 RETURN
      END
      SUBROUTINE RECIP(CP,RCP)
      DIMENSION CP(6),RCP(6)
      AR=CP(1)
      BR=CP(2)
      CR=CP(3)
      COSAR=CP(4)
      COSBR=CP(5)
      COSGR=CP(6)
      VR=AR*BR*CR*SQRT(1.+2.*COSAR*COSBR*COSGR-
     1                       COSAR**2-COSBR**2-COSGR**2)
      SINAR=SQRT(1.-COSAR**2)
      SINBR=SQRT(1.-COSBR**2)
      SINGR=SQRT(1.-COSGR**2)
      A=BR*CR*SINAR/VR
      B=CR*AR*SINBR/VR
      C=AR*BR*SINGR/VR
      COSA=(COSBR*COSGR-COSAR)/(SINBR*SINGR)
      COSB=(COSGR*COSAR-COSBR)/(SINGR*SINAR)
      COSC=(COSAR*COSBR-COSGR)/(SINAR*SINBR)
      RCP(1)=A
      RCP(2)=B
      RCP(3)=C
      RCP(4)=COSA
      RCP(5)=COSB
      RCP(6)=COSC
      RETURN
      END
      SUBROUTINE RESETX(X)
      COMMON/UNITS/ W,M40,M50
      INTEGER W
      DIMENSION X(3,1)
      DATA IFIRST/0/
      IF(IFIRST.EQ.0) WRITE(W,100)
  100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETX')
      IFIRST=1
      RETURN
      END
      SUBROUTINE RESETA(A0,AX,AY)
      COMMON/UNITS/ W,M40,M50
      INTEGER W
      DIMENSION A0(1),AX(4,1),AY(4,1)
      DATA IFIRST/0/
      IF(IFIRST.EQ.0) WRITE(W,100)
  100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETA')
      IFIRST=1
      RETURN
      SUBROUTINE RESETC(BX,BY)
      COMMON/UNITS/ W,M40,M50
      INTEGER W
      DIMENSION BX(6,6,1),BY(6,4,1)
      DATA IFIRST/0/
      IF(IFIRST.EQ.0) WRITE(W,100)
  100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETC')
      IFIRST=1
      RETURN
      END
      SUBROUTINE RESETB(BETA)
      COMMON/UNITS/ W,M40,M50
      INTEGER W
      DIMENSION BETA(6,1)
      DATA IFIRST/0/
      IF(IFIRST.EQ.0) WRITE(W,100)
  100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETB')
      IFIRST=1
      RETURN
      END
      SUBROUTINE RESETD(UX,UY)
      COMMON/UNITS/ W,M40,M50
      INTEGER W
      DIMENSION UX(3,4,1),UY(3,4,1)
      DATA IFIRST/0/
      IF(IFIRST.EQ.0) WRITE(W,100)
  100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETD')
      IFIRST=1
      RETURN
      END
      SUBROUTINE RESETM(UTX,UTY,URX,URY)
      COMMON/UNITS/ W,M40,M50
      INTEGER W
      DIMENSION UTX(3,4,1),UTY(3,4,1),URX(3,4,1),URY(3,4,1)
      DATA IFIRST/0/
      IF(IFIRST.EQ.0) WRITE(W,100)
  100 FORMAT('0YOU USED A DUMMY SUBROUTINE RESETM')
      IFIRST=1
      RETURN
      END
      SUBROUTINE DSETA(DER,NA)
      DIMENSION DER(92,NA)
C
C     SUBROUTINE FOR SETTING OF DERIVATIONS OF ATOMIC PARAMETERS
C
C     FOR EACH ATOM WE HAVE 92 PARAMETERS :
C
C     #1       AI              - OCCUPATIONAL PARAMETER
C     #2-#4    X(I),I=1,3      - POSITIONAL PARAMETERS
C     #5-#10   BETA(I),I=1,6   - ANISOTROPIC TEMPERATURE PARAMETERS
C                                IF ISOTROPIC ATOM ONLY #5 IS USED
C                                AND #6-#10 ARE VOID
C           *********************************** 
C           ****   ATOMIC PARAMETERS FOR   ****
C           ****  OCCUPATIONAL MODULATION  ****
C           ***********************************
C     #11      A0              - ABSOLUTE TERM
C     #12      AX(1)           - SIN WAVE OF 1ST HARMONIC
C     #13      AY(1)           - COS WAVE OF 1ST HARMINIC 
C     #14      AX(2)           - SIN WAVE OF 2ND HARMONIC 
C     #15      AY(2)           - COS WAVE OF 2ND HARMINIC 
C     #16      AX(3)           - SIN WAVE OF 3RD HARMONIC 
C     #17      AY(3)           - COS WAVE OF 3RD HARMINIC 
C     #18      AX(4)           - SIN WAVE OF 4TH HARMONIC 
C     #19      AY(4)           - COS WAVE OF 4TH HARMINIC
C           *********************************** 
C           ****   ATOMIC PARAMETERS FOR   ****
C           ****  POSITIONAL  MODULATION   ****
C           ***********************************
C     #20-#22  UX(I,1,1),I=1,3 - SIN WAVE OF 1ST HARMONIC 
C     #23-#25  UY(I,1,1),I=1,3 - COS WAVE OF 1ST HARMINIC 
C     #26-#28  UX(I,2,1),I=1,3 - SIN WAVE OF 2ND HARMONIC 
C     #29-#31  UY(I,2,1),I=1,3 - COS WAVE OF 2ND HARMINIC 
C     #32-#34  UX(I,3,1),I=1,3 - SIN WAVE OF 3RD HARMONIC 
C     #35-#37  UY(I,3,1),I=1,3 - COS WAVE OF 3RD HARMINIC 
C     #38-#40  UX(I,4,1),I=1,3 - SIN WAVE OF 4TH HARMONIC 
C     #41-#43  UY(I,4,1),I=1,3 - COS WAVE OF 4TH HARMINIC 
C           *********************************** 
C           ****   ATOMIC PARAMETERS FOR   ****
C           ****  TEMP.FACTOR MODULATION   ****
C           ***********************************
C     #44-#49  BX(I,4,1),I=1,6 - SIN WAVE OF 1ST HARMONIC
C     #50-#55  BY(I,4,1),I=1,6 - COS WAVE OF 1ST HARMONIC
C     #56-#61  BX(I,4,2),I=1,6 - SIN WAVE OF 2ND HARMONIC
C     #62-#67  BY(I,4,2),I=1,6 - COS WAVE OF 2ND HARMONIC
C     #68-#73  BX(I,4,3),I=1,6 - SIN WAVE OF 3RD HARMONIC
C     #74-#79  BY(I,4,3),I=1,6 - COS WAVE OF 3RD HARMONIC
C     #80-#85  BX(I,4,4),I=1,6 - SIN WAVE OF 4TH HARMONIC
C     #86-#91  BY(I,4,4),I=1,6 - COS WAVE OF 4TH HARMONIC
C     #92      PH              - PHASON FACTOR
      RETURN
      END
      SUBROUTINE DSETM(DER,NMOL)
      DIMENSION DER(49,NMOL)
C
C     SUBROUTINE FOR SETTING OF DERIVATIONS OF MOLECULAR PARAMETERS
C
C     FOR EACH MOLECULE WE HAVE 49 PARAMETERS :
C
C     #1-#3    UTX(I,1,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 1ST HARMONIC
C     #4-#6    UTY(I,1,1),I=1,3 - TRANSLATIONAL COS WAVE OF 1ST HARMONIC
C     #7-#10   UTX(I,2,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 2ND HARMONIC
C     #11-#12  UTY(I,2,1),I=1,3 - TRANSLATIONAL COS WAVE OF 2ND HARMONIC
C     #13-#15  UTX(I,3,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 3RD HARMONIC
C     #16-#18  UTY(I,3,1),I=1,3 - TRANSLATIONAL COS WAVE OF 3RD HARMONIC
C     #19-#21  UTX(I,4,1),I=1,3 - TRANSLATIONAL SIN WAVE OF 4TH HARMONIC
C     #25-#27  URX(I,1,1),I=1,3 - ROTATIONAL SIN WAVE OF 1ST HARMONIC
C     #28-#30  URY(I,1,1),I=1,3 - ROTATIONAL COS WAVE OF 1ST HARMONIC
C     #31-#33  URX(I,2,1),I=1,3 - ROTATIONAL SIN WAVE OF 2ND HARMONIC
C     #34-#36  URY(I,2,1),I=1,3 - ROTATIONAL COS WAVE OF 2ND HARMONIC
C     #37-#39  URX(I,3,1),I=1,3 - ROTATIONAL SIN WAVE OF 3RD HARMONIC
C     #40-#42  URY(I,3,1),I=1,3 - ROTATIONAL COS WAVE OF 3RD HARMONIC
C     #43-#45  URX(I,4,1),I=1,3 - ROTATIONAL SIN WAVE OF 4TH HARMONIC
C     #46-#48  URY(I,4,1),I=1,3 - ROTATIONAL COS WAVE OF 4TH HARMONIC
C     #49      PH               - MOLECULAR PHASON FACTOR
      RETURN
      END
      SUBROUTINE TRMAT(A,B,N)
C
C  .... SUBROUTINE FOR CALCULATION OF TRANSPOSE MATRIX
C
      DIMENSION A(N,N),B(N,N)
      DO 2000J=1,N
      DO 1000I=1,N
 1000 B(J,I)=A(I,J)
 2000 CONTINUE
      RETURN
      END
      SUBROUTINE QL(A,Z,D,N)
      DIMENSION E(20),D(1),A(1),Z(1)
      DATA TOL,EPS/.0000001,.0000001/
      IMD=1-2*MOD(N,2)
      NB=N*(N-1)+1
      K=0
      DO 20I=1,N
      DO 10J=1,N
      K=K+1
      IF(I-J) 8,9,8
    9 Z(K)=1.
      GO TO 10
    8 Z(K)=0.
   10 CONTINUE
   20 CONTINUE
      IZ=N*(N+1)/2
      L=N
      DO  900IP=3,N
      I=L
      L=L-1
      H=0.
      IZ=IZ-I
      DO  100K=1,L
      IND=IZ+K
      F=A(IND)
      D(K)=F
  100 H=H+F**2
      IF(H.GT.TOL) GO TO 150
      E(I)=0.
      H=0.
      GO TO 890
  150 G=-SQRT(H)
      IF(F.LT.0.) G=-G
      E(I)=G
      H=H-F*G
      RH=1./H
      DL=F-G
      D(L)=DL
      IND=IZ+L
      A(IND)=DL
      F=0.
      IND=1
      DO  300J=1,L
      G=0.
      JK=IND
      DO  200K=1,L
      G=G+A(JK )*D(K)
      ID=1
      IF(J.LE.K) ID=K
      JK=JK +ID
  200 CONTINUE
      G=G*RH
      E(J)=G
      F=F+G*D(J)
      IND=IND+J
  300 CONTINUE
      HH=F*.5*RH
      JK=0
      DO  500J=1,L
      F=D(J)
      G=E(J)-HH*F
      E(J)=G
      DO  400K=1,J
      JK=JK+1
  400 A(JK)=A(JK)-F*E(K)-G*D(K)
  500 CONTINUE
  890 JK=IZ+I
      D(I)=A(JK)
      A(JK)=H
  900 CONTINUE
      D(1)=A(1)
      E(2)=A(2)
      D(2)=A(3)
      B=0.
      F=0.
      DO 2000L=1,N
      LN=L*N-N
      IF(L.EQ.N) GO TO 1800
      L1=L+1
      H=EPS*(ABS(D(L))+ABS(E(L1)))
      IF(B.LT.H) B=H
      DO 1100M=L1,N
      IF(ABS(E(M)).LE.B) GO TO 1200
 1100 CONTINUE
      M=N+1
 1200 M=M-1
      IF(M.EQ.L) GO TO 1800
      ML=M-L
      MN=M*N
 1250 G=D(L)
      EL=E(L1)
      DL=(D(L1)-G)/(2.*EL)
      H=ABS(DL)+SQRT(DL**2+1.)
      IF(DL.LT.0.) H=-H
      DL=EL/H
      D(L)=DL
      H=G-DL
      IF(L1.GT.N) GO TO 1350
      DO 1300I=L1,N
 1300 D(I)=D(I)-H
 1350 F=F+H
      DL=D(M)
      C=1.
      S=0.
      IND=MN
      DO 1600JK=1,ML
      IND=IND-N
      I=M-JK
      I1=I+1
      EL=E(I1)
      G=C*EL
      H=C*DL
      IF(ABS(DL).LT.ABS(EL)) GO TO 1400
      C=EL/DL
      R=SQRT(C**2+1.)
      E(I+2)=S*DL*R
      S=C/R
      C=1./R
      GO TO 1500
 1400 C=DL/EL
      R=SQRT(C**2+1.)
      E(I+2)=S*EL*R
      S=1./R
      C=C*S
 1500 DL=C*D(I)-S*G
      D(I1 )=H+S*(C*G+S*D(I))
      DO 1550K=1,N
      IZ=K+IND
      IZN=IZ-N
      H=Z(IZN)
      G=Z(IZ)
      Z(IZN) =C*H-S*G
 1550 Z(IZ  )=S*H+C*G
 1600 CONTINUE
      EL=S*DL
      E(L1 )=EL
      D(L)=C*DL
      IF(ABS(EL).GT.B)GO TO 1250
 1800 DL=D(L)+F
      JK=L-1
      IF(JK.LT.1) GO TO 1855
      IZ=LN
      DO 1850K=1,JK
      I=L-K
      IZ=IZ-N
      IF(DL.GE.D(I)) GO TO 1860
      IMD=-IMD
      DO 1820J=1,N
      IND=IZ+J
      INZ=IND+N
      P=Z(IND)
      Z(IND)=Z(INZ)
 1820 Z(INZ)=P
 1850 D(I+1)=D(I)
 1855 I=0
 1860 D(I+1)=DL
 2000 CONTINUE
      IF(IMD) 2010,2100,2100
 2010 DO 2020I=1,N
 2020 Z(I)=-Z(I)
 2100 IZ=1
      DO 2500I=3,N
      L=I-1
      IZ=IZ+L
      IND=IZ+I
      H=A(IND)
      IF(H.EQ.0.) GO TO 2500
      DO 2400J=1,NB,N
      JK=J-1
      S=0.
      DO 2200K=1,L
      M=IZ+K
      IND=JK+K
 2200 S=S+A(M)*Z(IND)
      S=S/H
      DO 2300K=1,L
      IND=JK+K
      M=IZ+K
 2300 Z(IND)=Z(IND)-S*A(M)
 2400 CONTINUE
 2500 CONTINUE
      RETURN
      END
      SUBROUTINE BMOD(F,DF1,DF2,DF3,DB,D1,D2,M,ITYP)
      COMMON/STALY/ PI2
      GO TO (1000,2000,3000),ITYP
 1000 POM=EXP(DB)
      IF(M.EQ.0) THEN
      F=D1*POM+D2*(POM-1.)/DB-D1-D2+1.
      DF1=D1*POM+D2*(POM*DB-POM+1.)/DB**2
      DF2=POM-1.
      DF3=(POM-1.)/DB-1.
      ELSE
      ARG=PI2*FLOAT(M)*.5
      ALFA1=ARG*D1
      CS1=COS(ALFA1)
      SN1=SIN(ALFA1)
      IF(D2.NE.0.) THEN
      ALFA2=ARG*D2
      AQ=ALFA2**2
      CS12=COS(ALFA1+ALFA2)
      SN12=SIN(ALFA1+ALFA2)
      BQ=DB**2
      AB=ALFA2*DB
      FCIT=AB*POM*CS1-AB*CS12+BQ*POM*SN1-BQ*SN12
      FJMEN=ALFA2*(AQ+BQ)
      FJQ=FJMEN**2
      Z=FCIT/FJMEN
      F=D2*Z
      DF1=D2*((POM*CS1*(AB+ALFA2)-ALFA2*CS12+2.*DB*(POM*SN1-SN12))*FJMEN
     1        -2.*AB*FCIT)/FJQ
      DF2=(-AB*ALFA2*POM*SN1+AB*ALFA2*SN12+BQ*POM*ALFA2*CS1
     1     -BQ*ALFA2*CS12)/FJMEN
      DF3=Z+((AB*POM*CS1-AB*CS12+ALFA2*AB*SN12-BQ*ALFA2*CS12)*FJMEN
     1       -(3.*AQ+BQ)*ALFA2*FCIT)/FJQ
      ELSE
      FCIT=SN1/ARG
      F=(POM-1.)*FCIT
      DF1=POM*FCIT
      DF2=(POM-1.)*CS1
      DF3=0.
      ENDIF
      ENDIF
      RETURN
 2000 IF(DB.NE.0.) THEN
      D=SQRT(ABS(DB))
      POM=SIN(D)/D
      POMD=-(COS(D)*D-SIN(D))/D**3*.5
      ELSE
      POM=1.
      POMD=0.
      ENDIF
      IF(M.EQ.0) THEN
      F=D1*(POM-1.)+1.
      DF1=D1*POMD
      DF2=POM-1.
      ELSE
      ARG=PI2*FLOAT(M)*.5
      ALFA1=ARG*D1
      CS1=COS(ALFA1)
      SN1=SIN(ALFA1)
      FCIT=SN1/ARG
      F=(POM-1.)*FCIT
      DF1=POMD*FCIT
      DF2=(POM-1.)*CS1
      DF3=0.
      ENDIF
      RETURN
 3000 IF(M.EQ.0) THEN
      F=D1
      DF1=0.
      DF2=1.
      DF3=0.
      ELSE
      ARG=PI2*FLOAT(M)*.5
      ALFA1=ARG*D1
      CS1=COS(ALFA1)
      SN1=SIN(ALFA1)
      F=SN1/ARG
      DF1=0.
      DF2=CS1
      DF3=0.
      ENDIF
      RETURN
      END
      SUBROUTINE XMOD(F,DF1,DF2,DF3,U,U1,U2,M,ITYP)
      COMMON/STALY/ PI2
      POM=PI2*FLOAT(M)*.5*U1
      ARG=POM+U
      IF(ARG.NE.0.) THEN
      ARG2=ARG**2
      SN=SIN(ARG)
      CS=COS(ARG)
      F=U1*SN/ARG
      DF1=U1*(CS*ARG-SN)/ARG2
      DF2=((SN+POM*CS)*ARG-POM*SN)/ARG2
      ELSE
      F=U1
      DF1=0.
      DF2=1.
      ENDIF
      DF3=0.
      RETURN
      END
/EXEC $TSOSLNK
PROG PXXX,FILENAM=P.JANA4S.POM,IDA=N
INCLUDE *
RESOLVE ,$X.USER.FOR1
END
/COPY P.JANA4S.POM,P.JANA4S
/ER P.JANA4S.POM
/SYSFILE SYSLST=(PRIMARY)
/EXEC $X.P.SRCXREF
C.JANA4S.POM
C.JANA4S
END
/ER C.JANA4S
/CAT C.JANA4S.POM,C.JANA4S,STATE=U
/LOGOFF NOSPOOL

