Sources Esope | Cast3M

distrt
C DISTRT    SOURCE    CHAT      05/01/12    22:51:42     5004
       SUBROUTINE DISTRT(MMTRA,NNT,DT, DPI, TE,ID,ERR, RES,IOK,INEW)
      IMPLICIT INTEGER(I-N)
       IMPLICIT REAL*8(A-H,O-Z)
C=======================================================================
C  =  CALCUL DU MAXIMUM DE DISTRIBUTION DANS LE CAS                    =
C  =  NON STATIONNAIRE (ET STATIONAIRE)                                =
C  =                                                                   =
C  =  ENTREE XLTIME(1,NNT)=XL0,XL1,XL2 MOMENT STATIQUES AU COURS DU    =
C  =                       TEMPS                                       =
C  =         ID=1 NEWMARK/GUMBEL 1 + SUBSTITUTION                      =
C  =         ID=2 NEWMARK/GUMBEL 1 + DICHOTOMIE                        =
C  =                                                                   =
C  =  ERREUR                                                           =
C  =         IOK=1    NUMERO INEXISTANT                                =
C  =  ID=1   IOK=2    NON CONVERGENCE                                  =
C  =  ID=2   IOK=2    PAS D'ENCADREMENT INFERIEUR                      =
C  =         IOK=3    PAS D'ENCADREMENT SUPERIEUR                      =
C  =         IOK=100  NON CONVERGENCE INFERIEUR                        =
C  =         IOK=101  NON CONVERGENCE SUPERIEURE                       =
C=======================================================================
      REAL*8 L1,L2,N1,N2
      LOGICAL LMOD2
C======+++++++++++ bon fontionnement erreur PSRS/PRNS +++++++++++ ======
 
-INC PPARAM
-INC CCOPTIO
      SEGMENT,MMTRA
        IMPLIED XLTIME(NNT,3)
      ENDSEGMENT
C======+++++++++++ bon fontionnement erreur PSRS/PRNS +++++++++++ ======
C
      INEW=0
C
      IF(ID.LE.0.OR.ID.GT.2)THEN
        IOK=1
        CALL ERREUR(571)
        RETURN
      ENDIF
C
      GOTO (10,20),ID
C
C     1)METHODE DE NEWMARK/GUMBEL 1 (METHODE DE SUBSTITUTION)
C
 10   CONTINUE
C
C     1.1) FACTEUR DE GUMBEL 1
C
      A=LOG(-LOG(0.05D0))
      B=A-LOG(-LOG(0.95D0))
C
C     1.2) CALCUL DE LAM0 MOYEN ET DE LAM0 MAX
C
      XL0MA=XLTIME(1,1)
          IF(NNT.NE.1)THEN
      DO 100 IE1=2,NNT-1
        IF(XLTIME(IE1,1).GT.XL0MA)XL0MA=XLTIME(IE1,1)
 100    CONTINUE
      IF(XLTIME(NNT,1).GT.XL0MA)XL0MA=XLTIME(NNT,1)
          ENDIF
C
C     1.3) INITIALISATION DE L1 ET L2
C
      L2=SQRT(XL0MA)*3
      L1=L2/4
C
C     1.4) BOUCLE DE DETERMINATION DE L1 ET L2
C
      DO 101 IE1=1,20
C
C     1.4.1) CALCUL DES DISTRIBUTION DE NEWMARK
C
        CALL NEWMAR(L1,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >              DPI,NNT,N1)
        CALL NEWMAR(L2,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >              DPI,NNT,N2)
C       WRITE(6,'(I2,3X,2H1:,1PD13.6,1X,1PD13.6,3X,
C    >                  2H2:,1PD13.6,1X,1PD13.6,1X)')IE1,L1,N1,L2,N2
        INEW=INEW+2
C
C     1.4.2) CALCUL DES PARAMETRE DE GUMBEL
C
        AGUMB=LOG(LOG(N1)/LOG(N2))/(L2-L1)
        UGUMB=L1+LOG(-LOG(N1))/AGUMB
C
C     1.4.3) TEST DE CONVERGENCE
C
        IF((ABS(N1-.05D0)+ABS(N2-.95D0)).LT.ERR) GOTO 102
C
C     1.4.4) CALCUL DES NOUVEAUX L
C
        L1 = UGUMB - A / AGUMB
        L2 = L1    + B / AGUMB
C
 101    CONTINUE
C
C     1.5) NON CONVERGENCE
C
      IOK=2
      CALL ERREUR(572)
      RETURN
C
C     1.6) CONVERGENCE
C
 102  IOK=0
      RES=UGUMB + 0.57722/AGUMB
      RETURN
C
C     3) METHODE DE NEWMARK/GUMBEL 1 (METHODE DE DICHOTOMIE)
C
 20   CONTINUE
      IOK=0
C
C     2.1) CALCUL DE LAM0 MOYEN ET DE LAM0 MAX
C
      XL0MA=XLTIME(1,1)
           IF(NNT.NE.1)THEN
      DO 200 IE1=2,NNT-1
        IF(XLTIME(IE1,1).GT.XL0MA)XL0MA=XLTIME(IE1,1)
 200    CONTINUE
      IF(XLTIME(NNT,1).GT.XL0MA)XL0MA=XLTIME(NNT,1)
           ENDIF
C
C     2.2) INITIALISATION DE L1 ET L2
C
      L2=SQRT(XL0MA)*3
      L1=L2/4
C
C     2.3) VERIFICATION DES INITIALISATIONS
C
C     2.3.1) INITIALISATION
C
      CALL NEWMAR(L1,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >            DPI,NNT,N1)
      BL2=L1
      BN2=N1
C
      CALL NEWMAR(L2,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >            DPI,NNT,N2)
      BL1=L2
      BN1=N2
