Télécharger la plaquette Cast3M


fr en

Télécharger strate.eso

Retour à la liste

Numérotation des lignes : 

strate

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

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