C PHAPRO    SOURCE    CHAT      11/03/16    21:28:55     6902           
      subroutine phapro
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
*
*  Calcule la nouvelle proportions de phase
* en entree 
*           : le modele et son chamelem de materiau associé
*           : un champ de temperature proposée avec les mult de
*             lagrange
*  en sortie chamelem de proportion de phase mis a jour
*
-INC CCOPTIO
-INC SMMODEL
-INC SMCHAML
-INC SMCHPOI
-INC SMELEME
-INC SMRIGID
-INC SMCOORD
-INC CCREEL

      segment icp
        integer icpr1(npt)
      endsegment
      pointeur icpr.icp
      segment icon(imax)
      segment ipos(ino)
      segment ipo2(ino)
      segment imoy(npart)
      segment idej(npart)
      segment qrest(ino)
      segment qtot(ino)
      segment ivoi(nvoi)
      segment  qmoye
          real*8 qmoy(nbelem)
      endsegment
      xpre=xpetit
      umxpre=1.d0 - xpre
* lecture du chamelem de source de proportions de phase
      call lirobj('MCHAML',mchelm,1,iretou)
      if(ierr.ne.0) return
* lecture du chamelem de source de chaleur latente
      call lirobj('MCHAML',mchel1,1,iretou)
      if(ierr.ne.0) return
* lecture des reactions
      call lirobj('CHPOINT',mchpo1,1,iretou)
      if(ierr.ne.0) return
* recopie du mcheml entré
      call copie8(mchelm,ipoin2)
      mchelm=ipoin2
* ON regarde d'abord le chpoint des reactions. s'il est vide il n'y a
* rien a faire 
      segact mchpo1
      if(mchpo1.ipchp(/1).eq.0) then    
        call ecrobj('MCHAML',ipoin2)
        segdes mchpo1
        return
      endif
*  test du chpoint en entree
      if(mchpo1.ipchp(/1) .ne.1) then
       call erreur(180)
       return
      endif
      msoupo=mchpo1.ipchp(1)
      segact msoupo
      meleme=igeoc
      segact meleme
      mpoval=ipoval
      segact mpoval
* creation du tableau des connexions a partir du chamelem de chaleur
*  latente
      segact mchel1
      npt = xcoor(/1)/(idim+1)
      segini icp
      ngrand=0
      do 520 ima=1,mchel1.imache(/1)
         meleme=mchel1.imache(ima)
         segact meleme
          ngrand = max(ngrand,num(/2)) 
         do 521 nel=1,num(/2)
         do 521 nbn=1,num(/1)
         ia = num(nbn,nel)
         icpr1(ia)=icpr1(ia)+1
  521    continue
  520 continue
      imax=0
      ino=1
      ngrand=ngrand+1
      do 522 ic=1,icpr1(/1)
       if(icpr1(ic).ne.0) then
        ino=ino+1
        imax=imax+icpr1(ic)
       endif
  522 continue
      segini icon,ipos,ipo2
      ia=1
      nvoi=0
      do 525 ic=1,icpr1(/1)
       if(icpr1(ic).ne.0) then
       ia = ia + 1
       nvoi = max(icpr1(ic),nvoi)
       ipos(ia)=ipos(ia-1) + icpr1(ic)
       icpr1(ic)=ia-1
       endif
  525 continue
      do 523 ima = 1,mchel1.imache(/1)
      meleme= mchel1.imache(ima)
      nplu= ngrand*(ima-1)
         do 524 nel=1,num(/2)
         do 524 nbn=1,num(/1)
         ia = num(nbn,nel)
         ib=icpr1(ia)
         ipo2(ib)=ipo2(ib)+1
         id=ipos(ib)+ipo2(ib)
         icon(id)=nel+nplu
  524 continue
  523 continue
      nvoi=nvoi+1
      segini ivoi
      segsup ipo2  
* activation du chamelem de chaleur latente et de proportion de phase
       segact mchelm
      do 560 n1=1,imache(/1)
      meleme=imache(n1)
      mchaml=ichaml(n1)
      mcham1=mchel1.ichaml(n1)
      segact meleme,mchaml,mcham1
      do 561 n2=1,ielval(/1)
      melval=ielval(n2)
      melva1=mcham1.ielval(n2)
      segact melval*MOD,melva1*MOD
  561 continue
  560 continue
* fabrication de la proportion de phase deja mangee par element
      ifois=0
      segini,icpr=icp
 1000 continue
