Sources Esope | Cast3M

gflex1
C GFLEX1    SOURCE    CHAT      05/01/13    00:18:02     5004
      SUBROUTINE GFLEX1(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    CREATION    : 21/09/87
C    PROGRAMMEUR : GUILBAUD
C  =====================================================================
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
      DIMENSION AB(10,*)
      WRITE(IOIMP,*) ' DEBUT DE GFLEX1 '
      PIS4=XPI*0.25D0
      RPI=1.D0/SQRT(XPI)
      XP=1.D0
      XG=20.D0
      XTG=200.D0
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--------BOUCLE SUR LES EXTREMITES----------------------------
C
      DO 20 NE=1,2
      DL=0.D0
      IF(NE.EQ.2) DL=DLL
      DKSI=DL/RF
      CSRF=CTC/RF
      TETA=CSRF*T
      RTETA=SQRT(TETA)
      STK=0.5D0/RTETA
      U2=DKSI*DKSI/(4.D0*TETA)
      SDK=RPI*SIN(U2-PIS4)
      CDK=RPI*COS(U2-PIS4)
      IF(NE.EQ.1) THEN
C
C  PREMIERE EXTREMITE
C
        G4=-STK*CDK
        G3=0.5D0
        G2=RTETA*SDK
        G1=0.D0
        G0=(2.D0/3.D0)*RTETA*TETA*CDK
        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
        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
        IF(U2.GT.XTG) THEN
          TK=2.D0*TETA/DKSI
          TK1=STK*TK
          TK2=TK1*TK
          TK3=TK2*TK
          TK4=TK3*TK
          TK5=TK4*TK
          TK6=TK5*TK
          G3=-TK1*SDK
          G2=TK2*CDK
          G1=TK3*SDK
          G0=-TK4*CDK
          GG3=-TK3*CDK
          GG2=-TK4*SDK
          GG1= TK5*CDK
        ELSE
          IF(T.EQ.DELTAT) THEN
            U2P=XG+1.D0
          ELSE
            U2P=DKSI*DKSI/(4.D0*CSRF*(T-DELTAT))
          ENDIF
          IF(U2.GT.XP.AND.U2P.LT.XG) THEN
            CALL INTFR2(DCC,DSS,U2P,U2)
            SS=SS-DSS
            CC=CC-DCC
          ELSE
            CALL INTFRE(CC,SS,U2)
          ENDIF
          SC=SS-CC
          SC1=1.D0-SS-CC
          U=SQRT(U2)
          U3=U*U2
          U4=U*U3
          GG4=-RTETA*(U*SC+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*GG4+U*G1)
          GG3=TETA*(0.5D0*SC1-U*CDK-U2*SC)
          GG2=(2.D0/3.D0)*TETA*RTETA*(SDK+1.5D0*U*SC1-U2*CDK-U3*SC)
          GG1=TETA*RTETA*(1.5D0*SC+U*SDK+2.D0*(U3*CDK+U4*SC))/6.D0
          GG1=RTETA*(GG1+U*GG2)
        ENDIF
        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
          G4=-STK*CDK
          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=0.D0
      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
      WRITE(IOIMP,*)' FIN DE GFLEX1 '
      RETURN
      END