Télécharger ottocp.eso

Retour à la liste

Numérotation des lignes :

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

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