Télécharger la plaquette Cast3M


fr en

Télécharger ottocp.eso

Retour à la liste

Numérotation des lignes : 

ottocp

  1. C OTTOCP    SOURCE    CHAT      05/01/13    02:07:14     5004
  2.       SUBROUTINE OTTOCP(SIGMA,FCRIT,XLTR,DF,DG,H,
  3.      &    PRECIE,PRECIZ,XCOMP,XLC,KERRE)
  4. C=========================================================================
  5. C
  6. C  ENTREES :
  7. C    SIGMA,
  8. C
  9. C  SORTIES :
  10. C    FCRIT,
  11. C
  12. C==========================================================================
  13. C
  14.       IMPLICIT INTEGER(I-N)
  15.       IMPLICIT REAL*8(A-H,O-Z)
  16.  
  17. -INC PPARAM
  18. -INC CCOPTIO
  19. -INC CCREEL
  20. C
  21.       PARAMETER (XZER=0.D0,UNDEMI=.5D0,UN=1.D0,DEUX=2.D0,TROIS=3.D0)
  22. C
  23.       DIMENSION SIGMA(6),S(6),DF(6),DG(6),AUX(6),XCOMP(*)
  24. *
  25.       DIMENSION XLTR(3)
  26. C
  27. C     INITIALISATIONS
  28. C
  29.       KERRE=0
  30.       RAC3=SQRT(TROIS)
  31.  
  32. ******************************
  33.  
  34.       IZOB=0
  35.  
  36. ******************************
  37. *          CAS 1
  38. ******************************
  39.  
  40.       IF(IZOB.EQ.0) THEN
  41.  
  42. *
  43. *--------------------------------
  44. *     ON RECUPERE LES CONSTANTES
  45. *--------------------------------
  46. *
  47.       EPCMAX=XCOMP(1)
  48.       EPCULT=XCOMP(2)
  49.       RCBI  =XCOMP(3)
  50.       XK2   =XCOMP(4)
  51.       XGB   =XCOMP(5)
  52.       XPA   =XCOMP(6)
  53.       RCMAX =XCOMP(7)
  54.       XK1   =XCOMP(8)
  55.       XGA   =XCOMP(9)
  56.       AL    =XCOMP(10)
  57.       BE    =XCOMP(11)
  58.       GAMA  =XCOMP(12)
  59.       DELT  =XCOMP(13)
  60.       XMU1  =XCOMP(14)
  61.       XLCMAX=XCOMP(15)
  62.       A2    =XCOMP(16)
  63.       XMU2  =XCOMP(17)
  64.       XLCULT=XCOMP(18)
  65.  
  66.  
  67. *     IF(IIMPI.EQ.42) THEN
  68. *       WRITE(IOIMP,89010) (XCOMP(I),I=1,18)
  69. *89010   FORMAT( 2X, ' OTTOCP - XCOMP'/6(1X,1PE12.5)/)
  70. *       WRITE(IOIMP,80911) (SIGMA(I),I=1,6)
  71. *80911   FORMAT( 2X, ' ENTREE OTTOCP  - SIGMA '/6(1X,1PE12.5)/)
  72. *       WRITE(IOIMP,80912) XLC,XLCMAX,XLCULT
  73. *80912   FORMAT( 2X, ' XLC=',1PE12.5,2X,' XLCMAX=',1PE12.5,2X,
  74. *    &     ' XLCULT=',1PE12.5/)
  75. *     ENDIF
  76.  
  77. *
  78. *     CALCUL DES DIFFERENTES EXPRESSIONS
  79. *
  80.       XLAMC=RAC3*(UN+XGB-XGA/TROIS)
  81.       DELTAP=EPCULT-EPCMAX
  82.       IF(XLC.LE.XLCMAX) THEN
  83.         X=(GAMA-BE)*(EXP(XMU1*XLC)-UN)/(UN+DELT*EXP(XMU1*XLC))
  84.         SBAR=RCMAX*(UN+4.D0*X-DEUX*X*X)/TROIS
  85.         FMU=XGA*SBAR*SBAR/TROIS/RCMAX + XLAMC*SBAR/RAC3
  86.         FMU=FMU/(XGB*SBAR+RCMAX)
  87.         RR=RCMAX
  88.         XGBP=XGB*(XPA-(UN+XPA)*X)
  89.  
  90. *
  91.       ELSE IF(XLC.GT.XLCMAX.AND.XLC.LE.XLCULT) THEN
  92.         DLC=XLC-XLCMAX
  93.         TRA=EXP(XMU2*DLC)
  94.         X=A2*(TRA-UN)/(TRA+UN)
  95.         SBAR=RCMAX*(UN-X*X)
  96.         RR=XGA*SBAR*SBAR/TROIS/RCMAX +(UN-XGA/TROIS)*SBAR
  97.         XGBP=-XGB
  98.         FMU=1.D0
  99. *
  100.       ELSE
  101.         RR=0.D0
  102.         XGBP=-XGB
  103.         FMU=1.D0
  104.       ENDIF
  105.  
  106. *
  107.       XI1=TRACE(SIGMA)
  108.       PP=XI1/TROIS
  109. *
  110.       DO I=1,3
  111.          S(I)=SIGMA(I)-PP
  112.       ENDDO
  113.       DO I=4,6
  114.          S(I)=SIGMA(I)
  115.       ENDDO
  116. *
  117.       XJ2=UNDEMI*(S(1)*S(1) + S(2)*S(2) + S(3)*S(3))
  118.      &         +  S(4)*S(4) + S(5)*S(5) + S(6)*S(6)
  119.       TAU=SQRT(XJ2)
  120. *
  121.       XJ3= (S(1)**3 + S(2)**3 + S(3)**3)/TROIS
  122.      &    + S(4)*S(4)*(S(1)+S(2)) + S(5)*S(5)*(S(1)+S(3))
  123.      &    + S(6)*S(6)*(S(2)+S(3)) + DEUX*S(4)*S(5)*S(6)
  124. *
  125.       if(tau.NE.0.)THEN
  126.         C3T=1.5D0*RAC3*XJ3/(TAU**3)
  127.         C3T = MAX(-1.D0,C3T)
  128.         C3T = MIN( 1.D0,C3T)
  129.       else
  130.         C3T=1.
  131.       endif
  132.  
  133. *
  134.       IF(C3T.GE.0.D0) THEN
  135.          PHI=ACOS(XK2*C3T)/TROIS
  136.       ELSE
  137.          PHI=(XPI-ACOS(-XK2*C3T))/TROIS
  138.       ENDIF
  139.       XLAM=XK1*COS(PHI)
  140.  
  141. *
  142.       FCRIT=XGA*XJ2/RCMAX + XLAM*TAU +FMU*(XGB*XI1-RR)
  143.  
  144. *
  145. *
  146. *     IF(IIMPI.EQ.42) THEN
  147. *       WRITE(IOIMP,79000) XI1,TAU,XJ3
  148. *79000   FORMAT( 2X, ' OTTOCP - I1=',1PE12.5,2X,
  149. *    &       'TAU=',1PE12.5,2X,'J3=',1PE12.5/)
  150. *       WRITE(IOIMP,79001) ZOZOB,XLAM,PHI
  151. *79001   FORMAT( 2X, ' OTTOCP - TETA=',1PE12.5,2X,
  152. *    &       'XLAM=',1PE12.5,2X,'PHI=',1PE12.5/)
  153. *     ENDIF
  154.  
  155. *
  156. *  ECOULEMENT
  157. *  ----------
  158. *
  159.       DO I=1,6
  160.          DF(I)=XZER
  161.          DG(I)=XZER
  162.       ENDDO
  163.       H=XZER
  164. *
  165.       IF(TAU.GE.PRECIZ) THEN
  166. *
  167.         DLA=(XK1*XK2*RAC3/DEUX)*SIN(PHI)/SQRT(UN-(XK2*C3T)**2)
  168.  
  169.  
  170. *
  171.         AUX(1)=S(1)**2 + S(4)**2 + S(5)**2
  172.         AUX(2)=S(2)**2 + S(4)**2 + S(6)**2
  173.         AUX(3)=S(3)**2 + S(5)**2 + S(6)**2
  174.         AUX(4)=S(1)*S(4) + S(2)*S(4) + S(5)*S(6)
  175.         AUX(5)=S(1)*S(5) + S(4)*S(6) + S(5)*S(3)
  176.         AUX(6)=S(4)*S(5) + S(2)*S(6) + S(3)*S(6)
  177. *
  178.         TAU4=XJ2*XJ2
  179.         DEUTAU=DEUX*TAU
  180.         DEUTIE=DEUX/TROIS
  181. *
  182.         DO I=1,3
  183.            TRA =(XGA/RCMAX)*S(I)
  184.      &       + DLA*(AUX(I)/XJ2 -DEUTIE -1.5D0*XJ3*S(I)/TAU4 )
  185.      &       + XLAM*S(I)/DEUTAU
  186.  
  187.            DF(I)=TRA + XGB*FMU
  188.            DG(I)=TRA - XGBP
  189.         ENDDO
  190. *
  191.         DO I=4,6
  192.            TRA =(XGA/RCMAX)*S(I)
  193.      &       + DLA*(AUX(I)/XJ2  -1.5D0*XJ3*S(I)/TAU4 )
  194.      &       + XLAM*S(I)/DEUTAU
  195.            DF(I)=TRA
  196.            DG(I)=TRA
  197.         ENDDO
  198. *
  199.       ELSE
  200.  
  201.         DO I=1,3
  202.            DF(I)= XGB*FMU
  203.            DG(I)= - XGBP
  204.         ENDDO
  205. *
  206.       ENDIF
  207.  
  208. *
  209. *  ECROUISSAGE
  210. *  -----------
  211. *
  212.  
  213.       IF(XLC.LE.XLCMAX) THEN
  214.         H=XGB*XI1-RCMAX
  215.         H=H*(XGA*XGB*SBAR*SBAR/TROIS/RCMAX +DEUX*XGA*SBAR/TROIS
  216.      &       +XLAMC*RCMAX/RAC3)/((XGB*SBAR+RCMAX)**2)
  217.         H=H*4.D0*RCMAX*(UN-X)/TROIS
  218.         H=H*(GAMA-BE)*XMU1/(UN+DELT*EXP(XMU1*XLC))
  219.  
  220.       ELSE IF(XLC.GT.XLCMAX.AND.XLC.LE.XLCULT) THEN
  221.  
  222.         TRA=EXP(XMU2*DLC)
  223.         H= DEUX*XGA*SBAR/TROIS/RCMAX+UN-XGA/TROIS
  224.         H= H*DEUX*X*RCMAX*DEUX*A2*XMU2*TRA/((UN+TRA)**2)
  225.  
  226.       ELSE
  227.         H=0.D0
  228.       ENDIF
  229.  
  230.  
  231. ******************************
  232. *      FIN DU CAS 1
  233. ******************************
  234.  
  235.       ENDIF
  236. *
  237.  
  238. ******************************
  239. *          CAS 2
  240. ******************************
  241.  
  242.  
  243.       IF(IZOB.EQ.1) THEN
  244. *
  245. *--------------------------------
  246. *     ON RECUPERE LES CONSTANTES
  247. *--------------------------------
  248. *
  249.       EPCMAX=XCOMP(1)
  250.       EPCULT=XCOMP(2)
  251.       RCBI  =XCOMP(3)
  252.       REPS  =XCOMP(4)
  253.       A0    =XCOMP(6)
  254.       RCMAX =XCOMP(7)
  255.       AA    =XCOMP(8)
  256.       C0    =XCOMP(9)
  257.       FAC   =XCOMP(10)
  258.       XLCMAX=XCOMP(1)*FAC
  259.       XLCULT=XCOMP(2)*FAC
  260.  
  261.  
  262. *     IF(IIMPI.EQ.42) THEN
  263. *       WRITE(IOIMP,89010) (XCOMP(I),I=1,18)
  264. *89010   FORMAT( 2X, ' OTTOCP - XCOMP'/6(1X,1PE12.5)/)
  265. *       WRITE(IOIMP,80911) (SIGMA(I),I=1,6)
  266. *80911   FORMAT( 2X, ' ENTREE OTTOCP  - SIGMA '/6(1X,1PE12.5)/)
  267. *       WRITE(IOIMP,80912) XLC,XLCMAX,XLCULT
  268. *80912   FORMAT( 2X, ' XLC=',1PE12.5,2X,' XLCMAX=',1PE12.5,2X,
  269. *    &     ' XLCULT=',1PE12.5/)
  270. *     ENDIF
  271.  
  272. *
  273. *     CALCUL DES DIFFERENTES EXPRESSIONS
  274. *
  275.  
  276.       IF(XLC.LE.XLCMAX) THEN
  277.         X=XLC/XLCMAX
  278.         SBAR=RCMAX*(UN+4.D0*X-DEUX*X*X)/TROIS
  279.  
  280.       ELSE
  281.         X=(XLC-XLCMAX)/(XLCULT-XLCMAX)
  282.         SBAR=RCMAX*(UN-X*X)
  283.  
  284.       ENDIF
  285.  
  286.       TC= SBAR/AA
  287.       TC=MAX(1.D-6,TC)
  288.  
  289. *
  290.       XI1=TRACE(SIGMA)
  291.       PP=-XI1/TROIS
  292. *
  293.       DO I=1,3
  294.          S(I)=SIGMA(I)+PP
  295.       ENDDO
  296.       DO I=4,6
  297.          S(I)=SIGMA(I)
  298.       ENDDO
  299. *
  300.       XJ2=UNDEMI*(S(1)*S(1) + S(2)*S(2) + S(3)*S(3))
  301.      &         +  S(4)*S(4) + S(5)*S(5) + S(6)*S(6)
  302.       TAU=SQRT(XJ2)
  303.       QQ=SQRT(3.D0)*TAU
  304. *
  305.       ZEZE = 3.D0*QQ*QQ-A0*TC*PP
  306.       ZEZE = MAX(ZEZE,1.D-7)
  307.  
  308.       FCRIT=SQRT(ZEZE) - TC
  309.  
  310. *
  311. *     IF(IIMPI.EQ.42) THEN
  312. *       WRITE(IOIMP,79000) XI1,TAU,XJ3
  313. *79000   FORMAT( 2X, ' OTTOCP - I1=',1PE12.5,2X,
  314. *    &       'TAU=',1PE12.5,2X,'J3=',1PE12.5/)
  315. *       WRITE(IOIMP,79001) ZOZOB,XLAM,PHI
  316. *79001   FORMAT( 2X, ' OTTOCP - TETA=',1PE12.5,2X,
  317. *    &       'XLAM=',1PE12.5,2X,'PHI=',1PE12.5/)
  318. *     ENDIF
  319.  
  320. *
  321. *  ECOULEMENT
  322. *  ----------
  323. *
  324.       DO I=1,6
  325.          DF(I)=XZER
  326.          DG(I)=XZER
  327.       ENDDO
  328.       H=XZER
  329. *
  330.       FF = 1.D0+C0*(PP/TC)**2
  331.       DO I=1,3
  332.          DF(I)=(9.D0*S(I)+A0*TC)  /2.D0/TC
  333.          DG(I)=FF*DF(I)
  334.       ENDDO
  335. *
  336.       DO I=4,6
  337.          DF(I)=(9.D0*S(I))  /2.D0/TC
  338.          DG(I)=FF*DF(I)
  339.       ENDDO
  340. *
  341. *
  342. *  ECROUISSAGE
  343. *  -----------
  344. *
  345.  
  346.       IF(XLC.LE.XLCMAX) THEN
  347.         H=-4.D0*RCMAX*(UN-X)/TROIS
  348.         H=H*(1.D0+ A0*PP/2.D0/SQRT(ZEZE) )
  349.         H=H/XLCMAX/AA
  350.  
  351.       ELSE IF(XLC.GT.XLCMAX.AND.XLC.LE.XLCULT) THEN
  352.  
  353.         H= 2.D0*RCMAX*X
  354.         H=H*(1.D0+ A0*PP/2.D0/SQRT(ZEZE) )
  355.         H= H/(XLCULT-XLCMAX)/AA
  356.  
  357.       ELSE
  358.         H=0.D0
  359.       ENDIF
  360.  
  361. ******************************
  362. *      FIN DU CAS 2
  363. ******************************
  364.  
  365.       ENDIF
  366. *
  367.  
  368.  
  369. *     IF(IIMPI.EQ.42) THEN
  370. *      WRITE(IOIMP,77000) FCRIT
  371. *77000   FORMAT( 2X, ' OTTOCP - FCRIT '/(1X,1PE12.5)/)
  372. *      WRITE(IOIMP,77001) (DF(I),I=1,6)
  373. *77001   FORMAT( 2X, ' OTTOCP - DF '/6(1X,1PE12.5)/)
  374. *      WRITE(IOIMP,77002) (DG(I),I=1,6)
  375. *77002   FORMAT( 2X, ' OTTOCP - DG '/6(1X,1PE12.5)/)
  376. *      WRITE(IOIMP,77003) H
  377. *77003   FORMAT( 2X, ' OTTOCP - H     '/(1X,1PE12.5)/)
  378. *     ENDIF
  379. *
  380.       RETURN
  381.       END
  382.  
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  
  389.  
  390.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales