C EXCPHA    SOURCE    CHAT      11/03/16    21:21:36     6902           
      subroutine excpha
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z) 
-INC CCOPTIO
-INC SMCHAML
-INC SMCHPOI
-INC SMMODEL
-INC SMRIGID
-INC CCREEL
-INC SMCOORD
-INC SMELEME
      segment icpr(xcoor(/1)/(idim+1))
      segment icpr1(xcoor(/1)/(idim+1))
      segment icpr2(xcoor(/1)/(idim+1))
      segment tdeb(node)
      segment tprop(node)
      segment qrest(node)
      segment qtot(node)
      segment qsecon(node)
      segment timpo(node)
      segment iact(node)
      segment prop(node)
* lecture du modele
      call lirobj('MMODEL',mmodel,1,iretou)
      if(ierr.ne.0) return
* lecture du materiau ( pour les TPHAS)
      call lirobj('MCHAML',mchelm,1,iretou)
      if(ierr.ne.0) return
* lecture des proportions de phase
      call lirobj('MCHAML',mchel1,1,iretou)
      if(ierr.ne.0) return
* lecture des chaleur latentes reduies aux noeuds
      call lirobj('MCHAML',mchel2,1,iretou)
      if(ierr.ne.0) return
* lecture du chpoint de temperatures initiales
      call lirobj('CHPOINT',mchpoi,1,iretou)
      if(ierr.ne.0) return
* lecture du chpoint de temperatures proposées
      call lirobj('CHPOINT',mchpo1,1,iretou)
*      write(6,*) ' entree dans excpha temp proposee'
*      call ecchpo(mchpo1)
      if(ierr.ne.0) return
* lecture des rigidites de blocages de phases.
      call lirobj('RIGIDITE',mrigid,1,iretou)
      segact mrigid
      mpt = irigel(1,1)
*      write(6,*) ' maillage sous tendant les rigidites'
*      call ecmail ( mpt)
      if(ierr.ne.0) return
* fin des lectures
      xpre=1.d-8
      unxpre=1.D0-xpre
*
* algorithmes on boucle sur les rigidite s de blocages ou sur les
* modeles ( ils sont de structure identiques par construction des
* rigidites)
* puis si la rigidite avait ete prise en compte on teste le sens 
* et la valeur de la chaleur latente proposee
* sinon on teste la non traverse de la temperature de changement de
* phase et au besoin on active le blocage ( temperature imposee et 
*flux de chaleur latente limite)
*
* Au préalable repérage des champs de temperatures initiales et
* proposees d'abord champs de temperatures debut on melange les
* multiplicateurs

      segini icpr
      segact mchpoi
      ib = 0
      do 80 I=1,ipchp(/1)
      msoupo=ipchp(i)
      segact msoupo
      ipt1=igeoc
      segact ipt1
      do 81 j=1,ipt1.num(/2)
      ia = ipt1.num(1,j)
      if(icpr(ia).eq.0) then
      ib = ib +1
      icpr(ia)=ib
      endif
   81 continue
   80 continue
      node=ib
      segini tdeb
      do 82 I = 1 , ipchp(/1)
      msoupo=ipchp(i)
      ipt1=igeoc
      mpoval=ipoval
      segact mpoval 
      j = 1
      if(nocomp(/2) .ne.1) then
        do 83 j=1,nocomp(/2)
        if( nocomp(J).eq.'T') go to 84
   83   continue
        moterr(1:4)='   T'
        call erreur (181)
        return
   84   continue
      endif
      do 85 k=1,ipt1.num(/2)
      ia = ipt1.num(1,k)
      ib = icpr(ia)
      tdeb(ib)=vpocha(k,j)
   85 continue
   82 continue
* maintenant champs de temperatures proposees
      segact mchpo1
      ib = 0
      segini icpr1
      do 90 I=1,mchpo1.ipchp(/1)
      msoupo=mchpo1.ipchp(i)
      segact msoupo
      ipt1=igeoc
      segact ipt1
      do 91 j=1,ipt1.num(/2)
      ia = ipt1.num(1,j)
      if(icpr1(ia).eq.0) then
      ib = ib +1
      icpr1(ia)=ib
      endif
   91 continue
   90 continue
      node=max(node,ib)
      mnode=node
      segini tprop
      do 92 I = 1 , mchpo1.ipchp(/1)
      msoupo=mchpo1.ipchp(i)
      ipt1=igeoc
      mpoval=ipoval
      segact mpoval 
      j = 1
      if(nocomp(/2) .ne.1) then
      do 93 j=1,nocomp(/2)
      if( nocomp(J).eq.'T') go to 94
   93 continue
      moterr(1:4)='   T'
      call erreur (181)
      return
   94 continue
      endif
      do 95 k=1,ipt1.num(/2)
      ia = ipt1.num(1,k)
      ib = icpr1(ia)
      tprop(ib)=vpocha(k,j)
   95 continue
   92 continue
