1 : C STRATE SOURCE BP208322 13/04/11 21:16:13 7754 2 : subroutine strate(ipsolu,nmodma,nbmode,frshif,nmod,ifini 3 : $ ,shifti,idist,nbonZ,icycle,Northo 4 : $ ,IPRIGI,IPMASS,IPKW2M,nvp0,IFLU,incvrai) 5 : * 6 : ********************************************************************** 7 : * STRATE 8 : * 9 : * fonction : petit systeme "expert" pour definir la meilleure strategie 10 : * pour poursuivre la recherche 11 : * 12 : * creation : chat 11.2009 13 : * modifs : bp 02.2010 : ajout du tri des modes par freq croissante (si 14 : * numero mode errone) 15 : * bp 11.2010 : modif stratégie alternatives, verification de la 16 : * completude du spectre, completion de la base si 17 : * mode multiple a cheval sur le spectre 18 : * bp 12.2010 : ajout tri des modes selon le shift (ordshi.eso) 19 : * bp 05.2011 : ajout cas K non-definie positive (nvp0) 20 : * 21 : * shifti = shift initial (on le range dans shifin du segment idist) 22 : * frshif = shift courant (en entrée) 23 : * = shift proposé pour le prochain cycle (en sortie) 24 : * 25 : ********************************************************************** 26 : * 27 : IMPLICIT REAL*8(A-H,O-Z) 28 : IMPLICIT INTEGER (I-N) 29 : * 30 : -INC CCREEL 31 : -INC CCOPTIO 32 : -INC SMSOLUT 33 : *-INC SMLREEL 34 : -INC SMLENTI 35 : -INC SMELEME 36 : * 37 : * idist = reunion des ipsolu des differents cycles 38 : segment idist 39 : integer iexis(jg),ipomod(jg),immode(jg),ipve(jg),imil,ienc 40 : integer ialter 41 : real*8 dist(jg),frequ(jg),shifin,dismin 42 : endsegment 43 : pointeur mmod2.mmode 44 : pointeur mso3.msolut,msole3.msolen,mmod3.mmode,mmod4.mmode 45 : c 46 : logical test1,test2 47 : c 48 : REAL*8 FRSHIF,W2SHIF,DEUXPI 49 : PARAMETER (DEUXPI = (2.D0*XPI)) 50 : c 51 : c 52 : c quelques initialisations et informations utiles 53 : ifini=0 54 : inouv=0 55 : * nvp1 = numero de la 1ere valeur propre non negative 56 : nvp1=nvp0+1 57 : c IPKW2M a deja été factorisé par Lanczos ou simul1 57 : CALL DIAGN1(IPKW2M,imoshi) 58 : if(iimpi.ge.5) write(IOIMP,*) 59 : & 'strate: IPKW2M,frshif,imoshi=',IPKW2M,frshif,imoshi 60 : *ajout bp 06/01/2012: 60 : CALL DIAGN1(IPMASS,nvp0M) 61 : * correction pour elements fluides (inconnue PI mise à 0 via INITFL) 62 : c nvp0M = nvp0M - IFLU 63 : * bp 10/01/2012: nvp0M et NEMSM ne semblent pas bien calculés ... 64 : * (resultats dependant machine -> cf. dyna7.dgibi) 65 : * on propose la solution qui suppose que nvp0M si M est LIQUIDE 66 : if(IFLU.gt.0) nvp0M=0 67 : if(iimpi.ge.5) write(IOIMP,*) 68 : & 'strate: IPMASS, nvp0M',IPMASS,' ',nvp0M 69 : if(nvp0M.ne.0) then 70 : if (frshif.gt.0.D0) then 71 : imoshi=nvp0M+imoshi 72 : elseif (frshif.lt.0.D0) then 73 : imoshi=nvp0M-imoshi 74 : else 75 : imoshi=nvp0M 76 : endif 77 : endif 78 : if(iimpi.ge.5) write(IOIMP,*) 79 : & 'strate: ... => imoshi=',imoshi 80 : * 81 : *---- recup ou creation de idist -------------------------------------* 82 : * 83 : if (idist.ne.0) then 84 : segact idist*mod 85 : jg=iexis(/1) 86 : c iplus=jg 87 : c bp (10/12/2010) : idist existait deja, mais les numeros de mode 88 : c attribués sont-ils coherents avec le shift actuel ? 88 : call ordshi(idist,imoshi,frshif) 89 : IF(IERR.NE.0) RETURN 90 : 91 : else 92 : jg=0 93 : c iplus=0 94 : segini idist 95 : shifin=shifti 96 : ialter=0 97 : c bp (13/05/2011) : on initialise imil avec imoshi du 1er cycle (=shift initial) 98 : imil=imoshi+1 99 : endif 100 : * 101 : *---- travail preliminaire sur ipsolu --------------------------------* 102 : * 103 : mso3=ipsolu 104 : if(mso3.eq.0) goto 7 105 : * rem bp : ligne ci dessus ne doit jamais arriver? (cf. simul3 et 7) 106 : jg=0 107 : segact mso3 108 : imodgr=0 109 : msole3=mso3.msolis(4) 110 : msole2=mso3.msolis(5) 111 : ipt3=mso3.msolis(3) 112 : segact msole3,ipt3,msole2 113 : do iou=1,msole3.isolen(/1) 114 : mmode=msole3.isolen(iou) 115 : segact mmode 116 : ia=immodd(1) 117 : fr=fmmodd(1) 118 : fr2 = fr * (abs(fr)) 119 : shifin2=shifin * (abs(shifin)) 120 : if (dismin.eq.0.d0) then 121 : dismin=abs(shifin2 - fr2)*10000.d0 122 : endif 123 : if(ia.gt.jg) jg=ia 124 : enddo 125 : if(jg.gt.iexis(/1)) segadj,idist 126 : * on a parfois besoin d'aller chercher dans (imoshi+1) 127 : jg=iexis(/1) 128 : if(jg.lt.(imoshi+1)) then 129 : jg=imoshi+1 130 : segadj,idist 131 : endif 132 : * 133 : *---- remplissage de idist par la nouvelle solution ------------------* 134 : * 135 : ai=frshif*(abs(frshif)) 136 : c ----boucle sur les modes de ce cycle 137 : c bp: avec simul3, ils sont tous nouveaux, reste a les placer 138 : do iou=1,msole3.isolen(/1) 139 : 140 : mmode=msole3.isolen(iou) 141 : imo=immodd(1) 142 : fre=fmmodd(1) 143 : c -on suggere que ce mode de frequence fre a le numero imo 144 : if (iexis(imo).eq.0) then 145 : c imo est effectivement libre -> c'est un bon candidat pour jmo 146 : jmo = imo 147 : else 148 : c imo contient deja un mode -> pb dans la numerotation 149 : c -> on cherche une place jmo libre a proximite 150 : c (avec imo+idec1 en priorite sur imo+idec2) 151 : if (fre.lt.frequ(imo)) then 152 : idec1 = -1 153 : idec2 = 1 154 : else 155 : idec1 = 1 156 : idec2 = -1 157 : endif 158 : jmo = 0 159 : jmo1 = imo 160 : jmo2 = imo 161 : ndec = max((jg-imo),imo) 162 : do idec=1,ndec 163 : jmo1 = jmo1 + idec1 164 : if (jmo1.gt.0.and.jmo1.le.jg) then 165 : if (iexis(jmo1).eq.0) then 166 : jmo = jmo1 167 : goto 1 168 : endif 169 : c write(IOIMP,*) ' ',jmo1,' occupée',iexis(jmo1) 170 : endif 171 : jmo2 = jmo2 + idec2 172 : if (jmo2.gt.0.and.jmo2.le.jg) then 173 : if (iexis(jmo2).eq.0) then 174 : jmo = jmo2 175 : goto 1 176 : endif 177 : c write(IOIMP,*) ' ',jmo1,' occupée',iexis(jmo1) 178 : endif 179 : enddo 180 : c on n a pas trouve de place libre 181 : jg = jg + 1 182 : segadj,idist 183 : jmo = jg 184 : endif 185 : 1 continue 186 : c => la place jmo est libre. 187 : c write(IOIMP,*) ' iou,imo,jmo est libre,jg=',iou,imo,jmo,jg 188 : 189 : c -avant de remplir idist(jmo), 190 : c il reste a verifier que les freq sont bien ordonnées : 191 : c par construction elles le sont toutes, 192 : c sauf fre qu on cherche a insérer 193 : 5 continue 194 : c il peut arriver que finalement il faille augmenter la taille de idist 195 : if (jmo.gt.jg) then 196 : jg = jmo 197 : segadj,idist 198 : endif 199 : if (jmo.le.0.or.jmo.gt.jg) then 200 : write(IOIMP,*)'ERREUR dans STRATE : ' 201 : write(IOIMP,*)'plage des numeros de mode incorrecte',jmo,jg 202 : write(IOIMP,*)'Contactez votre support' 203 : stop 204 : endif 205 : jmo1 = jmo 206 : jmo2 = jmo 207 : c recherche des plus proches voisins 208 : 2 continue 209 : jmo1 = jmo1 - 1 210 : if(jmo1.le.0) goto 3 211 : if(iexis(jmo1).eq.0) goto 2 212 : 3 continue 213 : jmo2 = jmo2 + 1 214 : if(jmo2.gt.jg) goto 4 215 : if(iexis(jmo2).eq.0) goto 3 216 : 4 continue 217 : c write(IOIMP,*) 'jmo1,jmo2=',jmo1,jmo2 218 : c il faut que frequ(jmo1) < fre < frequ(jmo2) 219 : c pour attribuer frequ(jmo) = fre 220 : if (jmo1.le.0) then 221 : test1 = .true. 222 : else 223 : test1 = (frequ(jmo1).le.fre) 224 : endif 225 : if (jmo2.gt.jg) then 226 : test2 = .true. 227 : else 228 : test2 = (fre.le.frequ(jmo2)) 229 : endif 230 : if (.not.test1) then 231 : c on met jmo1 -> jmo 232 : c write(IOIMP,*) ' on met ',jmo1,' -> ',jmo 233 : c iexis(jmo) = 1 234 : iexis(jmo) = iexis(jmo1) 235 : frequ(jmo) = frequ(jmo1) 236 : dist (jmo) = dist (jmo1) 237 : ipve (jmo) = ipve (jmo1) 238 : ipomod(jmo)= ipomod(jmo1) 239 : immode(jmo)= immode(jmo1) 240 : iexis(jmo1) = 0 241 : frequ(jmo1) = 0.D0 242 : c jmo = jmo1 les 2 doivent etre corrects (et svt equivalentes) 243 : jmo = jmo - 1 244 : * rem (bp): on pourrait aussi décider en fonction de dist(jmo) et de 245 : * dist(jmo1) si on ne met pas plutot jmo1 -> jmo1+1 et jmo=jmo1 246 : * mais on risque de créer des faux trous ... 247 : goto 5 248 : elseif (.not.test2) then 249 : c on met jmo2 -> jmo 250 : c write(IOIMP,*) ' on met ',jmo2,' -> ',jmo 251 : c iexis(jmo) = 1 252 : iexis(jmo) = iexis(jmo2) 253 : frequ(jmo) = frequ(jmo2) 254 : dist (jmo) = dist (jmo2) 255 : ipve (jmo) = ipve (jmo2) 256 : ipomod(jmo)= ipomod(jmo2) 257 : immode(jmo)= immode(jmo2) 258 : iexis(jmo2) = 0 259 : frequ(jmo2) = 0.D0 260 : c jmo = jmo2 les 2 doivent etre corrects (et svt equivalentes) 261 : jmo = jmo + 1 262 : goto 5 263 : endif 264 : 265 : c -on peut enfin remplir idist(jmo) 266 : c write(IOIMP,*) 'on remplit idist (',jmo 267 : inouv=inouv+1 268 : iexis(jmo) = icycle 269 : frequ(jmo)=fre 270 : fr2 = fre * (abs(fre)) 271 : dist(jmo)=abs(ai-fr2) 272 : if(iimpi.ge.6) write(IOIMP,*) 273 : & '--> iou,fre,imo,jmo',iou,fre,imo,jmo 274 : if (dismin.gt.dist(jmo)) then 275 : c if(iplus.eq.0) imil=jmo 276 : c cette ligne ci dessus conduit a des erreurs si modes multiples 277 : dismin=dist(jmo) 278 : endif 279 : ipve(jmo)=msole2.isolen(iou) 280 : ipomod(jmo)=ipt3.num(1,iou) 281 : immode(jmo)=mmode 282 : 283 : segdes mmode 284 : 285 : enddo 286 : c ----fin de la boucle sur les modes de ce cycle 287 : segdes msole2,msole3,mso3 288 : 289 : c on trie les modes de idist en fonction du shift actuel 289 : call ordshi(idist,imoshi,frshif) 290 : IF(IERR.NE.0) RETURN 291 : 292 : 7 continue 293 : * 294 : * as-t-on progressé? -------------------* 295 : if(iimpi.ge.6) 296 : & write(IOIMP,*) ' nombre de nouveaux modes' ,inouv,'/',nbonZ 297 : nbonZ=inouv 298 : 299 : if (inouv.eq.0) then 300 : * -on n a pas de nouveaux modes 301 : ialter=ialter+1 302 : write(IOIMP,*) ' strategie alternative : ',ialter 303 : c autorisées 4 fois de suite 304 : c if (ialter.ge.3) then 305 : if (ialter.ge.5) then 306 : write(IOIMP,*) 'ERREUR dans STRATE :', 307 : & ' nombre d alternatives maxi atteinte' 308 : write(IOIMP,*) ' Changez de décalage initial,', 309 : & ' Utilisez une machine plus puissante', 310 : & ' ou Contactez votre support...' 310 : call erreur(6) 311 : return 312 : endif 313 : c -0eme chance (bp 13/05/2011): on a jamais rien trouvé 314 : c on choisit arbitrairement un shift 315 : if (jg.eq.0) then 316 : shipro = 10.D0**(icycle-2) 317 : write(IOIMP,*) ' On choisit arbitrairement le shift : ', 318 : & frshif,' devient ',shipro 319 : goto 2000 320 : endif 321 : c -1ere chance : on est allé trop loin et on a rien trouvé : 322 : c on shifte en direction du shift initial 323 : c if (((imoshi+1).gt.jg).or.(imoshi.eq.0)) then 324 : if ((imoshi.ge.incvrai).or.(imoshi.eq.0)) then 325 : if(iimpi.ge.6) write(IOIMP,*) 'jg,incvrai,nbmode,imoshi=' 326 : & ,jg,incvrai,nbmode,imoshi 327 : shipro = shifti + ((frshif - shifti) / 2.D0) 328 : write(IOIMP,*) ' on se rapproche du shift initial : ', 329 : & frshif,' devient ',shipro 330 : goto 2000 331 : endif 332 : c -2eme chance : on augmente le nombre de vecteurs sans shifter 333 : if (ienc.lt.nmodma) then 334 : c nmodma = int(1.5 * ienc) + 4 335 : c ienc = ienc + max(4,(nmodma/2)) 336 : ienc = ienc + max(4,nmodma) 337 : ienc = min(ienc,2*nmodma) 338 : Northo=ienc 339 : shipro=frshif 340 : write(IOIMP,*)' augmentation du nombre de mode cherché',ienc 341 : c else 342 : endif 343 : if (ialter.ge.2) then 344 : c -3eme chance : on shifte en direction du mode qui semble mal placé 345 : if (iexis(imoshi).ne.0.and.iexis(imoshi+1).eq.0) then 346 : shipro = (2.*frequ(imoshi)) - frshif 347 : Northo=ienc 348 : write(IOIMP,*) ' on shifte en direction du mode ',imoshi, 349 : & ' mal placé : ',frshif,' devient ',shipro 350 : c goto 2000 351 : elseif (iexis(imoshi).eq.0.and.iexis(imoshi+1).ne.0) then 352 : shipro = (2.*frequ(imoshi+1)) - frshif 353 : Northo=ienc 354 : write(IOIMP,*)' on shifte en direction du mode ',(imoshi+1), 355 : & ' mal placé : ',frshif,' devient ',shipro 356 : c goto 2000 357 : else 358 : if(iimpi.ge.6) 359 : & write(IOIMP,*) ' jg,nbmode,imoshi=',jg,nbmode,imoshi 360 : Northo=ienc 361 : if ((frshif.ne.shifti).and.(ialter.lt.4)) then 362 : shipro = shifti + ((frshif - shifti) / 2.D0) 363 : write(IOIMP,*) ' on se rapproche du shift initial : ', 364 : & frshif,' devient ',shipro 365 : else 366 : ifin=jg+1 367 : 8 ifin=ifin-1 368 : if (ifin.ne.0) then 369 : if(iexis(ifin).eq.0) goto 8 370 : shipro = frequ(ifin) 371 : else 372 : shipro = shifti 373 : endif 374 : shipro = max(shipro,frshif) 375 : shipro = 4.D0 * shipro 376 : shipro = max(shipro,100.D0) 377 : write(IOIMP,*) ' on perturbe arbitrairement le shift : ', 378 : & frshif,' devient ',shipro 379 : endif 380 : c goto 2000 381 : endif 382 : endif 383 : goto 2000 384 : 385 : else 386 : c on a eu des nouveaux modes => on reinitialise ialter 387 : ialter=0 388 : endif 389 : * ----fin du cas inouv=0 (pas de nouveau modes) -----* 390 : 391 : * 392 : * on cherche la longueur entourant imil 393 : ibas=0 394 : if(imil.eq.1) goto 9 395 : do ibas=imil-1,1,-1 396 : if( iexis(ibas).eq.0) goto 9 397 : enddo 398 : ibas= 0 399 : 9 continue 400 : ibas=ibas+1 401 : c bp : mise en commentaire le 06/06/2011 402 : c if (imil.lt.nvp1) then 403 : c write(IOIMP,*) 'Le shift proposé ne permet pas ', 404 : c & 'd obtenir des valeurs propres positives !' 405 : c write(IOIMP,*) 'ibas,imil,nvp1=',ibas,imil,nvp1 406 : c write(IOIMP,*) 'iexis=',(iexis(iou),iou=1,jg) 406 : c call erreur(6) 407 : c return 408 : c endif 409 : c bp (13/05/2011): eventuellement (si .not.limage cf. simul1), on limite 410 : c ibas pour ne pas recuperer les valeurs propres negatives 411 : ibas0=ibas 412 : ibas=max(ibas,nvp1) 413 : *bp: on rajoute cette ligne pour le cas ou l on doit trouver ibas > imil 414 : if(iexis(ibas).eq.0) goto 9 415 : imil1=max(imil,ibas) 416 : do ihaut=imil1,jg 417 : if(iexis(ihaut).eq.0) goto 10 418 : enddo 419 : ihaut=jg+1 420 : 10 continue 421 : ihaut = ihaut-1 422 : if(iexis(ihaut).eq.0) goto 10 423 : ilong= ihaut - ibas + 1 424 : if(iimpi.ge.6) write(IOIMP,*) 'numéros trouvés: ' 425 : & ,'[ibas(imil)ihaut]ilong,jg=',ibas,imil,ihaut,ilong,jg 426 : * 427 : * 428 : * --- on a, peut etre, suffisament de mode -------------------------* 429 : if (ilong.ge.nbmode) then 430 : * 431 : * -on partait de zero (ou presque) : pas besoin de réfléchir 432 : if (ibas.eq.1.and.shifin.le.frequ(1)) then 433 : idep=1 434 : ifin=nbmode 435 : goto 1000 436 : endif 437 : 438 : * -on ne peut descendre plus bas : 439 : if (ibas.eq.1) then 440 : dbas= dist(ibas) 441 : dhaut=dist(ihaut) 442 : if(iimpi.ge.6) write(IOIMP,*) 'ibas,dbas,ihaut,dhaut' 443 : & ,ibas,dbas,ihaut,dhaut 444 : 445 : * on recupere la plage [idep:ifin] la mieux centree 446 : if (dhaut.ge.dbas) then 447 : idep=ibas 448 : ifin=ihaut 449 : 13 ilong=ifin-idep+1 450 : if(ilong.eq.nbmode) goto 1000 451 : if (dist(idep).gt.dist(ifin)) then 452 : idep=idep+1 453 : else 454 : ifin=ifin-1 455 : endif 456 : goto 13 457 : 458 : * a-t-on suffisament de frequence vers le haut? 459 : else 460 : * on regarde combien de frequence sont bien centrées sur imil 461 : idep=ibas-1 462 : ifin=ihaut 463 : 11 idep=idep+1 464 : if(dist(idep).gt.dhaut) goto 11 465 : * bp2011: en theorie il faudrait garder idep mais pour les petits sytemes, 466 : * idep-1 est mieux adapté sinon on enleve trop de modes, d'ou : 467 : idep=idep-1 468 : ilong=ifin-idep+1 469 : if (ilong.ge.nbmode) then 470 : * c'est ok il faut selectionner la plage [idep:ifin] la mieux centree 471 : 14 ilong=ifin-idep+1 472 : if(ilong.eq.nbmode) goto 1000 473 : if (dist(idep).gt.dist(ifin)) then 474 : idep=idep+1 475 : else 476 : ifin=ifin-1 477 : endif 478 : goto 14 479 : else 480 : * combien au max doit-on trouver de modes supplementaires 481 : iveut=nbmode-ilong 482 : cbp ienc=min(int(iveut*1.1),nmodma) 483 : ienc=min(int(iveut*1.05),(nmodma-2)) 484 : ienc=max(ienc,1) 485 : icle=1 485 : call queshi(shifti,ienc,ifin,frequ,nmodma,shipro,iexis 486 : $ ,icle,iexis(/1),Northo) 487 : goto 2000 488 : endif 489 : endif 490 : 491 : * ici on a ibas.ne.1 et ilong > nbmode 492 : else 493 : * on regarde combien de frequence sont bien centrees sur imil 494 : * bp 2012-09-25 : ajout cas particuliers nbmode = 1 ou 2 495 : if(nbmode.eq.1) then 496 : idep=imil 497 : ifin=imil 498 : if(iimpi.ge.6) write(IOIMP,*)'nbmode=1 et imil=',imil 499 : goto 1000 500 : else 501 : if (dist(imil-1).lt.dist(imil+1).and.imil.gt.nvp1) then 502 : idep=imil-1 503 : ifin=imil 504 : else 505 : idep=imil 506 : ifin=imil+1 507 : endif 508 : if(nbmode.eq.2) then 509 : if(iimpi.ge.6) write(IOIMP,*)'nbmode=2 et imil=',imil 510 : & ,' plage: [idep - ifin]=',idep,' - ',ifin 511 : if(idep.lt.ibas.or.ifin.gt.ihaut) goto 12 512 : goto 1000 513 : endif 514 : endif 515 : do iou=1,nbmode-2 516 : if (dist(idep).le.dist(ifin)) then 517 : * sans depasser la limite fixee eventuellement par nvp1 518 : if (idep.gt.nvp1) then 519 : idep=idep-1 520 : else 521 : ifin=ifin+1 522 : endif 523 : else 524 : ifin=ifin+1 525 : endif 526 : if(idep.lt.ibas.or.ifin.gt.ihaut) goto 12 527 : enddo 528 : * selection de la plage [idep:ifin] de dim nbmode effectuee! 529 : if(iimpi.ge.6) 530 : &write(IOIMP,*) 'plage: [idep - ifin]=',idep,' - ',ifin 531 : goto 1000 532 : 12 continue 533 : * on a atteint la limite de [ibas:ihaut] sans obtenir une plage [idep:ifin] de dim nbmode 534 : if(iimpi.ge.6) 535 : &write(IOIMP,*) 'numéros centrés: [idep - ifin]=',idep,' - ',ifin 536 : * dans quel sens faut-il chercher? 537 : icle=1 538 : icher=ihaut 539 : if (idep.lt.ibas) then 540 : icle=2 541 : icher=ibas 542 : endif 543 : iveut= nbmode-ifin+idep 544 : cbp ienc=min(int(iveut*1.1),nmodma) 545 : ienc=min(int(iveut*1.05),(nmodma-2)) 546 : ienc=max(ienc,1) 547 : c write(IOIMP,*) 'iveut,ienc,icher,icle=', 548 : c $ iveut,ienc,icher,icle 548 : call queshi(shifti,ienc,icher,frequ,nmodma,shipro,iexis, 549 : $ icle,iexis(/1),Northo) 550 : goto 2000 551 : endif 552 : * 553 : * --- on a pas assez de mode de maniere certaine -------------------* 554 : else 555 : * maintenant ilong n'est plus assez grand 556 : * dans quelle direction faut-il chercher ? 557 : icle=1 558 : icher=ihaut 559 : if (dist(ibas).lt.dist(ihaut).and.ibas.ne.nvp1) then 560 : icle=2 561 : icher=ibas 562 : endif 563 : iveut= nbmode-ilong 564 : cbp ienc=min(int(iveut*1.1),nmodma) 565 : ienc=min(int(iveut*1.05),(nmodma-2)) 566 : ienc=max(ienc,1) 567 : c write(IOIMP,*) 'iveut,ienc,icher,icle=', 568 : c $ iveut,ienc,icher,icle 568 : call queshi(shifti,ienc,icher,frequ,nmodma,shipro,iexis, 569 : $ icle,iexis(/1),Northo) 570 : ccc $ icle,iexis(/1),fpet,fgra) 571 : goto 2000 572 : endif 573 : * ------------------------------------------------------------------* 574 : 575 : 576 : 1000 continue 577 : *---- on pense avoir gagné ! 578 : 579 : *---- on commence par rajouter les eventuels modes doubles 580 : c pour ne pas gener lors du controle du spectre 581 : c en bas 582 : if (idep.gt.ibas) then 583 : idep0=idep 584 : fref = max(1.D0,abs(frequ(idep0))) 585 : fref = max(fref,abs(shifti)) 586 : do ia=(idep0-1),ibas,-1 587 : c dfrq = (frequ(idep0) - frequ(ia)) / max(1.D0,abs(frequ(idep0))) 588 : * bp 13/01/2012 : si on est dans la plage des modes de corps rigides 589 : * on les inclut par defaut, sinon on va avoir des souci lors de la 590 : * verif de completude du spectre 591 : dfrq = (frequ(idep0) - frequ(ia)) / fref 592 : if((abs(dfrq)).ge.1.D-3) goto 1001 593 : idep=ia 594 : enddo 595 : endif 596 : 1001 continue 597 : c en haut 598 : if (ifin.lt.ihaut) then 599 : ifin0=ifin 600 : fref = max(1.D0,abs(frequ(ifin0))) 601 : fref = max(fref,abs(shifti)) 602 : do ia=(ifin0+1),ihaut,1 603 : c dfrq = (frequ(ifin0) - frequ(ia)) / max(1.D0,abs(frequ(ifin0))) 604 : dfrq = (frequ(ifin0) - frequ(ia)) / fref 605 : if((abs(dfrq)).ge.1.D-3) goto 1002 606 : ifin=ia 607 : enddo 608 : endif 609 : 1002 continue 610 : nbmode2 = ifin - idep + 1 611 : if(iimpi.ge.6) write(IOIMP,*) 612 : & 'ibas,ihaut,idep,ifin,nbmode2,nbmode', 613 : & ibas,ihaut,idep,ifin,nbmode2,nbmode 614 : if(nbmode2.ne.nbmode) write(IOIMP,*) 615 : & 'strate: Presence de modes multiples: on les ajoute' 616 : 617 : * -- ultime verif : controle de la completude du spectre -- 618 : c en bas 619 : if (idep.gt.1) then 620 : frshi1 = (1.D0-2.D-4) * frequ(idep) 621 : W2SHI1= DEUXPI * frshi1 622 : W2SHI1= W2SHI1 * W2SHI1 622 : CALL DECALE(IPRIGI,IPMASS,W2SHI1, IPKW21) 622 : CALL DIAGN1(IPKW21,nvpdep) 622 : CALL DTRIGI(IPKW21) 623 : else 624 : nvpdep=0 625 : endif 626 : c en haut 627 : frshi2 = (1.D0+2.D-4) * frequ(ifin) 628 : W2SHIF= DEUXPI * frshi2 629 : W2SHIF= W2SHIF * W2SHIF 629 : CALL DTRIGI(IPKW2M) 629 : CALL DECALE(IPRIGI,IPMASS,W2SHIF, IPKW2M) 629 : CALL DIAGN1(IPKW2M,nvpfin) 630 : nvpneg = nvpfin-nvpdep 631 : if(iimpi.ge.6) 632 : & write(IOIMP,*) 'nvpdep,nvpfin,nvpneg=',nvpdep,nvpfin,nvpneg 633 : * il manque des modes => on redemarre avec ce shift en orthogonalisant 634 : if (nvpneg.gt.nbmode2) then 635 : c ienc0=nvpneg-nbmode2 636 : c cbp: on en demande 2 * nbre vraiment manquant a cause des modes doubles 637 : c c + on avait calculé beaucoup, + on a de chance d'en avoir raté beaucoup 638 : c ienc = (2*ienc0) + (nbonZ/8) 639 : c bp : on orthogonalise par % aux modes deja convergé 640 : ienc=nvpneg-nbmode2 641 : Northo=ienc 642 : shipro = frshi2 643 : cbp: pour l instant on suppose qu il manque seulement de modes vers le haut 644 : if(iimpi.ge.6) write(IOIMP,*) 645 : & 'DIAGN1 trouve',ienc,' mode manquant' 646 : goto 2000 647 : elseif(nbmode2.ne.nbmode) then 648 : write(IOIMP,*) '...definitivement ',nbmode,'devient',nbmode2 649 : nbmode=nbmode2 650 : endif 651 : * il ne manque pas de modes => on termine 652 : 653 : 654 : *---- on refabrique un ipsolu qui contient juste les bonnes valeurs ---* 655 : segini,mso1=mso3 656 : nbnn=1 657 : nbelem=nbmode 658 : nbsous=0 659 : nbref=0 660 : segini ipt1 661 : n= nbelem 662 : segini msole1,msole2 663 : mso1.msolis(3)=ipt1 664 : mso1.msolis(4)=msole1 665 : mso1.msolis(5)=msole2 666 : idep1=idep-1 667 : if(iimpi.ge.6) write(ioimp,*) 'strate: idep1=',idep1 668 : do ia= 1,msole1.isolen(/1) 669 : if(iimpi.ge.6) 670 : & write(ioimp,*) 'strate: msole1(',ia,'/',n,immode(idep1+ia) 671 : msole1.isolen(ia)=immode(idep1+ia) 672 : ipt1.num(1,ia)=ipomod(idep1+ia) 673 : msole2.isolen(ia)=ipve(idep1+ia) 674 : enddo 675 : segdes msole1,msole2 676 : segdes ipt1 677 : segdes mso1 678 : ipsolu=mso1 679 : ifini=1 680 : if(iimpi.ge.7) 681 : & write(IOIMP,*) 'iexis=',(iexis(iio),iio=1,jg) 682 : if(iimpi.ge.6) 683 : & write(IOIMP,*) 'frequ=',(frequ(iio),iio=1,jg) 684 : return 685 : 686 : 687 : 2000 continue 688 : *---- on n'a pas fini : menage et nouveau cycle ----------------------* 689 : c ifini=0 690 : c bp: fait au debut 691 : nmod=min ( ienc+2, nmodma) 692 : frshif=shipro 693 : if(mso3.ne.0) segsup,mso3 694 : if(iimpi.ge.7) 695 : & write(IOIMP,*) 'iexis=',(iexis(iio),iio=1,jg) 696 : if(iimpi.ge.6) 697 : & write(IOIMP,*) 'frequ=',(frequ(iio),iio=1,jg) 698 : segdes,idist 699 : return 700 : 701 : end 702 : 703 : 704 : 705 : 706 : 707 : 708 : 709 : 710 : 711 : 712 : 713 : 714 :
© Cast3M 2003 - All rights reserved.
Disclaimer