Sources Esope | Cast3M

queshi
C QUESHI    SOURCE    BP208322  13/11/13    21:15:21     7865
      subroutine queshi(shifti,ienc,ifin,frequ,nmodma,shipro,iexis,
     $icle,nn,Northo)
ccc     $icle,nn,fpet,fgra)
*
      implicit real*8(a-h,o-z)
      implicit integer(i-n)
C-INC CCREEL
 
-INC PPARAM
-INC CCOPTIO
*
* routine pour calculer le shift a appliquer
* appelée par strate
*
* creation : chat 11.2009
* modifs : bp 01.2010 (traitements de qq cas particuliers)
*
* icle=1 on cherche à partir de ifin au dessus
* icle=2 on cherche au dessous
*
*
      dimension frequ(*),iexis(*)
*
      ideca= 4
      imul=1
      ifin2=nn
      if (icle.eq.2) then
        imul = -1
        ifin2=1
      endif
*
*---- on commence par regarder si on est pas en presence d'un trou ----*
*
      itrou=0
      do iou=ifin+imul,ifin2,imul
         if (iexis(iou).ne.0) then
cbp           itrou=(iou-ifin)*imul
           itrou=((iou-ifin)*imul) - 1
           go to 2
         endif
      enddo
    2 continue
c       write(6,*) 'iou,ifin,ifin2,itrou,ienc=',iou,ifin,ifin2,itrou,ienc
 
 
*---- on est en presence d'un trou -----------------------------------*
      if (itrou.ne.0) then
 
cbp reporté + loin        if (itrou.lt.ienc) then
 
   22      continue
           if (icle.eq.1) then
              fpet=frequ(ifin)
              fgra=frequ(iou)
           else
              fpet=frequ(iou)
              fgra=frequ(ifin)
           endif
cbp : ajout pour le cas ou delta =0 car on a raté un mode et mal numeroté
c     et on cherche entre 2 frequences identiques trouvées avec 2 cycles
           df1 = abs(fpet - fgra) / (max(1.D-6,(abs(fpet)+abs(fgra))))
           if (df1.lt.1.D-4) then
               if (iou.ne.ifin2.and.iexis(iou+imul).ne.0) then
                 iou=iou+imul
                 goto 22
               elseif(iexis(ifin-imul).ne.0) then
                 ifin=ifin-imul
                 goto 22
               else
                 write(6,*) 'ERREUR QUESHI : Contactez support '
                 stop
               endif
           endif
 
         if (itrou.lt.ienc) then
           fgra2=fgra*fgra
           fpet2=fpet*fpet
c            b2= 4.d0*fpet2
c            quatac=-2.d0*(fgra2-fpet2)
c            delta=b2 - quatac
c            if ( delta.lt.0.D0) then
c              call erreur(5)
c              return
c            endif
c            aa = sqrt(delta)
c            x2= (-2.d0 * fpet + aa )/ 2.d0
c            shipro=x2+fpet
cbp : simplification
           fmoy2 = (fgra2 + fpet2) / 2.D0
           shipro = sqrt(fmoy2)
cbp : cas d un petit trou, hyp. modes de multiplicite=2 => ienc*2
c           ienc=itrou + 1
c           ienc = 2 * (itrou + 1)
c puisqu il y a dans strate:  nmod=min ( ienc+2, nmodma), on peut reduire:
           ienc = 2 * itrou
           Northo = ienc - itrou
c
      if (iimpi.ge.6) then
        WRITE(6,*) '----------------- Queshi -----------------'
        WRITE(IOIMP,1991) itrou,ienc,shipro
1991    FORMAT(2X,'PETIT TROU DE ',I5,' MODES'
     C  ,'=> ON EN REDEMANDE ',I5,' AUTOUR DE ',F10.3)
      endif
c
           return
 
         else
* on est en presence d'un trou et itrou>ienc donc on avait visé beaucoup
* trop loin
           fpet=frequ(ifin)
           friou=frequ(iou)
           fgra= fpet + ( friou -fpet)/itrou*ienc
           Northo = 0
           go to 3
         endif
 
*---- pas de trou ----------------------------------------------------*
      else