* au prealable encore, fabrication et reperage de Qsecon qui est le 
* champ de LX passé a resou( remplace les chaleurs latentes)
      segact mmodel,mrigid,mchelm,mchel1,mchel2
      segini icpr2
      ib=0
      do 200 meri = 1,irigel(/2)
      ipt3=irigel(1,meri)
      segact ipt3
      do 201 nel=1,ipt3.num(/2)
      do 201 nod = 1,1
      ia = ipt3.num(nod,nel)
      if(icpr2(ia).eq.0) then
      ib=ib+1
      icpr2(ia)=ib
      endif
  201 continue
  200 continue
      node = ib
      segini qsecon,timpo,iact
*
* debut du travail on boucle sur les sous zones du modeles
*
      do 100 mo=1,kmodel(/1)
      imodel = kmodel(mo)
      segact imodel
      meleme = imamod
      segact meleme
* recherche des chamelems correspondants 
* d'abord le champ de materiau recherche de TPHA
      do51 mchm=1,imache(/1)
       if( imache(mchm) . eq. meleme) then
         mchaml=ichaml(mchm)
         go to 52
       endif
   51  continue
       call erreur ( 472)
       return
   52  continue
       segact mchaml    
      do 56 n2=1,nomche(/2)
      if( nomche(n2).eq.'TPHA') then
        melval=ielval(n2)
        go to 57
      endif
   56 continue
      moterr(1:8) = 'TPHA'
      call erreur ( 677)
      return
   57 continue
* recherche de la proportion de phase PPHA
      do61 mchm=1,mchel1.imache(/1)
       if(mchel1. imache(mchm) . eq. meleme) then
         mcham1=mchel1.ichaml(mchm)
         go to 62
       endif
   61  continue
       call erreur ( 472)
       return
   62  continue
       segact mcham1    
      do 66 n2=1,mcham1.nomche(/2)
      if(mcham1. nomche(n2).eq.'PPHA') then
        melva1=mcham1.ielval(n2)
        go to 67
      endif
   66 continue
      moterr(1:8) = 'PPHA'
      call erreur ( 677)
      return
   67 continue
* recherche des chaleurs latentes reduites aux noeuds 
      do71 mchm=1,mchel2.imache(/1)
       if(mchel2. imache(mchm) . eq. meleme) then
         mcham2=mchel2.ichaml(mchm)
         go to 72
       endif
   71  continue
       call erreur ( 472)
       return
   72  continue
       segact mcham2    
      do 76 n2=1,mcham2.nomche(/2)
      if(mcham2. nomche(n2).eq.'Q   ') then
        melva2=mcham2.ielval(n2)
        go to 77
      endif
   76 continue
      moterr(1:8) = 'QLAT'
      call erreur ( 677)
      return
   77 continue
      segact melval,melva1,melva2
* par sous-modele on suppose TPHA constant.
* on extrait  la Temperature de changement de phase.
* fabrication de qtot : chaleur latente additionnée par point
* A l'aide de PPHA et de QLAR on fabrique un QREST qui est la chaleur
* latente qui reste à consommer pour ce sous modele.
*
      node=mnode
      segini qrest,qtot,prop
      do 110 mel=1,num(/2)
      do 110 nod = 1,num(/1)
      ia = num(nod,mel)
      ib = icpr(ia)
      pp = melva1.velche(nod,mel)
      if( pp.gt.1.d0) then
*       write(6,*) ' excp : erreur 110 1'
       pp = 1.d0
      endif
       if( pp. lt.0.d0) then
*       write(6,*) ' excp : erreur 110 2'
       pp = 0.D0
      endif
      prop(ib)=pp
      ql=melva2.velche(nod,mel)
      qr=ql* ( 1.d0 - pp)
      qrest(ib) = qrest(ib) + qr
      qtot(ib) = qtot(ib) + ql
      prop(ib)=1.d0 - qrest(ib)/qtot(ib)
  110 continue
      tt = velche(1,1) 
      ipt2= irigel(1,mo)
      segact ipt2
      do 120 meri=1,ipt2.num(/2)
      ip = ipt2.num(2,meri)
      il1=ipt2.num(1,meri)
