      PROGRAM LFFFT
CCC   140630
CCC   LE=20, LM=20, LE<=LM
CCC   FRPIF IS INTENSITY FACTOR ANALOGOUS TO I0 FOR FOCAL REGION POINTS;
CCC   JFT <= 0 --> OUTPUTS ON TAPE FORCE*(c/FRPIF) IN m**2; 
CCC   JFT >= 0 --> OUTPUTS ON TAPE TORQUE*(c/FRPIF)*(VK*EXRI)**3;
CCC   (TORQUE OUTPUTS)/(EXDC*VK*VK) ARE IN m**2
CCC   JSS = 1 USES TWS FOR INPUT
CCC   JTW = 1 FOR OUTPUTS ON PAPER
      IMPLICIT REAL*8(A-H,O-Z)
CCC   DIMENSION ZPV(LE,3,2,2),AC(NLEMT),WS(NLEMT)
      DIMENSION ZPV(20,3,2,2),AC(880),WS(880),
     1FFFE(3),FFFS(3),FFTE(3),FFTS(3)
      COMPLEX*16 AC,WS
CCC   DIMENSION XV(NXV),YV(NYV),ZV(NZV),WSL(NLMMT)
      DIMENSION XV(51),YV(51),ZV(51),WSL(160)
      COMPLEX*16 WSL
CCC   DIMENSION AM0M(NLEMT,NLEMT)
      DIMENSION AM0M(880,880)
      COMPLEX*16 AM0M
CCC   DIMENSION AMD(NVAM,4),INDAM(LE,LE)
CCC   NVAM=LE*LE+(LE*LEPO*LETPO)/3
      DIMENSION AMD(6140,4),INDAM(20,20)
      COMPLEX*16 AMD
CCC   DIMENSION TMSM(LE),TMSE(LE)
      DIMENSION TMSM(20),TMSE(20)
      COMPLEX*16 TMSM,TMSE