C
      INEW=INEW+2
C
C     2.3.2) L1 (AVEC RAFINEMENT EVENTUEL DES BORNES) ...
C
      DO 201 IE1=1,10
C       WRITE(6,'(1X,A5,2(1X,1PD14.7))')'L1,N1',L1,N1
        IF(N1.LT.0.05D0)GOTO 202
        BL1=L1
        BN1=N1
        L1=L1/2
        CALL NEWMAR(L1,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >              DPI,NNT,N1)
        INEW=INEW+1
 201    CONTINUE
C
C     2.3.3) ... ET ERREUR L1
C
      IOK=2
      MOTERR='inferieur'
      CALL ERREUR(573)
      RETURN
C
C     2.3.4) L2 (AVEC RAFINEMENT EVENTUEL DES BORNES) ...
C
 202  DO 203 IE1=1,10
C       WRITE(6,'(1X,A5,2(1X,1PD14.7))')'L2,N2',L2,N2
        IF(N2.GT.0.95D0)GOTO 204
        BL2=L2
        BN2=N2
        L2=L2*2
        CALL NEWMAR(L2,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >              DPI,NNT,N2)
        INEW=INEW+1
 203    CONTINUE
C
C     2.3.5) ... ET ERREUR L2
C
      IOK=3
      MOTERR='superieur'
      CALL ERREUR(573)
      RETURN
C
C     2.4) DETERMINATION DE L1
C
C     2.4.1) UNE DES BORNES EST-ELLE CORRECTE ???
C
 204  IF(ABS(N1-0.05D0).LT.ERR)GOTO 210
      IF(ABS(BN1-0.05D0).LT.ERR)THEN
        L1=BL1
        N1=BN1
        GOTO 210
      ENDIF
C
C     2.4.2) PEUT-ON AMELIORER LES ESTIMATIONS POUR L2 ?
C
      IF(BN1.GE.0.95D0)THEN
        L2=BL1
        N2=BN1
        LMOD2=.TRUE.
      ELSE
        LMOD2=.FALSE.
      ENDIF
C
C     2.4.3) BOUCLE DE DICHOTOMIE
C
      DO 205 IE1=1,25
C
C     2.4.3.1) SOLUTION AU MILIEU DE l'INNTERVALLE
C
        TL1=(L1+BL1)/2
        CALL NEWMAR(TL1,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >              DPI,NNT,TN1)
        INEW=INEW+1
C
C     2.4.3.2) AMELIORATION ESTIMATION L2 (SI POSSIBLE)
C
        IF (LMOD2)THEN
          IF(TN1.GE.0.95D0)THEN
            L2=TL1
            N2=TN1
          ELSE
            LMOD2=.FALSE.
            BL2=TL1
            BN2=TN1
          ENDIF
        ENDIF
C
C     2.4.3.3) TEST DE CONVERGENCE
C
        IF(ABS(TN1-0.05D0).LT.ERR)THEN
          L1=TL1
          N1=TN1
          GOTO 210
        ENDIF
C
C     2.4.3.4) RESTRICTION DE L'INTERVALLE DE RECHERCHE
C
        IF(TN1.LT.0.05D0)THEN
          L1=TL1
          N1=TN1
        ELSE
          BL1=TL1
          BN1=TN1
        ENDIF
C
 205    CONTINUE
C
C     2.4.4) NON CONVERGENCE
C
      IOK=100
      MOTERR='inferieure'
      CALL ERREUR(574)
C     RETURN (L'ERREUR 100 N'EST PAS ELIMINATOIRE)
C
C     2.5) DETERMINATION DE L2
C
C     2.5.1) UNE DES BORNES EST-ELLE CORRECTE ???
C
 210  IF(ABS(N2-0.95D0).LT.ERR)GOTO 220
      IF(ABS(BN2-0.95D0).LT.ERR)THEN
        L2=BL2
        N2=BN2
        GOTO 220
      ENDIF
C
C     2.5.2) BOUCLE DE DICHOTOMIE
C
      DO 215 IE1=1,25
C
C     2.5.2.1) SOLUTION AU MILIEU DE l'INTERVALLE
C
        TL2=(L2+BL2)/2
        CALL NEWMAR(TL2,DT,XLTIME(1,1),XLTIME(1,2),XLTIME(1,3),
     >              DPI,NNT,TN2)
        INEW=INEW+1
C
C     2.5.2.2) TEST DE CONVERGENCE
C
        IF(ABS(TN2-0.95D0).LT.ERR)THEN
          L2=TL2
          N2=TN2
          GOTO 220
        ENDIF
C
C     2.5.2.3) RESTRICTION DE L'INTERVALLE DE RECHERCHE
C
        IF(TN2.GT.0.95D0)THEN
          L2=TL2
          N2=TN2
        ELSE
          BL2=TL2
          BN2=TN2
        ENDIF
C
 215    CONTINUE
C
C     2.5.3) NON CONVERGENCE
C
      IOK=101 + IOK
      MOTERR='superieure'
      CALL ERREUR(574)
C     RETURN (L'ERREUR 101 N'EST PAS ELIMINATOIRE)
C            (ET SE COMBINE AVEC 100             )
C
C     2.6) CALCUL DES PARAMETRE DE GUMBEL
C
 220  AGUMB=LOG(LOG(N1)/LOG(N2))/(L2-L1)
      UGUMB=L1+LOG(-LOG(N1))/AGUMB
C
C     2.6) CONVERGENCE
C
      RES=UGUMB + 0.57722/AGUMB
C
C     WRITE(6,'(1X,A5,2(1X,1PD14.7))')'L1,N1',L1,N1
C     WRITE(6,'(1X,A5,2(1X,1PD14.7))')'L2,N2',L2,N2
C
      RETURN
C
      END