* ce blocage existe -t-il?
      iex= icpr1(il1)
      tde = tdeb(icpr(ip))
*   tpro est la valeur proposée
      tpro=tprop(icpr1(ip))  
      if(iex.ne.0) then
*            if (iimpi.eq.1947) then
*         write(6,*) 'ip iex qtot(icpr(ip)',ip,iex,qtot(icpr(ip))
*         write(6,*) 'ip qres prop' ,qrest(icpr(ip)),prop(icpr(ip))
*         write(6,*) 'tde tpro ' , tde,tpro
*         write(6,*) 'il1  tprop(iex) ', il1,tprop(iex)
*        write(6,*) ' tt - tde ', tt-tde
*      write(6,*) 'tprop(iex)*qtot(icpr(ip)',tprop(iex)*qtot(icpr(ip))
*          endif
* le blocage existe on regarde si le point etait a la temperature de 
* changement de phase dans ce cas il faut garder le blocage si pas tout 
* deja dans un etat
        if( abs(tt-tde)/max(abs(tt),1.d0).le.1.e-8 ) then
*        if( iimpi.eq.1947) then 
*          write(6,*) ' detection d un point a la temp de changement'
*         endif
           if( (tprop(iex).gt.0.d0.and . prop(icpr(ip)).lT.unxpre).OR.
     $   (tprop(iex).lt.0.D0. and. prop(icpr(ip)). gt.xpre)) then
*            write(6,*) ' on garde la condition '
           ipl1=icpr2(il1)           
           iact(ipl1)=il1
           if(tprop(iex)*qtot(icpr(ip)).gt.0.D0) then
* la tempe cherche à monter
             qsecon(ipl1)=-qrest(icpr(ip))+qsecon(ipl1)
           else 
             qsecon(ipl1)=qtot(icpr(ip))-qrest(icpr(ip))+qsecon(ipl1)
           endif
         endif
* maintenant il faut regarder si le blocage est dans le bon sens
        elseif(tprop(iex)*qtot(icpr(ip)).gt.0.D0) then 
* la temperature cherche à monter 
* la temperature initiale etait-elle superieure a tt ?
* si non il faut prendre -qrest si oui relacher la condition
          if(tde . le . tt) then
*          if( iimpi.eq.1947) then
*           write(6,*) 'on garde le blocage'
*           endif
           ipl1=icpr2(il1)           
           qsecon(ipl1)=-qrest(icpr(ip))+qsecon(ipl1)
* on impose aussi la temperature 
           iact(ipl1)=il1
          endif
        else
* la temperature cherche a descendre
* la temperature initiale etait-elle inferieure a tt?
* si oui relacher la condition sinon prendre qtot - qrest
          if(tde.ge. tt) then
*          if(iimpi.eq.1947) then
*           write(6,*) ' temperature chute on solidifie'
*           endif
           ipl1=icpr2(il1)           
           qsecon(ipl1)=qtot(icpr(ip))-qrest(icpr(ip))+qsecon(ipl1)
* on impose aussi la temperature 
           iact(ipl1)=il1
          endif
       endif
      else
* le blocage n'existe pas
* verification que tdeb et tprop son du meme cote que tt
* si c'est vrai rien a faire sinon impose le blocage et  +ou- qtot
        aa = (tt - tde)* ( tt - tpro)
        tsens= tpro-tde
        if( abs((tt-tde)/(tt+1.d-30)).le.1.e-9 ) aa=0.d0
        if(aa.eq.0.D0) then
*       if(aa.le.0.D0) then
*         write(6,*) 'ip iex qtot(icpr(ip)',ip,iex,qtot(icpr(ip))
*         write(6,*) ' tt - tde ', tt-tde
*         write(6,*) 'tde tpro aa ', tde,tpro, aa
*        write(6,*)  'tt - tde' , 'tt - tpro',tt - tde,tt - tpro
*        uuut=tprop(iex)
*        write(6,*)' tprop(iex)' , uuut 
*       endif     
         if( (tsens.gt.0.d0.and . prop(icpr(ip)).lT.unxpre).OR.
     $   (tsens.lt.0.D0. and. prop(icpr(ip)). gt.xpre)) then
*          write(6,*) ' on passe ici'
*         if( iimpi.eq.1947) then
*           write(6,*) ' un point changeait de phase'
*         endif 
           xde=1.d0
          if( tpro.gt.tt*1.0000001) xde=-1.D0
          ipl1=icpr2(il1) 
                  