CCC   DIMENSION TMS(LE,3)
      DIMENSION TMS(20,3)
      COMPLEX*16 TMS
      COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM
      COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I
      COMPLEX*16 CC0,UIM,SQ2ITI
      IR=5
      IT=7
      OPEN(IR,FILE='../../test_data/trapping/DLFFFT',STATUS='OLD')
      READ(IR,*)JFT,JSS,JTW
      CLOSE(IR)
      OPEN(IT,FILE='TTMS',FORM='UNFORMATTED',STATUS='OLD')
      READ(IT)IS,LE
      READ(IT)VKS,EXRIS
      NLEM=LE*(LE+2)
      NLEMT=NLEM+NLEM
      IF(IS.GE.2222)GO TO 120
      IF(IS.GE.1111)GO TO 125
      IF(IS.GE.0)GO TO 135
      IF(IS.LT.0)GO TO 145
  120 CONTINUE
      READ(IT)(TMS(I,1),TMS(I,2),TMS(I,3),I=1,LM)
      GO TO 150
  125 CONTINUE
      READ(IT)(TMSM(I),TMSE(I),I=1,LE)
      GO TO 150
  135 CONTINUE
      READ(IT)((AM0M(I,J),J=1,NLEMT),I=1,NLEMT)
      GO TO 150
  145 CONTINUE
      NVAM=LE*LE+(LE*(LE+1)*(LE*2+1))/3
      READ(IT)((AMD(I,J),J=1,4),I=1,NVAM)
      READ(IT)((INDAM(I,J),J=1,LE),I=1,LE)
      READ(IT)MXMPO
      MXIM=MXMPO*2-1
  150 CLOSE(IT)
      IF(JSS.NE.1)
     1OPEN(IT,FILE='TFRFME',FORM='UNFORMATTED',STATUS='OLD')
      IF(JSS.EQ.1)
     1OPEN(IT,FILE='TWS',FORM='UNFORMATTED',STATUS='OLD')
      READ(IT)LMODE,LM,NKV,NXV,NYV,NZV
      IF(LM.LT.LE)GO TO 500
      READ(IT)VK,EXRI,AN,FF,TRA
      IF((VK.NE.VKS).OR.(EXRI.NE.EXRIS))GO TO 500
      READ(IT)SPD,FRSH,EXRIL
      READ(IT)(XV(I),I=1,NXV)
      READ(IT)(YV(I),I=1,NYV)
      READ(IT)(ZV(I),I=1,NZV)
      CC0=(0.0D0,0.0D0)
      UIM=(0.0D0,1.0D0)
      IF(JFT.GT.0)GO TO 155
      CALL THDPS(LE,ZPV)
      EXDC=EXRI*EXRI
      SQK=VK*VK*EXDC
      COF=1.0D0/SQK
      CIMU=COF/DSQRT(2.0D0)
      IF(JSS.EQ.1)GO TO 155
      ITF=8
      OPEN(ITF,FILE='TLFFF',FORM='UNFORMATTED',STATUS='UNKNOWN')
      WRITE(ITF)LMODE,LE,NKV,NXV,NYV,NZV
      WRITE(ITF)VK,EXRI,AN,FF,TRA
      WRITE(ITF)SPD,FRSH,EXRIL
      WRITE(ITF)(XV(I),I=1,NXV)
      WRITE(ITF)(YV(I),I=1,NYV)
      WRITE(ITF)(ZV(I),I=1,NZV)
      IF(JFT.LT.0)GO TO 160
  155 SQ2I=1.0D0/DSQRT(2.0D0)
      SQ2ITI=SQ2I*UIM
      IF(JSS.EQ.1)GO TO 160
      ITT=9
      OPEN(ITT,FILE='TLFFT',FORM='UNFORMATTED',STATUS='UNKNOWN')
      WRITE(ITT)LMODE,LE,NKV,NXV,NYV,NZV
      WRITE(ITT)VK,EXRI,AN,FF,TRA
      WRITE(ITT)SPD,FRSH,EXRIL
      WRITE(ITT)(XV(I),I=1,NXV)
      WRITE(ITT)(YV(I),I=1,NYV)
      WRITE(ITT)(ZV(I),I=1,NZV)
  160 NLMM=LM*(LM+2)
      NLMMT=NLMM+NLMM
      DO 475 IZ=1,NZV
      DO 475 IY=1,NYV
      DO 475 IX=1,NXV
      IF(LM.GT.LE)GO TO 170
      READ(IT)(WS(I),I=1,NLMMT)
      GO TO 180
  170 READ(IT)(WSL(I),I=1,NLMMT)
      DO 175 I=1,NLEM
      IE=I+NLEM
      IEL=I+NLMM
      WS(I)=WSL(I)
      WS(IE)=WSL(IEL)
  175 CONTINUE
  180 IF(IS.EQ.2222)GO TO 200
      IF(IS.EQ.1111)GO TO 205
      IF(IS.GT.0)GO TO 305
      IF(IS.LT.0)GO TO 405
  200 CONTINUE
      CALL SAMPOA(AC,TMS,WS)
      GO TO 445
  205 CONTINUE
      CALL SAMP(AC,TMSM,TMSE,WS)
      GO TO 445
  305 CONTINUE
      CALL CAMP(AC,AM0M,WS)
      GO TO 445
  405 CONTINUE
      CALL CZAMP(AC,AMD,INDAM,WS)
  445 CONTINUE
      IF(JFT.GT.0)GO TO 460
      CALL FFRF(ZPV,AC,WS,FFFE,FFFS)
      IF(JSS.NE.1)GO TO 450
      WRITE(66,*)FFFE(1),FFFS(1),FFFE(1)-FFFS(1)
      WRITE(66,*)FFFE(2),FFFS(2),FFFE(2)-FFFS(2)
      WRITE(66,*)FFFE(3),FFFS(3),FFFE(3)-FFFS(3)
      GO TO 455
  450 WRITE(ITF)(FFFE(I)-FFFS(I),I=1,3)
      IF(JTW.EQ.1)WRITE(66,6666)
     1IX,IY,IZ,FFFE(1)-FFFS(1),FFFE(2)-FFFS(2),FFFE(3)-FFFS(3)
  455 IF(JFT.LT.0)GO TO 475
  460 CALL FFRT(AC,WS,FFTE,FFTS)
      IF(JSS.NE.1)GO TO 470
      WRITE(67,*)FFTE(1),FFTS(1),FFTE(1)-FFTS(1)
      WRITE(67,*)FFTE(2),FFTS(2),FFTE(2)-FFTS(2)
      WRITE(67,*)FFTE(3),FFTS(3),FFTE(3)-FFTS(3)
      GO TO 475
  470 WRITE(ITT)(FFTE(I)-FFTS(I),I=1,3)
      IF(JTW.EQ.1)WRITE(67,6666)
     1IX,IY,IZ,FFTE(1)-FFTS(1),FFTE(2)-FFTS(2),FFTE(3)-FFTS(3)
  475 CONTINUE
      IF(JSS.EQ.1)GO TO 500
      IF(JFT.LE.0)CLOSE(ITF)
      IF(JFT.GE.0)CLOSE(ITT)
  500 CONTINUE
      CLOSE(IT)
      STOP
 6666 FORMAT(I5,I4,I4,1PD15.4,1PD12.4,1PD12.4)
      END
      SUBROUTINE SAMPOA(AC,TMS,WS)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM
