1 : C TRJPAR SOURCE CHAT 11/03/16 21:32:44 6902 2 : SUBROUTINE TRJPAR(MELEME,IZVIT,IZCENT,IELTFA,IFACEL,MTAB2, 3 : * IICAL,IZPOR,IZDIFF,IZDISP,MCHELM,MMODEL) 4 : C 5 : C 6 : C 7 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 8 : C 9 : C calcul de trajectoires appelé par TRAJEC decode la table de 10 : C lacher et pilote le calcul 11 : C 12 : C sous programmes appelés: 13 : C TRJCN2 controle des données 14 : C TRJINI premiere localisation et initialisation des 15 : C segments de travail 16 : C TRJAVA fait avancer une particule 17 : C TRJRAL en fin de calcul passe des coordonnées de 18 : C références aux coordonnées réelles 19 : C Les tables résultats sont chargées par ce sous programme 20 : C 21 : C ARGUMENTS 22 : C E MELEME pointeur du maillage de calcul ( DOMAINE.MAILLAGE) 23 : C E IZVIT pointeur du segment décrivant les vitesses . 24 : C Ce segment est créé dans TRJVIT et TRJFLU 25 : C E IZCENT pointeur de DOMAINE.CENTRE 26 : C E IELTFA pointeur de DOMAINE.ELTFA 27 : C E IFACEL pointeur de DOMAINE.FACEL 28 : C E MTAB2 pointeur de la table des lachers 29 : C 30 : C S MCHELM CHAMELEM resultat 31 : C S MMODEL modele sur lequel s'appuie MCHELM 32 : C 33 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 34 : C 35 : C 36 : IMPLICIT INTEGER(I-N) 37 : IMPLICIT REAL*8 (A-H,O-Z) 38 : PARAMETER(MPNO9=9) 39 : C MPNO9 NOMBRE MAX DE NOEUDS POUR LES ELEMENTS CONSIDERES 40 : C 41 : -INC SMELEME 42 : -INC SMTABLE 43 : -INC SMMODEL 44 : -INC SMCHAML 45 : -INC SMCHPOI 46 : -INC CCOPTIO 47 : -INC CCREEL 48 : -INC SMLREEL 49 : POINTEUR MLREE4.MLREEL,MLREE5.MLREEL,MLREE6.MLREEL 50 : *** pteurs sur les MAILLAGES 51 : POINTEUR IZCENT.MELEME,IELTFA.MELEME,IZLAC.MELEME,IFACEL.MELEME 52 : *** variable de travail pour sp 'ACCTAB' % table de lacher 'MTAB2' 53 : CHARACTER*20 MTYPI,CHARI,MTYPR,CHARR 54 : LOGICAL LOGRE 55 : *** 'IZPART' segment contenant position ref des particules au départ du lacher : 56 : *** 'NPART' nbre de particules lachées 57 : *** 'NLEPA' = n° global de l'element contenant ieme particule lachée, 58 : *** 'NUMPA' = n° ieme particule lachée, 'COORPA' = coords ieme particule lachée 59 : SEGMENT IZPART 60 : INTEGER NLEPA(NPART),NUMPA(NPART) 61 : REAL*8 COORPA(NDIM,NPART) 62 : ENDSEGMENT 63 : POINTEUR IZN.IZPART,IZL.IZPART 64 : *** 'IZCOU' segment content pas de tps des mailles traversées (selon critere convectif) 65 : SEGMENT IZCOU 66 : REAL*8 DTCO(NEL),COU 67 : ENDSEGMENT 68 : *** 'IZSH' segment de travail (pour passage réel -> réf) 69 : *** 'SHP' fonctions de forme et leurs dérivées 70 : *** 'SHY' fonctions de base et leurs dérivées 71 : *** 'XYZL(i, j)' = ieme coord reelle du jeme noeud d'1 elem donne 72 : SEGMENT IZSH 73 : REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) 74 : ENDSEGMENT 75 : *** 'IZVIT' segment décrivant les vitesses ou flux 76 : *** 'TEMTRA' = tableau des divers tps de palier (cas transitoire) 77 : *** 'NVIPT' = nbre de pas de tps (chgt) (cas transitoire) 78 : *** 'NBS' = nbre de sous-maillages 79 : *** 'IPUN' = pteurs sur les 'IZUN' 80 : *** 'IDUN(I)' = nbre d'éléments avant le sous-maillage I 81 : *** 'IPVPT(J)' = pteur sur les 'IZVPT' pour le pas de tps n° J (cas transitoire) 82 : SEGMENT IZVIT 83 : REAL*8 TEMTRA(NVIPT) 84 : INTEGER IPUN(NBS),IDUN(NBS),IPVPT(NVIPT),IFORML 85 : ENDSEGMENT 86 : *** 'IPARPO' segment permettant stockage resultats pdt 'avancee particule 87 : *** 'NPOS' = nbre d'elements traversés, 'NAPAR' = n° des élémts traversés 88 : *** 'NUMP' = inutile, 'CREF' = coord sauvegardées de chaque elemt traversee 89 : *** 'TPAR' = tps sauvegardé associé aux coords sauvegardées 90 : SEGMENT IPARPO 91 : INTEGER NAPAR(NPOS),NUMP(NPOS) 92 : REAL*8 CREF(NDIM,NPOS),TPAR(NPOS) 93 : ENDSEGMENT 94 : POINTEUR IZN3.IPARPO 95 : *** tableau contenant tps et coordonnées de sortie (et non sauvegarde) 96 : DIMENSION SORTIE(4) 97 : 98 : **************************************************************** 99 : C 100 : C 101 : SEGACT MTAB2 102 : SEGACT IZCENT 103 : NEL=IZCENT.NUM(/2) 104 : SEGINI IZCOU 104 : CALL INITD(DTCO,NEL,0.D0) 105 : MNO9=MPNO9 106 : SEGINI IZSH 107 : C 108 : C RECUPERATION D'ELEMENTS PROVENANT DE LA TABLE DE LACHER 109 : IVALI=1 110 : XVALI=0.D0 111 : IRETI=0 112 : IVALR=0 113 : XVALR=0.D0 114 : IRETR=0 115 : MTYPI='MOT ' 116 : C CHARI='CFL' 117 : C MTYPR='FLOTTANT' 118 : MTYPR=' ' 119 : CHARR=' ' 120 : COU=0.05D0 120 : CALL ACCTAB(MTAB2,MTYPI,IVALI,XVALI,'CFL',.TRUE.,IRETI, 121 : * MTYPR,IVALR,XVALR,CHARR,.TRUE.,IRETR) 122 : IF(IERR.NE.0)RETURN 123 : IF(MTYPR.NE.' ')THEN 124 : IF(MTYPR.NE.'FLOTTANT')THEN 125 : MOTERR(1:11)='CFL ' 126 : MOTERR(12:20)='FLOTTANT' 126 : CALL ERREUR(627) 127 : RETURN 128 : ENDIF 129 : COU=XVALR 130 : IF(COU.LT.1.D-8)THEN 131 : MOTERR(1:8)='CFL ' 132 : REAERR(1)=REAL(COU) 133 : REAERR(2)=1.E-8 133 : CALL ERREUR(41) 134 : RETURN 135 : ENDIF 136 : ENDIF 137 : IVALI=1 138 : XVALI=0.D0 139 : IRETI=0 140 : IVALR=0 141 : XVALR=0.D0 142 : IRETR=0 143 : MTYPI='MOT ' 144 : C MTYPR='FLOTTANT' 145 : MTYPR=' ' 146 : CHARR=' ' 146 : CALL ACCTAB(MTAB2,MTYPI,IVALI,XVALI,'DELTAT_SAUVE',.TRUE., 147 : * IRETI,MTYPR,IVALR,XVALR,CHARR,.TRUE.,IRETR) 148 : IF(MTYPR.NE.' ')THEN 149 : IF(MTYPR.NE.'FLOTTANT')THEN 150 : MOTERR(1:11)='DELTAT_SAUV' 151 : MOTERR(12:20)='FLOTTANT' 151 : CALL ERREUR(627) 152 : RETURN 153 : ENDIF 154 : ENDIF 155 : IF(IERR.NE.0)RETURN 156 : DELTAT=XVALR 157 : IF ((IICAL.EQ.3).AND.(DELTAT.GT.0.D0)) THEN 158 : JG=1 159 : SEGINI MLREE6 160 : MLREE6.PROG(1)=DELTAT 161 : NSAUV=0 162 : ENDIF 163 : IVALI=1 164 : XVALI=0.D0 165 : IRETI=0 166 : IVALR=0 167 : XVALR=0.D0 168 : IRETR=0 169 : MTYPI='MOT ' 170 : MTYPR='FLOTTANT' 171 : CHARR=' ' 171 : CALL ACCTAB(MTAB2,MTYPI,IVALI,XVALI,'TEMPS_LIMITE',.TRUE., 172 : * IRETI,MTYPR,IVALR,XVALR,CHARR,.TRUE.,IRETR) 173 : IF(IERR.NE.0)RETURN 174 : TMAX=XVALR 174 : CALL TRJCN7(TMAX,IZVIT) 175 : IF(IERR.NE.0)RETURN 176 : 177 : C LECTURE DU TABLEAU DES TEMPS DE LACHER 178 : IVALI=1 179 : XVALI=0.D0 180 : IRETI=0 181 : IVALR=0 182 : XVALR=0.D0 183 : IRETR=0 184 : MTYPI='MOT ' 185 : C CHARI='BORNES' 186 : MTYPR='LISTREEL' 187 : CHARR=' ' 187 : CALL ACCTAB(MTAB2,MTYPI,IVALI,XVALI,'TEMPS_LACHER',.TRUE.,IRETI, 188 : * MTYPR,IVALR,XVALR,CHARR,.TRUE.,IRETR) 189 : IF(IERR.NE.0)RETURN 190 : MLREE3=IRETR 191 : SEGACT MLREE3 192 : NLAC=MLREE3.PROG(/1) 193 : 194 : C DONNEES SUPLEMENTAIRES POUR CONVECTION_DIFFUSION 195 : IVALI=1 196 : XVALI=0.D0 197 : IRETI=0 198 : IVALR=0 199 : XVALR=0.D0 200 : IRETR=0 201 : MTYPI='MOT ' 202 : MTYPR=' ' 203 : CHARR=' ' 203 : CALL ACCTAB(MTAB2,MTYPI,IVALI,XVALI,'IMPERMEABLE',.TRUE.,IRETI, 204 : * MTYPR,IVALR,XVALR,CHARR,LOGRE,IRETR) 205 : IF(MTYPR.NE.' ')THEN 206 : IF(MTYPR.NE.'MAILLAGE')THEN 207 : MOTERR(1:11)='IMPERMEABLE' 208 : MOTERR(12:20)='MAILLAGE' 208 : CALL ERREUR(627) 209 : RETURN 210 : ENDIF 211 : ENDIF 212 : IF(IERR.NE.0)RETURN 213 : IZFACE=IRETR 214 : IVALI=1 215 : XVALI=0.D0 216 : IRETI=0 217 : IVALR=0 218 : XVALR=0.D0 219 : IRETR=0 220 : MTYPI='MOT ' 221 : MTYPR=' ' 222 : CHARR=' ' 223 : *** recuperation de la liste des tps de sauvegarde des resultats 223 : CALL ACCTAB(MTAB2,MTYPI,IVALI,XVALI,'TEMPS_SAUVES',.TRUE., 224 : $ IRETI,MTYPR,IVALR,XVALR,CHARR,.TRUE.,IRETR) 225 : IF (MTYPR.NE.' ') THEN 226 : IF (MTYPR.NE.'LISTREEL') THEN 227 : MOTERR(1:12)='TEMPS_SAUVES' 228 : MOTERR(13:20)='LISTREEL' 228 : CALL ERREUR(627) 229 : RETURN 230 : ENDIF 231 : ENDIF 232 : IF (IERR.NE.0) RETURN 233 : IZSAUV=IRETR 234 : IF (IZSAUV.GT.0) THEN 235 : IF (DELTAT.GT.0.D0) THEN 236 : MOTERR(1:8)= 'DELTAT_S' 237 : MOTERR(9:16)='TEMPS_SA' 237 : CALL ERREUR(135) 238 : RETURN 239 : ELSE 240 : MLREE6=IZSAUV 241 : SEGACT MLREE6 242 : NSAUV=MLREE6.PROG(/1) 243 : ENDIF 244 : ELSE 245 : IF (DELTAT.LE.0.D0) NSAUV=-1 246 : ENDIF 247 : 248 : C ON ACTIVE LES TABLEAUX 249 : *** activation du chamelem de dispersion si nécessaire 250 : IF (IZDISP.GT.0) THEN 251 : MCHEL1=IZDISP 252 : SEGACT MCHEL1 253 : MCHAML=MCHEL1.ICHAML(1) 254 : SEGACT MCHAML 255 : MELVA1=IELVAL(1) 256 : SEGACT MELVA1 257 : MELVA2=IELVAL(2) 258 : SEGACT MELVA2 259 : SEGDES MCHEL1,MCHAML 260 : ENDIF 261 : *** activation du chamelem de diffusivité si nécessaire 262 : IF (IZDIFF.GT.0) THEN 263 : MCHEL1=IZDIFF 264 : SEGACT MCHEL1 265 : MCHAML=MCHEL1.ICHAML(1) 266 : SEGACT MCHAML 267 : MELVA3=IELVAL(1) 268 : SEGACT MELVA3 269 : SEGDES MCHEL1,MCHAML 270 : ENDIF 271 : *** activation du maillage des faces impermeables si nécessaire 272 : IF (IZFACE.GT.0) THEN 273 : IPT9=IZFACE 274 : SEGACT IPT9 275 : ENDIF 276 : 277 : C PRISE EN COMPTE DE LA POROSITE 278 : C 279 : IF(IZPOR.NE.0)THEN 279 : CALL TRJPOR(IZPOR,IZVIT,MELEME) 280 : IF(IERR.NE.0)RETURN 281 : ENDIF 282 : C 283 : C STRUCTURES DE RESULTATS 284 : C 285 : L1=5 286 : N1=1 287 : N3=6 288 : SEGINI MCHELM 289 : TITCHE(1:5)='NOEUD' 290 : SEGINI MMODEL 291 : 292 : C ------------------------------ 293 : C BOUCLE SUR LES TEMPS DE LACHER 294 : C ------------------------------ 295 : 296 : EPSILO=10*SQRT(XZPREC) 297 : NPART0=0 298 : DO 80 KLAC=1,NLAC 299 : 300 : TMIN=MLREE3.PROG(KLAC) 300 : CALL TRJCN7(TMIN,IZVIT) 300 : IF(ABS(DELTAT).GT.1.D-15)CALL TRJCN4(TMIN,TMAX,DELTAT) 301 : 302 : C LECTURE D UN MAILLAGE DE LACHER 303 : IVALI=KLAC 304 : XVALI=0.D0 305 : IRETI=0 306 : IVALR=0 307 : XVALR=0.D0 308 : IRETR=0 309 : MTYPI='ENTIER ' 310 : CHARI=' ' 311 : MTYPR='MAILLAGE' 312 : CHARR='MAILLAGE' 312 : CALL ACCTAB(MTAB2,MTYPI,IVALI,XVALI,CHARI,.TRUE.,IRETI, 313 : * MTYPR,IVALR,XVALR,CHARR,.TRUE.,IRETR) 314 : IF(IERR.NE.0)RETURN 315 : IZLAC=IRETR 315 : CALL TRJCN2(IZLAC) 316 : IF(IERR.NE.0)RETURN 317 : C 317 : CALL TRJINI(IZPART,IZCOU,IZLAC,MELEME,IZVIT,IZCENT,IELTFA,IZSH, 318 : * TMIN) 319 : NPARTI=NUMPA(/1) 320 : N1=NPART0+NPARTI 321 : SEGADJ MCHELM 322 : SEGADJ MMODEL 323 : 324 : C ........................................ 325 : C 326 : C BOUCLE SUR LES PARTICULES DU LACHER KLAC 327 : C ........................................ 328 : 329 : DO 50 IPART=1,NPARTI 330 : IIPOS=0 331 : IREP=0 332 : NDIM=COORPA(/1) 333 : TMINI=TMIN 334 : DT1=0.D0 335 : IDD1=1 336 : 53 CONTINUE 337 : 338 : IF(IICAL.EQ.1) THEN 338 : CALL TRJAVA(IZVIT,IZPART,IZN3,IPART,TMINI,TMAX,IZCOU,MELEME, 339 : * IELTFA,IZCENT,IFACEL,DELTAT,IZSH,IEROR,DT1) 340 : IF(IERR.NE.0)RETURN 340 : CALL TRJRAL(IZN3,IPARPO,MELEME,IZSH) 341 : SEGSUP IZN3 342 : ELSEIF(IICAL.EQ.2) THEN 343 : C convection seule avec méthode analytique 343 : CALL TRJAV2(IZVIT,IZPART,IZN3,IPART,TMINI,TMAX,MELEME, 344 : * IELTFA,IZCENT,IFACEL,DELTAT,IZSH) 345 : IF(IERR.NE.0)RETURN 345 : CALL TRJRAL(IZN3,IPARPO,MELEME,IZSH) 346 : SEGSUP IZN3 347 : ELSE 348 : C** convection + dispersion + diffusion 348 : CALL TRJAV3(EPSILO,NSAUV,MLREE6,MELVA1,MELVA2, 349 : $ MELVA3,IPT9,IZVIT,IZPART,IPART,TMINI,TMAX,MELEME, 350 : $ IELTFA,IZCENT,IFACEL,IPARPO,IZSH,SORTIE) 351 : IF (IERR.NE.0) RETURN 352 : C** cas ou particule i est eliminee du probleme 353 : IF (SORTIE(1).LT.0.D0) THEN 354 : NPART0=NPART0-1 355 : N1=NPART0+NPARTI 356 : SEGADJ MCHELM 357 : SEGADJ MMODEL 358 : GOTO 50 359 : ENDIF 360 : ENDIF 361 : 362 : NPOS=NAPAR(/1) 363 : IF(IIMPI.GT.0)THEN 364 : WRITE(6,*)' PARTICULE ',IPART,NPOS 365 : WRITE(6,*)NAPAR(1),NUMP(1),(CREF(J,1),J=1,NDIM),TPAR(1) 366 : I2=NPOS 367 : IF(IIMPI.GT.1)I2=2 368 : DO 10 I=I2,NPOS 369 : WRITE(6,*)NAPAR(I),NUMP(I),(CREF(J,I),J=1,NDIM),TPAR(I) 370 : 10 CONTINUE 371 : ENDIF 372 : 373 : C CAS OU TOUS LES ELEMENTS N'ONT PAS LA MEME ORIENTATION 374 : C ON MET LES RESULTATS DANS LE MEME SEGMENT 375 : NPOS1=NPOS 376 : IF(IREP.NE.0)THEN 377 : IZN3=IREP 378 : NPOS=IIPOS+NPOS+1-IDD1 379 : SEGADJ IZN3 380 : DO 54 II=IDD1,NPOS1 381 : II1=II+1-IDD1 382 : IZN3.CREF(1,IIPOS+II1)=CREF(1,II) 383 : IZN3.CREF(2,IIPOS+II1)=CREF(2,II) 384 : IZN3.TPAR(IIPOS+II1)=TPAR(II) 385 : 54 CONTINUE 386 : IF(NDIM.EQ.3)THEN 387 : DO 55 II=1,NPOS1 388 : IZN3.CREF(3,IIPOS+II+1-IDD1)=CREF(3,II) 389 : 55 CONTINUE 390 : ENDIF 391 : SEGSUP IPARPO 392 : IPARPO=IZN3 393 : ENDIF 394 : C 395 : IF(IEROR.EQ.1)THEN 396 : TMINI=TPAR(NPOS) 397 : TTEMP=TPAR(NPOS) 398 : NPART=1 399 : SEGINI IZN,IZL 400 : DO 56 I=1,NDIM 401 : IZL.COORPA(I,1)=CREF(I,NPOS) 402 : 56 CONTINUE 403 : 404 : IZN.NLEPA(1)=0 405 : IZN.NUMPA(1)=1 406 : IZL.NLEPA(1)=0 407 : IZL.NUMPA(1)=1 408 : IIPOS=IIPOS+NPOS1-IDD1 408 : CALL TRJPEL(IZL,IZN,MELEME,IZVIT,IZCOU,IZCENT,IELTFA,IZSH, 409 : & TTEMP) 410 : IF(IZN.NLEPA(1).NE.0)THEN 411 : IREP=IPARPO 412 : DO 57 I=1,NDIM 413 : COORPA(I,IPART)=IZN.COORPA(I,1) 414 : 57 CONTINUE 415 : NLEPA(IPART)=IZN.NLEPA(1) 416 : DT1=0.D0 417 : IDD1=1 418 : C CAS D'UN PAS CONSTANT, ON RECALE LE CALCUL AVEC DT1 419 : IF((ABS(DELTAT).GT.1.D-15).AND. 420 : * ((TPAR(NPOS)-TPAR(NPOS-1)).NE.DELTAT))THEN 421 : DT1=TPAR(NPOS)-TPAR(NPOS-1) 422 : IDD1=2 423 : ENDIF 424 : SEGSUP IZL,IZN 425 : GO TO 53 426 : ENDIF 427 : SEGSUP IZL,IZN 428 : ENDIF 429 : C 430 : IVALI=IPART+NPART0 430 : CALL TRJCHA(IPARPO,MCHELM,MMODEL,IVALI) 431 : SEGSUP IPARPO 432 : 50 CONTINUE 433 : 434 : C ............................................ 435 : C 436 : C FIN BOUCLE SUR LES PARTICULES DU LACHER KLAC 437 : C ............................................ 438 : 439 : NPART0=NPART0+NPARTI 440 : SEGSUP IZPART 441 : 80 CONTINUE 442 : 443 : C ---------------------------------- 444 : C FIN BOUCLE SUR LES TEMPS DE LACHER 445 : C ---------------------------------- 446 : 447 : SEGDES MLREE3 448 : SEGDES MTAB2 449 : SEGSUP IZCOU 450 : SEGSUP IZSH 451 : C 452 : C 453 : RETURN 454 : END 455 : 456 : 457 : 458 :
© Cast3M 2003 - All rights reserved.
Disclaimer