Sources Esope | Cast3M

gflex0
C GFLEX0    SOURCE    CHAT      05/01/13    00:17:58     5004
      SUBROUTINE GFLEX0(AB,DLL,RF,CTC,DELTAT,LANBN)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C
C  =====================================================================
C    SOUS-PROGRAMME FORTRAN APPELE PAR GREEN1 POUR CALCULER LES FCTS DE
C    GREEN EN FLEXION
C
C    LES MODIFICATIONS APPORTEES AU PROGRAMME GFLEX1 SONT LES SUIVANTES:
C        LE CALCUL DES INTEGRALES DE FRESNEL EST REALISE PAR LA METHODE
C        DE LANCZOS (SUBROUTINE FRESNE)
C
C        LES VALEURS DE LA FONCTION DE GREEN ET DE SES DERIVEES EN L AU
C        PREMIER PAS DE TEMPS NE SONT PAS NULLES
C
C    CREATION    : 21/09/87
C    PROGRAMMEUR : VACELET (10/03/89)
C  =====================================================================
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
      DIMENSION AB(10,*),CC(1),SS(1),U2(1)
      WRITE(IOIMP,*) ' DEBUT DE GFLEX0 '
      PIS4=XPI*0.25D0
      RPI=1.D0/SQRT(XPI)
C
C--------BOUCLE SUR LES PAS DE TEMPS -------------------------
C
      GP4=0.D0
      GP3=0.D0
      GP2=0.D0
      GP1=0.D0
      GP0=0.D0
      GPL0=0.D0
      G4P=0.D0
      G3P=0.D0
      G2P=0.D0
      G1P=0.D0
      G0P=0.D0
      DTM=DLL*DLL/(3.D0*CTC*RF*XPI)
      LANBN1=LANBN+1
      DO 40 L=1,LANBN1
      T=L*DELTAT
C
C--------ETUDE DES DEUX EXTREMITES----------------------------
C
      DO 20 NE=1,2
C
C  PREMIERE EXTREMITE
C
C
      IF(NE.EQ.1) THEN
C
      DL=0.D0
      DKSI=DL/RF
      CSRF=CTC/RF
      TETA=CSRF*T
      RTETA=SQRT(TETA)
      STK=0.5D0/RTETA
      U2(1)=DKSI*DKSI/(4.D0*TETA)
      SDK=RPI*SIN(U2(1)-PIS4)
      CDK=RPI*COS(U2(1)-PIS4)
C
        G4=-STK*CDK
        G3=0.5D0
        G2=RTETA*SDK
        G1=0.D0
        G0=(2.D0/3.D0)*RTETA*TETA*CDK
C
        AB4=G4-GP4
        GP4=G4
        AB3=G3-GP3
        GP3=G3
        AB2=G2-GP2
        GP2=G2
        AB1=G1-GP1
        GP1=G1
        AB0=G0-GP0
        GP0=G0
C
        AB(1,L)=AB0
        AB(2,L)=AB1
        AB(3,L)=AB2
        AB(4,L)=AB3
        AB(5,L)=AB4
      ELSE
C
C  SECONDE EXTREMITE
C
        DL=DLL
        DKSI=DL/RF
        CSRF=CTC/RF
        TETA=CSRF*T
        RTETA=SQRT(TETA)
        STK=0.5D0/RTETA
        U2(1)=DKSI*DKSI/(4.D0*TETA)
        SDK=RPI*SIN(U2(1)-PIS4)
        CDK=RPI*COS(U2(1)-PIS4)
C
        CALL FRESNE(CC,SS,U2,1)
C
        U=SQRT(U2(1))
        U3=U*U2(1)
        U4=U*U3
        SC=SS(1)-CC(1)
        SC1=1.D0-SS(1)-CC(1)
C
        GG=RTETA*(CDK+U*SC)
        G4=-STK*CDK
        G3=0.5D0*SC1
        G2=RTETA*(SDK+U*SC1)
        G1=RTETA*(0.5D0*RTETA*SC+U*G2)
        G0=(2.D0/3.D0)*RTETA*(RTETA*GG+U*G1)
C
        GG3=TETA*(0.5D0*SC1-U*CDK-U2(1)*SC)
        GG2=(2.D0/3.D0)*TETA*RTETA*(SDK+1.5D0*U*SC1-U2(1)*CDK-U3*SC)
        GG1=TETA*RTETA*(1.5D0*SC+U*SDK+2.D0*(U3*CDK+U4*SC))/6.D0
        GG1=RTETA*(GG1+U*GG2)
C
        IF(DELTAT.LE.DTM) THEN
          AB0=G0-GPL0
          GPL0=G0
          DELT=CSRF*DELTAT
          DELT2=DELT*DELT
          AB(6,L)=AB0
          AB(7,L)=GG1/DELT
          AB(8,L)=GG2/DELT
          AB(9,L)=GG3/DELT
          AB(10,L)=-AB0/DELT2
        ELSE
          AB(10,L)=G4-G4P
          G4P=G4
          AB(9,L)=G3-G3P
          G3P=G3
          AB(8,L)=G2-G2P
          G2P=G2
          AB(7,L)=G1-G1P
          G1P=G1
          AB(6,L)=G0-G0P
          G0P=G0
        ENDIF
      ENDIF
   20 CONTINUE
   40 CONTINUE
C
      IF(DELTAT.GT.DTM) RETURN
      DO 70 I=7,10
      GP=AB(I,2)-2.D0*AB(I,1)
      DO 50 L=2,LANBN
      GA=AB(I,L+1)-2.D0*AB(I,L)+AB(I,L-1)
      AB(I,L-1)=GP
      GP=GA
   50 CONTINUE
      AB(I,LANBN)=GP
      GP=AB(I,1)
      DO 60 L=2,LANBN
      GA=0.5D0*(AB(I,L-1)+AB(I,L))
      AB(I,L-1)=GP
      GP=GA
   60 CONTINUE
      AB(I,LANBN)=GP
   70 CONTINUE
C
C
C
      WRITE(IOIMP,*)' FIN DE GFLEX0 '
      RETURN
      END