CCC   DIMENSION AC(NLEMT),TMS(LE,3),WS(NLEMT)
      DIMENSION AC(880),TMS(20,3),WS(880)
      COMPLEX*16 AC,TMS,WS
      DIMENSION TM(2,2)
      COMPLEX*16 TM
      I=0
      DO 20 L=1,LE
      TM(1,1)=TMS(L,1)
      TM(1,2)=TMS(L,2)
      TM(2,2)=TMS(L,3)
      TM(2,1)=TM(1,2)
      LTPO=L+L+1
      DO 20 IM=1,LTPO
      I=I+1
      IE=I+NLEM
      AC(I)=TM(1,1)*WS(I)+TM(1,2)*WS(IE)
      AC(IE)=TM(2,1)*WS(I)+TM(2,2)*WS(IE)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE SAMP(AC,TMSM,TMSE,WS)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM
CCC   DIMENSION AC(NLEMT),TMSM(LE),TMSE(LE),WS(NLEMT)
      DIMENSION AC(880),TMSM(20),TMSE(20),WS(880)
      COMPLEX*16 AC,TMSM,TMSE,WS
      I=0
      DO 20 L=1,LE
      LTPO=L+L+1
      DO 20 IM=1,LTPO
      I=I+1
      IE=I+NLEM
      AC(I)=TMSM(L)*WS(I)
      AC(IE)=TMSE(L)*WS(IE)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE CAMP(AC,AM0M,WS)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM
      COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I
      COMPLEX*16 CC0,UIM,SQ2ITI
CCC   DIMENSION AC(NLEMT),AM0M(NLEMT,NLEMT),WS(NLEMT)
      DIMENSION AC(880),AM0M(880,880),WS(880)
      COMPLEX*16 AC,AM0M,WS
      DO 20 J=1,NLEMT
      AC(J)=CC0
      DO 20 I=1,NLEMT
      AC(J)=AC(J)+AM0M(J,I)*WS(I)
   20 CONTINUE
      RETURN
      END
      SUBROUTINE CZAMP(AC,AMD,INDAM,WS)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM
      COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I
      COMPLEX*16 CC0,UIM,SQ2ITI
CCC   DIMENSION AC(NLEMT),AMD(NVAM,4),INDAM(LE,LE),WS(NLEMT)
      DIMENSION AC(880),AMD(6140,4),INDAM(20,20),WS(880)
      COMPLEX*16 AC,AMD,WS
      COMPLEX*16 SUMM,SUME
      DO 10 I=1,NLEMT
   10 AC(I)=CC0
      DO 20 IM=1,MXIM
      M=IM-MXMPO
      LMN=MAX0(IABS(M),1)
      DO 20 L=LMN,LE
      I=M+L*(L+1)
      IE=I+NLEM
      SUMM=CC0
      SUME=CC0
      DO 15 LS=LMN,LE
      IS=M+LS*(LS+1)
      ISE=IS+NLEM
      NUM=INDAM(L,LS)+M
      SUMM=SUMM+AMD(NUM,1)*WS(IS)+AMD(NUM,2)*WS(ISE)
      SUME=SUME+AMD(NUM,3)*WS(IS)+AMD(NUM,4)*WS(ISE)
   15 CONTINUE
      AC(I)=SUMM
      AC(IE)=SUME
   20 CONTINUE
      RETURN
      END
      SUBROUTINE FFRF(ZPV,AC,WS,FFFE,FFFS)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM
      COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I
      COMPLEX*16 CC0,UIM,SQ2ITI