*      write(6,*) ' icpr   1 '
*      write(6,*) ( icpr.icpr1(iop),iop=1,10)
      ifois=ifois+1
      npart=imache(/1)
      if(ifois.eq.2) then
         ino=qrest(/1)
         segsup qrest,qtot
         do 1005 IMO=1,imoy(/1)
         qmoye = imoy(imo)
         nbelem=qmoy(/1)
         segsup qmoye
 1005    continue
         segsup imoy
      endif
      segini imoy,qrest,qtot
      do 570 n1=1,imache(/1)
      meleme=imache(n1)
      do 580 ico=1,mchel1.imache(/1)
      if(meleme.eq.mchel1.imache(ico)) go to 581
  580 continue
      call erreur(21)
  581 continue
      mcham1=mchel1.ichaml(ico)
      mchaml=ichaml(n1) 
      nbelem = num(/2)
      segini qmoye
      imoy(n1)=qmoye
      do 571 inom=1,nomche(/2)
      if( nomche(inom).eq.'PPHA') go to 573
  571 continue
      moterr(1:4)='PPHA'
      call erreur(236)
      return
  573 continue
      do 583  inom1=1,mcham1.nomche(/2)
      if(mcham1.nomche(inom1).eq.'Q') go to 584
  583 continue
      moterr(1:4)='Q'
      call erreur(236)
      return
  584 continue
      melva1=mcham1.ielval(inom1)  
      melval=ielval(inom)
      nbnn=num(/1)      
      do 575 nbe=1,num(/2)
        qs=0.d0
        nbnn=num(/1)
*        ipr=1
        do 574 nbn=1,num(/1)
        is=icpr1(num(nbn,nbe))
        qtot(is) =melva1.velche(nbn,nbe)+qtot(is)
        qrest(is)=qrest(is)+
     $             melva1.velche(nbn,nbe)*(1.d0-velche(nbn,nbe))
        qs=qs + velche(nbn,nbe) 
  574   continue
        qmoy(nbe)=qs/nbnn
  575 continue
*      write(6,*) ' qmoy'
*      write(6,*) ( qmoy(ip),ip=1,qmoy(/1))
  570 continue             
* boucle sur les noeuds concernés on y passe deux fois, la premiere
* pour traiter les noeuds qui basculent completement
       do 600 ipc=1,mchpo1.ipchp(/1)
         msoupo=mchpo1.ipchp(ipc)
         segact msoupo
         meleme=igeoc
         segact meleme
         mpoval=ipoval
         segact mpoval
         do 602 noc=1,nocomp(/2)
           if( nocomp(noc).eq.'Q') go to 603
  602    continue
  603    continue
*        write(6,*)'num'
*      write(6,*) ( num(1,iop),iop=1,num(/2))
         do 601 nel=1,num(/2)
* on commence par comparer qreac avec qlat pour voir si tout est mangé
           qreac=vpocha(nel,noc)
           ia = num(1,nel)
           ib = icpr1(ia)
*           write(6,*)' nel,ia,ib,icpr.icpr1(ia)',
*     $     nel,ia,ib,icpr.icpr1(ia)
           if(icpr.icpr1(ia).EQ.0) go to 601
           qres= qrest(ib)
           qto=qtot(ib) 
           ic=ipos(ib)+1
           id=ipos(ib+1)            
*      write(6,*)'ia,ib,qreac,qres,qto',
*     $   ia,ib,qreac,qres,qto,icpr.icpr1(ia)
*      write(6,*) 'ifois' , ifois 
           if(ifois.eq.2) go to 620
           if( qreac.lt.0.D0) then
                 if(-qreac.lt.qres*umxpre) go to 601
                 xval=1.d0
          else 
                if(qreac .lt.(qto-qres)*umxpre) go to 601
                xval=0.d0
         endif
* tous les elements sont manges  
          ipt2= imache(1)
*          write(6,987) ((ipt2.num(iu,ip),iu=1,4),ip=1,92)
*  987 format(8I8) 
          icpr.icpr1(ia)=0  
*        write(6,*) 'traitement complet de  ia ib', ia,ib       
              do 610 ie=ic,id
              il=icon(ie)
              ilv=mod(il,ngrand)
              ieme= (il - ilv)/ngrand  + 1
              ipt1= imache(ieme)
              mchaml=ichaml(ieme)
              melval=ielval(1)                       
              do 611 ino=1,ipt1.num(/1)
                 if( ipt1.num(ino,ilv).eq.ia) then
*       write(6,*)'ie,il,ilv,ieme,mchaml,melval',ie,il,ilv,
*     $    ieme,mchaml,melval    
*       write(6,*) ' ia,ie,ipt1.num(ino,ilv)',ia,ie,ipt1.num(ino,ilv) 
                  velche(ino,ilv)=xval 
                 go to 610
                 endif
  611         continue
  610        continue
             go to 601
  620     continue
* il va falloir en tartiner sur certains elements
* on prend d'abord l'élement dont la plus gande proportion est dans le
* sens qui nous interesse
* si ifon =1 on fonds sinon on solidifie
       do 631 ivo=ic,id
       ivoi(ivo-ic+1)=0
 631   continue
 630   continue
       ifon=0 
       iiem=-1
       qma= -1.5
       qmi=1.5           
       if(qreac.lt.0.D0)  then
         ifon=1
         qreac=-qreac   
       endif   
