Sources Esope | Cast3M

dezoma
C DEZOMA    SOURCE    CHAT      05/01/12    22:47:58     5004
CCC
C **********************************************************************
CCC
      SUBROUTINE DETERZONAMAC(SIGTRI,NDIMS,VARINT,IPLCON,IPLCAP,IPLAPEX,
     .                        nescri,ues)
      IMPLICIT INTEGER(I-N)
      integer ndims,iplcon,iplcap,iplapex,nescri,ues
      real*8    sigtri(ndims),varint(2),siginvari(3),Eauxp,Eauxq,
     .          Yicon,Yicap,p,q,t,xco,xca,psepara,aux3,qcor,Nr_qfunc,
     .          vecp,vecq,pcap,prodesc,wwfunc,etaconfunc,pcapfunc
      real*8    exm,qaa,pcc,eww,xk1,xk2,etb,exv,xep,paa,cco,
     .          exl,xnn,aaa,bbb,pc0,exr,cca,phi,alp,pcc2,paa2,qaa2,
     .          exm2,eww2,cpr,cpm,cfr,you,xnu
      common /elasdata/ you,xnu
      common /MRScone/ exm,qaa,pcc,eww,xk1,xk2,etb,exv,xep,paa,cco,
     .                 exl,xnn,aaa,bbb
      common /MRScapp/ pc0,exr,cca,phi,alp,pcc2,paa2,qaa2,exm2,eww2,
     .                 cpr,cpm,cfr
      iplcon =0
      iplcap =0
      iplapex=0
      do i=1,3
          siginvari(i)=0.D0
        enddo
      call InvariantesPQT(sigtri,ndims,siginvari)
      p=siginvari(1)
      q=siginvari(2)
      t=siginvari(3)
      xco=varint(1)
      xca=varint(2)
      Eauxp=you/(1.D0-2.d0*xnu)/9.D0
      Eauxq=you/(1.d0+xnu)/2.D0
      call yieldcoMAC(siginvari,xco,Yicon)
      call yieldcaMAC(siginvari,xca,xco,Yicap)
      pcap=pcapfunc(xca)
      psepara=alp*pcap
      aux3=etaconfunc(xco)*(psepara-pcc)/wwfunc(t,eww)
      qcor=NR_qfunc(aux3,qaa,exm,aux3,1.D-10)
      if(nescri.eq.1) then
         write(ues,999) '   Invari : ',p,q,t
         write(ues,999) '   Var int: ',xco,xca
         write(ues,999) '    Yicon : ',yicon
         write(ues,999) '    Yicap : ',yicap
         write(ues,999) '   Sep.p_q: ',psepara,qcor
 999     format(a12,3(E15.9,1x))
      endif
c zona del cone
      If (p.le.psepara) then
         If (Yicon.le.0.D0) return
         iplcon=1
c mirar si es apex
         vecp=-xnn*etaconfunc(xco)*(pcc-alp*pcap)/(pcc+alp*pcap)
         vecq=wwfunc(t,eww)*Eauxq/Eauxp
         prodesc=((p-pcc)*(-vecq)+q*(vecp))
         if (nescri.eq.1) then
            write(ues,999)' CONE GrApex: ', vecp,vecq,xnn
            write(ues,999)'      PtApex: ', pcc,0.D0
            write(ues,999)'      Prodes: ', prodesc
         endif
         if (prodesc.gt.0.D0) then
            iplapex=1
            return
         endif
c zona del cap
      else
         If (Yicap.le.0.D0) return
         iplcap=1
      endif
      return
      end