CCC   DIMENSION ZPV(LE,3,2,2),AC(NLEMT),WS(NLEMT)
      DIMENSION ZPV(20,3,2,2),AC(880),WS(880),FFFE(3),FFFS(3)
      COMPLEX*16 AC,WS
      DIMENSION GAP(3)
      COMPLEX*16 GAP
      COMPLEX*16 UIMMP,SUMM,SUME,SUEM,SUEE
      DO 50 IMU=1,3
      MU=IMU-2
      GAP(IMU)=CC0
      DO 40 L=1,LE
      LPO=L+1
      LTPO=LPO+L
      IMM=L*LPO
      DO 40 ILMP=1,3
      IF(((L.EQ.1).AND.(ILMP.EQ.1)).OR.((L.EQ.LE).AND.(ILMP.EQ.3)))
     1GO TO 40
      LMPML=ILMP-2
      LMP=L+LMPML
      UIMMP=-LMPML*UIM
      IMPMMMP=LMP*(LMP+1)
      DO 30 IM=1,LTPO
      M=IM-LPO
      MMP=M-MU
      IF(IABS(MMP).GT.LMP)GO TO 30
      I=IMM+M
      IE=I+NLEM
      IMP=IMPMMMP+MMP
      IMPE=IMP+NLEM
      CGC=CG1(LMPML,MU,L,M)
      SUMM=DCONJG(WS(I))*AC(IMP)
      SUME=DCONJG(WS(I))*AC(IMPE)
      SUEM=DCONJG(WS(IE))*AC(IMP)
      SUEE=DCONJG(WS(IE))*AC(IMPE)
      IF(LMPML.EQ.0)GO TO 25
      SUMM=SUMM*UIMMP
      SUME=SUME*UIMMP
      SUEM=SUEM*UIMMP
      SUEE=SUEE*UIMMP
   25 GAP(IMU)=GAP(IMU)+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))
   30 CONTINUE
   40 CONTINUE
   50 CONTINUE
      SUME=-GAP(1)*CIMU
      SUEE=-GAP(2)*COF
      SUEM=-GAP(3)*CIMU
      FFFE(1)=DREAL(SUME-SUEM)
      FFFE(2)=DREAL((SUME+SUEM)*UIM)
      FFFE(3)=DREAL(SUEE)
      DO 90 IMU=1,3
      MU=IMU-2
      GAP(IMU)=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)
      SUMM=DCONJG(AC(I))*AC(IMP)
      SUME=DCONJG(AC(I))*AC(IMPE)
      SUEM=DCONJG(AC(IE))*AC(IMP)
      SUEE=DCONJG(AC(IE))*AC(IMPE)
      IF(LMPML.EQ.0)GO TO 65
      SUMM=SUMM*UIMMP
      SUME=SUME*UIMMP
      SUEM=SUEM*UIMMP
      SUEE=SUEE*UIMMP
   65 GAP(IMU)=GAP(IMU)+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))
   70 CONTINUE
   80 CONTINUE
   90 CONTINUE
      SUME=GAP(1)*CIMU
      SUEE=GAP(2)*COF
      SUEM=GAP(3)*CIMU
      FFFS(1)=DREAL(SUME-SUEM)
      FFFS(2)=DREAL((SUME+SUEM)*UIM)
      FFFS(3)=DREAL(SUEE)
      RETURN
      END
      SUBROUTINE THDPS(LM,ZPV)
      IMPLICIT REAL*8(A-H,O-Z)
CCC   DIMENSION ZPV(LM,3,2,2)
      DIMENSION ZPV(20,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 FFRT(AC,WS,FFTE,FFTS)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/CIL/LE,NLEM,NLEMT,MXMPO,MXIM
      COMMON/CCR/CC0,UIM,SQ2ITI,COF,CIMU,SQ2I
      COMPLEX*16 CC0,UIM,SQ2ITI
CCC   DIMENSION AC(NLEMT),WS(NLEMT)
      DIMENSION AC(880),WS(880),FFTE(3),FFTS(3)
      COMPLEX*16 AC,WS
      DIMENSION CTQCE(3),CTQCS(3)
      COMPLEX*16 CTQCE,CTQCS
      COMPLEX*16 ACW,ACA
      CTQCE(1)=CC0
      CTQCE(2)=CC0
      CTQCE(3)=CC0
      CTQCS(1)=CC0
      CTQCS(2)=CC0
      CTQCS(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(AC(I))*WS(IMMU)+DCONJG(AC(IE))*WS(IMMUE)
      ACA=DCONJG(AC(I))*AC(IMMU)+DCONJG(AC(IE))*AC(IMMUE)
      CTQCE(1)=CTQCE(1)+ACW*RMU
      CTQCS(1)=CTQCS(1)+ACA*RMU
   30 RMU=-M
      ACW=DCONJG(AC(I))*WS(I)+DCONJG(AC(IE))*WS(IE)
      ACA=DCONJG(AC(I))*AC(I)+DCONJG(AC(IE))*AC(IE)
      CTQCE(2)=CTQCE(2)+ACW*RMU
      CTQCS(2)=CTQCS(2)+ACA*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(AC(I))*WS(IMMU)+DCONJG(AC(IE))*WS(IMMUE)
      ACA=DCONJG(AC(I))*AC(IMMU)+DCONJG(AC(IE))*AC(IMMUE)
      CTQCE(3)=CTQCE(3)+ACW*RMU
      CTQCS(3)=CTQCS(3)+ACA*RMU
   50 CONTINUE
   60 CONTINUE
      FFTE(1)=DREAL(CTQCE(1)-CTQCE(3))*SQ2I
      FFTE(2)=DREAL((CTQCE(1)+CTQCE(3))*SQ2ITI)
      FFTE(3)=DREAL(CTQCE(2))
      FFTS(1)=-DREAL(CTQCS(1)-CTQCS(3))*SQ2I
      FFTS(2)=-DREAL((CTQCS(1)+CTQCS(3))*SQ2ITI)
      FFTS(3)=-DREAL(CTQCS(2))
      RETURN
      END
CCC
