Upload New File

cluster/clu.f

+      PROGRAM CLU
+CCC   111130
+CCC   SUPPORTS EXPERIMENTAL DIELECTRIC FUNCTIONS ONLY
+CCC   NSPH=6; LI=8, LE=8; NXI=200
+CCC
+CCC   FOR IDFC>=0 ONLY,
+CCC   WHEN JXI=JWTM WRITES THE JWTM-TH TRANSITION MATRIX
+CCC
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   COMMON/C1/VH(NCOU*LITPO),VJ0(NSPH*LMTPO),
+CCC  1VYHJ(NCOU*LITPOS),VYJ0(NSPH*LMTPOS),
+CCC  2RMI(LI,NSPH),REI(LI,NSPH),W(NLEMT,4),
+CCC  3AM0M(NLEMT,NLEMT),FSAS(NSPH),SAS(NSPH,2,2),VINTS(NSPH,16),
+CCC  4V3J0(NV3J),SSCS(NSPH),SEXS(NSPH),SABS(NSPH),
+CCC  5SQSCS(NSPH),SQEXS(NSPH),SQABS(NSPH),GCSV(NSPH),
+CCC  6RXX(NSPH),RYY(NSPH),RZZ(NSPH),ROS(NSPH),RC(NSPH,NSHL),
+CCC  7IND3J(LMPO,LM),IOG(NSPH),NSHL(NSPH)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+CCC   COMMON/C2/RIS(NHSPO),DLRI(NHSPO),DC0(NSHL-NTL+1),
+CCC  1VKT(NSPH),VSZ(NSPH)
+CCC   NHSPO=MAX0(NPNT,NPNTTS)*2-1
+      COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+      COMPLEX*16 RIS,DLRI,DC0,VKT
+      COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+      COMPLEX*16 TFSAS,TSAS
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+CCC   COMMON/C6/RAC3J(LMTPO)
+      COMMON/C6/RAC3J(17)
+CCC   COMMON/C9/GIS(NDI,NLEM),GLS(NDI,NLEM),SAM(NDIT,NLEMT)
+      COMMON/C9/GIS(480,80),GLS(480,80),SAM(960,160)
+      COMPLEX*16 GIS,GLS,SAM
+CCC   DIMENSION AM(MXNDM,MXNDM)
+CCC   NDIT.LE.MXNDM
+      DIMENSION AM(960,960)
+      COMPLEX*16 AM
+CCC   DIMENSION RCF(NSPH,NSHL)
+      DIMENSION RCF(6,8),CMULLR(4,4),CMUL(4,4),CEXTLR(4,4),CEXT(4,4),
+     1U(3),US(3),UN(3),UNS(3),UP(3),UPS(3),
+     2UNMP(3),UNSMP(3),UPMP(3),UPSMP(3),DUK(3),ARGI(1),ARGS(1)
+CCC   DIMENSION ZPV(LM,3,2,2),GAPS(NSPH)
+      DIMENSION ZPV(8,3,2,2),GAPS(6),GAPV(3),
+     1GAP(3,2),GAPP(3,2),GAPM(3,2),GAPPM(3,2)
+      COMPLEX*16 GAPP,GAPPM
+CCC   DIMENSION TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
+      DIMENSION TQSE(2,6),TQSPE(2,6),TQSS(2,6),TQSPS(2,6),
+     1TQEV(3),TQSV(3),TQCE(2,3),TQCPE(2,3),TQCS(2,3),TQCPS(2,3)
+      COMPLEX*16 TQSPE,TQSPS,TQCPE,TQCPS
+CCC   DIMENSION DC0M(NSPH,NSHL-NTL),XIV(NXI)
+      DIMENSION DC0M(6,4),XIV(200)
+      COMPLEX*16 DC0M
+      COMPLEX*16 ARG,S0,S0M
+ 8010 FORMAT(1H ,16I5)
+ 8015 FORMAT(1H ,
+     1'READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM')
+ 8050 FORMAT(3(1PD17.8))
+ 8060 FORMAT(1H ,'READ(IR,*)RXX(I),RYY(I),RZZ(I)')
+ 8068 FORMAT(1H ,'READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST')
+ 8069 FORMAT(1H ,'READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST')
+ 8170 FORMAT(1H ,7(1PD10.3))
+ 8171 FORMAT(1H ,'READ(IR,*)JWTM')
+ 6000 FORMAT(1H )
+ 6050 FORMAT(1H ,' SIZE=',1PD15.7,', REFRACTIVE INDEX=',2(1PD15.7))
+ 6055 FORMAT(1H ,' SIZE=',1PD15.7)
+ 6060 FORMAT(1H ,' CIRC ',I2)
+ 6061 FORMAT(1H ,' CIRC')
+ 6065 FORMAT(1H ,'  LIN ',I2)
+ 6066 FORMAT(1H ,'  LIN')
+ 6068 FORMAT(1H ,4X,'SPHERES; LMX=LI')
+ 6069 FORMAT(1H ,4X,'SPHERES; LMX=MIN0(LI,LE)')
+ 6070 FORMAT(1H ,4X,'SPHERE ',I2)
+ 6080 FORMAT(1H ,4X,'CLUSTER')
+ 6100 FORMAT(1H ,'----- SCS ----- ABS ----- EXS ----- ALBEDS --')
+ 6150 FORMAT(1H ,1PD14.7,5(1PD15.7))
+ 6165 FORMAT(1H ,' UN=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
+ 6166 FORMAT(1H ,' UNI=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
+ 6167 FORMAT(1H ,' UNS=(',1PD12.5,',',1PD12.5,',',1PD12.5,')')
+ 6200 FORMAT(1H ,'---- SCS/GS -- ABS/GS -- EXS/GS ---')
+ 6240 FORMAT(1H ,1X,'FSAS=',2(1PD15.7))
+ 6242 FORMAT(1H ,1X,'SAS(1,1)=',2(1PD15.7),', SAS(2,1)=',2(1PD15.7))
+ 6243 FORMAT(1H ,1X,'SAS(1,2)=',2(1PD15.7),', SAS(2,2)=',2(1PD15.7))
+ 6250 FORMAT(1H ,1X,'FSAT=',2(1PD15.7))
+ 6252 FORMAT(1H ,1X,'SAT(1,1)=',2(1PD15.7),', SAT(2,1)=',2(1PD15.7))
+ 6253 FORMAT(1H ,1X,'SAT(1,2)=',2(1PD15.7),', SAT(2,2)=',2(1PD15.7))
+ 6300 FORMAT(1H ,' STOP IN LUCIN: MATRIX SINGULAR')
+ 6331 FORMAT(1H ,1X,'MULS') 
+ 6332 FORMAT(1H ,1X,'MULC') 
+ 6336 FORMAT(1H ,1X,'MULSLR') 
+ 6337 FORMAT(1H ,1X,'MULCLR') 
+ 6350 FORMAT((1H ,7X,4(1PD15.7))) 
+ 6460 FORMAT(1H ,4X,'CLUSTER (ENSEMBLE AVERAGE, MODE',I2,')')
+ 6500 FORMAT(1H ,'----- SCC ----- ABC ----- EXC ----- ALBEDC --')
+ 6550 FORMAT
+     1(1H ,' AT ',I1,' LCALC=',I3,' TOO SMALL WITH ARG=',2(1PD15.7))
+ 6600 FORMAT(1H ,'--- SCC/TGS - ABC/TGS - EXC/TGS ---')
+ 6650 FORMAT(1H ,' STOP IN HJV')
+ 6700 FORMAT(1H ,'---- SCCRT --- ABCRT --- EXCRT ----')
+ 6710 FORMAT(1H ,' REFR. INDEX OF EXTERNAL MEDIUM=',1PD15.7)
+ 6714 FORMAT(1H ,' VK=',1PD15.7,', XI IS SCALE FACTOR FOR LENGTHS')
+ 6715 FORMAT(1H ,' VK=',1PD15.7,', XI=',1PD15.7)
+ 6716 FORMAT(1H ,' XI=',1PD15.7)
+ 6720 FORMAT(1H ,' TIDG=',1PD10.3,', PIDG=',1PD10.3,
+     1', TSDG=',1PD10.3,', PSDG=',1PD10.3)
+ 6721 FORMAT(1H ,' CFMP=',1PD15.7,', SFMP=',1PD15.7)
+ 6722 FORMAT(1H ,' CFSP=',1PD15.7,', SFSP=',1PD15.7)
+ 6740 FORMAT(1H ,' SCAND=',1PD10.3)
+ 6750 FORMAT(1H ,' STOP IN DME')
+ 6800 FORMAT(1H ,1X,2('FSAC(',I1,',',I1,')=',2(1PD15.7),3X))
+ 6810 FORMAT(1H ,1X,2(' SAC(',I1,',',I1,')=',2(1PD15.7),3X))
+ 6850 FORMAT(1H ,1X,'RE(FSAC(',I1,',',I1,'))/RE(TFSAS)=',1PD15.7,
+     1', IM(FSAC(',I1,',',I1,'))/IM(TFSAS)=',1PD15.7)
+ 6930 FORMAT(1H ,1X,'QSCHU=',1PD15.7,', PSCHU=',1PD15.7,
+     1', S0MAG=',1PD15.7)
+ 6940 FORMAT(1H ,1X,'(RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=',
+     11PD15.7)
+ 6945 FORMAT(1H ,1X,'(IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=',
+     11PD15.7)
+ 6970 FORMAT(1H ,' COSAV=',1PD15.7,', RAPRS=',1PD15.7)
+ 6971 FORMAT(1H ,' Fx=',1PD15.7,', Fy=',1PD15.7,', Fz=',1PD15.7)
+ 6973 FORMAT(1H ,' Fl=',1PD15.7,', Fr=',1PD15.7,', Fk=',1PD15.7)
+ 6974 FORMAT(1H ,' Fk=',1PD15.7)
+ 6977 FORMAT(1H ,' IPO=',I2,', TQEk=',1PD15.7,', TQSk=',1PD15.7)
+ 6981 FORMAT
+     1(1H ,'  TQEx=',1PD15.7,',  TQEy=',1PD15.7,',  TQEz=',1PD15.7)
+ 6982 FORMAT
+     1(1H ,'  TQSx=',1PD15.7,',  TQSy=',1PD15.7,',  TQSz=',1PD15.7)
+ 6983 FORMAT
+     1(1H ,'  TQEl=',1PD15.7,',  TQEr=',1PD15.7,',  TQEk=',1PD15.7)
+ 6984 FORMAT
+     1(1H ,'  TQSl=',1PD15.7,',  TQSr=',1PD15.7,',  TQSk=',1PD15.7)
+ 6991 FORMAT('========== JXI =',I3,' ====================')
+ 6992 FORMAT('********** JTH =',I3,', JPH =',I3,
+     1', JTHS =',I3,', JPHS =',I3,' ********************')
+      IR=5
+      IW=6
+      IT=7
+      ITIN=17
+CCC
+CCC   OTHER COMMENTS IN FILE GLOBALCOMS
+CCC
+CCC
+CCC   CONSTANTS APPEAR AS (...).X(...)D0, OR (...).X(...)D+X(X),
+CCC   OR (...).X(...)D-X(X)
+CCC
+CCC
+CCC   MXNDM IS FIRST NUMERICAL DECLARED DIMENSION OF AM;
+CCC   NDIT=2*NDI; NDI=NSPH*NLIM
+CCC
+CCC
+CCC   VK*ROS(I) IS SIZE PARAMETER in vacuo
+CCC
+CCC
+CCC   OUTPUTS (NORMALIZED SCATTERING AMPLITUDE)*(2*PIG/(VK*VK)) IN m**3;
+CCC   OUTPUTS FORCE*(c/I0) IN m**2; 
+CCC   OUTPUTS (EXT_TO_TORQUE) AND (SCA_TO_TORQUE) CONTRIBUTIONS 
+CCC   *(c/I0)*(VK*EXRI)**3; (TORQUE OUTPUTS)/(EXDC*VK*VK) ARE IN m**2
+CCC
+CCC
+CCC   ON A LINEAR BASIS,
+CCC   POLARIZATION (1), (1,ni), (ns,1) OR -1 IS PARALLE(l),
+CCC   POLARIZATION (2), (2,ni), (ns,2) OR  1 IS PERPENDICULA(r);
+CCC   ON A CIRCULAR BASIS,
+CCC   POLARIZATION (1), (1,ni), (ns,1) OR -1 IS (R)ight,
+CCC   POLARIZATION (2), (2,ni), (ns,2) OR  1 IS (L)eft;
+CCC   ns=SCATWAVE POLARIZATION, ni=INCWAVE POLARIZATION
+CCC
+CCC
+CCC   NV3J=(LM*LMPO*(LMT+7))/6;
+CCC   LM=MAX0(LI,LE)
+CCC
+CCC
+CCC   IAVM=0 --> RANDOM AVERAGES FOR AMPLITUDE AND CROSS SECTIONS;
+CCC   IAVM=1 --> SAME AS IAVM=0, PLUS RANDOM AVERAGES FOR INTENSITY
+CCC
+CCC
+CCC   ISAM<0 --> ARBITRARILY VARYING SCATTERING PLANE,
+CCC              AS IT IS UNDEFINED FOR FORWARD OR BACKWARD SCATTERING,
+CCC              MERIDIONAL PLANE FOR INCIDENCE ANGLES IS ASSUMED;
+CCC   ISAM=0 --> MERIDIONAL PLANES
+CCC              (I. E. REFERENCE PLANES ACROSS Z AXIS)
+CCC              AS DETERMINED BY INCIDENCE AND OBSERVATION ANGLES;
+CCC   ISAM=1 --> SCATTERING PLANE ACROSS THE Z AXIS; SCATTERING
+CCC              AS A FUNCTION OF INCIDENCE AND OF OBSERVATION ANGLES;
+CCC              NEEDS ONLY TH, PH, THS; THSLST GT 180 IF NECESSARY;
+CCC   ISAM>1 --> SCATTERING PLANE ACROSS THE Z AXIS; SCATTERING
+CCC              AS A FUNCTION OF INCIDENCE ANGLES
+CCC              FOR A FIXED SCATTERING ANGLE THSCA=THS-TH;
+CCC              NEEDS ONLY TH, PH, THS
+CCC
+      OPEN(IR,FILE='DCLU',STATUS='OLD')
+      READ(IR,*)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM
+      DO 101 I=1,NSPH
+  101 READ(IR,*)RXX(I),RYY(I),RZZ(I)
+      READ(IR,*)TH,THSTP,THLST,THS,THSSTP,THSLST
+      READ(IR,*)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
+      READ(IR,*)JWTM
+      CLOSE(IR)
+      OPEN(IW,FILE='OCLU',STATUS='UNKNOWN')
+      WRITE(IW,8015)
+      WRITE(IW,8010)NSPH,LI,LE,MXNDM,INPOL,NPNT,NPNTTS,IAVM,ISAM
+      WRITE(IW,8060)
+      WRITE(IW,8050)(RXX(I),RYY(I),RZZ(I),I=1,NSPH)
+      WRITE(IW,8068)
+      WRITE(IW,8170)TH,THSTP,THLST,THS,THSSTP,THSLST
+      WRITE(IW,8069)
+      WRITE(IW,8170)PH,PHSTP,PHLST,PHS,PHSSTP,PHSLST
+      WRITE(IW,8171)
+      WRITE(IW,8010)JWTM
+      OPEN(ITIN,FILE='TEDF',FORM='UNFORMATTED',STATUS='OLD')
+      READ(ITIN)NSPHT
+      IF(NSPHT.NE.NSPH)GO TO 497
+      READ(ITIN)(IOG(I),I=1,NSPH)
+      READ(ITIN)EXDC,WP,XIP,IDFC,NXI
+      READ(ITIN)(XIV(I),I=1,NXI)
+      DO 103 I=1,NSPH
+      IF(IOG(I).LT.I)GO TO 103
+      READ(ITIN)NSHL(I),ROS(I)
+      NSH=NSHL(I)
+      READ(ITIN)(RCF(I,NS),NS=1,NSH)
+  103 CONTINUE
+      WRITE(IW,*)' READ(ITIN)NSPHT'
+      WRITE(IW,*)' READ(ITIN)(IOG(I),I=1,NSPH)'
+      WRITE(IW,*)' READ(ITIN)EXDC,WP,XIP,IDFC,NXI'
+      WRITE(IW,*)' READ(ITIN)(XIV(I),I=1,NXI)'
+      WRITE(IW,*)' READ(ITIN)NSHL(I),ROS(I)'
+      WRITE(IW,*)' READ(ITIN)(RCF(I,NS),NS=1,NSH)'
+      WRITE(IW,6000)
+      C0=0.0D0
+      SML=1.0D-3
+      IF(THSTP.EQ.C0)NTH=0
+      IF(THSTP.NE.C0)NTH=IDINT((THLST-TH)/THSTP+SML)
+      NTH=NTH+1
+      IF(PHSTP.EQ.C0)NPH=0
+      IF(PHSTP.NE.C0)NPH=IDINT((PHLST-PH)/PHSTP+SML)
+      NPH=NPH+1
+      IF(ISAM.LE.1)GO TO 110 
+      NTHS=1
+      THSCA=THS-TH
+      GO TO 111
+  110 IF(THSSTP.EQ.C0)NTHS=0
+      IF(THSSTP.NE.C0)NTHS=IDINT((THSLST-THS)/THSSTP+SML)
+      NTHS=NTHS+1
+      IF(ISAM.LT.1)GO TO 112
+  111 NPHS=1
+      GO TO 113
+  112 IF(PHSSTP.EQ.C0)NPHS=0
+      IF(PHSSTP.NE.C0)NPHS=IDINT((PHSLST-PHS)/PHSSTP+SML)
+      NPHS=NPHS+1
+  113 NK=NTH*NPH
+      NKS=NTHS*NPHS
+      NKKS=NK*NKS
+      TH1=TH
+      PH1=PH
+      THS1=THS
+      PHS1=PHS
+      CALL STR(RCF)
+      LM=MAX0(LI,LE)
+      CALL THDPS(LM,ZPV)
+      EXRI=DSQRT(EXDC)
+      WRITE(IW,6710)EXRI
+      OPEN(IT,FILE='TPPOAN',FORM='UNFORMATTED',STATUS='UNKNOWN')
+      WRITE(IT)IAVM,ISAM,INPOL
+      WRITE(IT)NXI,NTH,NPH,NTHS,NPHS
+
+C
+      NLEMT=NLEM+NLEM
+C
+
+      WN=WP/3.0D08
+      SQSFI=1.0D0
+      IF(IDFC.GE.0)GO TO 118
+      VK=XIP*WN
+      WRITE(IW,6714)VK
+      WRITE(IW,6000)
+  118 DO 488 JXI=1,NXI
+      JAW=1
+      WRITE(IW,6991)JXI
+      XI=XIV(JXI)
+      IF(IDFC.LT.0)GO TO 119
+      VK=XI*WN
+      VKARG=VK
+      WRITE(IW,6715)VK,XI
+      GO TO 120
+  119 VKARG=XI*VK
+      SQSFI=1.0D0/(XI*XI)
+      WRITE(IW,6716)XI
+  120 CALL HJV(EXRI,VKARG,JER,LCALC,ARG)
+      IF(JER.EQ.0)GO TO 122
+      WRITE(IW,6650)
+      GO TO 490
+  122 DO 132 I=1,NSPH
+      IOGI=IOG(I)
+      IF(IOGI.EQ.I)GO TO 125
+      DO 123 L=1,LI
+      RMI(L,I)=RMI(L,IOGI)
+  123 REI(L,I)=REI(L,IOGI)
+      GO TO 132
+  125 NSH=NSHL(I)
+      ICI=(NSH+1)/2
+      IF(IDFC.NE.0)GO TO 127
+      READ(ITIN)(DC0(IC),IC=1,ICI)
+      GO TO 130
+  127 IF(JXI.EQ.1)READ(ITIN)(DC0M(I,IC),IC=1,ICI)
+      DO 128 IC=1,ICI
+  128 DC0(IC)=DC0M(I,IC)
+  130 IF(MOD(NSH,2).EQ.0)DC0(ICI+1)=EXDC
+      CALL DME(LI,I,NPNT,NPNTTS,VKARG,EXDC,EXRI,JER,LCALC,ARG)
+      IF(JER.EQ.0)GO TO 132
+      WRITE(IW,6750)
+      GO TO 490
+  132 CONTINUE
+      CALL CMS(AM)
+      NDIT=2*NSPH*NLIM
+      CALL LUCIN(AM,MXNDM,NDIT,JER)
+      IF(JER.EQ.1)GO TO 495
+      CALL ZTM(AM)
+
+C
+      IF(IDFC.LT.0)GO TO 156
+      IF(JXI.NE.JWTM)GO TO 156
+      IS=1
+      ITS=8
+      OPEN(ITS,FILE='TTMS',FORM='UNFORMATTED',STATUS='UNKNOWN')
+      WRITE(ITS)IS,LM
+      WRITE(ITS)VK,EXRI
+      WRITE(ITS)((AM0M(I,J),J=1,NLEMT),I=1,NLEMT)
+      CLOSE(ITS)
+  156 CONTINUE
+C
+
+      IF(INPOL.NE.0)GO TO 158
+      WRITE(IW,6066)
+      GO TO 160
+  158 WRITE(IW,6061)
+  160 CS0=0.25D0*VK*VK*VK/DACOS(C0)
+      CALL SCR0(VK,EXRI)
+      SQK=VK*VK*EXDC
+      CALL APS(ZPV,LI,NSPH,IOG,RMI,REI,SQK,GAPS)
+      CALL RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
+      IF(LI.NE.LE)WRITE(IW,6068)
+      DO 170 I=1,NSPH
+      IF(IOG(I).LT.I)GO TO 170
+      ALBEDS=SSCS(I)/SEXS(I)
+      SQSCS(I)=SQSCS(I)*SQSFI
+      SQABS(I)=SQABS(I)*SQSFI
+      SQEXS(I)=SQEXS(I)*SQSFI
+      WRITE(IW,6070)I
+      IF(NSHL(I).EQ.1)GO TO 162
+      WRITE(IW,6055)VSZ(I)
+      GO TO 164
+  162 WRITE(IW,6050)VSZ(I),VKT(I)
+  164 WRITE(IW,6100)
+      WRITE(IW,6150)SSCS(I),SABS(I),SEXS(I),ALBEDS
+      WRITE(IW,6200)
+      WRITE(IW,6150)SQSCS(I),SQABS(I),SQEXS(I)
+      WRITE(IW,6240)FSAS(I)
+      CSCH=2.0D0*VK*SQSFI/GCSV(I)
+      S0=FSAS(I)*EXRI
+      QSCHU=DIMAG(S0)*CSCH
+      PSCHU=DREAL(S0)*CSCH
+      S0MAG=CDABS(S0)*CS0
+      WRITE(IW,6930)QSCHU,PSCHU,S0MAG
+      RAPR=SEXS(I)-GAPS(I)
+      COSAV=GAPS(I)/SSCS(I)
+      WRITE(IW,6970)COSAV,RAPR
+      WRITE(IW,6977)1,TQSE(1,I),TQSS(1,I)
+      WRITE(IW,6977)2,TQSE(2,I),TQSS(2,I)
+  170 CONTINUE
+      WRITE(IW,6250)TFSAS
+      CSCH=2.0D0*VK*SQSFI/GCS
+      S0=TFSAS*EXRI
+      QSCHU=DIMAG(S0)*CSCH
+      PSCHU=DREAL(S0)*CSCH
+      S0MAG=CDABS(S0)*CS0
+      WRITE(IW,6930)QSCHU,PSCHU,S0MAG
+      WRITE(IT)VK
+      CALL PCRSM0(VK,EXRI,INPOL)
+      CALL APCRA(ZPV,LE,AM0M,INPOL,SQK,GAPM,GAPPM)
+      TH=TH1
+      DO 486 JTH=1,NTH
+      PH=PH1
+      DO 484 JPH=1,NPH
+      IF((NK.EQ.1).AND.(JXI.GT.1))GO TO 180
+      CALL UPVMP(TH,PH,0,COST,SINT,COSP,SINP,U,UPMP,UNMP)
+      IF(ISAM.LT.0)GO TO 184
+      CALL WMAMP
+     1(0,COST,SINT,COSP,SINP,INPOL,LE,0,NSPH,ARGI,U,UPMP,UNMP)
+      GO TO 182
+  180 IF(ISAM.LT.0)GO TO 184
+  182 CALL APC(ZPV,LE,AM0M,W,SQK,GAP,GAPP)
+      CALL RABA(LE,AM0M,W,TQCE,TQCPE,TQCS,TQCPS)
+      JW=1
+  184 THSL=THS1
+      DO 482 JTHS=1,NTHS
+      THS=THSL
+      ICSPNV=0
+      IF(ISAM.GT.1)THS=TH+THSCA
+      IF(ISAM.LT.1)GO TO 186
+      PHSPH=C0
+      IF((THS.LT.C0).OR.(THS.GT.180.0D0))PHSPH=180.0D0
+      IF(THS.LT.C0)THS=-THS
+      IF(THS.GT.180.0D0)THS=360.0D0-THS
+      IF(PHSPH.NE.C0)ICSPNV=1
+  186 PHS=PHS1
+      DO 480 JPHS=1,NPHS
+      IF(ISAM.LT.1)GO TO 188
+      PHS=PH+PHSPH
+      IF(PHS.GE.360.0D0)PHS=PHS-360.0D0
+  188 IF((NKS.EQ.1).AND.((JXI.GT.1).OR.(JTH.GT.1).OR.(JPH.GT.1)))
+     1GO TO 190
+      CALL UPVMP(THS,PHS,ICSPNV,COSTS,SINTS,COSPS,SINPS,US,UPSMP,UNSMP)
+      IF(ISAM.GE.0)CALL WMAMP
+     1(2,COSTS,SINTS,COSPS,SINPS,INPOL,LE,0,NSPH,ARGS,US,UPSMP,UNSMP)
+  190 IF((NKKS.EQ.1).AND.(JXI.GT.1))GO TO 194
+      CALL UPVSP(U,UPMP,UNMP,US,UPSMP,UNSMP,UP,UN,UPS,UNS,DUK,
+     1ISQ,IBF,SCAN,CFMP,SFMP,CFSP,SFSP)
+      IF(ISAM.GE.0)GO TO 192
+      CALL WMASP(COST,SINT,COSP,SINP,COSTS,SINTS,COSPS,SINPS,
+     1U,UP,UN,US,UPS,UNS,ISQ,IBF,INPOL,LE,0,NSPH,ARGI,ARGS)
+      GO TO 194
+  192 DO 193 I=1,3
+      UP(I)=UPMP(I)
+      UN(I)=UNMP(I)
+      UPS(I)=UPSMP(I)
+  193 UNS(I)=UNSMP(I)
+  194 IF(IAVM.EQ.1)CALL CRSM1(VK,EXRI)
+      IF(ISAM.GE.0)GO TO 196
+      CALL APC(ZPV,LE,AM0M,W,SQK,GAP,GAPP)
+      CALL RABA(LE,AM0M,W,TQCE,TQCPE,TQCS,TQCPS)
+      JW=1
+  196 WRITE(IT)TH,PH,THS,PHS
+      WRITE(IT)SCAN
+      IF(JAW.EQ.0)GO TO 212
+      JAW=0
+      CALL MEXTC(VK,EXRI,FSACM,CEXTLR,CEXT)
+      WRITE(IT)((CEXT(I,J),J=1,4),I=1,4)
+      WRITE(IT)(SCSCM(I),SCSCPM(I),ECSCM(I),ECSCPM(I),I=1,2)
+      WRITE(IT)((GAPM(I,J),GAPPM(I,J),J=1,2),I=1,3)
+      WRITE(IW,6460)IAVM
+      JLR=2
+      DO 210 ILR=1,2
+      IPOL=(-1)**ILR
+      IF(ILR.EQ.2)JLR=1
+      EXTSM=ECSCM(ILR)
+      QEXTM=EXTSM*SQSFI/GCS
+      EXTRM=EXTSM/ECS
+      SCASM=SCSCM(ILR)
+      ALBDM=SCASM/EXTSM
+      QSCAM=SCASM*SQSFI/GCS
+      SCARM=SCASM/SCS
+      ABSSM=EXTSM-SCASM
+      QABSM=ABSSM*SQSFI/GCS
+      IF(DABS(ACS/ECS).GT.1.0D-06)GO TO 202
+      ABSRM=1.0D0
+      GO TO 204
+  202 ABSRM=ABSSM/ACS
+  204 S0M=FSACM(ILR,ILR)*EXRI
+      QSCHUM=DIMAG(S0M)*CSCH
+      PSCHUM=DREAL(S0M)*CSCH
+      S0MAGM=CDABS(S0M)*CS0
+      RFINRM=DREAL(FSACM(ILR,ILR))/DREAL(TFSAS)
+      EXTCRM=DIMAG(FSACM(ILR,ILR))/DIMAG(TFSAS)
+      IF(INPOL.NE.0)GO TO 206
+      WRITE(IW,6065)IPOL
+      GO TO 208
+  206 WRITE(IW,6060)IPOL
+  208 WRITE(IW,6500)
+      WRITE(IW,6150)SCASM,ABSSM,EXTSM,ALBDM
+      WRITE(IW,6600)
+      WRITE(IW,6150)QSCAM,QABSM,QEXTM
+      WRITE(IW,6700)
+      WRITE(IW,6150)SCARM,ABSRM,EXTRM
+      WRITE(IW,6800)ILR,ILR,FSACM(ILR,ILR),JLR,ILR,FSACM(JLR,ILR)
+      WRITE(IW,6850)ILR,ILR,RFINRM,ILR,ILR,EXTCRM
+      WRITE(IW,6930)QSCHUM,PSCHUM,S0MAGM
+      RAPR=ECSCM(ILR)-GAPM(3,ILR)
+      COSAV=GAPM(3,ILR)/SCSCM(ILR)
+      FZ=RAPR
+      WRITE(IW,6970)COSAV,RAPR
+      WRITE(IW,6974)FZ
+  210 CONTINUE
+      RMBRIF=(DREAL(FSACM(1,1))-DREAL(FSACM(2,2)))/DREAL(FSACM(1,1))
+      RMDCHR=(DIMAG(FSACM(1,1))-DIMAG(FSACM(2,2)))/DIMAG(FSACM(1,1))
+      WRITE(IW,6940)RMBRIF
+      WRITE(IW,6945)RMDCHR
+  212 WRITE(IW,6992)JTH,JPH,JTHS,JPHS
+      WRITE(IW,6720)TH,PH,THS,PHS
+      WRITE(IW,6740)SCAN
+      WRITE(IW,6721)CFMP,SFMP
+      WRITE(IW,6722)CFSP,SFSP
+      IF(ISAM.LT.0)GO TO 214
+      WRITE(IW,6166)UN(1),UN(2),UN(3)
+      WRITE(IW,6167)UNS(1),UNS(2),UNS(3)
+      GO TO 220
+  214 WRITE(IW,6165)UN(1),UN(2),UN(3)
+  220 IF(INPOL.NE.0)GO TO 222
+      WRITE(IW,6066)
+      GO TO 224
+  222 WRITE(IW,6061)
+  224 CALL SCR2(VK,VKARG,EXRI,DUK)
+      IF(LI.NE.LE)WRITE(IW,6069)
+      DO 226 I=1,NSPH
+      IF(IOG(I).LT.I)GO TO 226
+      WRITE(IW,6070)I
+      WRITE(IW,6242)SAS(I,1,1),SAS(I,2,1)
+      WRITE(IW,6243)SAS(I,1,2),SAS(I,2,2)
+      DO 225 J=1,16
+  225 VINT(J)=VINTS(I,J)
+      CALL MMULC(VINT,CMULLR,CMUL)
+      WRITE(IW,6331)
+      WRITE(IW,6350)((CMUL(I1,J),J=1,4),I1=1,4)
+      WRITE(IW,6336)
+      WRITE(IW,6350)((CMULLR(I1,J),J=1,4),I1=1,4)
+  226 CONTINUE
+      WRITE(IW,6252)TSAS(1,1),TSAS(2,1)
+      WRITE(IW,6253)TSAS(1,2),TSAS(2,2)
+      WRITE(IW,6080)
+      CALL PCROS(VK,EXRI)
+      CALL MEXTC(VK,EXRI,FSAC,CEXTLR,CEXT)
+      CALL MMULC(VINT,CMULLR,CMUL)
+      IF(JW.EQ.0)GO TO 254
+      JW=0
+      WRITE(IT)((CEXT(I,J),J=1,4),I=1,4)
+      WRITE(IT)(SCSC(I),SCSCP(I),ECSC(I),ECSCP(I),I=1,2)
+      WRITE(IT)((GAP(I,J),GAPP(I,J),J=1,2),I=1,3)
+      WRITE(IT)((TQCE(I,J),TQCPE(I,J),J=1,3),I=1,2)
+      WRITE(IT)((TQCS(I,J),TQCPS(I,J),J=1,3),I=1,2)
+      WRITE(IT)(U(I),UP(I),UN(I),I=1,3)
+  254 WRITE(IT)(VINT(I),I=1,16)
+      WRITE(IT)((CMUL(I,J),J=1,4),I=1,4)
+      JLR=2
+      DO 290 ILR=1,2
+      IPOL=(-1)**ILR
+      IF(ILR.EQ.2)JLR=1
+      EXTSEC=ECSC(ILR)
+      QEXT=EXTSEC*SQSFI/GCS
+      EXTRAT=EXTSEC/ECS
+      SCASEC=SCSC(ILR)
+      ALBEDC=SCASEC/EXTSEC
+      QSCA=SCASEC*SQSFI/GCS
+      SCARAT=SCASEC/SCS
+      ABSSEC=EXTSEC-SCASEC
+      QABS=ABSSEC*SQSFI/GCS
+      IF(DABS(ACS/ECS).GT.1.0D-06)GO TO 260
+      ABSRAT=1.0D0
+      GO TO 262
+  260 ABSRAT=ABSSEC/ACS
+  262 S0=FSAC(ILR,ILR)*EXRI
+      QSCHU=DIMAG(S0)*CSCH
+      PSCHU=DREAL(S0)*CSCH
+      S0MAG=CDABS(S0)*CS0
+      REFINR=DREAL(FSAC(ILR,ILR))/DREAL(TFSAS)
+      EXTCOR=DIMAG(FSAC(ILR,ILR))/DIMAG(TFSAS)
+      IF(INPOL.NE.0)GO TO 273
+      WRITE(IW,6065)IPOL
+      GO TO 275
+  273 WRITE(IW,6060)IPOL
+  275 WRITE(IW,6500)
+      WRITE(IW,6150)SCASEC,ABSSEC,EXTSEC,ALBEDC
+      WRITE(IW,6600)
+      WRITE(IW,6150)QSCA,QABS,QEXT
+      WRITE(IW,6700)
+      WRITE(IW,6150)SCARAT,ABSRAT,EXTRAT
+      WRITE(IW,6800)ILR,ILR,FSAC(ILR,ILR),JLR,ILR,FSAC(JLR,ILR)
+      WRITE(IW,6810)ILR,ILR,SAC(ILR,ILR),JLR,ILR,SAC(JLR,ILR)
+      WRITE(IW,6850)ILR,ILR,REFINR,ILR,ILR,EXTCOR
+      WRITE(IW,6930)QSCHU,PSCHU,S0MAG
+      IF((ISAM.GE.0).AND.((JTHS.GT.1).OR.(JPHS.GT.1)))GO TO 290
+      GAPV(1)=GAP(1,ILR)
+      GAPV(2)=GAP(2,ILR)
+      GAPV(3)=GAP(3,ILR)
+      EXTINS=ECSC(ILR)
+      SCATTS=SCSC(ILR)
+      CALL RFTR
+     1(U,UP,UN,GAPV,EXTINS,SCATTS,RAPR,COSAV,FP,FN,FK,FX,FY,FZ)
+      WRITE(IW,6970)COSAV,RAPR
+      WRITE(IW,6973)FP,FN,FK
+      WRITE(IW,6971)FX,FY,FZ
+      TQEV(1)=TQCE(ILR,1)
+      TQEV(2)=TQCE(ILR,2)
+      TQEV(3)=TQCE(ILR,3)
+      TQSV(1)=TQCS(ILR,1)
+      TQSV(2)=TQCS(ILR,2)
+      TQSV(3)=TQCS(ILR,3)
+      CALL TQR(U,UP,UN,TQEV,TQSV,TEP,TEN,TEK,TSP,TSN,TSK)
+      WRITE(IW,6983)TEP,TEN,TEK
+      WRITE(IW,6984)TSP,TSN,TSK
+      WRITE(IW,6981)TQCE(ILR,1),TQCE(ILR,2),TQCE(ILR,3)
+      WRITE(IW,6982)TQCS(ILR,1),TQCS(ILR,2),TQCS(ILR,3)
+  290 CONTINUE
+      RBIRIF=(DREAL(FSAC(1,1))-DREAL(FSAC(2,2)))/DREAL(FSAC(1,1))
+      RDICHR=(DIMAG(FSAC(1,1))-DIMAG(FSAC(2,2)))/DIMAG(FSAC(1,1))
+      WRITE(IW,6940)RBIRIF
+      WRITE(IW,6945)RDICHR
+      WRITE(IW,6332)
+      WRITE(IW,6350)((CMUL(I,J),J=1,4),I=1,4)
+      WRITE(IW,6337)
+      WRITE(IW,6350)((CMULLR(I,J),J=1,4),I=1,4)
+      IF(IAVM.EQ.0)GO TO 420
+      CALL MMULC(VINTM,CMULLR,CMUL)
+      WRITE(IT)(VINTM(I),I=1,16)
+      WRITE(IT)((CMUL(I,J),J=1,4),I=1,4)
+      WRITE(IW,6460)IAVM
+      IF(INPOL.NE.0)GO TO 316
+      WRITE(IW,6066)
+      GO TO 318
+  316 WRITE(IW,6061)
+  318 WRITE(IW,6332)
+      WRITE(IW,6350)((CMUL(I,J),J=1,4),I=1,4)
+      WRITE(IW,6337)
+      WRITE(IW,6350)((CMULLR(I,J),J=1,4),I=1,4)
+  420 CONTINUE
+  480 IF(ISAM.LT.1)PHS=PHS+PHSSTP
+  482 IF(ISAM.LE.1)THSL=THSL+THSSTP
+  484 PH=PH+PHSTP
+  486 TH=TH+THSTP
+  488 CONTINUE
+      GO TO 500
+  490 WRITE(IW,6550)JER,LCALC,ARG
+      GO TO 500  
+  495 WRITE(IW,6300)
+      GO TO 500
+  497 WRITE(IW,*)' WRONG INPUT TAPE'
+  500 CONTINUE
+      CLOSE(IW)
+      CLOSE(IT)
+      CLOSE(ITIN)
+      STOP
+      END
+      SUBROUTINE RABAS(INPOL,LI,NSPH,IOG,RMI,REI,TQSE,TQSPE,TQSS,TQSPS)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION IOG(NSPH),RMI(LI,NSPH),REI(LI,NSPH),
+CCC  1TQSE(2,NSPH),TQSPE(2,NSPH),TQSS(2,NSPH),TQSPS(2,NSPH)
+      DIMENSION IOG(6),RMI(8,6),REI(8,6),
+     1TQSE(2,6),TQSPE(2,6),TQSS(2,6),TQSPS(2,6)
+      COMPLEX*16 RMI,REI,TQSPE,TQSPS
+      COMPLEX*16 CC0,UIM,RM,RE,PCE,PCS
+      C0=0.0D0
+      PIG2=DACOS(C0)*4.0D0
+      CC0=(0.0D0,0.0D0)
+      UIM=(0.0D0,1.0D0)
+      DO 80 I=1,NSPH
+      IF(IOG(I).LT.I)GO TO 80
+      TQSE(1,I)=C0
+      TQSE(2,I)=C0
+      TQSPE(1,I)=CC0
+      TQSPE(2,I)=CC0
+      TQSS(1,I)=C0
+      TQSS(2,I)=C0
+      TQSPS(1,I)=CC0
+      TQSPS(2,I)=CC0
+      DO 70 L=1,LI
+      FL=L+L+1
+      RM=1.0D0/RMI(L,I)
+      RMM=DREAL(RM*DCONJG(RM))
+      RE=1.0D0/REI(L,I)
+      REM=DREAL(RE*DCONJG(RE))
+      IF(INPOL.NE.0)GO TO 40
+      PCE=((RM+RE)*UIM)*FL
+      PCS=((RMM+REM)*FL)*UIM
+      TQSPE(1,I)=TQSPE(1,I)-PCE
+      TQSPE(2,I)=TQSPE(2,I)+PCE
+      TQSPS(1,I)=TQSPS(1,I)-PCS
+      TQSPS(2,I)=TQSPS(2,I)+PCS
+      GO TO 70
+   40 CE=DREAL(RM+RE)*FL
+      CS=(RMM+REM)*FL
+      TQSE(1,I)=TQSE(1,I)-CE
+      TQSE(2,I)=TQSE(2,I)+CE
+      TQSS(1,I)=TQSS(1,I)-CS
+      TQSS(2,I)=TQSS(2,I)+CS
+   70 CONTINUE
+      IF(INPOL.NE.0)GO TO 75
+      TQSPE(1,I)=TQSPE(1,I)*PIG2
+      TQSPE(2,I)=TQSPE(2,I)*PIG2
+      TQSPS(1,I)=TQSPS(1,I)*PIG2
+      TQSPS(2,I)=TQSPS(2,I)*PIG2
+      GO TO 80
+   75 TQSE(1,I)=TQSE(1,I)*PIG2
+      TQSE(2,I)=TQSE(2,I)*PIG2
+      TQSS(1,I)=TQSS(1,I)*PIG2
+      TQSS(2,I)=TQSS(2,I)*PIG2
+   80 CONTINUE
+      RETURN
+      END
+      SUBROUTINE TQR(U,UP,UN,TQEV,TQSV,TEP,TEN,TEK,TSP,TSN,TSK)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION U(3),UP(3),UN(3),TQEV(3),TQSV(3)
+      TEP=UP(1)*TQEV(1)+UP(2)*TQEV(2)+UP(3)*TQEV(3)
+      TEN=UN(1)*TQEV(1)+UN(2)*TQEV(2)+UN(3)*TQEV(3)
+      TEK=U(1)*TQEV(1)+U(2)*TQEV(2)+U(3)*TQEV(3)
+      TSP=UP(1)*TQSV(1)+UP(2)*TQSV(2)+UP(3)*TQSV(3)
+      TSN=UN(1)*TQSV(1)+UN(2)*TQSV(2)+UN(3)*TQSV(3)
+      TSK=U(1)*TQSV(1)+U(2)*TQSV(2)+U(3)*TQSV(3)
+      RETURN
+      END
+      SUBROUTINE RABA(LE,AM0M,W,TQCE,TQCPE,TQCS,TQCPS)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION AM0M(NLEMT,NLEMT),W(NLEMT,4)
+      DIMENSION AM0M(160,160),W(160,4),
+     1TQCE(2,3),TQCPE(2,3),TQCS(2,3),TQCPS(2,3)
+      COMPLEX*16 AM0M,W,TQCPE,TQCPS
+CCC   DIMENSION A(NLEMT,2)
+      DIMENSION A(160,2),CTQCE(2,3),CTQCS(2,3)
+      COMPLEX*16 A,CTQCE,CTQCS
+      COMPLEX*16 CC0,UIM,ACW,ACWP,ACA,ACAP,C1,C2,C3
+      CC0=(0.0D0,0.0D0)
+      UIM=(0.0D0,1.0D0)
+      SQ2I=1.0D0/DSQRT(2.0D0)
+      NLEM=LE*(LE+2)
+      NLEMT=NLEM+NLEM
+      DO 20 I=1,NLEMT
+      C1=CC0
+      C2=CC0
+      DO 10 J=1,NLEMT
+      C1=C1+AM0M(I,J)*W(J,1)
+      C2=C2+AM0M(I,J)*W(J,2)
+   10 CONTINUE
+      A(I,1)=C1
+      A(I,2)=C2
+   20 CONTINUE
+      JPO=2
+      DO 70 IPO=1,2
+      IF(IPO.EQ.2)JPO=1
+      CTQCE(IPO,1)=CC0
+      CTQCE(IPO,2)=CC0
+      CTQCE(IPO,3)=CC0
+      TQCPE(IPO,1)=CC0
+      TQCPE(IPO,2)=CC0
+      TQCPE(IPO,3)=CC0
+      CTQCS(IPO,1)=CC0
+      CTQCS(IPO,2)=CC0
+      CTQCS(IPO,3)=CC0
+      TQCPS(IPO,1)=CC0
+      TQCPS(IPO,2)=CC0
+      TQCPS(IPO,3)=CC0   
+      DO 60 L=1,LE
+      LPO=L+1
+      IL=L*LPO
+      LTPO=L+LPO
+      DO 60 IM=1,LTPO
+      M=IM-LPO
+      I=M+IL
+      IE=I+NLEM      
+      MMMU=M+1
+      MMMMU=IABS(MMMU)
+      IF(MMMMU.GT.L)GO TO 30
+      IMMU=MMMU+IL
+      IMMUE=IMMU+NLEM
+      RMU=(L+MMMU)*(L-M)
+      RMU=-DSQRT(RMU)*SQ2I
+      ACW=DCONJG(A(I,IPO))*W(IMMU,IPO)+DCONJG(A(IE,IPO))*W(IMMUE,IPO)
+      ACWP=DCONJG(A(I,IPO))*W(IMMU,JPO)+DCONJG(A(IE,IPO))*W(IMMUE,JPO)
+      ACA=DCONJG(A(I,IPO))*A(IMMU,IPO)+DCONJG(A(IE,IPO))*A(IMMUE,IPO)
+      ACAP=DCONJG(A(I,IPO))*A(IMMU,JPO)+DCONJG(A(IE,IPO))*A(IMMUE,JPO)
+      CTQCE(IPO,1)=CTQCE(IPO,1)+ACW*RMU
+      TQCPE(IPO,1)=TQCPE(IPO,1)+ACWP*RMU
+      CTQCS(IPO,1)=CTQCS(IPO,1)+ACA*RMU
+      TQCPS(IPO,1)=TQCPS(IPO,1)+ACAP*RMU
+   30 RMU=-M
+      ACW=DCONJG(A(I,IPO))*W(I,IPO)+DCONJG(A(IE,IPO))*W(IE,IPO)
+      ACWP=DCONJG(A(I,IPO))*W(I,JPO)+DCONJG(A(IE,IPO))*W(IE,JPO)
+      ACA=DCONJG(A(I,IPO))*A(I,IPO)+DCONJG(A(IE,IPO))*A(IE,IPO)
+      ACAP=DCONJG(A(I,IPO))*A(I,JPO)+DCONJG(A(IE,IPO))*A(IE,JPO)
+      CTQCE(IPO,2)=CTQCE(IPO,2)+ACW*RMU
+      TQCPE(IPO,2)=TQCPE(IPO,2)+ACWP*RMU
+      CTQCS(IPO,2)=CTQCS(IPO,2)+ACA*RMU
+      TQCPS(IPO,2)=TQCPS(IPO,2)+ACAP*RMU
+      MMMU=M-1
+      MMMMU=IABS(MMMU)
+      IF(MMMMU.GT.L)GO TO 50
+      IMMU=MMMU+IL
+      IMMUE=IMMU+NLEM
+      RMU=(L-MMMU)*(L+M)
+      RMU=DSQRT(RMU)*SQ2I
+      ACW=DCONJG(A(I,IPO))*W(IMMU,IPO)+DCONJG(A(IE,IPO))*W(IMMUE,IPO)
+      ACWP=DCONJG(A(I,IPO))*W(IMMU,JPO)+DCONJG(A(IE,IPO))*W(IMMUE,JPO)
+      ACA=DCONJG(A(I,IPO))*A(IMMU,IPO)+DCONJG(A(IE,IPO))*A(IMMUE,IPO)
+      ACAP=DCONJG(A(I,IPO))*A(IMMU,JPO)+DCONJG(A(IE,IPO))*A(IMMUE,JPO)
+      CTQCE(IPO,3)=CTQCE(IPO,3)+ACW*RMU
+      TQCPE(IPO,3)=TQCPE(IPO,3)+ACWP*RMU
+      CTQCS(IPO,3)=CTQCS(IPO,3)+ACA*RMU
+      TQCPS(IPO,3)=TQCPS(IPO,3)+ACAP*RMU
+   50 CONTINUE
+   60 CONTINUE
+   70 CONTINUE
+      DO 78 IPO=1,2
+      TQCE(IPO,1)=DREAL(CTQCE(IPO,1)-CTQCE(IPO,3))*SQ2I
+      TQCE(IPO,2)=DREAL((CTQCE(IPO,1)+CTQCE(IPO,3))*UIM)*SQ2I
+      TQCE(IPO,3)=DREAL(CTQCE(IPO,2))
+      C1=TQCPE(IPO,1)
+      C2=TQCPE(IPO,2)
+      C3=TQCPE(IPO,3)
+      TQCPE(IPO,1)=(C1-C3)*SQ2I
+      TQCPE(IPO,2)=(C1+C3)*(UIM*SQ2I)
+      TQCPE(IPO,3)=C2
+      TQCS(IPO,1)=-DREAL(CTQCS(IPO,1)-CTQCS(IPO,3))*SQ2I
+      TQCS(IPO,2)=-DREAL((CTQCS(IPO,1)+CTQCS(IPO,3))*UIM)*SQ2I
+      TQCS(IPO,3)=-DREAL(CTQCS(IPO,2))
+      C1=TQCPS(IPO,1)
+      C2=TQCPS(IPO,2)
+      C3=TQCPS(IPO,3)
+      TQCPS(IPO,1)=-(C1-C3)*SQ2I
+      TQCPS(IPO,2)=-(C1+C3)*(UIM*SQ2I)
+      TQCPS(IPO,3)=-C2
+   78 CONTINUE
+      RETURN
+      END
+      SUBROUTINE RFTR
+     1(U,UP,UN,GAPV,EXTINS,SCATTS,RAPR,COSAV,FP,FN,FK,FX,FY,FZ)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION U(3),UP(3),UN(3),GAPV(3)
+      FK=U(1)*GAPV(1)+U(2)*GAPV(2)+U(3)*GAPV(3)
+      RAPR=EXTINS-FK
+      COSAV=FK/SCATTS
+      FP=-(UP(1)*GAPV(1)+UP(2)*GAPV(2)+UP(3)*GAPV(3))
+      FN=-(UN(1)*GAPV(1)+UN(2)*GAPV(2)+UN(3)*GAPV(3))
+      FK=RAPR
+      FX=U(1)*EXTINS-GAPV(1)
+      FY=U(2)*EXTINS-GAPV(2)
+      FZ=U(3)*EXTINS-GAPV(3)
+      RETURN
+      END
+      SUBROUTINE MEXTC(VK,EXRI,FSAC,CEXTLR,CEXT)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION FSAC(2,2),CEXTLR(4,4),CEXT(4,4)
+      COMPLEX*16 FSAC
+      FA11R=DREAL(FSAC(1,1))
+      FA11I=DIMAG(FSAC(1,1))
+      FA21R=DREAL(FSAC(2,1))
+      FA21I=DIMAG(FSAC(2,1))
+      FA12R=DREAL(FSAC(1,2))
+      FA12I=DIMAG(FSAC(1,2))
+      FA22R=DREAL(FSAC(2,2))
+      FA22I=DIMAG(FSAC(2,2))
+      CEXTLR(1,1)=FA11I*2.0D0
+      CEXTLR(1,2)=0.0D0
+      CEXTLR(1,3)=-FA12I
+      CEXTLR(1,4)=-FA12R
+      CEXTLR(2,1)=0.0D0
+      CEXTLR(2,2)=FA22I*2.0D0
+      CEXTLR(2,3)=-FA21I
+      CEXTLR(2,4)=FA21R
+      CEXTLR(3,1)=-FA21I*2.0D0
+      CEXTLR(3,2)=-FA12I*2.0D0
+      CEXTLR(3,3)=FA11I+FA22I
+      CEXTLR(3,4)=FA22R-FA11R
+      CEXTLR(4,1)=FA21R*2.0D0
+      CEXTLR(4,2)=-FA12R*2.0D0
+      CEXTLR(4,3)=FA11R-FA22R
+      CEXTLR(4,4)=CEXTLR(3,3)
+      CEXT(1,1)=CEXTLR(4,4)
+      CEXT(2,2)=CEXTLR(4,4)
+      CEXT(3,3)=CEXTLR(4,4)
+      CEXT(3,4)=CEXTLR(3,4)
+      CEXT(4,3)=CEXTLR(4,3)
+      CEXT(4,4)=CEXTLR(4,4)
+      CEXT(1,2)=FA11I-FA22I
+      CEXT(1,3)=-FA12I-FA21I
+      CEXT(1,4)=FA21R-FA12R
+      CEXT(2,1)=CEXT(1,2)
+      CEXT(2,3)=FA21I-FA12I
+      CEXT(4,2)=FA12R+FA21R
+      CEXT(2,4)=-CEXT(4,2)
+      CEXT(3,1)=CEXT(1,3)
+      CEXT(3,2)=-CEXT(2,3)
+      CEXT(4,1)=CEXT(1,4)
+      CKM=VK/EXRI
+      DO 10 I=1,4
+      DO 10 J=1,4
+      CEXTLR(I,J)=CEXTLR(I,J)*CKM
+   10 CEXT(I,J)=CEXT(I,J)*CKM
+      RETURN
+      END
+      SUBROUTINE MMULC(VINT,CMULLR,CMUL)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION VINT(16),CMULLR(4,4),CMUL(4,4)
+      COMPLEX*16 VINT
+      SM2=DREAL(VINT(1))
+      S24=DREAL(VINT(2))
+      D24=DIMAG(VINT(2))
+      SM4=DREAL(VINT(6))
+      S23=DREAL(VINT(9))
+      D32=DIMAG(VINT(9))
+      S34=DREAL(VINT(10))
+      D34=DIMAG(VINT(10))
+      SM3=DREAL(VINT(11))
+      S31=DREAL(VINT(12))
+      D31=DIMAG(VINT(12))
+      S21=DREAL(VINT(13))
+      D12=DIMAG(VINT(13))
+      S41=DREAL(VINT(14))
+      D14=DIMAG(VINT(14))
+      SM1=DREAL(VINT(16))
+      CMULLR(1,1)=SM2
+      CMULLR(1,2)=SM3
+      CMULLR(1,3)=-S23
+      CMULLR(1,4)=-D32
+      CMULLR(2,1)=SM4
+      CMULLR(2,2)=SM1
+      CMULLR(2,3)=-S41
+      CMULLR(2,4)=-D14
+      CMULLR(3,1)=-S24*2.0D0
+      CMULLR(3,2)=-S31*2.0D0
+      CMULLR(3,3)=S21+S34
+      CMULLR(3,4)=D34+D12
+      CMULLR(4,1)=-D24*2.0D0
+      CMULLR(4,2)=-D31*2.0D0
+      CMULLR(4,3)=D34-D12
+      CMULLR(4,4)=S21-S34
+      CMUL(1,1)=(SM2+SM3+SM4+SM1)*0.5D0
+      CMUL(1,2)=(SM2-SM3+SM4-SM1)*0.5D0
+      CMUL(1,3)=-S23-S41
+      CMUL(1,4)=-D32-D14
+      CMUL(2,1)=(SM2+SM3-SM4-SM1)*0.5D0
+      CMUL(2,2)=(SM2-SM3-SM4+SM1)*0.5D0
+      CMUL(2,3)=-S23+S41
+      CMUL(2,4)=-D32+D14
+      CMUL(3,1)=-S24-S31
+      CMUL(3,2)=-S24+S31
+      CMUL(3,3)=S21+S34
+      CMUL(3,4)=D34+D12
+      CMUL(4,1)=-D24-D31
+      CMUL(4,2)=-D24+D31
+      CMUL(4,3)=D34-D12
+      CMUL(4,4)=S21-S34
+      RETURN
+      END
+      SUBROUTINE SCR0(VK,EXRI)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+      COMPLEX*16 TFSAS,TSAS
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMPLEX*16 SUM21,RM,RE,CC0,CSAM
+      CC0=(0.0D0,0.0D0)
+      C0=0.0D0
+      EXDC=EXRI*EXRI
+      CCS=4.0D0*DACOS(C0)/(VK*VK)
+      CCCS=CCS/EXDC
+      CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+      SCS=C0
+      ECS=C0
+      ACS=C0
+      TFSAS=CC0
+      DO 14 I=1,NSPH
+      IOGI=IOG(I)
+      IF(IOGI.LT.I)GO TO 12
+      SUMS=C0
+      SUM21=CC0
+      DO 10 L=1,LI
+      FL=L+L+1
+      RM=1.0D0/RMI(L,I)
+      RE=1.0D0/REI(L,I)
+      SUMS=SUMS+DREAL(DCONJG(RM)*RM+DCONJG(RE)*RE)*FL
+   10 SUM21=SUM21+(RM+RE)*FL
+      SUM21=-SUM21
+      SCASEC=CCCS*SUMS
+      EXTSEC=-CCCS*DREAL(SUM21)
+      ABSSEC=EXTSEC-SCASEC
+      SSCS(I)=SCASEC
+      SEXS(I)=EXTSEC
+      SABS(I)=ABSSEC
+      GCSS=GCSV(I)
+      SQSCS(I)=SCASEC/GCSS
+      SQEXS(I)=EXTSEC/GCSS
+      SQABS(I)=ABSSEC/GCSS
+      FSAS(I)=SUM21*CSAM
+   12 SCS=SCS+SSCS(IOGI)
+      ECS=ECS+SEXS(IOGI)
+      ACS=ACS+SABS(IOGI)
+   14 TFSAS=TFSAS+FSAS(IOGI)
+      RETURN
+      END
+      SUBROUTINE SCR2(VK,VKARG,EXRI,DUK)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+      COMPLEX*16 TFSAS,TSAS
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      DIMENSION DUK(3)
+      COMPLEX*16 S11,S21,S12,S22,RM,RE,CC0,CSAM,CPH,PHAS,CC
+      CC0=(0.0D0,0.0D0)
+      C0=0.0D0
+      CCS=1.0D0/(VK*VK)
+      CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+      PIGFSQ=6.4D+1*DACOS(C0)**2
+      CFSQ=4.0D0/(PIGFSQ*CCS*CCS)
+      CPH=(0.0D0,1.0D0)*(EXRI*VKARG)
+      LS=MIN0(LI,LE)
+      TSAS(1,1)=CC0
+      TSAS(2,1)=CC0
+      TSAS(1,2)=CC0
+      TSAS(2,2)=CC0
+      DO 14 I=1,NSPH
+      IOGI=IOG(I)
+      IF(IOGI.LT.I)GO TO 12
+      K=0
+      S11=CC0
+      S21=CC0
+      S12=CC0
+      S22=CC0
+      DO 10 L=1,LS
+      RM=1.0D0/RMI(L,I)
+      RE=1.0D0/REI(L,I)
+      LTPO=L+L+1
+      DO 10 IM=1,LTPO
+      K=K+1
+      KE=K+NLEM
+      S11=S11-W(K,3)*W(K,1)*RM-W(KE,3)*W(KE,1)*RE
+      S21=S21-W(K,4)*W(K,1)*RM-W(KE,4)*W(KE,1)*RE
+      S12=S12-W(K,3)*W(K,2)*RM-W(KE,3)*W(KE,2)*RE
+   10 S22=S22-W(K,4)*W(K,2)*RM-W(KE,4)*W(KE,2)*RE
+      SAS(I,1,1)=S11*CSAM
+      SAS(I,2,1)=S21*CSAM
+      SAS(I,1,2)=S12*CSAM
+      SAS(I,2,2)=S22*CSAM
+   12 PHAS=CDEXP(CPH*(DUK(1)*RXX(I)+DUK(2)*RYY(I)+DUK(3)*RZZ(I)))
+      TSAS(1,1)=TSAS(1,1)+SAS(IOGI,1,1)*PHAS
+      TSAS(2,1)=TSAS(2,1)+SAS(IOGI,2,1)*PHAS
+      TSAS(1,2)=TSAS(1,2)+SAS(IOGI,1,2)*PHAS
+   14 TSAS(2,2)=TSAS(2,2)+SAS(IOGI,2,2)*PHAS
+      DO 24 I=1,NSPH
+      IOGI=IOG(I)
+      IF(IOGI.LT.I)GO TO 24
+      J=0
+      DO 22 IPO1=1,2
+      DO 22 JPO1=1,2
+      CC=DCONJG(SAS(I,JPO1,IPO1))
+      DO 22 IPO2=1,2
+      DO 22 JPO2=1,2
+      J=J+1
+   22 VINTS(I,J)=SAS(I,JPO2,IPO2)*CC*CFSQ
+   24 CONTINUE
+      J=0
+      DO 32 IPO1=1,2
+      DO 32 JPO1=1,2
+      CC=DCONJG(TSAS(JPO1,IPO1))
+      DO 32 IPO2=1,2
+      DO 32 JPO2=1,2
+      J=J+1
+   32 VINTT(J)=TSAS(JPO2,IPO2)*CC*CFSQ
+      RETURN
+      END
+      SUBROUTINE PCROS(VK,EXRI)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMPLEX*16 SUMP,SUM1,SUM2,SUM3,SUM4,AM,AMP,CC,CSAM,CC0
+      CC0=(0.0D0,0.0D0)
+      C0=0.0D0
+      EXDC=EXRI*EXRI
+      CCS=1.0D0/(VK*VK)
+      CCCS=CCS/EXDC
+      CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+      PIGFSQ=6.4D+1*DACOS(C0)**2
+      CFSQ=4.0D0/(PIGFSQ*CCS*CCS)
+      NLEMT=NLEM+NLEM
+      JPO=2
+      DO 18 IPO=1,2
+      IF(IPO.EQ.2)JPO=1
+      IPOPT=IPO+2
+      JPOPT=JPO+2
+      SUM=C0
+      SUMP=CC0
+      SUM1=CC0
+      SUM2=CC0
+      SUM3=CC0
+      SUM4=CC0
+      DO 12 I=1,NLEMT
+      AM=CC0
+      AMP=CC0
+      DO 10 J=1,NLEMT
+      AM=AM+AM0M(I,J)*W(J,IPO)
+   10 AMP=AMP+AM0M(I,J)*W(J,JPO)
+      SUM=SUM+DREAL(DCONJG(AM)*AM)
+      SUMP=SUMP+(DCONJG(AMP)*AM)
+      SUM1=SUM1+DCONJG(W(I,IPO))*AM
+      SUM2=SUM2+DCONJG(W(I,JPO))*AM
+      SUM3=SUM3+W(I,IPOPT)*AM
+   12 SUM4=SUM4+W(I,JPOPT)*AM
+      SCSC(IPO)=CCCS*SUM
+      SCSCP(IPO)=CCCS*SUMP
+      ECSC(IPO)=-CCCS*DREAL(SUM1)
+      ECSCP(IPO)=-CCCS*SUM2
+      FSAC(IPO,IPO)=CSAM*SUM1
+      FSAC(JPO,IPO)=CSAM*SUM2
+      SAC(IPO,IPO)=CSAM*SUM3
+   18 SAC(JPO,IPO)=CSAM*SUM4
+      I=0
+      DO 22 IPO1=1,2
+      DO 22 JPO1=1,2
+      CC=DCONJG(SAC(JPO1,IPO1))
+      DO 22 IPO2=1,2
+      DO 22 JPO2=1,2
+      I=I+1
+   22 VINT(I)=SAC(JPO2,IPO2)*CC*CFSQ
+      RETURN
+      END
+      SUBROUTINE PCRSM0(VK,EXRI,INPOL)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMPLEX*16 SUM1,SUM2,SUM3,SUM4,SUMPD,SUMS1,SUMS2,SUMS3,SUMS4,
+     1CSAM,CC0,UIM
+      CC0=(0.0D0,0.0D0)
+      UIM=(0.0D0,1.0D0)
+      C0=0.0D0
+      EXDC=EXRI*EXRI
+      CCS=4.0D0*DACOS(C0)/(VK*VK)
+      CCCS=CCS/EXDC
+      CSAM=-(CCS/(EXRI*VK))*(0.0D0,0.5D0)
+      SUM2=CC0
+      SUM3=CC0
+      DO 14 I=1,NLEM
+      IE=I+NLEM
+      SUM2=SUM2+(AM0M(I,I)+AM0M(IE,IE))
+   14 SUM3=SUM3+(AM0M(I,IE)+AM0M(IE,I))
+      SUMPI=C0
+      SUMPD=CC0
+      NLEMT=NLEM+NLEM
+      DO 16 I=1,NLEMT
+      DO 16 J=1,NLEM
+      JE=J+NLEM
+      SUMPI=SUMPI+DREAL
+     1(DCONJG(AM0M(I,J))*AM0M(I,J)+DCONJG(AM0M(I,JE))*AM0M(I,JE))
+   16 SUMPD=SUMPD+
+     1(DCONJG(AM0M(I,J))*AM0M(I,JE)+DCONJG(AM0M(I,JE))*AM0M(I,J))
+      IF(INPOL.NE.0)GO TO 18
+      SUM1=SUM2
+      SUM4=SUM3*UIM
+      SUM3=-SUM4
+      SUMS1=SUMPI
+      SUMS2=SUMPI
+      SUMS3=SUMPD*UIM
+      SUMS4=-SUMS3
+      GO TO 20
+   18 SUM1=SUM2+SUM3
+      SUM2=SUM2-SUM3
+      SUM3=CC0
+      SUM4=CC0
+      SUMS1=SUMPI-SUMPD
+      SUMS2=SUMPI+SUMPD
+      SUMS3=CC0
+      SUMS4=CC0
+   20 ECSCM(1)=-CCCS*DREAL(SUM2)
+      ECSCM(2)=-CCCS*DREAL(SUM1)
+      ECSCPM(1)=-CCCS*SUM4
+      ECSCPM(2)=-CCCS*SUM3
+      FSACM(1,1)=CSAM*SUM2
+      FSACM(2,1)=CSAM*SUM4
+      FSACM(2,2)=CSAM*SUM1
+      FSACM(1,2)=CSAM*SUM3
+      SCSCM(1)=CCCS*DREAL(SUMS1)
+      SCSCM(2)=CCCS*DREAL(SUMS2)
+      SCSCPM(1)=CCCS*SUMS3
+      SCSCPM(2)=CCCS*SUMS4
+      RETURN
+      END
+      SUBROUTINE CRSM1(VK,EXRI)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMMON/C6/RAC3J(17)
+CCC   DIMENSION SVF(LE4PO,LE4PO,4),SVW(LE4PO,4,4),SVS(LE4PO,4)
+CCC   LE4PO=LE*4+1
+      DIMENSION SVF(33,33,4),SVW(33,4,4),SVS(33,4)
+      COMPLEX*16 SVF,SVW,SVS
+      COMPLEX*16 CAM,CC0
+      CC0=(0.0D0,0.0D0)
+      C0=0.0D0
+      EXDC=EXRI*EXRI
+      CCS=1.0D0/(VK*VK)
+      PIGFSQ=6.4D+1*DACOS(C0)**2
+      CINT=CCS/(PIGFSQ*EXDC)
+      LETPO=LE+LE+1
+      DO 20 I=1,16
+   20 VINTM(I)=CC0
+      DO 40 LPO=1,LETPO
+      L=LPO-1
+      LTPO=LPO+L
+      IMMN=LETPO-L
+      IMMX=LETPO+L
+      DO 22 IMF=IMMN,IMMX
+      DO 22 IMS=IMMN,IMMX
+      DO 22 IPO=1,4
+   22 SVF(IMF,IMS,IPO)=CC0
+      DO 34 L1=1,LE
+      IL1=L1*(L1+1)
+      DO 34 L2=1,LE
+      IF((L.LT.IABS(L2-L1)).OR.(L.GT.(L2+L1)))GO TO 34
+      IL2=L2*(L2+1)
+      DO 24 IM=IMMN,IMMX
+      DO 24 IPA=1,4
+      SVS(IM,IPA)=CC0
+      DO 24 IPO=1,4
+   24 SVW(IM,IPA,IPO)=CC0
+      DO 30 IM=IMMN,IMMX
+      M=IM-LETPO
+      CALL R3JMR(L,L1,L2,M)
+      M1MNMO=MAX0(-L1,-L2-M)-1
+      NM1=MIN0(L1,L2-M)-M1MNMO
+      DO 30 IM1=1,NM1
+      M1=-IM1-M1MNMO
+      ISN=1
+      IF(MOD(M1,2).NE.0)ISN=-1
+      CG3J=RAC3J(IM1)*ISN
+CCC   CG=(LTPO**(1/2))*CG3J
+      ILM1=IL1+M1
+      ILM2=IL2+M1-M
+      IPA=0
+      DO 30 IPA1=1,2
+      I1=ILM1
+      IF(IPA1.EQ.2)I1=ILM1+NLEM
+      DO 30 IPA2=1,2
+      I2=ILM2
+      IF(IPA2.EQ.2)I2=ILM2+NLEM
+      IPA=IPA+1
+      SVS(IM,IPA)=SVS(IM,IPA)+AM0M(I1,I2)*CG3J
+      IPO=0
+      DO 30 IPO2=1,2
+      DO 30 IPO1=3,4
+      IPO=IPO+1
+   30 SVW(IM,IPA,IPO)=SVW(IM,IPA,IPO)+W(I1,IPO1)*W(I2,IPO2)*CG3J
+      DO 32 IMF=IMMN,IMMX
+      DO 32 IMS=IMMN,IMMX
+      DO 32 IPO=1,4
+      DO 32 IPA=1,4
+   32 SVF(IMF,IMS,IPO)=SVF(IMF,IMS,IPO)+SVW(IMF,IPA,IPO)*SVS(IMS,IPA)
+   34 CONTINUE
+      DO 40 IMF=IMMN,IMMX
+      DO 40 IMS=IMMN,IMMX
+      I=0
+      DO 40 IPO1=1,4
+      CAM=DCONJG(SVF(IMF,IMS,IPO1))
+      DO 40 IPO2=1,4
+      I=I+1
+   40 VINTM(I)=VINTM(I)+(SVF(IMF,IMS,IPO2)*CAM)*LTPO
+      DO 42 I=1,16
+   42 VINTM(I)=VINTM(I)*CINT
+      RETURN
+      END
+      SUBROUTINE APS(ZPV,LI,NSPH,IOG,RMI,REI,SQK,GAPS)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION ZPV(LM,3,2,2),IOG(NSPH),
+CCC  1RMI(LI,NSPH),REI(LI,NSPH),GAPS(NSPH)
+      DIMENSION ZPV(8,3,2,2),IOG(6),RMI(8,6),REI(8,6),GAPS(6)
+      COMPLEX*16 RMI,REI
+      COMPLEX*16 CC0,SUMM,SUME,SUEM,SUEE,SUM
+      CC0=(0.0D0,0.0D0)
+      C0=0.0D0
+      PIGH=DACOS(C0)
+      COFS=PIGH*2.0D0/SQK
+      DO 40 I=1,NSPH
+      IF(IOG(I).LT.I)GO TO 40
+      SUM=CC0
+      DO 30 L=1,LI
+      LTPO=L+L+1
+      DO 30 ILMP=1,3
+      IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LI).AND.(ILMP.EQ.3)))
+     1GO TO 30
+      LMPML=ILMP-2
+      LMP=L+LMPML
+      COFL=LTPO*(LMP+LMP+1)
+      COFL=DSQRT(COFL)
+      SUMM=ZPV(L,ILMP,1,1)/(DCONJG(RMI(L,I))*RMI(LMP,I))
+      SUME=ZPV(L,ILMP,1,2)/(DCONJG(RMI(L,I))*REI(LMP,I))
+      SUEM=ZPV(L,ILMP,2,1)/(DCONJG(REI(L,I))*RMI(LMP,I))
+      SUEE=ZPV(L,ILMP,2,2)/(DCONJG(REI(L,I))*REI(LMP,I))
+      SUM=SUM+(CG1(LMPML,0,L,-1)*(SUMM-SUME-SUEM+SUEE)+
+     1CG1(LMPML,0,L,1)*(SUMM+SUME+SUEM+SUEE))*COFL
+   30 CONTINUE
+      GAPS(I)=DREAL(SUM)*COFS
+   40 CONTINUE
+      RETURN
+      END
+      SUBROUTINE APC(ZPV,LE,AM0M,W,SQK,GAPR,GAPP)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION ZPV(LM,3,2,2),AM0M(NLEMT,NLEMT),W(NLEMT,4)
+      DIMENSION ZPV(8,3,2,2),AM0M(160,160),W(160,4),
+     1GAPR(3,2),GAPP(3,2)
+      COMPLEX*16 AM0M,W,GAPP
+CCC   DIMENSION AC(NLEMT,2)
+      DIMENSION AC(160,2),GAP(3,2)
+      COMPLEX*16 AC,GAP
+      COMPLEX*16 CC0,UIM,UIMMP,
+     1SUMM,SUME,SUEM,SUEE,SUMMP,SUMEP,SUEMP,SUEEP
+      CC0=(0.0D0,0.0D0)
+      UIM=(0.0D0,1.0D0)
+      COF=1.0D0/SQK
+      CIMU=COF/DSQRT(2.0D0)
+      NLEM=LE*(LE+2)
+      NLEMT=NLEM+NLEM
+      DO 45 J=1,NLEMT
+      AC(J,1)=CC0
+      AC(J,2)=CC0
+      DO 45 I=1,NLEMT
+      AC(J,1)=AC(J,1)+AM0M(J,I)*W(I,1)
+      AC(J,2)=AC(J,2)+AM0M(J,I)*W(I,2)
+   45 CONTINUE
+      DO 90 IMU=1,3
+      MU=IMU-2
+      GAP(IMU,1)=CC0
+      GAP(IMU,2)=CC0
+      GAPP(IMU,1)=CC0
+      GAPP(IMU,2)=CC0
+      DO 80 L=1,LE
+      LPO=L+1
+      LTPO=LPO+L
+      IMM=L*LPO
+      DO 80 ILMP=1,3
+      IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
+     1GO TO 80
+      LMPML=ILMP-2
+      LMP=L+LMPML
+      UIMMP=-LMPML*UIM
+      IMPMMMP=LMP*(LMP+1)
+      DO 70 IM=1,LTPO
+      M=IM-LPO
+      MMP=M-MU
+      IF(IABS(MMP).GT.LMP)GO TO 70
+      I=IMM+M
+      IE=I+NLEM
+      IMP=IMPMMMP+MMP
+      IMPE=IMP+NLEM
+      CGC=CG1(LMPML,MU,L,M)
+      JPO=2
+      DO 60 IPO=1,2
+      IF(IPO.EQ.2)JPO=1
+      SUMM=DCONJG(AC(I,IPO))*AC(IMP,IPO)
+      SUME=DCONJG(AC(I,IPO))*AC(IMPE,IPO)
+      SUEM=DCONJG(AC(IE,IPO))*AC(IMP,IPO)
+      SUEE=DCONJG(AC(IE,IPO))*AC(IMPE,IPO)
+      SUMMP=DCONJG(AC(I,JPO))*AC(IMP,IPO)
+      SUMEP=DCONJG(AC(I,JPO))*AC(IMPE,IPO)
+      SUEMP=DCONJG(AC(IE,JPO))*AC(IMP,IPO)
+      SUEEP=DCONJG(AC(IE,JPO))*AC(IMPE,IPO)
+      IF(LMPML.EQ.0)GO TO 55
+      SUMM=SUMM*UIMMP
+      SUME=SUME*UIMMP
+      SUEM=SUEM*UIMMP
+      SUEE=SUEE*UIMMP
+      SUMMP=SUMMP*UIMMP
+      SUMEP=SUMEP*UIMMP
+      SUEMP=SUEMP*UIMMP
+      SUEEP=SUEEP*UIMMP
+   55 GAP(IMU,IPO)=GAP(IMU,IPO)+CGC*
+     1(SUMM*ZPV(L,ILMP,1,1)+SUME*ZPV(L,ILMP,1,2)+
+     2SUEM*ZPV(L,ILMP,2,1)+SUEE*ZPV(L,ILMP,2,2))
+      GAPP(IMU,IPO)=GAPP(IMU,IPO)+CGC*
+     1(SUMMP*ZPV(L,ILMP,1,1)+SUMEP*ZPV(L,ILMP,1,2)+
+     2SUEMP*ZPV(L,ILMP,2,1)+SUEEP*ZPV(L,ILMP,2,2))
+   60 CONTINUE
+   70 CONTINUE
+   80 CONTINUE
+   90 CONTINUE
+      DO 95 IPO=1,2
+      SUME=GAP(1,IPO)*CIMU
+      SUEE=GAP(2,IPO)*COF
+      SUEM=GAP(3,IPO)*CIMU
+      GAPR(1,IPO)=DREAL(SUME-SUEM)
+      GAPR(2,IPO)=DREAL((SUME+SUEM)*UIM)
+      GAPR(3,IPO)=DREAL(SUEE)
+      SUMEP=GAPP(1,IPO)*CIMU
+      SUEEP=GAPP(2,IPO)*COF
+      SUEMP=GAPP(3,IPO)*CIMU
+      GAPP(1,IPO)=SUMEP-SUEMP
+      GAPP(2,IPO)=(SUMEP+SUEMP)*UIM
+      GAPP(3,IPO)=SUEEP
+   95 CONTINUE
+      RETURN
+      END
+      SUBROUTINE APCRA(ZPV,LE,AM0M,INPOL,SQK,GAPRM,GAPPM)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION ZPV(LM,3,2,2),AM0M(NLEMT,NLEMT)
+      DIMENSION ZPV(8,3,2,2),AM0M(160,160),GAPRM(3,2),GAPPM(3,2)
+      COMPLEX*16 AM0M,GAPPM
+CCC   DIMENSION SVW(LE,3,2,2),SVS(LE,3,2,2)
+      DIMENSION SVW(8,3,2,2),SVS(8,3,2,2)
+      COMPLEX*16 SVS
+      COMPLEX*16 CC0,UIM,UIMTL,UIMTLS,FC,
+     1CA11,CA12,CA21,CA22,A11,A12,A21,A22,SUM1,SUM2
+      CC0=(0.0D0,0.0D0)
+      UIM=(0.0D0,1.0D0)
+      NLEM=LE*(LE+2)
+      DO 28 L=1,LE
+      LPO=L+1
+      LTPO=LPO+L
+      FL=LTPO
+      FL=DSQRT(FL)
+      DO 26 ILMP=1,3
+      IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
+     1GO TO 26
+      LMPML=ILMP-2
+      LMP=L+LMPML
+      FLMP=LMP+LMP+1
+      FLMP=DSQRT(FLMP)
+      FLLMP=FLMP/FL
+      CGMMO=FLLMP*CG1(LMPML,0,L,1)
+      CGMPO=FLLMP*CG1(LMPML,0,L,-1)
+      IF(INPOL.NE.0)GO TO 22
+      CGS=CGMPO+CGMMO
+      CGD=CGMPO-CGMMO
+      SVW(L,ILMP,1,1)=CGS
+      SVW(L,ILMP,1,2)=CGD
+      SVW(L,ILMP,2,1)=CGD
+      SVW(L,ILMP,2,2)=CGS
+      GO TO 26
+   22 SVW(L,ILMP,1,1)=CGMPO
+      SVW(L,ILMP,2,1)=CGMPO
+      SVW(L,ILMP,1,2)=-CGMMO
+      SVW(L,ILMP,2,2)=CGMMO
+   26 CONTINUE
+   28 CONTINUE
+      DO 30 L=1,LE
+      DO 30 ILMP=1,3
+      DO 30 IPA=1,2
+      DO 30 IPAMP=1,2
+      SVS(L,ILMP,IPA,IPAMP)=CC0
+   30 CONTINUE
+      DO 58 L=1,LE
+      LPO=L+1
+      LTPO=L+LPO
+      IMM=L*LPO
+      DO 56 ILMP=1,3
+      IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
+     1GO TO 56
+      LMPML=ILMP-2
+      LMP=L+LMPML
+      IMPMM=LMP*(LMP+1)
+      UIMTL=LMPML*UIM
+      IF(LMPML.EQ.0)UIMTL=(1.0D0,0.0D0)
+      DO 54 IM=1,LTPO
+      M=IM-LPO
+      I=IMM+M
+      IE=I+NLEM
+      DO 52 IMU=1,3
+      MU=IMU-2
+      MMP=M-MU
+      IF(IABS(MMP).GT.LMP)GO TO 52
+      IMP=IMPMM+MMP
+      IMPE=IMP+NLEM
+      CGC=CG1(LMPML,-MU,L,-M)
+      DO 48 LS=1,LE
+      LSPO=LS+1
+      LSTPO=LS+LSPO
+      ISMM=LS*LSPO
+      DO 46 ILSMP=1,3
+      IF(((LS.EQ.1).AND.(ILSMP.EQ.1)).OR.((LS.EQ.LE).AND.(ILSMP.EQ.3)))
+     1GO TO 46
+      LSMPML=ILSMP-2
+      LSMP=LS+LSMPML
+      ISMPMM=LSMP*(LSMP+1)
+      UIMTLS=-LSMPML*UIM
+      IF(LSMPML.EQ.0)UIMTLS=(1.0D0,0.0D0)
+      DO 44 IMS=1,LSTPO
+      MS=IMS-LSPO
+      MSMP=MS-MU
+      IF(IABS(MSMP).GT.LSMP)GO TO 44
+      IS=ISMM+MS
+      ISE=IS+NLEM
+      ISMP=ISMPMM+MSMP
+      ISMPE=ISMP+NLEM
+      CGCS=CG1(LSMPML,MU,LS,MS)
+      FC=(UIMTL*UIMTLS)*(CGC*CGCS)
+      CA11=DCONJG(AM0M(IS,I))
+      CA12=DCONJG(AM0M(IS,IE))
+      CA21=DCONJG(AM0M(ISE,I))
+      CA22=DCONJG(AM0M(ISE,IE))
+      A11=AM0M(ISMP,IMP)
+      A12=AM0M(ISMP,IMPE)
+      A21=AM0M(ISMPE,IMP)
+      A22=AM0M(ISMPE,IMPE)
+      Z11=ZPV(LS,ILSMP,1,1)
+      Z12=ZPV(LS,ILSMP,1,2)
+      Z21=ZPV(LS,ILSMP,2,1)
+      Z22=ZPV(LS,ILSMP,2,2)
+      SVS(L,ILMP,1,1)=SVS(L,ILMP,1,1)+
+     1(CA11*A11*Z11+CA11*A21*Z12+CA21*A11*Z21+CA21*A21*Z22)*FC
+      SVS(L,ILMP,1,2)=SVS(L,ILMP,1,2)+
+     1(CA11*A12*Z11+CA11*A22*Z12+CA21*A12*Z21+CA21*A22*Z22)*FC
+      SVS(L,ILMP,2,1)=SVS(L,ILMP,2,1)+
+     1(CA12*A11*Z11+CA12*A21*Z12+CA22*A11*Z21+CA22*A21*Z22)*FC
+      SVS(L,ILMP,2,2)=SVS(L,ILMP,2,2)+
+     1(CA12*A12*Z11+CA12*A22*Z12+CA22*A12*Z21+CA22*A22*Z22)*FC
+   44 CONTINUE
+   46 CONTINUE
+   48 CONTINUE
+   52 CONTINUE
+   54 CONTINUE
+   56 CONTINUE
+   58 CONTINUE
+      SUM1=CC0
+      SUM2=CC0
+      DO 68 L=1,LE
+      LPO=L+1
+      LTPO=L+LPO
+      IMM=L*LPO
+      DO 66 ILMP=1,3
+      IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
+     1GO TO 66
+      IF(INPOL.NE.0)GO TO 62
+      SUM1=SUM1+
+     1SVW(L,ILMP,1,1)*SVS(L,ILMP,1,1)+SVW(L,ILMP,2,1)*SVS(L,ILMP,1,2)+
+     2SVW(L,ILMP,2,1)*SVS(L,ILMP,2,1)+SVW(L,ILMP,1,1)*SVS(L,ILMP,2,2)
+      SUM2=SUM2+
+     1SVW(L,ILMP,1,2)*SVS(L,ILMP,1,1)+SVW(L,ILMP,2,2)*SVS(L,ILMP,1,2)+
+     2SVW(L,ILMP,2,2)*SVS(L,ILMP,2,1)+SVW(L,ILMP,1,2)*SVS(L,ILMP,2,2)
+      GO TO 66
+   62 SUM1=SUM1+
+     1SVW(L,ILMP,2,1)*SVS(L,ILMP,1,1)+SVW(L,ILMP,1,1)*SVS(L,ILMP,1,2)+
+     2SVW(L,ILMP,1,1)*SVS(L,ILMP,2,1)+SVW(L,ILMP,2,1)*SVS(L,ILMP,2,2)
+      SUM2=SUM2+
+     1SVW(L,ILMP,2,2)*SVS(L,ILMP,1,1)+SVW(L,ILMP,1,2)*SVS(L,ILMP,1,2)+
+     2SVW(L,ILMP,1,2)*SVS(L,ILMP,2,1)+SVW(L,ILMP,2,2)*SVS(L,ILMP,2,2)
+   66 CONTINUE
+   68 CONTINUE
+      C0=0.0D0
+      PIGH=DACOS(C0)
+      COFS=PIGH*2.0D0/SQK
+      GAPRM(1,1)=C0
+      GAPRM(1,2)=C0
+      GAPRM(2,1)=C0
+      GAPRM(2,2)=C0
+      GAPPM(1,1)=CC0
+      GAPPM(1,2)=CC0
+      GAPPM(2,1)=CC0
+      GAPPM(2,2)=CC0
+      IF(INPOL.NE.0)GO TO 72
+      SUM1=SUM1*COFS
+      SUM2=SUM2*COFS
+      GAPRM(3,1)=DREAL(SUM1)
+      GAPRM(3,2)=DREAL(SUM1)
+      GAPPM(3,1)=SUM2*UIM
+      GAPPM(3,2)=-GAPPM(3,1)
+      RETURN
+   72 COFS=COFS*2.0D0
+      GAPRM(3,1)=DREAL(SUM1)*COFS
+      GAPRM(3,2)=DREAL(SUM2)*COFS
+      GAPPM(3,1)=CC0
+      GAPPM(3,2)=CC0
+      RETURN
+      END
+      SUBROUTINE THDPS(LM,ZPV)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION ZPV(LM,3,2,2)
+      DIMENSION ZPV(8,3,2,2)
+      C0=0.0D0
+      DO 10 L=1,LM
+      DO 10 ILMP=1,3
+      ZPV(L,ILMP,1,1)=C0
+      ZPV(L,ILMP,1,2)=C0
+      ZPV(L,ILMP,2,1)=C0
+      ZPV(L,ILMP,2,2)=C0
+   10 CONTINUE
+      DO 15 L=1,LM
+      XD=L*(L+1)
+      ZP=-1.0D0/DSQRT(XD)
+      ZPV(L,2,1,2)=ZP
+      ZPV(L,2,2,1)=ZP
+   15 CONTINUE
+      IF(LM.EQ.1)GO TO 30
+      DO 20 L=2,LM
+      XN=(L-1)*(L+1)
+      XD=L*(L+L+1)
+      ZP=DSQRT(XN/XD)
+      ZPV(L,1,1,1)=ZP
+      ZPV(L,1,2,2)=ZP
+   20 CONTINUE
+      LMMO=LM-1
+      DO 25 L=1,LMMO
+      XN=L*(L+2)
+      XD=(L+1)*(L+L+1)
+      ZP=-DSQRT(XN/XD)
+      ZPV(L,3,1,1)=ZP
+      ZPV(L,3,2,2)=ZP
+   25 CONTINUE
+   30 CONTINUE
+      RETURN
+      END
+      REAL*8 FUNCTION CG1(LMPML,MU,L,M)
+CCC   CG1(LMPML,MU,L,M)=CLGO(1,LMP,L;MU,M-MU,M)
+      REAL*8 XD,XN
+      IF(LMPML)30,5,60
+    5 IF((M.NE.0).OR.(MU.NE.0))GO TO 10
+      CG1=0.0D0
+      RETURN
+   10 IF(MU.EQ.0)GO TO 20
+      XD=(L+1)*L*2
+      IF(MU.GT.0)GO TO 15
+      XN=(L-M)*(L+M+1)
+      CG1=-DSQRT(XN/XD)
+      RETURN
+   15 XN=(L+M)*(L-M+1)
+      CG1=DSQRT(XN/XD)
+      RETURN
+   20 XD=(L+1)*L
+      XN=-M
+      CG1=XN/DSQRT(XD)
+      RETURN
+   30 XD=(L*2-1)*L*2 
+      IF(MU)35,40,45
+   35 XN=(L-1-M)*(L-M)
+      CG1=DSQRT(XN/XD)
+      RETURN
+   40 XN=(L-M)*(L+M)*2
+      CG1=DSQRT(XN/XD)
+      RETURN
+   45 XN=(L-1+M)*(L+M)
+      CG1=DSQRT(XN/XD)
+      RETURN
+   60 XD=(L*2+3)*(L+1)*2 
+      IF(MU)65,70,75
+   65 XN=(L+1+M)*(L+2+M)
+      CG1=DSQRT(XN/XD)
+      RETURN
+   70 XN=(L+1-M)*(L+1+M)*2
+      CG1=-DSQRT(XN/XD)
+      RETURN
+   75 XN=(L+1-M)*(L+2-M)
+      CG1=DSQRT(XN/XD)
+      RETURN
+      END
+      SUBROUTINE UPVMP(THD,PHD,ICSPNV,COST,SINT,COSP,SINP,U,UP,UN)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION U(3),UP(3),UN(3)
+      C0=0.0D0
+      PIGH=DACOS(C0)
+      RDR=PIGH/9.0D+1
+      TH=THD*RDR
+      PH=PHD*RDR
+      COST=DCOS(TH)
+      SINT=DSIN(TH)
+      COSP=DCOS(PH)
+      SINP=DSIN(PH)
+      U(1)=COSP*SINT
+      U(2)=SINP*SINT
+      U(3)=COST
+      UP(1)=COSP*COST
+      UP(2)=SINP*COST
+      UP(3)=-SINT
+      UN(1)=-SINP
+      UN(2)=COSP
+      UN(3)=C0
+      IF(ICSPNV.EQ.0)RETURN
+      UP(1)=-UP(1)
+      UP(2)=-UP(2)
+      UP(3)=-UP(3)
+      UN(1)=-UN(1)
+      UN(2)=-UN(2)
+      RETURN
+      END
+      SUBROUTINE WMAMP
+     1(IIS,COST,SINT,COSP,SINP,INPOL,LM,IDOT,NSPH,ARG,U,UP,UN)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION YLM(NLMM+2),ARG(NSPEF)
+CCC   NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
+      DIMENSION YLM(82),U(3),UP(3),UN(3),ARG(1)
+      COMPLEX*16 YLM
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      NLMP=LM*(LM+2)+2
+      YLM(NLMP)=(0.0D0,0.0D0)
+      IF(IDOT.EQ.0)GO TO 65
+      IF(IDOT.EQ.1)GO TO 45
+      DO 40 N=1,NSPH
+   40 ARG(N)=U(1)*RXX(N)+U(2)*RYY(N)+U(3)*RZZ(N)
+      GO TO 55
+   45 DO 50 N=1,NSPH
+   50 ARG(N)=COST*RZZ(N)
+   55 IF(IIS.NE.2)GO TO 65
+      DO 60 N=1,NSPH
+   60 ARG(N)=-ARG(N)
+   65 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+      CALL PWMA(UP,UN,YLM,INPOL,LM,IIS)
+      RETURN
+      END
+      SUBROUTINE UPVSP(U,UPMP,UNMP,US,UPSMP,UNSMP,UP,UN,UPS,UNS,DUK,
+     1ISQ,IBF,SCAND,CFMP,SFMP,CFSP,SFSP)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION U(3),UPMP(3),UNMP(3),US(3),UPSMP(3),UNSMP(3),
+     1UP(3),UN(3),UPS(3),UNS(3),DUK(3)
+      C0=0.0D0
+      RDR=DACOS(C0)/9.0D+1
+      SML=1.0D-6
+      ISQ=0
+      SCAND=U(1)*US(1)+U(2)*US(2)+U(3)*US(3)
+      IF(DABS(SCAND-1.0D0).LT.SML)GO TO 10
+      IF(DABS(SCAND+1.0D0).LT.SML)GO TO 15
+      SCAND=DACOS(SCAND)/RDR
+      DUK(1)=U(1)-US(1)
+      DUK(2)=U(2)-US(2)
+      DUK(3)=U(3)-US(3)
+      IBF=0
+      GO TO 25
+   10 SCAND=C0
+      DUK(1)=C0
+      DUK(2)=C0
+      DUK(3)=C0
+      IBF=-1
+      ISQ=-1
+      UPS(1)=UPSMP(1)
+      UPS(2)=UPSMP(2)
+      UPS(3)=UPSMP(3)
+      UNS(1)=UNSMP(1)
+      UNS(2)=UNSMP(2)
+      UNS(3)=UNSMP(3)
+      GO TO 20
+   15 SCAND=180.0D0
+      DUK(1)=U(1)*2.0D0
+      DUK(2)=U(2)*2.0D0
+      DUK(3)=U(3)*2.0D0
+      IBF=1
+      UPS(1)=-UPSMP(1)
+      UPS(2)=-UPSMP(2)
+      UPS(3)=-UPSMP(3)
+      UNS(1)=-UNSMP(1)
+      UNS(2)=-UNSMP(2)
+      UNS(3)=-UNSMP(3)
+   20 UP(1)=UPMP(1)
+      UP(2)=UPMP(2)
+      UP(3)=UPMP(3)
+      UN(1)=UNMP(1)
+      UN(2)=UNMP(2)
+      UN(3)=UNMP(3)
+      GO TO 85
+   25 CALL ORUNVE(U,US,UN,-1,SML)
+      UNS(1)=UN(1)
+      UNS(2)=UN(2)
+      UNS(3)=UN(3)
+      CALL ORUNVE(UN,U,UP,1,SML)
+      CALL ORUNVE(UNS,US,UPS,1,SML)
+   85 CFMP=UPMP(1)*UP(1)+UPMP(2)*UP(2)+UPMP(3)*UP(3)
+      SFMP=UNMP(1)*UP(1)+UNMP(2)*UP(2)+UNMP(3)*UP(3)
+      CFSP=UPS(1)*UPSMP(1)+UPS(2)*UPSMP(2)+UPS(3)*UPSMP(3)
+      SFSP=UNS(1)*UPSMP(1)+UNS(2)*UPSMP(2)+UNS(3)*UPSMP(3)
+      RETURN
+      END      
+      SUBROUTINE WMASP(COST,SINT,COSP,SINP,COSTS,SINTS,COSPS,SINPS,
+     1U,UP,UN,US,UPS,UNS,ISQ,IBF,INPOL,LM,IDOT,NSPH,ARGI,ARGS)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION YLM(NLMM+2),ARGI(NSPEF),ARGS(NSPEF)
+CCC   NSPEF=1 CALLING WITH IDOT=0, NSPEF=NSPH OTHERWISE
+      DIMENSION YLM(82),ARGI(1),ARGS(1),
+     1U(3),UP(3),UN(3),US(3),UPS(3),UNS(3)
+      COMPLEX*16 YLM
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      NLMP=LM*(LM+2)+2
+      YLM(NLMP)=(0.0D0,0.0D0)
+      IF(IDOT.EQ.0)GO TO 75 
+      IF(IDOT.EQ.1)GO TO 50 
+      DO 40 N=1,NSPH
+      ARGI(N)=U(1)*RXX(N)+U(2)*RYY(N)+U(3)*RZZ(N)
+      IF(IBF.EQ.0)GO TO 35
+      ARGS(N)=ARGI(N)*IBF
+      GO TO 40
+   35 ARGS(N)=-(US(1)*RXX(N)+US(2)*RYY(N)+US(3)*RZZ(N))
+   40 CONTINUE
+      GO TO 75
+   50 DO 60 N=1,NSPH
+      ARGI(N)=COST*RZZ(N)
+      IF(IBF.EQ.0)GO TO 55
+      ARGS(N)=ARGI(N)*IBF
+      GO TO 60
+   55 ARGS(N)=-COSTS*RZZ(N)
+   60 CONTINUE
+   75 CALL SPHAR(COST,SINT,COSP,SINP,LM,YLM)
+      CALL PWMA(UP,UN,YLM,INPOL,LM,ISQ)
+      IF(IBF.LT.0)RETURN
+      CALL SPHAR(COSTS,SINTS,COSPS,SINPS,LM,YLM)
+      CALL PWMA(UPS,UNS,YLM,INPOL,LM,2)
+      RETURN
+      END
+      SUBROUTINE STR(RCF)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION RCF(NSPH,NSHL),YLM(MAX0(LITPOS,LMTPOS))
+      DIMENSION RCF(6,8),YLM(289)
+      COMPLEX*16 YLM
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C3/TFSAS,TSAS(2,2),GCS,SCS,ECS,ACS
+      COMPLEX*16 TFSAS,TSAS
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMMON/C6/RAC3J(17)
+      C0=0.0D0
+      PIGH=DACOS(C0)
+      PIG=PIGH*2.0D0
+      GCS=C0
+      DO 18 I=1,NSPH
+      IOGI=IOG(I)
+      IF(IOGI.LT.I)GO TO 18
+      GCSS=PIG*ROS(I)*ROS(I)
+      GCSV(I)=GCSS
+      NSH=NSHL(I)
+      DO 16 J=1,NSH
+   16 RC(I,J)=RCF(I,J)*ROS(I)
+   18 GCS=GCS+GCSV(IOGI)
+      LITPO=LI+LI+1
+      LITPOS=LITPO*LITPO
+      LMTPO=LI+LE+1
+      LMTPOS=LMTPO*LMTPO
+      NLIM=LI*(LI+2)
+      NLEM=LE*(LE+2)
+      LM=MAX0(LI,LE)
+      LMPO=LM+1
+      I=0
+      DO 28 L1PO=1,LMPO
+      L1=L1PO-1
+      DO 28 L2=1,LM
+      CALL R3J000(L1,L2)
+      IND3J(L1PO,L2)=I
+      LMNPO=IABS(L2-L1)+1
+      LMXPO=L2+L1+1
+      IL=0
+      DO 28 LPO=LMNPO,LMXPO,2
+      I=I+1
+      IL=IL+1
+   28 V3J0(I)=RAC3J(IL)
+      NSPHMO=NSPH-1
+      LIT=LI+LI
+      IVY=0
+      DO 40 NF=1,NSPHMO
+      NFPO=NF+1
+      DO 40 NS=NFPO,NSPH
+      RX=RXX(NF)-RXX(NS)
+      RY=RYY(NF)-RYY(NS)
+      RZ=RZZ(NF)-RZZ(NS)
+      CALL POLAR(RX,RY,RZ,RR,CRTH,SRTH,CRPH,SRPH)
+      CALL SPHAR(CRTH,SRTH,CRPH,SRPH,LIT,YLM)
+      DO 38 IV=1,LITPOS
+   38 VYHJ(IV+IVY)=DCONJG(YLM(IV))
+   40 IVY=IVY+LITPOS
+      LMT=LI+LE
+      IVY=0
+      DO 50 NF=1,NSPH
+      RX=RXX(NF)
+      RY=RYY(NF)
+      RZ=RZZ(NF)
+      IF(RX.EQ.C0.AND.RY.EQ.C0.AND.RZ.EQ.C0)GO TO 50
+      CALL POLAR(RX,RY,RZ,RR,CRTH,SRTH,CRPH,SRPH)
+      CALL SPHAR(CRTH,SRTH,CRPH,SRPH,LMT,YLM)
+      DO 48 IV=1,LMTPOS
+   48 VYJ0(IV+IVY)=DCONJG(YLM(IV))
+   50 IVY=IVY+LMTPOS
+      RETURN
+      END
+      SUBROUTINE HJV(EXRI,VK,JER,LCALC,ARG)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION RFJ(MAX0(LIT,LMT)),RFN(LITPO)
+      DIMENSION RFJ(17),RFN(17)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMPLEX*16 ARG,UIM
+      JER=0
+      C0=0.0D0
+      UIM=(0.0D0,1.0D0)
+      NSPHMO=NSPH-1
+      LIT=LI+LI
+      IVHB=0
+      DO 40 NF=1,NSPHMO
+      NFPO=NF+1
+      DO 40 NS=NFPO,NSPH
+      RX=RXX(NF)-RXX(NS)
+      RY=RYY(NF)-RYY(NS)
+      RZ=RZZ(NF)-RZZ(NS)
+      RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
+      RARG=RR*VK*EXRI
+      ARG=RARG
+      CALL RBF(LIT,RARG,LCALC,RFJ)
+      IF(LCALC.GE.LIT)GO TO 35
+      JER=1
+      RETURN
+   35 CALL RNF(LIT,RARG,LCALC,RFN)
+      IF(LCALC.GE.LIT)GO TO 37
+      JER=2
+      RETURN
+   37 DO 38 LPO=1,LITPO
+   38 VH(LPO+IVHB)=RFJ(LPO)+UIM*RFN(LPO)
+   40 IVHB=IVHB+LITPO
+      LMT=LI+LE
+      IVHB=0
+      DO 50 NF=1,NSPH
+      RX=RXX(NF)
+      RY=RYY(NF)
+      RZ=RZZ(NF)
+      IF(RX.EQ.C0.AND.RY.EQ.C0.AND.RZ.EQ.C0)GO TO 50
+      RR=DSQRT(RX*RX+RY*RY+RZ*RZ)
+      RARG=RR*VK*EXRI
+      CALL RBF(LMT,RARG,LCALC,RFJ)
+      IF(LCALC.GE.LMT)GO TO 45
+      JER=3
+      ARG=RARG
+      RETURN
+   45 DO 47 LPO=1,LMTPO
+   47 VJ0(LPO+IVHB)=RFJ(LPO)
+   50 IVHB=IVHB+LMTPO
+      RETURN
+      END
+      SUBROUTINE PWMA(UP,UN,YLM,INPOL,LW,ISQ)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION YLM(NLWM+2)
+      DIMENSION YLM(82),UP(3),UN(3)
+      COMPLEX*16 YLM
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMPLEX*16 CP1,CM1,C1,CP2,CM2,C2,UIM,CL
+      PI4=DACOS(0.0D0)*8.0D0
+      IS=ISQ
+      IF(ISQ.EQ.-1)IS=0
+      ISPO=IS+1
+      ISPT=IS+2
+      NLWM=LW*(LW+2)
+      NLWMT=NLWM+NLWM
+      SQRTWI=1.0D0/DSQRT(2.0D0)
+      UIM=(0.0D0,1.0D0)
+      CM1=.5D0*DCMPLX(UP(1),UP(2))
+      CP1=.5D0*DCMPLX(UP(1),-UP(2))
+      CZ1=UP(3)
+      CM2=.5D0*DCMPLX(UN(1),UN(2))
+      CP2=.5D0*DCMPLX(UN(1),-UN(2))
+      CZ2=UN(3)
+      DO 20 L=1,LW
+      LF=L+1
+      LFTL=LF*L
+      X=LFTL
+      CL=(PI4/DSQRT(X))*UIM**L
+      MV=L+LF
+      M=-LF
+      DO 20 MF=1,MV
+      M=M+1
+      K=LFTL+M
+      X=LFTL-M*(M+1)
+      CP=DSQRT(X)
+      X=LFTL-M*(M-1)
+      CM=DSQRT(X)
+      CZ=M
+      W(K,ISPO)=DCONJG(CP1*CP*YLM(K+2)+CM1*CM*YLM(K)+CZ1*CZ*YLM(K+1))*CL
+   20 W(K,ISPT)=DCONJG(CP2*CP*YLM(K+2)+CM2*CM*YLM(K)+CZ2*CZ*YLM(K+1))*CL
+      DO 30 K=1,NLWM
+      I=K+NLWM
+      W(I,ISPO)=UIM*W(K,ISPT)
+   30 W(I,ISPT)=-UIM*W(K,ISPO)
+      IF(INPOL.EQ.0)GO TO 42
+      DO 40 K=1,NLWM
+      I=K+NLWM
+      C1=(W(K,ISPO)+UIM*W(K,ISPT))*SQRTWI
+      C2=(W(K,ISPO)-UIM*W(K,ISPT))*SQRTWI
+      W(K,ISPO)=C2
+      W(I,ISPO)=-C2
+      W(K,ISPT)=C1
+   40 W(I,ISPT)=C1
+   42 IF(ISQ.EQ.0)RETURN
+      DO 50 I=1,2
+      IPT=I+2
+      IPIS=I+IS
+      DO 50 K=1,NLWMT
+   50 W(K,IPT)=DCONJG(W(K,IPIS))
+      RETURN
+      END
+      SUBROUTINE ORUNVE(U1,U2,U3,IORTH,TORTH)
+      IMPLICIT REAL*8(A-H,O-Z)
+      DIMENSION U1(3),U2(3),U3(3)
+      IF(IORTH.GT.0)GO TO 10
+      CP=U1(1)*U2(1)+U1(2)*U2(2)+U1(3)*U2(3)
+      IF((IORTH.LT.0).AND.(DABS(CP).LT.TORTH))GO TO 10
+      FN=1.0D0/DSQRT(1.0D0-CP*CP)
+      U3(1)=(U1(2)*U2(3)-U1(3)*U2(2))*FN
+      U3(2)=(U1(3)*U2(1)-U1(1)*U2(3))*FN
+      U3(3)=(U1(1)*U2(2)-U1(2)*U2(1))*FN
+      RETURN
+   10 U3(1)=U1(2)*U2(3)-U1(3)*U2(2)
+      U3(2)=U1(3)*U2(1)-U1(1)*U2(3)
+      U3(3)=U1(1)*U2(2)-U1(2)*U2(1)
+      RETURN
+      END      
+      SUBROUTINE CMS(AM)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+CCC   DIMENSION AM(MXNDM,MXNDM)
+CCC   NDIT.LE.MXNDM
+      DIMENSION AM(960,960)
+      COMPLEX*16 AM
+      COMPLEX*16 DM,DE,CGH,CGK,GHIT,CC0
+      CC0=(0.0D0,0.0D0)
+      NDI=NSPH*NLIM
+      NBL=0
+      NSPHMO=NSPH-1
+      DO 26 N1=1,NSPHMO
+      IN1=(N1-1)*NLIM
+      N1PO=N1+1
+      DO 26 N2=N1PO,NSPH
+      IN2=(N2-1)*NLIM
+      NBL=NBL+1
+      DO 24 L1=1,LI
+      L1PO=L1+1
+      IL1=L1PO*L1
+      L1TPO=L1PO+L1
+      DO 24 IM1=1,L1TPO
+      M1=IM1-L1PO
+      ILM1=IL1+M1
+      ILM1E=ILM1+NDI
+      I1=IN1+ILM1
+      I1E=IN1+ILM1E
+      J1=IN2+ILM1
+      J1E=IN2+ILM1E
+      DO 24 L2=1,LI
+      L2PO=L2+1
+      IL2=L2PO*L2
+      L2TPO=L2PO+L2
+      ISH=(-1)**(L2+L1)
+      ISK=-ISH
+      DO 24 IM2=1,L2TPO
+      M2=IM2-L2PO
+      ILM2=IL2+M2
+      ILM2E=ILM2+NDI
+      I2=IN2+ILM2
+      I2E=IN2+ILM2E
+      J2=IN1+ILM2
+      J2E=IN1+ILM2E
+      CGH=GHIT(0,0,NBL,L1,M1,L2,M2)
+      CGK=GHIT(0,1,NBL,L1,M1,L2,M2)
+      AM(I1,I2)=CGH
+      AM(I1,I2E)=CGK
+      AM(I1E,I2)=CGK
+      AM(I1E,I2E)=CGH
+      AM(J1,J2)=CGH*ISH
+      AM(J1,J2E)=CGK*ISK
+      AM(J1E,J2)=CGK*ISK
+      AM(J1E,J2E)=CGH*ISH
+   24 CONTINUE
+   26 CONTINUE
+      DO 30 N1=1,NSPH
+      IN1=(N1-1)*NLIM
+      DO 30 L1=1,LI
+      DM=RMI(L1,N1)
+      DE=REI(L1,N1)
+      L1PO=L1+1
+      IL1=L1PO*L1
+      L1TPO=L1PO+L1
+      DO 30 IM1=1,L1TPO
+      M1=IM1-L1PO
+      ILM1=IL1+M1
+      I1=IN1+ILM1
+      I1E=I1+NDI
+      DO 28 ILM2=1,NLIM
+      I2=IN1+ILM2
+      I2E=I2+NDI
+      AM(I1,I2)=CC0
+      AM(I1,I2E)=CC0
+      AM(I1E,I2)=CC0
+   28 AM(I1E,I2E)=CC0
+      AM(I1,I1)=DM
+   30 AM(I1E,I1E)=DE
+      RETURN
+      END
+      SUBROUTINE ZTM(AM)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMMON/C9/GIS(480,80),GLS(480,80),SAM(960,160)
+      COMPLEX*16 GIS,GLS,SAM
+CCC   DIMENSION AM(MXNDM,MXNDM)
+CCC   NDIT.LE.MXNDM
+      DIMENSION AM(960,960)
+      COMPLEX*16 AM
+      COMPLEX*16 GIE,GLE,GHIT,A1,A2,A3,A4,
+     1SUM1,SUM2,SUM3,SUM4,CC0
+      CC0=(0.0D0,0.0D0)
+      NDI=NSPH*NLIM
+      I2=0
+      DO 15 N2=1,NSPH
+      DO 15 L2=1,LI
+      L2TPO=L2+L2+1
+      M2=-L2-1
+      DO 15 IM2=1,L2TPO
+      M2=M2+1
+      I2=I2+1
+      I3=0
+      DO 15 L3=1,LE
+      L3TPO=L3+L3+1
+      M3=-L3-1
+      DO 15 IM3=1,L3TPO
+      M3=M3+1
+      I3=I3+1
+      GIS(I2,I3)=GHIT(2,0,N2,L2,M2,L3,M3)
+   15 GLS(I2,I3)=GHIT(2,1,N2,L2,M2,L3,M3)
+      DO 25 I1=1,NDI
+      I1E=I1+NDI
+      DO 25 I3=1,NLEM
+      I3E=I3+NLEM
+      SUM1=CC0
+      SUM2=CC0
+      SUM3=CC0
+      SUM4=CC0
+      DO 20 I2=1,NDI
+      I2E=I2+NDI
+      GIE=GIS(I2,I3)
+      GLE=GLS(I2,I3)
+      A1=AM(I1,I2)
+      A2=AM(I1,I2E)
+      A3=AM(I1E,I2)
+      A4=AM(I1E,I2E)
+      SUM1=SUM1+A1*GIE+A2*GLE
+      SUM2=SUM2+A1*GLE+A2*GIE
+      SUM3=SUM3+A3*GIE+A4*GLE
+   20 SUM4=SUM4+A3*GLE+A4*GIE
+      SAM(I1,I3)=SUM1
+      SAM(I1,I3E)=SUM2
+      SAM(I1E,I3)=SUM3
+   25 SAM(I1E,I3E)=SUM4
+      DO 35 I1=1,NDI
+      DO 35 I0=1,NLEM
+      GIS(I1,I0)=DCONJG(GIS(I1,I0))
+   35 GLS(I1,I0)=DCONJG(GLS(I1,I0))
+      NLEMT=NLEM+NLEM
+      DO 45 I0=1,NLEM
+      I0E=I0+NLEM
+      DO 45 I3=1,NLEMT
+      SUM1=CC0
+      SUM2=CC0
+      DO 40 I1=1,NDI
+      I1E=I1+NDI
+      A1=SAM(I1,I3)
+      A2=SAM(I1E,I3)
+      GIE=GIS(I1,I0)
+      GLE=GLS(I1,I0)
+      SUM1=SUM1+A1*GIE+A2*GLE
+   40 SUM2=SUM2+A1*GLE+A2*GIE
+      AM0M(I0,I3)=-SUM1
+   45 AM0M(I0E,I3)=-SUM2
+      RETURN
+      END
+      COMPLEX*16 FUNCTION GHIT(IHI,IPAMO,NBL,L1,M1,L2,M2)
+CCC   NBL INDIVIDUATES TRANSFER VECTOR GOING FROM N2 TO N1;
+CCC   IHI=0 FOR HANKEL, IHI=1 FOR BESSEL, IHI=2 FOR BESSEL FROM ORIGIN;
+CCC   DEPENDING ON IHI, IPAMO=0 GIVES H OR I, IPAMO=1 GIVES K OR L
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C4/LITPO,LITPOS,LMTPO,LMTPOS,LI,NLIM,LE,NLEM,NSPH
+      COMMON/C6/RAC3J(17)
+      COMPLEX*16 CSUM,CC0,UIM,CFUN
+      CC0=(0.0D0,0.0D0)
+      GHIT=CC0
+      IF(IHI.NE.2)GO TO 10
+      C0=0.0D0
+      IF(RXX(NBL).NE.C0.OR.RYY(NBL).NE.C0.OR.RZZ(NBL).NE.C0)GO TO 10
+      IF(IPAMO.NE.0)RETURN
+      IF(L1.EQ.L2.AND.M1.EQ.M2)GHIT=(1.0D0,0.0D0)
+      RETURN
+   10 L1MP=L1-IPAMO
+      L1PO=L1+1
+      M1MM2=M1-M2
+      M1MM2M=IABS(M1MM2)+1
+      LMINPO=IABS(L2-L1MP)+1
+      LMAXPO=L2+L1MP+1
+      I3J0IN=IND3J(L1MP+1,L2)
+      ILIN=-1
+      IF((M1MM2M.GT.LMINPO).AND.(MOD(M1MM2M-LMINPO,2).NE.0))ILIN=0
+      ISN=1
+      IF(MOD(M1,2).NE.0)ISN=-ISN
+      IF(MOD(LMINPO,2).NE.0)GO TO 12
+      ISN=-ISN
+      IF(L2.GT.L1MP)ISN=-ISN
+   12 NBLMO=NBL-1
+      IF(IHI.EQ.2)GO TO 50
+      NBHJ=NBLMO*LITPO
+      NBY=NBLMO*LITPOS
+      IF(IHI.EQ.1)GO TO 30
+      DO 24 JM=1,3
+      CSUM=CC0
+      MU=JM-2
+      MUPM1=MU+M1
+      MUPM2=MU+M2
+      IF(MUPM1.LT.-L1MP.OR.MUPM1.GT.L1MP.OR.MUPM2.LT.-L2.OR.MUPM2.GT.L2)
+     1GO TO 24
+      JSN=-ISN
+      IF(MU.EQ.0)JSN=ISN
+      CR=CGEV(IPAMO,MU,L1,M1)*CGEV(0,MU,L2,M2)
+      I3J0=I3J0IN
+      IF((MUPM1.NE.0).OR.(MUPM2.NE.0))GO TO 16
+      DO 14 LT=LMINPO,LMAXPO,2
+      I3J0=I3J0+1
+      L3=LT-1
+      NY=L3*L3+LT
+      AORS=L3+LT
+      F3J=(V3J0(I3J0)*V3J0(I3J0)*DSQRT(AORS))*JSN
+      CFUN=(VH(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+      CSUM=CSUM+CFUN
+   14 JSN=-JSN
+      GO TO 22
+   16 CALL R3JJR(L1MP,L2,-MUPM1,MUPM2)
+      IL=ILIN
+      DO 20 LT=LMINPO,LMAXPO,2
+      I3J0=I3J0+1
+      IF(M1MM2M.GT.LT)GO TO 20
+      IL=IL+2
+      L3=LT-1
+      NY=L3*L3+LT+M1MM2
+      AORS=L3+LT
+      F3J=(RAC3J(IL)*V3J0(I3J0)*DSQRT(AORS))*JSN
+      CFUN=(VH(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+      CSUM=CSUM+CFUN
+   20 JSN=-JSN
+   22 CSUM=CSUM*CR
+      GHIT=GHIT+CSUM
+   24 CONTINUE
+      GO TO 70
+   30 DO 44 JM=1,3
+      CSUM=CC0
+      MU=JM-2
+      MUPM1=MU+M1
+      MUPM2=MU+M2
+      IF(MUPM1.LT.-L1MP.OR.MUPM1.GT.L1MP.OR.MUPM2.LT.-L2.OR.MUPM2.GT.L2)
+     1GO TO 44
+      JSN=-ISN
+      IF(MU.EQ.0)JSN=ISN
+      CR=CGEV(IPAMO,MU,L1,M1)*CGEV(0,MU,L2,M2)
+      I3J0=I3J0IN
+      IF((MUPM1.NE.0).OR.(MUPM2.NE.0))GO TO 36
+      DO 34 LT=LMINPO,LMAXPO,2
+      I3J0=I3J0+1
+      L3=LT-1
+      NY=L3*L3+LT
+      AORS=L3+LT
+      F3J=(V3J0(I3J0)*V3J0(I3J0)*DSQRT(AORS))*JSN
+      CFUN=(VJ(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+      CSUM=CSUM+CFUN
+   34 JSN=-JSN
+      GO TO 42
+   36 CALL R3JJR(L1MP,L2,-MUPM1,MUPM2)
+      IL=ILIN
+      DO 40 LT=LMINPO,LMAXPO,2
+      I3J0=I3J0+1
+      IF(M1MM2M.GT.LT)GO TO 40
+      IL=IL+2
+      L3=LT-1
+      NY=L3*L3+LT+M1MM2
+      AORS=L3+LT
+      F3J=(RAC3J(IL)*V3J0(I3J0)*DSQRT(AORS))*JSN
+      CFUN=(VJ(NBHJ+LT)*VYHJ(NBY+NY))*F3J
+      CSUM=CSUM+CFUN
+   40 JSN=-JSN
+   42 CSUM=CSUM*CR
+      GHIT=GHIT+CSUM
+   44 CONTINUE
+      GO TO 70
+   50 NBHJ=NBLMO*LMTPO
+      NBY=NBLMO*LMTPOS
+      DO 64 JM=1,3
+      CSUM=CC0
+      MU=JM-2
+      MUPM1=MU+M1
+      MUPM2=MU+M2
+      IF(MUPM1.LT.-L1MP.OR.MUPM1.GT.L1MP.OR.MUPM2.LT.-L2.OR.MUPM2.GT.L2)
+     1GO TO 64
+      JSN=-ISN
+      IF(MU.EQ.0)JSN=ISN
+      CR=CGEV(IPAMO,MU,L1,M1)*CGEV(0,MU,L2,M2)
+      I3J0=I3J0IN
+      IF((MUPM1.NE.0).OR.(MUPM2.NE.0))GO TO 56
+      DO 54 LT=LMINPO,LMAXPO,2
+      I3J0=I3J0+1
+      L3=LT-1
+      NY=L3*L3+LT
+      AORS=L3+LT
+      F3J=(V3J0(I3J0)*V3J0(I3J0)*DSQRT(AORS))*JSN
+      CFUN=(VJ0(NBHJ+LT)*VYJ0(NBY+NY))*F3J
+      CSUM=CSUM+CFUN
+   54 JSN=-JSN
+      GO TO 62
+   56 CALL R3JJR(L1MP,L2,-MUPM1,MUPM2)
+      IL=ILIN
+      DO 60 LT=LMINPO,LMAXPO,2
+      I3J0=I3J0+1
+      IF(M1MM2M.GT.LT)GO TO 60
+      IL=IL+2
+      L3=LT-1
+      NY=L3*L3+LT+M1MM2
+      AORS=L3+LT
+      F3J=(RAC3J(IL)*V3J0(I3J0)*DSQRT(AORS))*JSN
+      CFUN=(VJ0(NBHJ+LT)*VYJ0(NBY+NY))*F3J
+      CSUM=CSUM+CFUN
+   60 JSN=-JSN
+   62 CSUM=CSUM*CR
+      GHIT=GHIT+CSUM
+   64 CONTINUE
+   70 PI4=DACOS(0.0D0)*8.0D0
+      IF(IPAMO.EQ.1)GO TO 90
+      CR=PI4*(L1+L1PO)*(L2+L2+1)
+      CR=DSQRT(CR)
+      GHIT=GHIT*CR
+      RETURN
+   90 UIM=(0.0D0,1.0D0)
+      CR=L1PO
+      CR=(PI4*(L1+L1MP)*(L1+L1PO)*(L2+L2+1))/CR
+      CR=DSQRT(CR)
+      GHIT=GHIT*(CR*UIM)
+      RETURN
+      END
+      REAL*8 FUNCTION CGEV(IPAMO,MU,L,M)
+CCC   CGEV(IPAMO,MU,L,M)=CLGO(1,L-IPAMO,L;-MU,MU+M,M)      
+      REAL*8 XD,XN
+      IF(IPAMO.NE.0)GO TO 30
+      IF((M.NE.0).OR.(MU.NE.0))GO TO 10
+      CGEV=0.0D0
+      RETURN
+   10 IF(MU.EQ.0)GO TO 20
+      XD=(L+1)*L*2
+      IF(MU.GT.0)GO TO 15
+      XN=(L+M)*(L-M+1)
+      CGEV=DSQRT(XN/XD)
+      RETURN
+   15 XN=(L-M)*(L+M+1)
+      CGEV=-DSQRT(XN/XD)
+      RETURN
+   20 XD=(L+1)*L
+      XN=-M
+      CGEV=XN/DSQRT(XD)
+      RETURN
+   30 XD=(L*2-1)*L*2 
+      IF(MU)35,40,45
+   35 XN=(L-1+M)*(L+M)
+      GO TO 50
+   40 XN=(L-M)*(L+M)*2
+      GO TO 50
+   45 XN=(L-1-M)*(L-M)
+   50 CGEV=DSQRT(XN/XD)
+      RETURN
+      END
+      SUBROUTINE R3J000(J2,J3)
+CCC   3j(J,J2,J3;0,0,0)      
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C6/RAC3J(17)
+      ONE=1.0D0
+      JMX=J3+J2
+      IF(JMX.GT.0)GO TO 10
+      RAC3J(1)=ONE
+      RETURN
+   10 JMN=IABS(J3-J2)
+      NJMO=(JMX-JMN)/2
+      JF=JMX+JMX+1
+      ISN=1
+      IF(MOD(JMN,2).NE.0)ISN=-1
+      IF(NJMO.GT.0)GO TO 15
+      SJ=JF
+      CNR=(ONE/DSQRT(SJ))*ISN
+      RAC3J(1)=CNR
+      RETURN
+   15 SJR=JF
+      JMXPOS=(JMX+1)*(JMX+1)
+      JMNS=JMN*JMN
+      J1MO=JMX-1
+      J1S=(J1MO+1)*(J1MO+1)
+      CJ=(JMXPOS-J1S)*(J1S-JMNS)
+      CJ=DSQRT(CJ)
+      J1MOS=J1MO*J1MO
+      CJMO=(JMXPOS-J1MOS)*(J1MOS-JMNS)
+      CJMO=DSQRT(CJMO)
+      IF(NJMO.GT.1)GO TO 20
+      RAC3J(1)=-CJ/CJMO
+      SJ=SJR+(RAC3J(1)*RAC3J(1))*(JF-4)
+      CNR=(ONE/DSQRT(SJ))*ISN
+      RAC3J(2)=CNR
+      RAC3J(1)=RAC3J(1)*CNR
+      RETURN
+   20 NJ=NJMO+1
+      NMAT=(NJ+1)/2
+      RAC3J(NJ)=ONE
+      RAC3J(NJMO)=-CJ/CJMO
+      IF(NMAT.EQ.NJMO)GO TO 50
+      NBR=NJMO-NMAT
+      DO 45 IBR=1,NBR
+      IRR=NJ-IBR
+      JF=JF-4
+      J1MO=J1MO-2
+      J1S=(J1MO+1)*(J1MO+1)
+      CJ=(JMXPOS-J1S)*(J1S-JMNS)
+      CJ=DSQRT(CJ)
+      J1MOS=J1MO*J1MO
+      CJMO=(JMXPOS-J1MOS)*(J1MOS-JMNS)
+      CJMO=DSQRT(CJMO)
+      RAC3J(IRR-1)=RAC3J(IRR)*(-CJ/CJMO)
+   45 SJR=SJR+(RAC3J(IRR)*RAC3J(IRR))*JF
+   50 RACMAT=RAC3J(NMAT)
+      SJR=SJR+(RACMAT*RACMAT)*(JF-4)
+      RAC3J(1)=ONE
+      JF=JMN+JMN+1
+      SJL=JF
+      J1PT=JMN+2
+      J1POS=(J1PT-1)*(J1PT-1)
+      CJPO=(JMXPOS-J1POS)*(J1POS-JMNS)
+      CJPO=DSQRT(CJPO)
+      J1PTS=J1PT*J1PT
+      CJPT=(JMXPOS-J1PTS)*(J1PTS-JMNS)
+      CJPT=DSQRT(CJPT)
+      RAC3J(2)=-CJPO/CJPT
+      NMATMO=NMAT-1
+      IF(NMATMO.LT.2)GO TO 75
+      DO 70 IRL=2,NMATMO
+      JF=JF+4
+      J1PT=J1PT+2
+      J1POS=(J1PT-1)*(J1PT-1)
+      CJPO=(JMXPOS-J1POS)*(J1POS-JMNS)
+      CJPO=DSQRT(CJPO)
+      J1PTS=J1PT*J1PT
+      CJPT=(JMXPOS-J1PTS)*(J1PTS-JMNS)
+      CJPT=DSQRT(CJPT)
+      RAC3J(IRL+1)=RAC3J(IRL)*(-CJPO/CJPT)
+   70 SJL=SJL+(RAC3J(IRL)*RAC3J(IRL))*JF
+   75 RATRAC=RACMAT/RAC3J(NMAT)
+      RATS=RATRAC*RATRAC
+      SJ=SJR+SJL*RATS
+      RAC3J(NMAT)=RACMAT
+      CNR=(ONE/DSQRT(SJ))*ISN
+      DO 80 IRR=NMAT,NJ
+   80 RAC3J(IRR)=RAC3J(IRR)*CNR
+      CNL=CNR*RATRAC
+      DO 85 IRL=1,NMATMO
+   85 RAC3J(IRL)=RAC3J(IRL)*CNL
+      RETURN
+      END
+      SUBROUTINE R3JJR(J2,J3,M2,M3)
+CCC   3j(J,J2,J3;-M2-M3,M2,M3)      
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C6/RAC3J(17)
+      ONE=1.0D0
+      JMX=J3+J2
+      JDF=J3-J2
+      M1=-M2-M3
+      JMN=MAX0(IABS(JDF),IABS(M1))
+      NJMO=JMX-JMN
+      JF=JMX+JMX+1
+      ISN=1
+      IF(MOD(JDF+M1,2).NE.0)ISN=-1
+      IF(NJMO.GT.0)GO TO 15
+      SJ=JF
+      CNR=(ONE/DSQRT(SJ))*ISN
+      RAC3J(1)=CNR
+      RETURN
+   15 SJT=ONE
+      SJR=JF
+      JSMPOS=(JMX+1)*(JMX+1)
+      JDFS=JDF*JDF
+      M1S=M1*M1
+      MDF=M3-M2
+      IDJC=M1*(J3*(J3+1)-J2*(J2+1))
+      J1=JMX
+      J1S=J1*J1
+      J1PO=J1+1
+      CCJ=(J1S-JDFS)*(J1S-M1S)
+      CJ=CCJ*(JSMPOS-J1S)
+      CJ=DSQRT(CJ)
+      DJ=JF*(J1*J1PO*MDF+IDJC)
+      IF(NJMO.GT.1)GO TO 20
+      RAC3J(1)=-DJ/(CJ*J1PO)
+      SJ=SJR+(RAC3J(1)*RAC3J(1))*(JF-2)
+      CNR=(ONE/DSQRT(SJ))*ISN
+      RAC3J(2)=CNR
+      RAC3J(1)=RAC3J(1)*CNR
+      RETURN
+   20 NJ=NJMO+1
+      NMAT=(NJ+1)/2
+      RAC3J(NJ)=ONE
+      RAC3J(NJMO)=-DJ/(CJ*J1PO)
+      IF(NMAT.EQ.NJMO)GO TO 50
+      NBR=NJMO-NMAT
+      DO 45 IBR=1,NBR
+      IRR=NJ-IBR
+      JF=JF-2
+      J1=J1-1
+      J1S=J1*J1
+      J1PO=J1+1
+      CJP=CJ
+      CCJ=(J1S-JDFS)*(J1S-M1S)
+      CJ=CCJ*(JSMPOS-J1S)
+      CJ=DSQRT(CJ)
+      SJT=RAC3J(IRR)*RAC3J(IRR)
+      DJ=JF*(J1*J1PO*MDF+IDJC)
+      RAC3J(IRR-1)=-(RAC3J(IRR)*DJ+RAC3J(IRR+1)*(CJP*J1))/(CJ*J1PO)
+   45 SJR=SJR+SJT*JF
+   50 OSJT=SJT
+      SJT=RAC3J(NMAT)*RAC3J(NMAT)
+      IF(SJT.LT.OSJT)GO TO 55
+      SJR=SJR+SJT*(JF-2)
+      GO TO 60
+   55 NMAT=NMAT+1
+   60 RACMAT=RAC3J(NMAT)
+      RAC3J(1)=ONE
+      JF=JMN+JMN+1
+      SJL=JF
+      J1=JMN
+      IF(J1.EQ.0)GO TO 62
+      J1PO=J1+1
+      J1POS=J1PO*J1PO
+      CCJP=(J1POS-JDFS)*(J1POS-M1S)
+      CJP=CCJP*(JSMPOS-J1POS)
+      CJP=DSQRT(CJP)
+      DJ=JF*(J1*J1PO*MDF+IDJC)
+      RAC3J(2)=-DJ/(CJP*J1)
+      GO TO 63
+   62 CJP=JSMPOS-1
+      CJP=DSQRT(CJP)
+      DJ=MDF
+      RAC3J(2)=-DJ/CJP
+   63 NMATMO=NMAT-1
+      IF(NMATMO.LT.2)GO TO 75
+      DO 70 IRL=2,NMATMO
+      JF=JF+2
+      J1=J1+1
+      J1PO=J1+1
+      J1POS=J1PO*J1PO
+      CJ=CJP
+      CCJP=(J1POS-JDFS)*(J1POS-M1S)
+      CJP=CCJP*(JSMPOS-J1POS)
+      CJP=DSQRT(CJP)
+      SJT=RAC3J(IRL)*RAC3J(IRL)
+      DJ=JF*(J1*J1PO*MDF+IDJC)
+      RAC3J(IRL+1)=-(RAC3J(IRL)*DJ+RAC3J(IRL-1)*(CJ*J1PO))/(CJP*J1)
+   70 SJL=SJL+SJT*JF
+   75 RATRAC=RACMAT/RAC3J(NMAT)
+      RATS=RATRAC*RATRAC
+      SJ=SJR+SJL*RATS
+      RAC3J(NMAT)=RACMAT
+      CNR=(ONE/DSQRT(SJ))*ISN
+      DO 80 IRR=NMAT,NJ
+   80 RAC3J(IRR)=RAC3J(IRR)*CNR
+      CNL=CNR*RATRAC
+      DO 85 IRL=1,NMATMO
+   85 RAC3J(IRL)=RAC3J(IRL)*CNL
+      RETURN
+      END
+      SUBROUTINE R3JMR(J1,J2,J3,M1)
+CCC   3j(J1,J2,J3;M1,M,-M1-M)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C6/RAC3J(17)
+      ONE=1.0D0
+      MMX=MIN0(J2,J3-M1)
+      MMN=MAX0(-J2,-J3-M1)
+      NMMO=MMX-MMN
+      J1PO=J1+1
+      J1TPO=J1PO+J1
+      ISN=1
+      IF(MOD(J2-J3-M1,2).NE.0)ISN=-1
+      IF(NMMO.GT.0)GO TO 15
+      SJ=J1TPO
+      CNR=(ONE/DSQRT(SJ))*ISN
+      RAC3J(1)=CNR
+      RETURN
+   15 J1S=J1*J1PO
+      J2PO=J2+1
+      J2S=J2*J2PO
+      J3PO=J3+1
+      J3S=J3*J3PO
+      ID=J1S-J2S-J3S
+      M2=MMX
+      M3=M1+M2
+      CM=(J2PO-M2)*(J2+M2)*(J3PO-M3)*(J3+M3)
+      CM=DSQRT(CM)
+      DM=ID+M2*M3*2
+      IF(NMMO.GT.1)GO TO 20
+      RAC3J(1)=DM/CM
+      SJ=(ONE+RAC3J(1)*RAC3J(1))*J1TPO
+      CNR=(ONE/DSQRT(SJ))*ISN
+      RAC3J(2)=CNR
+      RAC3J(1)=RAC3J(1)*CNR
+      RETURN
+   20 NM=NMMO+1
+      NMAT=(NM+1)/2
+      RAC3J(NM)=ONE
+      RAC3J(NMMO)=DM/CM
+      SJT=ONE
+      SJR=ONE
+      IF(NMAT.EQ.NMMO)GO TO 50
+      NBR=NMMO-NMAT
+      DO 45 IBR=1,NBR
+      IRR=NM-IBR
+      M2=M2-1
+      M3=M1+M2
+      CMP=CM
+      CM=(J2PO-M2)*(J2+M2)*(J3PO-M3)*(J3+M3)
+      CM=DSQRT(CM)
+      SJT=RAC3J(IRR)*RAC3J(IRR)
+      DM=ID+M2*M3*2
+      RAC3J(IRR-1)=(RAC3J(IRR)*DM-RAC3J(IRR+1)*CMP)/CM
+   45 SJR=SJR+SJT
+   50 OSJT=SJT
+      SJT=RAC3J(NMAT)*RAC3J(NMAT)
+      IF(SJT.LT.OSJT)GO TO 55
+      SJR=SJR+SJT
+      GO TO 60
+   55 NMAT=NMAT+1
+   60 RACMAT=RAC3J(NMAT)
+      RAC3J(1)=ONE
+      M2=MMN
+      M3=M1+M2
+      CMP=(J2-M2)*(J2PO+M2)*(J3-M3)*(J3PO+M3)
+      CMP=DSQRT(CMP)
+      DM=ID+M2*M3*2
+      RAC3J(2)=DM/CMP
+      SJL=ONE
+      NMATMO=NMAT-1
+      IF(NMATMO.LT.2)GO TO 75
+      DO 70 IRL=2,NMATMO
+      M2=M2+1
+      M3=M1+M2
+      CM=CMP
+      CMP=(J2-M2)*(J2PO+M2)*(J3-M3)*(J3PO+M3)
+      CMP=DSQRT(CMP)
+      SJT=RAC3J(IRL)*RAC3J(IRL)
+      DM=ID+M2*M3*2
+      RAC3J(IRL+1)=(RAC3J(IRL)*DM-RAC3J(IRL-1)*CM)/CMP
+   70 SJL=SJL+SJT
+   75 RATRAC=RACMAT/RAC3J(NMAT)
+      RATS=RATRAC*RATRAC
+      SJ=(SJR+SJL*RATS)*J1TPO
+      RAC3J(NMAT)=RACMAT
+      CNR=(ONE/DSQRT(SJ))*ISN
+      DO 80 IRR=NMAT,NM
+   80 RAC3J(IRR)=RAC3J(IRR)*CNR
+      CNL=CNR*RATRAC
+      DO 85 IRL=1,NMATMO
+   85 RAC3J(IRL)=RAC3J(IRL)*CNL
+      RETURN
+      END
+      SUBROUTINE DME(LI,I,NPNT,NPNTTS,VK,EXDC,EXRI,JER,LCALC,ARG)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION CFJ(LIPT),RFJ(LIPT),RFN(LIPT),
+CCC  1FBI(LIPT),FB(LIPT),FN(LIPT),
+CCC  2RMF(LI),DRMF(LI),REF(LI),DREF(LI)
+      DIMENSION CFJ(10),RFJ(10),RFN(10),
+     1FBI(10),FB(10),FN(10),RMF(8),DRMF(8),REF(8),DREF(8)
+      COMPLEX*16 CFJ,FBI,FB,FN,RMF,DRMF,REF,DREF
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+      COMPLEX*16 RIS,DLRI,DC0,VKT
+      COMPLEX*16 DFBI,DFB,DFN,CCNA,CCNB,CCNC,CCND,
+     1Y1,DY1,Y2,DY2,ARIN,ARG,CRI,UIM
+      JER=0
+      UIM=(0.0D0,1.0D0)
+      NSTP=NPNT-1
+      NSTPTS=NPNTTS-1
+      LIPO=LI+1
+      LIPT=LI+2
+      SZ=VK*ROS(I)
+      VSZ(I)=SZ
+      VKR1=VK*RC(I,1)
+      NSH=NSHL(I)
+      VKT(I)=CDSQRT(DC0(1))
+      ARG=VKR1*VKT(I)
+      ARIN=ARG
+      IF(DIMAG(ARG).EQ.0.0D0)GO TO 26
+      CALL CBF(LIPO,ARG,LCALC,CFJ)
+      IF(LCALC.GE.LIPO)GO TO 22
+      JER=5
+      RETURN
+   22 DO 24 J=1,LIPT
+   24 FBI(J)=CFJ(J)
+      GO TO 32
+   26 RARG=DREAL(ARG)
+      CALL RBF(LIPO,RARG,LCALC,RFJ)
+      IF(LCALC.GE.LIPO)GO TO 28
+      JER=5
+      RETURN
+   28 DO 30 J=1,LIPT
+   30 FBI(J)=RFJ(J)
+   32 AREX=SZ*EXRI
+      ARG=AREX
+      RARG=AREX
+      CALL RBF(LIPO,RARG,LCALC,RFJ)
+      IF(LCALC.GE.LIPO)GO TO 40
+      JER=7
+      RETURN
+   40 CALL RNF(LIPO,RARG,LCALC,RFN)
+      IF(LCALC.GE.LIPO)GO TO 42
+      JER=8
+      RETURN
+   42 DO 43 J=1,LIPT
+      FB(J)=RFJ(J)
+   43 FN(J)=RFN(J)
+      IF(NSH.GT.1)GO TO 65
+      CRI=DC0(1)/EXDC
+      DO 60 L=1,LI
+      LPO=L+1
+      LTPO=LPO+L
+      LPT=LPO+1
+      DFBI=(L*FBI(L)-LPO*FBI(LPT))*ARIN+FBI(LPO)*LTPO
+      DFB=(L*FB(L)-LPO*FB(LPT))*AREX+FB(LPO)*LTPO
+      DFN=(L*FN(L)-LPO*FN(LPT))*AREX+FN(LPO)*LTPO
+      CCNA=FBI(LPO)*DFN
+      CCNB=FN(LPO)*DFBI
+      CCNC=FBI(LPO)*DFB
+      CCND=FB(LPO)*DFBI
+      RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
+   60 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+      RETURN
+   65 DO 80 L=1,LI
+      LPO=L+1
+      LTPO=LPO+L
+      LPT=LPO+1
+      DLTPO=LTPO
+cccccccccccccccccccccc
+      Y1=FBI(LPO)
+      DY1=(L*FBI(L)-LPO*FBI(LPT))*VKT(I)/DLTPO
+cccccccccccccccccccccc
+      Y2=Y1
+      DY2=DY1
+      IC=1
+      DO 76 NS=2,NSH
+      NSMO=NS-1
+      VKR=VK*RC(I,NSMO)
+      IF(MOD(NS,2).EQ.0)GO TO 70
+      IC=IC+1
+      STEP=NSTP
+      STEP=VK*(RC(I,NS)-RC(I,NSMO))/STEP
+      ARG=DC0(IC)
+      CALL RKC(NSTP,STEP,ARG,VKR,LPO,Y1,Y2,DY1,DY2)
+      GO TO 76
+   70 CALL DIEL(NSTPTS,NSMO,I,IC,VK)
+      STEPTS=NSTPTS
+      STEPTS=VK*(RC(I,NS)-RC(I,NSMO))/STEPTS
+      CALL RKT(NSTPTS,STEPTS,VKR,LPO,Y1,Y2,DY1,DY2)
+   76 CONTINUE
+      RMF(L)=Y1*SZ
+      DRMF(L)=DY1*SZ+Y1
+      REF(L)=Y2*SZ
+   80 DREF(L)=DY2*SZ+Y2
+      CRI=(1.0D0,0.0D0)
+      IF(MOD(NSH,2).NE.0)CRI=DC0(IC)/EXDC
+      DO 90 L=1,LI
+      LPO=L+1
+      LTPO=LPO+L
+      LPT=LPO+1
+      DFB=(L*FB(L)-LPO*FB(LPT))*AREX+FB(LPO)*LTPO
+      DFN=(L*FN(L)-LPO*FN(LPT))*AREX+FN(LPO)*LTPO
+      CCNA=RMF(L)*DFN
+      CCNB=DRMF(L)*FN(LPO)*SZ*LTPO
+      CCNC=RMF(L)*DFB
+      CCND=DRMF(L)*FB(LPO)*SZ*LTPO
+      RMI(L,I)=1.0D0+UIM*(CCNA-CCNB)/(CCNC-CCND)
+      CCNA=REF(L)*DFN
+      CCNB=DREF(L)*FN(LPO)*SZ*LTPO
+      CCNC=REF(L)*DFB
+      CCND=DREF(L)*FB(LPO)*SZ*LTPO
+   90 REI(L,I)=1.0D0+UIM*(CRI*CCNA-CCNB)/(CRI*CCNC-CCND)
+      RETURN
+      END
+      SUBROUTINE RKT(NPNTMO,STEP,X,LPO,Y1,Y2,DY1,DY2)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+      COMPLEX*16 RIS,DLRI,DC0,VKT
+      COMPLEX*16 Y1,DY1,Y2,DY2,CY1,CDY1,C11,
+     1CY23,CDY23,YC2,C12,C13,CY4,CDY4,YY,C14,C21,C22,C23,C24
+      HSTEP=.5D0*STEP
+      CL=LPO*(LPO-1)
+      DO 60 IPNT=1,NPNTMO
+      JPNT=IPNT+IPNT-1
+      CY1=CL/(X*X)-RIS(JPNT)
+      CDY1=-2.0D0/X
+      C11=(CY1*Y1+CDY1*DY1)*STEP
+      XH=X+HSTEP
+      JPNTPO=JPNT+1
+      CY23=CL/(XH*XH)-RIS(JPNTPO)
+      CDY23=-2.0D0/XH
+      YC2=Y1+DY1*HSTEP
+      C12=(CY23*YC2+CDY23*(DY1+.5D0*C11))*STEP
+      C13=(CY23*(YC2+.25D0*C11*STEP)+CDY23*(DY1+.5D0*C12))*STEP
+      XN=X+STEP
+      JPNTPT=JPNT+2
+      CY4=CL/(XN*XN)-RIS(JPNTPT)
+      CDY4=-2.0D0/XN
+      YY=Y1+DY1*STEP
+      C14=(CY4*(YY+.5D0*C12*STEP)+CDY4*(DY1+C13))*STEP
+      Y1=YY+(C11+C12+C13)*STEP/6.0D0
+      DY1=DY1+(.5D0*C11+C12+C13+.5D0*C14)/3.0D0
+      CY1=CY1-CDY1*DLRI(JPNT)
+      CDY1=CDY1+2.0D0*DLRI(JPNT)
+      C21=(CY1*Y2+CDY1*DY2)*STEP
+      CY23=CY23-CDY23*DLRI(JPNTPO)
+      CDY23=CDY23+2.0D0*DLRI(JPNTPO)
+      YC2=Y2+DY2*HSTEP
+      C22=(CY23*YC2+CDY23*(DY2+.5D0*C21))*STEP
+      C23=(CY23*(YC2+.25D0*C21*STEP)+CDY23*(DY2+.5D0*C22))*STEP
+      CY4=CY4-CDY4*DLRI(JPNTPT)
+      CDY4=CDY4+2.0D0*DLRI(JPNTPT)
+      YY=Y2+DY2*STEP
+      C24=(CY4*(YC2+.5D0*C22*STEP)+CDY4*(DY2+C23))*STEP
+      Y2=YY+(C21+C22+C23)*STEP/6.0D0
+      DY2=DY2+(.5D0*C21+C22+C23+.5D0*C24)/3.0D0
+   60 X=XN
+      RETURN
+      END
+      SUBROUTINE RKC(NPNTMO,STEP,DCC,X,LPO,Y1,Y2,DY1,DY2)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMPLEX*16 Y1,DY1,Y2,DY2,CY1,CDY1,C11,
+     1CY23,YC2,C12,C13,CY4,YY,C14,C21,C22,C23,C24,DCC
+      HSTEP=.5D0*STEP
+      CL=LPO*(LPO-1)
+      DO 60 IPNT=1,NPNTMO
+      CY1=CL/(X*X)-DCC
+      CDY1=-2.0D0/X
+      C11=(CY1*Y1+CDY1*DY1)*STEP
+      XH=X+HSTEP
+      CY23=CL/(XH*XH)-DCC
+      CDY23=-2.0D0/XH
+      YC2=Y1+DY1*HSTEP
+      C12=(CY23*YC2+CDY23*(DY1+.5D0*C11))*STEP
+      C13=(CY23*(YC2+.25D0*C11*STEP)+CDY23*(DY1+.5D0*C12))*STEP
+      XN=X+STEP
+      CY4=CL/(XN*XN)-DCC
+      CDY4=-2.0D0/XN
+      YY=Y1+DY1*STEP
+      C14=(CY4*(YY+.5D0*C12*STEP)+CDY4*(DY1+C13))*STEP
+      Y1=YY+(C11+C12+C13)*STEP/6.0D0
+      DY1=DY1+(.5D0*C11+C12+C13+.5D0*C14)/3.0D0
+      C21=(CY1*Y2+CDY1*DY2)*STEP
+      YC2=Y2+DY2*HSTEP
+      C22=(CY23*YC2+CDY23*(DY2+.5D0*C21))*STEP
+      C23=(CY23*(YC2+.25D0*C21*STEP)+CDY23*(DY2+.5D0*C22))*STEP
+      YY=Y2+DY2*STEP
+      C24=(CY4*(YC2+.5D0*C22*STEP)+CDY4*(DY2+C23))*STEP
+      Y2=YY+(C21+C22+C23)*STEP/6.0D0
+      DY2=DY2+(.5D0*C21+C22+C23+.5D0*C24)/3.0D0
+   60 X=XN
+      RETURN
+      END
+      SUBROUTINE DIEL(NPNTMO,NS,I,IC,VK)
+      IMPLICIT REAL*8(A-H,O-Z)
+      COMMON/C1/VH(255),VJ0(102),VYHJ(4335),VYJ0(1734),VJ(1),
+     1RMI(8,6),REI(8,6),W(160,4),AM0M(160,160),
+     2FSAS(6),SAS(6,2,2),VINTS(6,16),VINTT(16),
+     3FSAC(2,2),SAC(2,2),FSACM(2,2),
+     4VINT(16),VINTM(16),SCSCP(2),ECSCP(2),SCSCPM(2),ECSCPM(2),
+     5V3J0(276),SCSC(2),ECSC(2),SCSCM(2),ECSCM(2),
+     6SSCS(6),SEXS(6),SABS(6),SQSCS(6),SQEXS(6),SQABS(6),GCSV(6),
+     7RXX(6),RYY(6),RZZ(6),ROS(6),RC(6,8),
+     8IND3J(9,8),IOG(6),NSHL(6)
+      COMPLEX*16 VH,VJ0,VYHJ,VYJ0,VJ,RMI,REI,W,AM0M,FSAS,SAS,
+     1VINTS,VINTT,FSAC,SAC,FSACM,VINT,VINTM,SCSCP,ECSCP,SCSCPM,ECSCPM
+      COMMON/C2/RIS(999),DLRI(999),DC0(5),VKT(6),VSZ(6)
+      COMPLEX*16 RIS,DLRI,DC0,VKT
+      COMPLEX*16 DEL
+      DIF=RC(I,NS+1)-RC(I,NS)
+      HSTEP=NPNTMO
+      HSTEP=.5D0*DIF/HSTEP
+      RR=RC(I,NS)
+      DEL=DC0(IC+1)-DC0(IC)
+      KPNT=NPNTMO+NPNTMO
+      RIS(KPNT+1)=DC0(IC+1)
+      DLRI(KPNT+1)=(0.0D0,0.0D0)
+      DO 90 NP=1,KPNT
+      FF=(RR-RC(I,NS))/DIF
+      RIS(NP)=DEL*FF*FF*(-2.0D0*FF+3.0D0)+DC0(IC)
+      DLRI(NP)=3.0D0*DEL*FF*(1.0D0-FF)/(DIF*VK*RIS(NP))
+   90 RR=RR+HSTEP
+      RETURN
+      END
+      SUBROUTINE RBF(N,X,NM,SJ)
+C
+C     FROM SPHJ OF LIBRARY specfun 
+C
+C     =======================================================
+C       Purpose: Compute spherical Bessel functions j
+C       Input :  X --- Argument of j
+C                N --- Order of j  ( N = 0,1,2,... )
+C       Output:  SJ(N+1) --- j
+C                NM --- Highest order computed
+C       Routines called:
+C                MSTA1 and MSTA2 for computing the starting
+C                point for backward recurrence
+C     =======================================================
+C
+      IMPLICIT REAL*8 (A-H,O-Z)
+      DIMENSION SJ(N+1)
+      A0=DABS(X)
+      NM=N
+      IF(A0.LT.1.0D-60)THEN
+      DO  K=2,N+1
+      SJ(K)=0.0D0
+      ENDDO
+      SJ(1)=1.0D0
+      RETURN
+      ENDIF
+      SJ(1)=DSIN(X)/X
+      IF(N.EQ.0)RETURN
+      SJ(2)=(SJ(1)-DCOS(X))/X
+      IF(N.EQ.1)RETURN
+      SA=SJ(1)
+      SB=SJ(2)
+      M=MSTA1(A0,200)
+      IF(M.LT.N)THEN
+      NM=M
+      ELSE
+      M=MSTA2(A0,N,15)
+      ENDIF
+      F0=0.0D0
+      F1=1.0D0-100
+      DO  K=M,0,-1
+      F=(2.0D0*K+3.0D0)*F1/X-F0
+      IF(K.LE.NM)SJ(K+1)=F
+      F0=F1
+      F1=F
+      ENDDO
+      IF(DABS(SA).GT.DABS(SB))CS=SA/F
+      IF(DABS(SA).LE.DABS(SB))CS=SB/F0
+      DO K=1,NM+1
+      SJ(K)=CS*SJ(K)
+      ENDDO
+      END
+      SUBROUTINE CBF(N,Z,NM,CSJ)
+C
+C     FROM CSPHJY OF LIBRARY specfun 
+C
+C     ==========================================================
+C       Purpose: Compute spherical Bessel functions j
+C       Input :  Z --- Complex argument of j
+C                N --- Order of j ( N = 0,1,2,... )
+C       Output:  CSJ(N+1) --- j
+C                NM --- Highest order computed
+C       Routines called:
+C                MSTA1 and MSTA2 for computing the starting
+C                point for backward recurrence
+C     ==========================================================
+C
+      IMPLICIT COMPLEX*16 (C,Z)
+      REAL*8 A0
+      DIMENSION CSJ(N+1)
+      A0=CDABS(Z)
+      NM=N
+      IF(A0.LT.1.0D-60)THEN
+      DO K=2,N+1
+      CSJ(K)=0.0D0
+      ENDDO
+      CSJ(1)=(1.0D0,0.0D0)
+      RETURN
+      ENDIF
+      CSJ(1)=CDSIN(Z)/Z
+      IF(N.EQ.0)RETURN
+      CSJ(2)=(CSJ(1)-CDCOS(Z))/Z
+      IF(N.EQ.1)RETURN
+      CSA=CSJ(1)
+      CSB=CSJ(2)
+      M=MSTA1(A0,200)
+      IF(M.LT.N)THEN
+      NM=M
+      ELSE
+      M=MSTA2(A0,N,15)
+      ENDIF
+      CF0=0.0D0
+      CF1=1.0D0-100
+      DO K=M,0,-1
+      CF=(2.0D0*K+3.0D0)*CF1/Z-CF0
+      IF(K.LE.NM)CSJ(K+1)=CF
+      CF0=CF1
+      CF1=CF
+      ENDDO
+      IF(CDABS(CSA).GT.CDABS(CSB))CS=CSA/CF
+      IF(CDABS(CSA).LE.CDABS(CSB))CS=CSB/CF0
+      DO K=1,NM+1
+      CSJ(K)=CS*CSJ(K)
+      ENDDO
+      END
+      INTEGER FUNCTION MSTA1(X,MP)
+C
+C     ===================================================
+C       Purpose: Determine the starting point for backward  
+C                recurrence such that the magnitude of all    
+C                J or j at that point is about 10**(-MP)
+C       Input :  X     --- Abs of argument of J or j
+C                MP    --- Value of magnitude
+C       Output:  MSTA1 --- Starting point   
+C     ===================================================
+C
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      A0=DABS(X)
+      N0=INT(1.1*A0)+1
+      F0=ENVJ(N0,A0)-MP
+      N1=N0+5
+      F1=ENVJ(N1,A0)-MP
+      DO 10 IT=1,20             
+      NN=N1-IDINT((N1-N0)/(1.0D0-F0/F1))
+      F=ENVJ(NN,A0)-MP
+      IF(ABS(NN-N1).LT.1)GO TO 20
+      N0=N1
+      F0=F1
+      N1=NN
+   10 F1=F
+   20 MSTA1=NN
+      RETURN
+      END
+      INTEGER FUNCTION MSTA2(X,N,MP)
+C
+C     ===================================================
+C       Purpose: Determine the starting point for backward
+C                recurrence such that all J or j
+C                have MP significant digits
+C       Input :  X  --- Abs of argument of J or j
+C                N --- Order of J or j
+C                MP --- Significant digit
+C       Output:  MSTA2 --- Starting point
+C     ===================================================
+C
+      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+      A0=DABS(X)
+      HMP=0.5D0*MP
+      EJN=ENVJ(N,A0)
+      IF(EJN.LE.HMP)THEN
+      OBJ=MP
+C
+      N0=INT(1.1*A0)+1
+C
+      ELSE
+      OBJ=HMP+EJN
+      N0=N
+      ENDIF
+      F0=ENVJ(N0,A0)-OBJ
+      N1=N0+5
+      F1=ENVJ(N1,A0)-OBJ
+      DO 10 IT=1,20
+      NN=N1-IDINT((N1-N0)/(1.0D0-F0/F1))
+      F=ENVJ(NN,A0)-OBJ
+      IF(ABS(NN-N1).LT.1)GO TO 20
+      N0=N1
+      F0=F1
+      N1=NN
+   10 F1=F
+   20 MSTA2=NN+10
+      RETURN
+      END
+      REAL*8 FUNCTION ENVJ(N,X)
+      DOUBLE PRECISION X,XN
+C
+      IF(N.EQ.0)THEN
+      XN=1.0D-100
+      ENVJ=0.5D0*DLOG10(6.28D0*XN)-XN*DLOG10(1.36D0*X/XN)
+      ELSE
+C
+      ENVJ=0.5D0*DLOG10(6.28D0*N)-N*DLOG10(1.36D0*X/N)
+C
+      ENDIF
+C
+      RETURN
+      END
+      SUBROUTINE RNF(N,X,NM,SY)
+C
+C     FROM SPHJY OF LIBRARY specfun 
+C
+C
+C     ======================================================
+C       Purpose: Compute spherical Bessel functions y
+C       Input :  X --- Argument of y ( X > 0 )
+C                N --- Order of y ( N = 0,1,2,... )
+C       Output:  SY(N+1) --- y
+C                NM --- Highest order computed
+C     ======================================================
+C
+      IMPLICIT REAL*8 (A-H,O-Z)
+      DIMENSION SY(N+1)
+      IF(X.LT.1.0D-60)THEN
+      DO K=1,N+1
+      SY(K)=-1.0D+300
+      ENDDO
+      RETURN
+      ENDIF
+      SY(1)=-DCOS(X)/X
+      IF(N.EQ.0)RETURN
+      SY(2)=(SY(1)-DSIN(X))/X
+      IF(N.EQ.1)RETURN
+      F0=SY(1)
+      F1=SY(2)
+      DO K=2,N
+      F=(2.0D0*K-1.0D0)*F1/X-F0
+      SY(K+1)=F
+      IF(DABS(F).GE.1.0D+300)GO TO 20              
+      F0=F1
+      F1=F
+      ENDDO
+      RETURN
+   20 NM=K
+      RETURN
+      END
+      SUBROUTINE POLAR(X,Y,Z,R,CTH,STH,CPH,SPH)
+      IMPLICIT REAL*8(A-H,O-Z)
+      LOGICAL ONX,ONY,ONZ
+      ONE=1.0D0
+      C0=0.0D0
+      ONX=Y.EQ.C0
+      ONY=X.EQ.C0
+      ONZ=ONX.AND.ONY
+      IF(ONZ)GO TO 10
+      IF(ONX)GO TO 15
+      IF(ONY)GO TO 20
+      RHOS=X*X+Y*Y
+      RHO=DSQRT(RHOS)
+      CPH=X/RHO
+      SPH=Y/RHO
+      GO TO 25
+   10 CPH=ONE
+      SPH=C0
+      GO TO 25
+   15 RHOS=X*X
+      RHO=DABS(X)
+      CPH=DSIGN(ONE,X)
+      SPH=C0
+      GO TO 25
+   20 RHOS=Y*Y
+      RHO=DABS(Y)
+      CPH=C0
+      SPH=DSIGN(ONE,Y)
+   25 IF(Z.NE.C0)GO TO 35
+      IF(ONZ)GO TO 30
+      R=RHO
+      CTH=C0
+      STH=ONE
+      RETURN
+   30 R=C0
+      CTH=ONE
+      STH=C0
+      RETURN
+   35 IF(ONZ)GO TO 40
+      R=DSQRT(RHOS+Z*Z)
+      CTH=Z/R
+      STH=RHO/R
+      RETURN
+   40 R=DABS(Z)
+      CTH=DSIGN(ONE,Z)
+      STH=C0
+      RETURN
+      END
+      SUBROUTINE SPHAR(COSRTH,SINRTH,COSRPH,SINRPH,LL,YLM)
+      IMPLICIT REAL*8(A-H,O-Z)
+CCC   DIMENSION SINRMP(LL),COSRMP(LL),PLEGN((LLPO*LL)/2+LLPO),
+CCC  1YLM(LLPO*LLPO)
+CCC   LL=MAX0(LIT,LMT)
+      DIMENSION SINRMP(16),COSRMP(16),PLEGN(153),YLM(289)
+      COMPLEX*16 YLM
+      PI4=DACOS(0.0D0)*8.0D0
+      PI4IRS=1.0D0/DSQRT(PI4)
+      X=COSRTH
+      Y=DABS(SINRTH)
+      CLLMO=3.0D0
+      CLL=1.5D0
+      YTOL=Y
+      PLEGN(1)=1.0D0
+      PLEGN(2)=X*DSQRT(CLLMO)
+      PLEGN(3)=YTOL*DSQRT(CLL)
+      SINRMP(1)=SINRPH
+      COSRMP(1)=COSRPH
+      IF(LL.LT.2)GO TO 30
+      K=3
+      DO 20 L=2,LL
+      LMO=L-1
+      LTPO=L+L+1
+      LTMO=LTPO-2
+      LTS=LTPO*LTMO
+      CN=LTS
+      DO 10 MPO=1,LMO
+      M=MPO-1
+      MPOPK=MPO+K
+      LS=(L+M)*(L-M)
+      CD=LS
+      CNM=LTPO*(LMO+M)*(L-MPO)
+      CDM=(LTMO-2)*LS
+   10 PLEGN(MPOPK)=PLEGN(MPOPK-L)*X*DSQRT(CN/CD)-
+     1PLEGN(MPOPK-LTMO)*DSQRT(CNM/CDM)
+      LPK=L+K
+      CLTPO=LTPO
+      PLEGN(LPK)=PLEGN(K)*X*DSQRT(CLTPO)
+      K=LPK+1
+      CLT=LTPO-1
+      CLL=CLL*(CLTPO/CLT)
+      YTOL=YTOL*Y
+      PLEGN(K)=YTOL*DSQRT(CLL)
+      SINRMP(L)=SINRPH*COSRMP(LMO)+COSRPH*SINRMP(LMO)
+   20 COSRMP(L)=COSRPH*COSRMP(LMO)-SINRPH*SINRMP(LMO)
+   30 L=0
+   40 M=0
+      K=L*(L+1)
+      L0Y=K+1
+      L0P=K/2+1
+      YLM(L0Y)=PI4IRS*PLEGN(L0P)
+      GO TO 45
+   44 LMP=L0P+M
+      SAVE=PI4IRS*PLEGN(LMP)
+      LMY=L0Y+M
+      YLM(LMY)=SAVE*DCMPLX(COSRMP(M),SINRMP(M))
+      IF(MOD(M,2).NE.0)YLM(LMY)=-YLM(LMY)
+      LMY=L0Y-M
+      YLM(LMY)=SAVE*DCMPLX(COSRMP(M),-SINRMP(M))
+ 45   IF(M.GE.L)GO TO 47
+      M=M+1
+      GO TO 44
+ 47   IF(L.GE.LL) RETURN
+      L=L+1
+      GO TO 40
+      END
+      SUBROUTINE LUCIN(AM,NDDMST,N,IER)
+C
+C          NDDMST  FIRST DIMENSION OF AM AS DECLARED IN DIMENSION
+C                  STATEMENT.
+C          N       NUMBER OF ROWS IN AM.
+C          IER     IS REPLACED BY 1 FOR SINGULARITY.
+C
+C
+C     CALLS COMPLEX INNER PRODUCT FUNCTION CDTP.
+C
+      IMPLICIT REAL*8 (A-H,O-Z)
+CCC   DIMENSION AM(NDDMST,NDDMST),V(NDDMST)
+      DIMENSION AM(960,960),V(960)
+      COMPLEX*16 AM
+      COMPLEX*16 CTEMP,CDTP
+      IER=0
+      NMINUS=N-1
+C
+C  1.2 DEAL WITH THE SPECIAL CASE N=1.
+C
+C1200 IF(N.NE.1)GO TO 1300
+C     CTEMP=AM(1,1)
+C     IF(CTEMP.EQ.(0.0D0,0.0D0))GO TO 5200
+C     AM(1,1)=1.0D0/CTEMP
+C     RETURN
+C
+C  1.3 SET V(I)=1/R(I)**2, WHERE R(I) IS THE LENGTH OF ROW I.
+C
+C1300 DO 1309 I=1,N
+      DO 1309 I=1,N
+      SUM=0.0D0
+      DO 1319 J=1,N
+      SUM=SUM+DREAL(AM(I,J))**2+DIMAG(AM(I,J))**2
+ 1319 CONTINUE
+      V(I)=1.0D0/SUM
+ 1309 CONTINUE
+C
+C  2.  REPLACE AM BY TRIANGULAR MATRICES (L,U) WHERE AM=L*U.
+C      REPLACE L(I,I) BY 1/L(I,I), READY FOR SECTION 4.
+C      (ROW INTERCHANGES TAKE PLACE, AND THE INDICES OF THE PIVOTAL ROWS
+C      ARE PLACED IN V.)
+C
+      DO 2019 K=1,N
+      KPLUS=K+1
+      KMINUS=K-1
+      L=K
+      PSQMAX=0.0D0
+      DO 2029 I=K,N
+      CTEMP=-CDTP(-AM(I,K),AM,I,1,K,KMINUS,NDDMST)
+      AM(I,K)=CTEMP
+      PSQ=V(I)*(DREAL(CTEMP)**2+DIMAG(CTEMP)**2)
+      IF(PSQ.LE.PSQMAX)GO TO 2029
+      PSQMAX=PSQ
+      L=I
+ 2029 CONTINUE
+      IF(L.EQ.K)GO TO 2011
+      DO 2049 J=1,N
+      CTEMP=AM(K,J)
+      AM(K,J)=AM(L,J)
+      AM(L,J)=CTEMP
+ 2049 CONTINUE
+      V(L)=V(K)
+ 2011 V(K)=L
+      IF(PSQMAX.EQ.0.0D0)GO TO 5200
+      CTEMP=1.0D0/AM(K,K)
+      AM(K,K)=CTEMP
+      IF(KPLUS.GT.N)GO TO 2019
+      DO 2059 J=KPLUS,N
+      AM(K,J)=-CTEMP*CDTP(-AM(K,J),AM,K,1,J,KMINUS,NDDMST)
+ 2059 CONTINUE
+ 2019 CONTINUE
+C
+C  4.  REPLACE AM BY ITS INVERSE AMINV.
+C
+C
+C  4.1 REPLACE L AND U BY THEIR INVERSES LINV AND UINV.
+C
+      DO 4109 K=1,NMINUS
+      KPLUS=K+1
+      DO 4119 I=KPLUS,N
+      AM(I,K)=-AM(I,I)*CDTP((0.0D0,0.0D0),AM,I,K,K,I-K,
+     1NDDMST)
+      AM(K,I)=-CDTP(AM(K,I),AM,K,KPLUS,I,I-K-1,NDDMST)
+ 4119 CONTINUE
+ 4109 CONTINUE
+C
+C  4.2 FORM AMINV=UINV*LINV.
+C
+      DO 4209 K=1,N
+      DO 4219 I=1,N
+      IF(I.GE.K)GO TO 4212
+      AM(I,K)=CDTP((0.0D0,0.0D0),AM,I,K,K,N-K+1,NDDMST)
+      GO TO 4219
+ 4212 AM(I,K)=CDTP(AM(I,K),AM,I,I+1,K,N-I,NDDMST)
+ 4219 CONTINUE
+ 4209 CONTINUE
+C
+C  4.3 INTERCHANGE COLUMNS OF AMINV AS SPECIFIED BY V, BUT IN REVERSE
+C      ORDER.
+C
+      DO 4309 L=1,N
+      K=N-L+1
+      KCOL=IDINT(V(K))
+      IF(KCOL.EQ.K)GO TO 4309
+      DO 4319 I=1,N
+      CTEMP=AM(I,K)
+      AM(I,K)=AM(I,KCOL)
+      AM(I,KCOL)=CTEMP
+ 4319 CONTINUE
+ 4309 CONTINUE
+      RETURN
+ 5200 IER=1
+      RETURN
+      END
+      COMPLEX*16 FUNCTION CDTP(Z,AV,I,JF,K,NJ,ISTEP)
+CCC   DIMENSION AV(NDDMST*NDDMST)
+      DIMENSION AV(921600)
+      COMPLEX*16 AV,Z
+      CDTP=Z
+      IF(NJ.LE.0) RETURN
+      JL=JF+NJ-1
+      IR=(K-1)*ISTEP
+      DO 3 J=JF,JL
+ 3    CDTP=CDTP+AV((J-1)*ISTEP+I)*AV(J+IR)
+      RETURN
+      END
+CCC
\ No newline at end of file