cbp : on ne passe ici que pour un cas improbable
        if ( icle.eq.2) then
           if (ifin-ienc.le.0) then
        WRITE(6,*) '!!! on ne passe ici que pour un cas improbable !!!'
             ienc=ifin - 1
             shipro=0.d0
             Northo = ienc
             return
           endif
        endif
 
      endif
 
*---- calcul du shift a partir de la pente de y=freq(numero_mode) ----*
*     (cas sans trou)
      fpet= frequ(ifin)
      idec5  = ifin-  (ideca*imul)
      idec10 = ifin- 2*(ideca*imul)
      idec20 = ifin- 4*(ideca*imul)
      pentm=0.D0
      pent5=0.D0
      pent10=0.D0
      pent20=0.D0
cbp13/10/2010 : pentm = pente moyenne
      ndec = 0
cbp plutot que 3,2,1 on prend 1,2,3 pour une moyenne + equilibrée
      icoe1 = 1
      icoe2 = 2
      icoe3 = 3
      if (idec5.gt.0)  then
      if (iexis(idec5).ne.0)  then
         pent5 = fpet - frequ(idec5)
         pentm = pentm + (icoe1*(pent5/ideca))
         ndec = ndec + icoe1
      endif
      endif
      if (idec10.gt.0) then
      if (iexis(idec10).ne.0)  then
         pent10 = fpet - frequ(idec10)
         pentm = pentm + (icoe2*(pent10 /(2.*ideca)))
         ndec = ndec + icoe2
      endif
      endif
      if (idec20.gt.0) then
      if (iexis(idec20).ne.0)  then
         pent20 = fpet - frequ(idec20)
         pentm = pentm + (icoe3*pent20 / (4.*ideca))
         ndec = ndec + icoe3
      endif
      endif
cbp      if (ndec.eq.0) then
cbp        write(6,*) 'ERREUR dans queshi.eso : contactez le support'
cbp        return
cbp     endif
cbp        fgra=fpet + (pentm * ienc)
cbp 13/10/2010 :
      PENMIN = 3.E-4*fpet
      if ((ndec.ne.0.).and.(abs(pentm)).gt.PENMIN) then
         pentm = pentm/ndec
      else
c          pentm = imul * XPI * abs(fpet-shifti)
         pentm = imul * abs(fpet-shifti)
         write(6,*) 'DIFFICULTE dans QUESHI :',
     &              ' calcul de la pente majorée ',pentm
      endif
c     on change completement la strategie en se placant juste apres fpet
c       fgra=fpet + (pentm * (min(ienc,2.)))
c     on revient a une strategie moins aggressive
      fgra=fpet + (pentm*ienc)
      if (iimpi.ge.6) then
        WRITE(6,*) 'pent5,pent10,pent20,pentm,fpet,fgra='
        WRITE(6,*) pent5,pent10,pent20,pentm,fpet,fgra,imul
      endif
      Northo = max(4,(ienc/10))
 
*   calcul du shift (cas sans trou ou avec un grand trou)
    3 continue
      if (icle.eq.2) then
         az=fpet
         fpet=fgra
         fgra=az
      endif
c       fgra2= fgra*fgra
c       fpet2= fpet*fpet
c       b2= 4.d0*fpet2
c       quatac= -2.d0*( fgra2-fpet2)
c       delta=b2-quatac
c       if (delta.lt.0.D0) then
c          call erreur(5)
c          return
c       endif
c       aa=sqrt(delta)
c       x2= (-2.d0 * fpet + aa)/2.D0
c       shipro=x2 + fpet
cbp : simplification
*         fmoy2 = (fgra2 + fpet2) / 2.D0
*         shipro = sqrt(fmoy2)
         shipro = (fgra + fpet) / 2.D0
c       WRITE(6,*) 'fpet,fgra,shipro=',fpet,fgra,shipro
ccc         fgra = (shipro + fgra)/2.D0
c
      if (iimpi.ge.6) then
        WRITE(6,*) '----------------- Queshi -----------------'
        WRITE(IOIMP,1993) ienc,shipro
1993    FORMAT(2X,'=> ON DEMANDE ',I5,' MODES AUTOUR DE ',F12.3)
      endif
c
      return
c
      end