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 :
© Cast3M 2003 - All rights reserved.
Disclaimer