*         if( iimpi.eq.1947) then
*         write(6,*) ' on garde le point ip ',ip
*         write(6,*) '2 il1,il2,ipl1,ipl2 xde',il1,il2,ipl1,ipl2,xde
*         endif
          qsecon(ipl1)=qtot(icpr(ip))* xde+qsecon(ipl1)
          iact(ipl1)=il1
*          write(6,*) ' qsecon (ipl1), ip,qtot(icpr(ip))',
*     $       qsecon (ipl1), ip,qtot(icpr(ip)) 
         endif      
        elseif( aa . lt . 0.D0) then
* tous les deux de chaque cote impose le blocage et la temperature
* impose  
          xde=1.d0
          if((tt-tde*1.000000001). GT.0.) xde = -1.d0   
          ipl1=icpr2(il1)           
*          if( iimpi.eq.1947) then
*           write(6,*) 'il1,il2,ipl1,ipl2 xde',il1,il2,ipl1,ipl2,xde
*           endif
          qsecon(ipl1)=qtot(icpr(ip))* xde+qsecon(ipl1)
          iact(ipl1)=il1
*           write(6,*) 'qtot qsec',qtot(icpr(ip)),qsecon(ipl1)
       endif
      endif
* fin de la boucle sur les rigidites de blocages 
  120 continue
* fin de la boucle sur les sous paquets de rigidite
  100 continue
* on dispose maintenant d'un tableau qsecon et timpo et iact qui
* donnent
* les second membre de chaleur latentes a mettre en LX et les
* temperatures de fusion a mettre en FLX
* on cree le champoint
      nbelem=0
*      write(6,*) (iact(ij),ij=1,iact(/1))
      do 300 ia=1,iact(/1)
      if(iact(ia).ne.0) nbelem = nbelem + 1
  300 continue
      nbsous=0
      nbref=0
      nbnn=1
      if(nbelem.eq.0) nbnn=0
      segini ipt4
      ipt4.itypel=1
      nat=1
*       write(6,*) ' nbnn nbelem ' , nbnn,nbelem
      IF(nbelem.ne.0) then
        nsoupo=1
        segini mchpo3
        mchpo3.mochde(1:72)=' creation par excpha'
        mchpo3.mtypoi='contmult'
        nc=1
        segini msoupo
        mchpo3. jattri(1)=2
        mchpo3.ipchp(1)=msoupo
        nocomp(1) = 'LX'
        igeoc=ipt4
        n=nbelem
        segini mpoval
        ipoval = mpoval
        ipo=0
        do 301,ia=1,iact(/1)
        if (iact(ia).ne.0) then
          ipo=ipo+1
          ipt4.num(1,ipo) = iact(ia)
          vpocha(ipo,1)=qsecon(icpr2(iact(ia)))
        endif
 301    continue
*        write(6,*) ' nbelem ipo', nbelem ,ipo
         segdes mpoval,ipt4,msoupo,mchpo3
      else
        nsoupo=0
        segini mchpo3
           mchpo3.mochde(1:72)='creation par excpha '
           mchpo3.jattri(1)=2
           mchpo3.mtypoi='contmult'
        segdes mchpo3
      endif
      call ecrobj('CHPOINT',mchpo3)
*      write(6,*) ' ecriture du chpoint, dans excpha, aprs ecrobj'
*      call ecchpo ( mchpo3)
      segdes melvaL,melva1,melva2
      do 303 ir=1,irigel(/2)
      meleme  = irigel(1,ir)
      segdes meleme
      imodel = kmodel(ir)
      meleme =  imamod
      segdes meleme,imodel
  303 continue
      segdes mrigid,mmodel
      segdes mchelm,mchel1,mchel2
* desactivation champ de temperature initiale
      do 304 ia=1,ipchp(/1)
      msoupo=ipchp(ia)
      meleme = igeoc
      segdes meleme
      mpoval = ipoval
      segdes mpoval
      segdes msoupo
  304 continue
      segdes,mchpoi
* desactivation champ de temperature propose
      IF (mchpoi.NE.mchpo1) THEN
      mchpoi= mchpo1
      do 305 ia=1,ipchp(/1)
      msoupo=ipchp(ia)
      meleme = igeoc
      segdes meleme
      mpoval = ipoval
      segdes mpoval
      segdes msoupo
  305 continue
      segdes,mchpoi
      ENDIF
      segsup icpr,icpr1,icpr2,tdeb,qsecon,timpo,tprop,iact
      segsup qrest,qtot,prop
      return
      end