*       write(6,*) ' ifon ' , ifon
           do 502 ie=ic,id
           if(ivoi(ie-ic+1).eq.1) go to 502
           il=icon(ie)
           ilv = mod(il,ngrand)
           ieme= (il - ilv)/ngrand  + 1
           ipt1= imache(ieme)
           qmoye=imoy(ieme)
*            write(6,*) ' qmoy(ilv)  ilv ieme ',ilv, qmoy(ilv)          
           if( ifon.eq.1) then
             if(qmoy(ilv).lt.umxpre. and. qmoy(ilv).GT.qma) then
              iiem=ieme
              ieg=ie
              iilv=ilv
              qma=qmoy(ilv)
             endif
          else
             if(qmoy(ilv).ge.xpre.and.qmoy(ilv).lt.qmi) then
               iiem=ieme
               iilv=ilv
               ieg=ie
               qmi=qmoy(ilv)
             endif
          endif
  502  continue
*      write(6,*) ' qma qmi iiem iilv', qma, qmi,iiem,iilv
       if ( iiem.eq.-1) go to 601
       ivoi(ieg-ic+1)=1
*    mise a jour du point en question dans l'element trouvé
       mchaml=ichaml(iiem)
       mcham1=mchel1.ichaml(iiem)
       do 506 n2=1,mcham1.nomche(/2)
         if(mcham1.nomche(n2).eq.'Q') GO TO 507
  506 continue
  507 continue
      melva1=mcham1.ielval(n2) 
      melval=mchaml.ielval(1)
      ipt1=imache(iiem)
       do 505 ino=1,ipt1.num(/1)
         ir= ipt1.num(ino,iilv)
*         write(6,*) ' ir, ia ifon', ir,ia,ifon
         if( ir.eq.ia) then
          if(ifon.eq.1) then
            qmax=melva1.velche(ino,iilv)*(1.d0 - velche(ino,iilv))
*            write(6,*) ' qmax qreac melva1. ',qmax,qreac,
*     $  melva1.velche(ino,iilv) 
            if(qmax.le.qreac) then
               velche(ino,iilv)=1.d0
                if( (qreac-qmax)/melva1.velche(ino,iilv).le.xpre) 
     $           go to 601
                qreac=qreac-qmax
                qreac = - qreac
*               write(6,*) ' on repart en 630 qreac ', qreac
               go to 630
           else
               xpro= qreac / melva1.velche(ino,iilv)
               velche(ino,iilv)=velche(ino,iilv) + xpro
           endif
        else
           qmin=melva1.velche(ino,iilv)*velche(ino,iilv)
           if(qmin.le.qreac) then
             velche(ino,iilv)=0.d0
             if(abs((qmin-qreac)/melva1.velche(ino,iilv)).le.xpre)
     $          go to 601
             qreac = qreac - qmin
             go to 630
          else
             xpro = qreac /melva1.velche(ino,iilv)
             velche(ino,iilv)=velche(ino,iilv) - xpro   
          endif
       endif 
       go to 601
       endif
  505  continue
  601 continue  
  600 continue   
      if(ifois.eq.1) go to 1000      
*  travail terminé il faut desactiver
      segsup icon,ipos,ipo2,icp,ivoi,icpr,imoy,qrest,qtot,qmoye
* desactivation mchelm,mchel1
      do 961 ich=1,mchelm.imache(/1)
       meleme=imache(ich)
       segdes meleme 
      mchaml=ichaml(ich)
      do 962 ive=1,ielval(/1)      
      melval=ielval(ive)
* test de validité des valeur   0.d0 <=  val <= 1.d0
       do 976 iu=1,velche(/2)
        do 976 io=1,velche(/1)
        if( velche(Io,iu). lT.0.D0) then
*        write(6,*) ' phapro: proportion < 0. : ' ,velche(Io,iu) 
         velche(Io,iu)=0.d0
        elseif (velche(Io,iu). gT.1.D0) then
*        write(6,*) ' phapro: proportion > 1. : ',velche(Io,iu) 
       velche(Io,iu)=1.d0
      endif
  976 continue
      segdes melval
  962 continue
      segdes mchaml
  961 continue
      segdes mchelm
      call ecrobj ('MCHAML',mchelm)
      mchelm=mchel1
      do 963 ich=1,mchelm.imache(/1)
       meleme=imache(ich)
       segdes meleme 
      mchaml=ichaml(ich)
      do 964 ive=1,ielval(/1)      
      melval=ielval(ive)
      segdes melval
  964 continue
      segdes mchaml
  963 continue
      do 965 ipc=1,mchpo1.ipchp(/1)
      msoupo = mchpo1.ipchp(ipc)
      meleme=igeoc
      segdes meleme
      mpoval=ipoval
      segdes mpoval
      segdes msoupo
  965 continue
      segdes mchpo1      
      return
      end