text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2021 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine incplas(elconloc,plconloc,xstate,xstateini,
     &  elas,emec,ithermal,icmd,beta,stre,vj,kode,
     &  ielas,amat,t1l,dtime,time,ttime,iel,iint,nstate_,mi,
     &  eloc,pgauss,nmethod,pnewdt,depvisc)
!
!     calculates stiffness and stresses for the incremental plasticity
!     material law (Ref: J.C. Simo, A framework for finite strain
!     elastoplasticity, Comp. Meth. Appl. Mech. Engng., 66(1988)199-219
!     and 68(1988)1-31)
!
!     icmd=3: calculates stress at mechanical strain
!     else: calculates stress at mechanical strain and the stiffness
!           matrix
!
!     the stresses in the routine proposed by Simo are Kirchhoff
!     stresses. Since the stress in the hardening laws are Chauchy
!     stresses, they are converted into Kirchhoff stress by
!     multiplication with the Jacobian determinant
!
      implicit none
!
      character*80 amat
!
      integer ithermal(*),icmd,i,j,kode,ivisco,ielastic,
     &  niso,nkin,ielas,iel,iint,nstate_,mi(*),id,leximp,lend,layer,
     &  kspt,kstep,kinc,iloop,nmethod,user_hardening,user_creep
!
      real*8 elconloc(*),elas(21),emec(6),beta(6),stre(6),
     &  vj,plconloc(802),stbl(6),stril(6),xitril(6),
     &  ee,un,um,al,xk,cop,umb,umbb,dxitril,f0,d0,f1,d1,d2,xg(3,3),
     &  xs(3,3),xx(3,3),xn(3,3),xd(3,3),cpl(6),c(6),ci(6),
     &  c1,c2,c3,c4,c5,c6,c7,c8,c9,cplb(6),stblb(6),
     &  ftrial,xiso(200),yiso(200),xkin(200),ykin(200),
     &  fiso,dfiso,fkin,dfkin,fiso0,fkin0,ep,t1l,dtime,
     &  epini,a1,dsvm,xxa,xxn,vj2,vj23,sc(3,3),depvisc,
     &  cop1,cop2,fu1,fu2,fu,dcop,time,ttime,eloc(6),
     &  xstate(nstate_,mi(1),*),xstateini(nstate_,mi(1),*),
     &  g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12,g13,g14,g15,g16,
     &  g17,g18,g28,g29,g30,g31,g32,g33,decra(5),deswa(5),serd,
     &  esw(2),ec(2),p,qtild,predef(1),dpred(1),timeabq(2),pgauss(3),
     &  dtemp,pnewdt
!
c      data kk /1,1,1,1,1,1,2,2,2,2,2,2,1,1,3,3,2,2,3,3,3,3,3,3,
c     &  1,1,1,2,2,2,1,2,3,3,1,2,1,2,1,2,1,1,1,3,2,2,1,3,3,3,1,3,
c     &  1,2,1,3,1,3,1,3,1,1,2,3,2,2,2,3,3,3,2,3,1,2,2,3,1,3,2,3,
c     &  2,3,2,3/
!
c      pnewdt=-1.d0
      leximp=1
      lend=2
      user_creep=0
!
!     localizing the plastic fields
!
      do i=1,6
         cpl(i)=-2.d0*xstateini(1+i,iint,iel)
         stbl(i)=xstateini(7+i,iint,iel)
      enddo
      do i=1,3
         cpl(i)=cpl(i)+1.d0
      enddo
      epini=xstateini(1,iint,iel)
!
      ee=elconloc(1)
      un=elconloc(2)
      um=ee/(1.d0+un)
      al=um*un/(1.d0-2.d0*un)
      xk=al+um/3.d0
      um=um/2.d0
!
      ep=epini
!
!     right Cauchy-Green tensor (eloc contains the Lagrange strain,
!     including thermal strain)
!
      c(1)=2.d0*emec(1)+1.d0
      c(2)=2.d0*emec(2)+1.d0
      c(3)=2.d0*emec(3)+1.d0
      c(4)=2.d0*emec(4)
      c(5)=2.d0*emec(5)
      c(6)=2.d0*emec(6)
!
!     calculating the Jacobian
!
      vj=c(1)*(c(2)*c(3)-c(6)*c(6))
     &  -c(4)*(c(4)*c(3)-c(6)*c(5))
     &  +c(5)*(c(4)*c(6)-c(2)*c(5))
      if(vj.gt.1.d-30) then
         vj=dsqrt(vj)
      else
ccc         write(*,*) '*WARNING in incplas: deformation inside-out'
!
!        deformation is reset to zero in order to continue the
!        calculation. Alternatively, a flag could be set forcing
!        a reiteration of the increment with a smaller size (to
!        be done)
!
         c(1)=1.d0
         c(2)=1.d0
         c(3)=1.d0
         c(4)=0.d0
         c(5)=0.d0
         c(6)=0.d0
         vj=1.d0
      endif
!
!     check whether the user activated a viscous calculation
!     (*VISCO instead of *STATIC)
!
      if((nmethod.ne.1).or.(ithermal(1).eq.3)) then
         ivisco=1
      else
         ivisco=0
      endif
!
!     check for user subroutines
!
      if((plconloc(801).lt.0.8d0).and.(plconloc(802).lt.0.8d0)) then
         user_hardening=1
      else
         user_hardening=0
      endif
c      if(kode.eq.-52) then
      if((kode.eq.-52).and.(ivisco.eq.1)) then
         if(elconloc(3).lt.0.d0) then
            user_creep=1
         else
            user_creep=0
            xxa=elconloc(3)*(ttime+time)**elconloc(5)
            if(xxa.lt.1.d-20) xxa=1.d-20
            xxn=elconloc(4)
            a1=xxa*dtime
         endif
      endif
!
!        inversion of the right Cauchy-Green tensor
!
      vj2=vj*vj
      ci(1)=(c(2)*c(3)-c(6)*c(6))/vj2
      ci(2)=(c(1)*c(3)-c(5)*c(5))/vj2
      ci(3)=(c(1)*c(2)-c(4)*c(4))/vj2
      ci(4)=(c(5)*c(6)-c(4)*c(3))/vj2
      ci(5)=(c(4)*c(6)-c(2)*c(5))/vj2
      ci(6)=(c(4)*c(5)-c(1)*c(6))/vj2
!
!        reducing the plastic right Cauchy-Green tensor and
!        the back stress to "isochoric" quantities (b stands
!        for bar)
!
      vj23=vj**(2.d0/3.d0)
      do i=1,6
         cplb(i)=cpl(i)/vj23
         stblb(i)=stbl(i)/vj23
      enddo
!
!        calculating the (n+1) trace and the (n+1) deviation of
!        the (n) "isochoric" plastic right Cauchy-Green tensor
!
      umb=(c(1)*cplb(1)+c(2)*cplb(2)+c(3)*cplb(3)+
     &     2.d0*(c(4)*cplb(4)+c(5)*cplb(5)+c(6)*cplb(6)))/3.d0
      do i=1,6
         cplb(i)=cplb(i)-umb*ci(i)
      enddo
!
!        calculating the (n+1) trace and the (n+1) deviation of
!        the (n) "isochoric" back stress tensor
!
      umbb=(c(1)*stblb(1)+c(2)*stblb(2)+c(3)*stblb(3)+
     &     2.d0*(c(4)*stblb(4)+c(5)*stblb(5)+c(6)*stblb(6)))/3.d0
      do i=1,6
         stblb(i)=stblb(i)-umbb*ci(i)
      enddo
!
!        calculating the trial stress
!
      do i=1,6
         stril(i)=um*cplb(i)-beta(i)
      enddo
!
!        calculating the trial radius vector of the yield surface
!
      do i=1,6
         xitril(i)=stril(i)-stblb(i)
      enddo
!
      sc(1,1)=xitril(1)*c(1)+xitril(4)*c(4)+xitril(5)*c(5)
      sc(1,2)=xitril(1)*c(4)+xitril(4)*c(2)+xitril(5)*c(6)
      sc(1,3)=xitril(1)*c(5)+xitril(4)*c(6)+xitril(5)*c(3)
      sc(2,1)=xitril(4)*c(1)+xitril(2)*c(4)+xitril(6)*c(5)
      sc(2,2)=xitril(4)*c(4)+xitril(2)*c(2)+xitril(6)*c(6)
      sc(2,3)=xitril(4)*c(5)+xitril(2)*c(6)+xitril(6)*c(3)
      sc(3,1)=xitril(5)*c(1)+xitril(6)*c(4)+xitril(3)*c(5)
      sc(3,2)=xitril(5)*c(4)+xitril(6)*c(2)+xitril(3)*c(6)
      sc(3,3)=xitril(5)*c(5)+xitril(6)*c(6)+xitril(3)*c(3)
!
      dxitril=sc(1,1)*sc(1,1)+sc(2,2)*sc(2,2)+sc(3,3)*sc(3,3)+
     &  2.d0*(sc(1,2)*sc(2,1)+sc(1,3)*sc(3,1)+sc(2,3)*sc(3,2))
c      g1=c(6)
c      g2=xitril(6)
c      g3=xitril(3)
c      g4=xitril(2)
c      g5=c(5)
c      g6=xitril(5)
c      g7=xitril(4)
c      g8=c(4)
c      g9=c(3)
c      g10=c(2)
c      g11=c(1)
c      g12=xitril(1)
c      g13=g12*g11
c      g14=g10*g4
c      g15=g9*g3
c      g16=g8*g7
c      g17=g6*g5
c      g18=g2*g1
c      g28=4*(g16 + g15)
c      g29=4*g13
c      g30=4*g14
c      g31=4*g6*g1
c      g32=4*g8*g2
c      g33=4*g7*g5
c      dxitril=(g13*g13 + g14*g14 + g15*g15 + g16*(g30 + g29 + 2*
c     &     g16) + g17*(g29 + g28 + 2*g17) + g18*(g30 + g28 + 2*
c     &     g18 + 4*g17) + g11*g7*(g31 + 2*g10*g7) + g9*g6*(g32 +
c     &     2*g11*g6) + g10*g2*(g33 + 2*g9*g2) + g8*g4*(g31 + 2*
c     &     g12*g8) + g12*g5*(g32 + 2*g5*g3) + g3*g1*(g33 + 2*g4*
c     &     g1))
!
!     in principal dxitril is a norm and cannot be less than zero
!     however, due to round-off it can happen anyway.
!     In order to avoid this, one would have to use the deformation
!     gradient and calculate tau_ij*tau_kl*g_ik*g_jl instead of
!     S_IJ*S_KL*C_IK*C_JL, however, this would increase the
!     calculational effort substantially.
!
      if(dxitril.lt.0.d0) then
ccc         write(*,*) '*WARNING in incplas: dxitril < 0'
         dxitril=0.d0
      else
         dxitril=dsqrt(dxitril)
      endif
!
!        restoring the hardening curves for the actual temperature
!        plconloc contains the true stresses. By multiplying by
!        the Jacobian, yiso and ykin are Kirchhoff stresses, as
!        required by the hyperelastic theory (cf. Simo, 1988).
!
      niso=int(plconloc(801))
      nkin=int(plconloc(802))
      if(niso.ne.0) then
         do i=1,niso
            xiso(i)=plconloc(2*i-1)
            yiso(i)=vj*plconloc(2*i)
         enddo
      endif
      if(nkin.ne.0) then
         do i=1,nkin
            xkin(i)=plconloc(399+2*i)
            ykin(i)=vj*plconloc(400+2*i)
         enddo
      endif
!
!     if no viscous calculation is performed a pure creep calculation
!     (without plasticity) is reduced to an elastic calculation
!
      ielastic=0
      if(ivisco.eq.0) then
         if(niso.eq.2) then
            if((dabs(yiso(1)).lt.1.d-10).and.
     &         (dabs(yiso(2)).lt.1.d-10)) then
               ielastic=1
            endif
         endif
      endif
!
!     check for yielding
!
      if(user_hardening.eq.1) then
         call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
     &        fiso,dfiso,fkin,dfkin)
         fiso=fiso*vj
      else
         if(niso.ne.0) then
            call ident(xiso,ep,niso,id)
            if(id.eq.0) then
               fiso=yiso(1)
            elseif(id.eq.niso) then
               fiso=yiso(niso)
            else
               dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
               fiso=yiso(id)+dfiso*(ep-xiso(id))
            endif
         elseif(nkin.ne.0) then
            fiso=ykin(1)
         else
            fiso=0.d0
         endif
      endif
!
      ftrial=dxitril-dsqrt(2.d0/3.d0)*fiso
      if((ftrial.le.1.d-10).or.(ielas.eq.1).or.(ielastic.eq.1)
     &    .or.(dtime.lt.1.d-30)) then
!
!        no plastic deformation
!        beta contains the Cauchy residual stresses
!
c         write(*,*) 'no plastic deformation'
         c8=xk*(vj2-1.d0)/2.d0
!
!           residual stresses are de facto PK2 stresses
!           (Piola-Kirchhoff of the second kind)
!
         stre(1)=c8*ci(1)+stril(1)-beta(1)
         stre(2)=c8*ci(2)+stril(2)-beta(2)
         stre(3)=c8*ci(3)+stril(3)-beta(3)
         stre(4)=c8*ci(4)+stril(4)-beta(4)
         stre(5)=c8*ci(5)+stril(5)-beta(5)
         stre(6)=c8*ci(6)+stril(6)-beta(6)
!
!        updating the plastic fields
!
         do i=1,3
            cpl(i)=cpl(i)-1.d0
         enddo
         do i=1,6
            xstate(1+i,iint,iel)=-cpl(i)/2.d0
            xstate(7+i,iint,iel)=stbl(i)
         enddo
         xstate(1,iint,iel)=ep
!
         if(icmd.ne.3) then
!
            umb=um*umb
!
!           calculating the local stiffness matrix
!
            xg(1,1)=(c(2)*c(3)-c(6)*c(6))/vj2
            xg(2,2)=(c(1)*c(3)-c(5)*c(5))/vj2
            xg(3,3)=(c(1)*c(2)-c(4)*c(4))/vj2
            xg(1,2)=(c(5)*c(6)-c(4)*c(3))/vj2
            xg(1,3)=(c(4)*c(6)-c(2)*c(5))/vj2
            xg(2,3)=(c(4)*c(5)-c(1)*c(6))/vj2
            xg(2,1)=xg(1,2)
            xg(3,1)=xg(1,3)
            xg(3,2)=xg(2,3)
!
            xs(1,1)=stril(1)
            xs(2,2)=stril(2)
            xs(3,3)=stril(3)
            xs(1,2)=stril(4)
            xs(2,1)=stril(4)
            xs(1,3)=stril(5)
            xs(3,1)=stril(5)
            xs(2,3)=stril(6)
            xs(3,2)=stril(6)
!
            elas(1)=umb*(xg(1,1)*xg(1,1)+xg(1,1)*xg(1,1)-
     &           2.d0*xg(1,1)*xg(1,1)/3.d0)
     &           -2.d0*(xs(1,1)*xg(1,1)+xg(1,1)*xs(1,1))/3.d0
     &           +xk*vj2*xg(1,1)*xg(1,1)
     &           -xk*(vj2-1.d0)*(xg(1,1)*xg(1,1)
     &           +xg(1,1)*xg(1,1))/2.d0
            elas(2)=umb*(xg(1,2)*xg(1,2)+xg(1,2)*xg(1,2)-
     &           2.d0*xg(1,1)*xg(2,2)/3.d0)
     &           -2.d0*(xs(1,1)*xg(2,2)+xg(1,1)*xs(2,2))/3.d0
     &           +xk*vj2*xg(1,1)*xg(2,2)
     &           -xk*(vj2-1.d0)*(xg(1,2)*xg(1,2)
     &           +xg(1,2)*xg(1,2))/2.d0
            elas(3)=umb*(xg(2,2)*xg(2,2)+xg(2,2)*xg(2,2)-
     &           2.d0*xg(2,2)*xg(2,2)/3.d0)
     &           -2.d0*(xs(2,2)*xg(2,2)+xg(2,2)*xs(2,2))/3.d0
     &           +xk*vj2*xg(2,2)*xg(2,2)
     &           -xk*(vj2-1.d0)*(xg(2,2)*xg(2,2)
     &           +xg(2,2)*xg(2,2))/2.d0
            elas(4)=umb*(xg(1,3)*xg(1,3)+xg(1,3)*xg(1,3)-
     &           2.d0*xg(1,1)*xg(3,3)/3.d0)
     &           -2.d0*(xs(1,1)*xg(3,3)+xg(1,1)*xs(3,3))/3.d0
     &           +xk*vj2*xg(1,1)*xg(3,3)
     &           -xk*(vj2-1.d0)*(xg(1,3)*xg(1,3)
     &           +xg(1,3)*xg(1,3))/2.d0
            elas(5)=umb*(xg(2,3)*xg(2,3)+xg(2,3)*xg(2,3)-
     &           2.d0*xg(2,2)*xg(3,3)/3.d0)
     &           -2.d0*(xs(2,2)*xg(3,3)+xg(2,2)*xs(3,3))/3.d0
     &           +xk*vj2*xg(2,2)*xg(3,3)
     &           -xk*(vj2-1.d0)*(xg(2,3)*xg(2,3)
     &           +xg(2,3)*xg(2,3))/2.d0
            elas(6)=umb*(xg(3,3)*xg(3,3)+xg(3,3)*xg(3,3)-
     &           2.d0*xg(3,3)*xg(3,3)/3.d0)
     &           -2.d0*(xs(3,3)*xg(3,3)+xg(3,3)*xs(3,3))/3.d0
     &           +xk*vj2*xg(3,3)*xg(3,3)
     &           -xk*(vj2-1.d0)*(xg(3,3)*xg(3,3)
     &           +xg(3,3)*xg(3,3))/2.d0
            elas(7)=umb*(xg(1,1)*xg(1,2)+xg(1,2)*xg(1,1)-
     &           2.d0*xg(1,1)*xg(1,2)/3.d0)
     &           -2.d0*(xs(1,1)*xg(1,2)+xg(1,1)*xs(1,2))/3.d0
     &           +xk*vj2*xg(1,1)*xg(1,2)
     &           -xk*(vj2-1.d0)*(xg(1,1)*xg(1,2)
     &           +xg(1,2)*xg(1,1))/2.d0
            elas(8)=umb*(xg(2,1)*xg(2,2)+xg(2,2)*xg(2,1)-
     &           2.d0*xg(2,2)*xg(1,2)/3.d0)
     &           -2.d0*(xs(2,2)*xg(1,2)+xg(2,2)*xs(1,2))/3.d0
     &           +xk*vj2*xg(2,2)*xg(1,2)
     &           -xk*(vj2-1.d0)*(xg(2,1)*xg(2,2)
     &           +xg(2,2)*xg(2,1))/2.d0
            elas(9)=umb*(xg(3,1)*xg(3,2)+xg(3,2)*xg(3,1)-
     &           2.d0*xg(3,3)*xg(1,2)/3.d0)
     &           -2.d0*(xs(3,3)*xg(1,2)+xg(3,3)*xs(1,2))/3.d0
     &           +xk*vj2*xg(3,3)*xg(1,2)
     &           -xk*(vj2-1.d0)*(xg(3,1)*xg(3,2)
     &           +xg(3,2)*xg(3,1))/2.d0
            elas(10)=umb*(xg(1,1)*xg(2,2)+xg(1,2)*xg(2,1)-
     &           2.d0*xg(1,2)*xg(1,2)/3.d0)
     &           -2.d0*(xs(1,2)*xg(1,2)+xg(1,2)*xs(1,2))/3.d0
     &           +xk*vj2*xg(1,2)*xg(1,2)
     &           -xk*(vj2-1.d0)*(xg(1,1)*xg(2,2)
     &           +xg(1,2)*xg(2,1))/2.d0
            elas(11)=umb*(xg(1,1)*xg(1,3)+xg(1,3)*xg(1,1)-
     &           2.d0*xg(1,1)*xg(1,3)/3.d0)
     &           -2.d0*(xs(1,1)*xg(1,3)+xg(1,1)*xs(1,3))/3.d0
     &           +xk*vj2*xg(1,1)*xg(1,3)
     &           -xk*(vj2-1.d0)*(xg(1,1)*xg(1,3)
     &           +xg(1,3)*xg(1,1))/2.d0
            elas(12)=umb*(xg(2,1)*xg(2,3)+xg(2,3)*xg(2,1)-
     &           2.d0*xg(2,2)*xg(1,3)/3.d0)
     &           -2.d0*(xs(2,2)*xg(1,3)+xg(2,2)*xs(1,3))/3.d0
     &           +xk*vj2*xg(2,2)*xg(1,3)
     &           -xk*(vj2-1.d0)*(xg(2,1)*xg(2,3)
     &           +xg(2,3)*xg(2,1))/2.d0
            elas(13)=umb*(xg(3,1)*xg(3,3)+xg(3,3)*xg(3,1)-
     &           2.d0*xg(3,3)*xg(1,3)/3.d0)
     &           -2.d0*(xs(3,3)*xg(1,3)+xg(3,3)*xs(1,3))/3.d0
     &           +xk*vj2*xg(3,3)*xg(1,3)
     &           -xk*(vj2-1.d0)*(xg(3,1)*xg(3,3)
     &           +xg(3,3)*xg(3,1))/2.d0
            elas(14)=umb*(xg(1,1)*xg(2,3)+xg(1,3)*xg(2,1)-
     &           2.d0*xg(1,2)*xg(1,3)/3.d0)
     &           -2.d0*(xs(1,2)*xg(1,3)+xg(1,2)*xs(1,3))/3.d0
     &           +xk*vj2*xg(1,2)*xg(1,3)
     &           -xk*(vj2-1.d0)*(xg(1,1)*xg(2,3)
     &           +xg(1,3)*xg(2,1))/2.d0
            elas(15)=umb*(xg(1,1)*xg(3,3)+xg(1,3)*xg(3,1)-
     &           2.d0*xg(1,3)*xg(1,3)/3.d0)
     &           -2.d0*(xs(1,3)*xg(1,3)+xg(1,3)*xs(1,3))/3.d0
     &           +xk*vj2*xg(1,3)*xg(1,3)
     &           -xk*(vj2-1.d0)*(xg(1,1)*xg(3,3)
     &           +xg(1,3)*xg(3,1))/2.d0
            elas(16)=umb*(xg(1,2)*xg(1,3)+xg(1,3)*xg(1,2)-
     &           2.d0*xg(1,1)*xg(2,3)/3.d0)
     &           -2.d0*(xs(1,1)*xg(2,3)+xg(1,1)*xs(2,3))/3.d0
     &           +xk*vj2*xg(1,1)*xg(2,3)
     &           -xk*(vj2-1.d0)*(xg(1,2)*xg(1,3)
     &           +xg(1,3)*xg(1,2))/2.d0
            elas(17)=umb*(xg(2,2)*xg(2,3)+xg(2,3)*xg(2,2)-
     &           2.d0*xg(2,2)*xg(2,3)/3.d0)
     &           -2.d0*(xs(2,2)*xg(2,3)+xg(2,2)*xs(2,3))/3.d0
     &           +xk*vj2*xg(2,2)*xg(2,3)
     &           -xk*(vj2-1.d0)*(xg(2,2)*xg(2,3)
     &           +xg(2,3)*xg(2,2))/2.d0
            elas(18)=umb*(xg(3,2)*xg(3,3)+xg(3,3)*xg(3,2)-
     &           2.d0*xg(3,3)*xg(2,3)/3.d0)
     &           -2.d0*(xs(3,3)*xg(2,3)+xg(3,3)*xs(2,3))/3.d0
     &           +xk*vj2*xg(3,3)*xg(2,3)
     &           -xk*(vj2-1.d0)*(xg(3,2)*xg(3,3)
     &           +xg(3,3)*xg(3,2))/2.d0
            elas(19)=umb*(xg(1,2)*xg(2,3)+xg(1,3)*xg(2,2)-
     &           2.d0*xg(1,2)*xg(2,3)/3.d0)
     &           -2.d0*(xs(1,2)*xg(2,3)+xg(1,2)*xs(2,3))/3.d0
     &           +xk*vj2*xg(1,2)*xg(2,3)
     &           -xk*(vj2-1.d0)*(xg(1,2)*xg(2,3)
     &           +xg(1,3)*xg(2,2))/2.d0
            elas(20)=umb*(xg(1,2)*xg(3,3)+xg(1,3)*xg(3,2)-
     &           2.d0*xg(1,3)*xg(2,3)/3.d0)
     &           -2.d0*(xs(1,3)*xg(2,3)+xg(1,3)*xs(2,3))/3.d0
     &           +xk*vj2*xg(1,3)*xg(2,3)
     &           -xk*(vj2-1.d0)*(xg(1,2)*xg(3,3)
     &           +xg(1,3)*xg(3,2))/2.d0
            elas(21)=umb*(xg(2,2)*xg(3,3)+xg(2,3)*xg(3,2)-
     &           2.d0*xg(2,3)*xg(2,3)/3.d0)
     &           -2.d0*(xs(2,3)*xg(2,3)+xg(2,3)*xs(2,3))/3.d0
     &           +xk*vj2*xg(2,3)*xg(2,3)
     &           -xk*(vj2-1.d0)*(xg(2,2)*xg(3,3)
     &           +xg(2,3)*xg(3,2))/2.d0
!
         endif
!
         return
      endif
!
!        plastic deformation
!
      umb=um*umb
      umbb=umb-umbb
!
!        calculating the consistency parameter
!
      c1=2.d0/3.d0
      c2=dsqrt(c1)
      c3=c1/um
      c4=c2/um
!
      iloop=0
      cop=0.d0
!
      loop: do
         iloop=iloop+1
         ep=epini+c2*cop
!
         if(user_hardening.eq.1) then
            call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
     &           fiso,dfiso,fkin,dfkin)
            fiso=fiso*vj
            dfiso=dfiso*vj
            fkin=fkin*vj
            dfkin=dfkin*vj
         else
            if(niso.ne.0) then
               call ident(xiso,ep,niso,id)
               if(id.eq.0) then
                  fiso=yiso(1)
                  dfiso=0.d0
               elseif(id.eq.niso) then
                  fiso=yiso(niso)
                  dfiso=0.d0
               else
                  dfiso=(yiso(id+1)-yiso(id))/(xiso(id+1)-xiso(id))
                  fiso=yiso(id)+dfiso*(ep-xiso(id))
               endif
            elseif(nkin.ne.0) then
               fiso=ykin(1)
               dfiso=0.d0
            else
               fiso=0.d0
               dfiso=0.d0
            endif
!
            if(nkin.ne.0) then
               call ident(xkin,ep,nkin,id)
               if(id.eq.0) then
                  fkin=ykin(1)
                  dfkin=0.d0
               elseif(id.eq.nkin) then
                  fkin=ykin(nkin)
                  dfkin=0.d0
               else
                  dfkin=(ykin(id+1)-ykin(id))/(xkin(id+1)-xkin(id))
                  fkin=ykin(id)+dfkin*(ep-xkin(id))
               endif
            elseif(niso.ne.0) then
               fkin=yiso(1)
               dfkin=0.d0
            else
               fkin=0.d0
               dfkin=0.d0
            endif
         endif
!
         if(dabs(cop).lt.1.d-10) then
            fiso0=fiso
            fkin0=fkin
         endif
!
         if((kode.eq.-51).or.(ivisco.eq.0)) then
            dcop=(ftrial-c2*(fiso-fiso0)
     &           -umbb*(2.d0*cop+c4*(fkin-fkin0)))/
     &           (-c1*dfiso-umbb*(2.d0+c3*dfkin))
         else
            if(user_creep.eq.1) then
               if(ithermal(1).eq.0) then
ccc                  write(*,*) '*ERROR in incplas: no temperature defined'
                  call exit(201)
               endif
               timeabq(1)=time
               timeabq(2)=ttime+time
               qtild=(ftrial-c2*(fiso-fiso0)
     &              -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
!
!              the Von Mises stress must be positive
!
               if(qtild.lt.1.d-10) qtild=1.d-10
               ec(1)=epini
               call creep(decra,deswa,xstateini(1,iint,iel),serd,ec,
     &             esw,p,qtild,t1l,dtemp,predef,dpred,timeabq,dtime,
     &             amat,leximp,lend,pgauss,nstate_,iel,iint,layer,kspt,
     &             kstep,kinc)
               dsvm=1.d0/decra(5)
               dcop=-(decra(1)-c2*cop)/
     &                 (c2*(decra(5)*(dfiso+umbb*(3.d0+dfkin/um))+1.d0))
            else
               qtild=(ftrial-c2*(fiso-fiso0)
     &              -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
!
!              the Von Mises stress must be positive
!
               if(qtild.lt.1.d-10) qtild=1.d-10
               decra(1)=a1*qtild**xxn
               decra(5)=xxn*decra(1)/qtild
               dsvm=1.d0/decra(5)
               dcop=-(decra(1)-c2*cop)/
     &                 (c2*(decra(5)*(dfiso+umbb*(3.d0+dfkin/um))+1.d0))
            endif
         endif
         cop=cop-dcop
!
         if((dabs(dcop).lt.cop*1.d-4).or.
     &      (dabs(dcop).lt.1.d-10)) exit
!
!        check for endless loops or a negative consistency
!        parameter
!
         if((iloop.gt.15).or.(cop.le.0.d0)) then
            iloop=1
            cop=0.d0
            do
               if(iloop.gt.100) then
c                  NOTE: write statements cause problems for
c                        parallellized execution
c                  write(*,*) '*WARNING in incplas: material loop'
c                  write(*,*) '         did not converge in integration'
c                  write(*,*) '         point',iint,'in element',iel,';'
c                  write(*,*) '         the increment size is reduced'
c                  write(*,*)
                  pnewdt=0.25d0
                  return
               endif
               ep=epini+c2*cop
!
               if(user_hardening.eq.1) then
                  call uhardening(amat,iel,iint,t1l,epini,ep,dtime,
     &                 fiso,dfiso,fkin,dfkin)
                  fiso=fiso*vj
                  fkin=fkin*vj
               else
                  if(niso.ne.0) then
                     call ident(xiso,ep,niso,id)
                     if(id.eq.0) then
                        fiso=yiso(1)
                     elseif(id.eq.niso) then
                        fiso=yiso(niso)
                     else
                        dfiso=(yiso(id+1)-yiso(id))/
     &                        (xiso(id+1)-xiso(id))
                        fiso=yiso(id)+dfiso*(ep-xiso(id))
                     endif
                  elseif(nkin.ne.0) then
                     fiso=ykin(1)
                  else
                     fiso=0.d0
                  endif
!
                  if(nkin.ne.0) then
                     call ident(xkin,ep,nkin,id)
                     if(id.eq.0) then
                        fkin=ykin(1)
                     elseif(id.eq.nkin) then
                        fkin=ykin(nkin)
                     else
                        dfkin=(ykin(id+1)-ykin(id))/
     &                        (xkin(id+1)-xkin(id))
                        fkin=ykin(id)+dfkin*(ep-xkin(id))
                     endif
                  elseif(niso.ne.0) then
                     fkin=yiso(1)
                  else
                     fkin=0.d0
                  endif
               endif
!
               if(dabs(cop).lt.1.d-10) then
                  fiso0=fiso
                  fkin0=fkin
               endif
!
               if((kode.eq.-51).or.(ivisco.eq.0)) then
                  fu=(ftrial-c2*(fiso-fiso0)
     &                 -umbb*(2.d0*cop+c4*(fkin-fkin0)))
               else
                  if(user_creep.eq.1) then
                     timeabq(1)=time
                     timeabq(2)=ttime+time
                     qtild=(ftrial-c2*(fiso-fiso0)
     &                    -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
!
!                    the Von Mises stress must be positive
!
                     if(qtild.lt.1.d-10) qtild=1.d-10
                     ec(1)=epini
                     call creep(decra,deswa,xstateini(1,iint,iel),serd,
     &                    ec,esw,p,qtild,t1l,dtemp,predef,dpred,timeabq,
     &                    dtime,amat,leximp,lend,pgauss,nstate_,iel,
     &                    iint,layer,kspt,kstep,kinc)
                     dsvm=1.d0/decra(5)
                     fu=decra(1)-c2*cop
                  else
                     qtild=(ftrial-c2*(fiso-fiso0)
     &                    -umbb*(2.d0*cop+c4*(fkin-fkin0)))/(c2*vj)
!
!                    the Von Mises stress must be positive
!
                     if(qtild.lt.1.d-10) qtild=1.d-10
                     decra(1)=a1*qtild**xxn
                     decra(5)=xxn*decra(1)/qtild
                     dsvm=1.d0/decra(5)
                     fu=decra(1)-c2*cop
                  endif
               endif
!
               if(iloop.eq.1) then
c                  write(*,*) 'cop,fu ',cop,fu
                  cop1=0.d0
                  fu1=fu
                  iloop=2
                  cop=1.d-10
               elseif(iloop.eq.2) then
                  if(fu*fu1.le.0.d0) then
c                     write(*,*) cop,fu
                     iloop=3
                     fu2=fu
                     cop2=cop
                     cop=(cop1+cop2)/2.d0
                     dcop=(cop2-cop1)/2.d0
                  else
c                     write(*,*) cop,fu
                     cop=cop*10.d0
                     if(cop.gt.100.d0) then
                        write(*,*) '*ERROR: no convergence in incplas'
                        pnewdt=0.25d0
                        return
c                        call exit(201)
                     endif
                  endif
               else
c                     write(*,*) cop,fu
                  if(fu*fu1.ge.0.d0) then
                     cop1=cop
                     fu1=fu
                  else
                     cop2=cop
                     fu2=fu
                  endif
                  cop=(cop1+cop2)/2.d0
                  dcop=(cop2-cop1)/2.d0
                  if((dabs(dcop).lt.cop*1.d-4).or.
     &                 (dabs(dcop).lt.1.d-10)) exit loop
               endif
            enddo
         endif
!
      enddo loop
!
!        updating the equivalent plastic strain
!
      ep=epini+c2*cop
!
!        updating the back stress
!
      c5=2.d0*umbb*cop/dxitril
      c6=c5/(3.d0*um)
      c7=c6*dfkin*vj23
      do i=1,6
         stbl(i)=stbl(i)+c7*xitril(i)
      enddo
!
!        updating the stress
!        vj: Jacobian of the total deformation gradient
!
      c8=xk*(vj2-1.d0)/2.d0
!
      do i=1,6
         stre(i)=c8*ci(i)-beta(i)+stril(i)-c5*xitril(i)
      enddo
!
!        updating the plastic right Cauchy-Green tensor
!
      c9=c6*3.d0*vj23
      do i=1,6
         cpl(i)=cpl(i)-c9*xitril(i)
      enddo
!
      if(icmd.ne.3) then
!
!        calculating the local stiffness matrix
!
         xg(1,1)=(c(2)*c(3)-c(6)*c(6))/vj2
         xg(2,2)=(c(1)*c(3)-c(5)*c(5))/vj2
         xg(3,3)=(c(1)*c(2)-c(4)*c(4))/vj2
         xg(1,2)=(c(5)*c(6)-c(4)*c(3))/vj2
         xg(1,3)=(c(4)*c(6)-c(2)*c(5))/vj2
         xg(2,3)=(c(4)*c(5)-c(1)*c(6))/vj2
         xg(2,1)=xg(1,2)
         xg(3,1)=xg(1,3)
         xg(3,2)=xg(2,3)
!
         xs(1,1)=stril(1)
         xs(2,2)=stril(2)
         xs(3,3)=stril(3)
         xs(1,2)=stril(4)
         xs(2,1)=stril(4)
         xs(1,3)=stril(5)
         xs(3,1)=stril(5)
         xs(2,3)=stril(6)
         xs(3,2)=stril(6)
!
         f0=2.d0*umbb*cop/dxitril
         d0=1.d0+(dfkin/um+dfiso/umbb)/3.d0
!
!        creep contribution
!
         if((kode.eq.-52).and.(ivisco.eq.1)) then
            d0=d0+dsvm/(3.d0*umbb)
         endif
!
         f1=1.d0/d0-f0
         d1=2.d0*f1*umbb-((1.d0+dfkin/(3.d0*um))/d0-1.d0)*
     &      4.d0*cop*dxitril/3.d0
         d2=2d0*dxitril*f1
!
         xx(1,1)=xitril(1)
         xx(2,2)=xitril(2)
         xx(3,3)=xitril(3)
         xx(1,2)=xitril(4)
         xx(2,1)=xitril(4)
         xx(1,3)=xitril(5)
         xx(3,1)=xitril(5)
         xx(2,3)=xitril(6)
         xx(3,2)=xitril(6)
!
         xn(1,1)=xitril(1)/dxitril
         xn(2,2)=xitril(2)/dxitril
         xn(3,3)=xitril(3)/dxitril
         xn(1,2)=xitril(4)/dxitril
         xn(2,1)=xitril(4)/dxitril
         xn(1,3)=xitril(5)/dxitril
         xn(3,1)=xitril(5)/dxitril
         xn(2,3)=xitril(6)/dxitril
         xn(3,2)=xitril(6)/dxitril
!
         do i=1,3
            do j=i,3
               xd(i,j)=xn(i,1)*xn(1,j)*c(1)+xn(i,1)*xn(2,j)*c(4)+
     &                 xn(i,1)*xn(3,j)*c(5)+xn(i,2)*xn(1,j)*c(4)+
     &                 xn(i,2)*xn(2,j)*c(2)+xn(i,2)*xn(3,j)*c(6)+
     &                 xn(i,3)*xn(1,j)*c(5)+xn(i,3)*xn(2,j)*c(6)+
     &                 xn(i,3)*xn(3,j)*c(3)
            enddo
         enddo
         xd(2,1)=xd(1,2)
         xd(3,1)=xd(1,3)
         xd(3,2)=xd(2,3)
!
!        deviatoric part
!
         c1=(xd(1,1)*c(1)+xd(2,2)*c(2)+xd(3,3)*c(3)+
     &      2.d0*(xd(1,2)*c(4)+xd(1,3)*c(5)+xd(2,3)*c(6)))/3.d0
         do i=1,3
            do j=i,3
               xd(i,j)=xd(i,j)-c1*xg(i,j)
            enddo
         enddo
         xd(2,1)=xd(1,2)
         xd(3,1)=xd(1,3)
         xd(3,2)=xd(2,3)
!
         elas(1)=(umb-f0*umbb)*(xg(1,1)*xg(1,1)+xg(1,1)*xg(1,1)-
     &        2.d0*xg(1,1)*xg(1,1)/3.d0)
     &        -2.d0*(xs(1,1)*xg(1,1)+xg(1,1)*xs(1,1))/3.d0
     &        +f0*2.d0*(xx(1,1)*xg(1,1)+xg(1,1)*xx(1,1))/3.d0
     &        -d1*xn(1,1)*xn(1,1)-d2*(xn(1,1)*xd(1,1)+
     &        xd(1,1)*xn(1,1))/2.d0+xk*vj2*xg(1,1)*xg(1,1)
     &        -xk*(vj2-1.d0)*(xg(1,1)*xg(1,1)+xg(1,1)*xg(1,1))/2.d0
         elas(2)=(umb-f0*umbb)*(xg(1,2)*xg(1,2)+xg(1,2)*xg(1,2)-
     &        2.d0*xg(1,1)*xg(2,2)/3.d0)
     &        -2.d0*(xs(1,1)*xg(2,2)+xg(1,1)*xs(2,2))/3.d0
     &        +f0*2.d0*(xx(1,1)*xg(2,2)+xg(1,1)*xx(2,2))/3.d0
     &        -d1*xn(1,1)*xn(2,2)-d2*(xn(1,1)*xd(2,2)+
     &        xd(1,1)*xn(2,2))/2.d0+xk*vj2*xg(1,1)*xg(2,2)
     &        -xk*(vj2-1.d0)*(xg(1,2)*xg(1,2)+xg(1,2)*xg(1,2))/2.d0
         elas(3)=(umb-f0*umbb)*(xg(2,2)*xg(2,2)+xg(2,2)*xg(2,2)-
     &        2.d0*xg(2,2)*xg(2,2)/3.d0)
     &        -2.d0*(xs(2,2)*xg(2,2)+xg(2,2)*xs(2,2))/3.d0
     &        +f0*2.d0*(xx(2,2)*xg(2,2)+xg(2,2)*xx(2,2))/3.d0
     &        -d1*xn(2,2)*xn(2,2)-d2*(xn(2,2)*xd(2,2)+
     &        xd(2,2)*xn(2,2))/2.d0+xk*vj2*xg(2,2)*xg(2,2)
     &        -xk*(vj2-1.d0)*(xg(2,2)*xg(2,2)+xg(2,2)*xg(2,2))/2.d0
         elas(4)=(umb-f0*umbb)*(xg(1,3)*xg(1,3)+xg(1,3)*xg(1,3)-
     &        2.d0*xg(1,1)*xg(3,3)/3.d0)
     &        -2.d0*(xs(1,1)*xg(3,3)+xg(1,1)*xs(3,3))/3.d0
     &        +f0*2.d0*(xx(1,1)*xg(3,3)+xg(1,1)*xx(3,3))/3.d0
     &        -d1*xn(1,1)*xn(3,3)-d2*(xn(1,1)*xd(3,3)+
     &        xd(1,1)*xn(3,3))/2.d0+xk*vj2*xg(1,1)*xg(3,3)
     &        -xk*(vj2-1.d0)*(xg(1,3)*xg(1,3)+xg(1,3)*xg(1,3))/2.d0
         elas(5)=(umb-f0*umbb)*(xg(2,3)*xg(2,3)+xg(2,3)*xg(2,3)-
     &        2.d0*xg(2,2)*xg(3,3)/3.d0)
     &        -2.d0*(xs(2,2)*xg(3,3)+xg(2,2)*xs(3,3))/3.d0
     &        +f0*2.d0*(xx(2,2)*xg(3,3)+xg(2,2)*xx(3,3))/3.d0
     &        -d1*xn(2,2)*xn(3,3)-d2*(xn(2,2)*xd(3,3)+
     &        xd(2,2)*xn(3,3))/2.d0+xk*vj2*xg(2,2)*xg(3,3)
     &        -xk*(vj2-1.d0)*(xg(2,3)*xg(2,3)+xg(2,3)*xg(2,3))/2.d0
         elas(6)=(umb-f0*umbb)*(xg(3,3)*xg(3,3)+xg(3,3)*xg(3,3)-
     &        2.d0*xg(3,3)*xg(3,3)/3.d0)
     &        -2.d0*(xs(3,3)*xg(3,3)+xg(3,3)*xs(3,3))/3.d0
     &        +f0*2.d0*(xx(3,3)*xg(3,3)+xg(3,3)*xx(3,3))/3.d0
     &        -d1*xn(3,3)*xn(3,3)-d2*(xn(3,3)*xd(3,3)+
     &        xd(3,3)*xn(3,3))/2.d0+xk*vj2*xg(3,3)*xg(3,3)
     &        -xk*(vj2-1.d0)*(xg(3,3)*xg(3,3)+xg(3,3)*xg(3,3))/2.d0
         elas(7)=(umb-f0*umbb)*(xg(1,1)*xg(1,2)+xg(1,2)*xg(1,1)-
     &        2.d0*xg(1,1)*xg(1,2)/3.d0)
     &        -2.d0*(xs(1,1)*xg(1,2)+xg(1,1)*xs(1,2))/3.d0
     &        +f0*2.d0*(xx(1,1)*xg(1,2)+xg(1,1)*xx(1,2))/3.d0
     &        -d1*xn(1,1)*xn(1,2)-d2*(xn(1,1)*xd(1,2)+
     &        xd(1,1)*xn(1,2))/2.d0+xk*vj2*xg(1,1)*xg(1,2)
     &        -xk*(vj2-1.d0)*(xg(1,1)*xg(1,2)+xg(1,2)*xg(1,1))/2.d0
         elas(8)=(umb-f0*umbb)*(xg(2,1)*xg(2,2)+xg(2,2)*xg(2,1)-
     &        2.d0*xg(2,2)*xg(1,2)/3.d0)
     &        -2.d0*(xs(2,2)*xg(1,2)+xg(2,2)*xs(1,2))/3.d0
     &        +f0*2.d0*(xx(2,2)*xg(1,2)+xg(2,2)*xx(1,2))/3.d0
     &        -d1*xn(2,2)*xn(1,2)-d2*(xn(2,2)*xd(1,2)+
     &        xd(2,2)*xn(1,2))/2.d0+xk*vj2*xg(2,2)*xg(1,2)
     &        -xk*(vj2-1.d0)*(xg(2,1)*xg(2,2)+xg(2,2)*xg(2,1))/2.d0
         elas(9)=(umb-f0*umbb)*(xg(3,1)*xg(3,2)+xg(3,2)*xg(3,1)-
     &        2.d0*xg(3,3)*xg(1,2)/3.d0)
     &        -2.d0*(xs(3,3)*xg(1,2)+xg(3,3)*xs(1,2))/3.d0
     &        +f0*2.d0*(xx(3,3)*xg(1,2)+xg(3,3)*xx(1,2))/3.d0
     &        -d1*xn(3,3)*xn(1,2)-d2*(xn(3,3)*xd(1,2)+
     &        xd(3,3)*xn(1,2))/2.d0+xk*vj2*xg(3,3)*xg(1,2)
     &        -xk*(vj2-1.d0)*(xg(3,1)*xg(3,2)+xg(3,2)*xg(3,1))/2.d0
         elas(10)=(umb-f0*umbb)*(xg(1,1)*xg(2,2)+xg(1,2)*xg(2,1)-
     &        2.d0*xg(1,2)*xg(1,2)/3.d0)
     &        -2.d0*(xs(1,2)*xg(1,2)+xg(1,2)*xs(1,2))/3.d0
     &        +f0*2.d0*(xx(1,2)*xg(1,2)+xg(1,2)*xx(1,2))/3.d0
     &        -d1*xn(1,2)*xn(1,2)-d2*(xn(1,2)*xd(1,2)+
     &        xd(1,2)*xn(1,2))/2.d0+xk*vj2*xg(1,2)*xg(1,2)
     &        -xk*(vj2-1.d0)*(xg(1,1)*xg(2,2)+xg(1,2)*xg(2,1))/2.d0
         elas(11)=(umb-f0*umbb)*(xg(1,1)*xg(1,3)+xg(1,3)*xg(1,1)-
     &        2.d0*xg(1,1)*xg(1,3)/3.d0)
     &        -2.d0*(xs(1,1)*xg(1,3)+xg(1,1)*xs(1,3))/3.d0
     &        +f0*2.d0*(xx(1,1)*xg(1,3)+xg(1,1)*xx(1,3))/3.d0
     &        -d1*xn(1,1)*xn(1,3)-d2*(xn(1,1)*xd(1,3)+
     &        xd(1,1)*xn(1,3))/2.d0+xk*vj2*xg(1,1)*xg(1,3)
     &        -xk*(vj2-1.d0)*(xg(1,1)*xg(1,3)+xg(1,3)*xg(1,1))/2.d0
         elas(12)=(umb-f0*umbb)*(xg(2,1)*xg(2,3)+xg(2,3)*xg(2,1)-
     &        2.d0*xg(2,2)*xg(1,3)/3.d0)
     &        -2.d0*(xs(2,2)*xg(1,3)+xg(2,2)*xs(1,3))/3.d0
     &        +f0*2.d0*(xx(2,2)*xg(1,3)+xg(2,2)*xx(1,3))/3.d0
     &        -d1*xn(2,2)*xn(1,3)-d2*(xn(2,2)*xd(1,3)+
     &        xd(2,2)*xn(1,3))/2.d0+xk*vj2*xg(2,2)*xg(1,3)
     &        -xk*(vj2-1.d0)*(xg(2,1)*xg(2,3)+xg(2,3)*xg(2,1))/2.d0
         elas(13)=(umb-f0*umbb)*(xg(3,1)*xg(3,3)+xg(3,3)*xg(3,1)-
     &        2.d0*xg(3,3)*xg(1,3)/3.d0)
     &        -2.d0*(xs(3,3)*xg(1,3)+xg(3,3)*xs(1,3))/3.d0
     &        +f0*2.d0*(xx(3,3)*xg(1,3)+xg(3,3)*xx(1,3))/3.d0
     &        -d1*xn(3,3)*xn(1,3)-d2*(xn(3,3)*xd(1,3)+
     &        xd(3,3)*xn(1,3))/2.d0+xk*vj2*xg(3,3)*xg(1,3)
     &        -xk*(vj2-1.d0)*(xg(3,1)*xg(3,3)+xg(3,3)*xg(3,1))/2.d0
         elas(14)=(umb-f0*umbb)*(xg(1,1)*xg(2,3)+xg(1,3)*xg(2,1)-
     &        2.d0*xg(1,2)*xg(1,3)/3.d0)
     &        -2.d0*(xs(1,2)*xg(1,3)+xg(1,2)*xs(1,3))/3.d0
     &        +f0*2.d0*(xx(1,2)*xg(1,3)+xg(1,2)*xx(1,3))/3.d0
     &        -d1*xn(1,2)*xn(1,3)-d2*(xn(1,2)*xd(1,3)+
     &        xd(1,2)*xn(1,3))/2.d0+xk*vj2*xg(1,2)*xg(1,3)
     &        -xk*(vj2-1.d0)*(xg(1,1)*xg(2,3)+xg(1,3)*xg(2,1))/2.d0
         elas(15)=(umb-f0*umbb)*(xg(1,1)*xg(3,3)+xg(1,3)*xg(3,1)-
     &        2.d0*xg(1,3)*xg(1,3)/3.d0)
     &        -2.d0*(xs(1,3)*xg(1,3)+xg(1,3)*xs(1,3))/3.d0
     &        +f0*2.d0*(xx(1,3)*xg(1,3)+xg(1,3)*xx(1,3))/3.d0
     &        -d1*xn(1,3)*xn(1,3)-d2*(xn(1,3)*xd(1,3)+
     &        xd(1,3)*xn(1,3))/2.d0+xk*vj2*xg(1,3)*xg(1,3)
     &        -xk*(vj2-1.d0)*(xg(1,1)*xg(3,3)+xg(1,3)*xg(3,1))/2.d0
         elas(16)=(umb-f0*umbb)*(xg(1,2)*xg(1,3)+xg(1,3)*xg(1,2)-
     &        2.d0*xg(1,1)*xg(2,3)/3.d0)
     &        -2.d0*(xs(1,1)*xg(2,3)+xg(1,1)*xs(2,3))/3.d0
     &        +f0*2.d0*(xx(1,1)*xg(2,3)+xg(1,1)*xx(2,3))/3.d0
     &        -d1*xn(1,1)*xn(2,3)-d2*(xn(1,1)*xd(2,3)+
     &        xd(1,1)*xn(2,3))/2.d0+xk*vj2*xg(1,1)*xg(2,3)
     &        -xk*(vj2-1.d0)*(xg(1,2)*xg(1,3)+xg(1,3)*xg(1,2))/2.d0
         elas(17)=(umb-f0*umbb)*(xg(2,2)*xg(2,3)+xg(2,3)*xg(2,2)-
     &        2.d0*xg(2,2)*xg(2,3)/3.d0)
     &        -2.d0*(xs(2,2)*xg(2,3)+xg(2,2)*xs(2,3))/3.d0
     &        +f0*2.d0*(xx(2,2)*xg(2,3)+xg(2,2)*xx(2,3))/3.d0
     &        -d1*xn(2,2)*xn(2,3)-d2*(xn(2,2)*xd(2,3)+
     &        xd(2,2)*xn(2,3))/2.d0+xk*vj2*xg(2,2)*xg(2,3)
     &        -xk*(vj2-1.d0)*(xg(2,2)*xg(2,3)+xg(2,3)*xg(2,2))/2.d0
         elas(18)=(umb-f0*umbb)*(xg(3,2)*xg(3,3)+xg(3,3)*xg(3,2)-
     &        2.d0*xg(3,3)*xg(2,3)/3.d0)
     &        -2.d0*(xs(3,3)*xg(2,3)+xg(3,3)*xs(2,3))/3.d0
     &        +f0*2.d0*(xx(3,3)*xg(2,3)+xg(3,3)*xx(2,3))/3.d0
     &        -d1*xn(3,3)*xn(2,3)-d2*(xn(3,3)*xd(2,3)+
     &        xd(3,3)*xn(2,3))/2.d0+xk*vj2*xg(3,3)*xg(2,3)
     &        -xk*(vj2-1.d0)*(xg(3,2)*xg(3,3)+xg(3,3)*xg(3,2))/2.d0
         elas(19)=(umb-f0*umbb)*(xg(1,2)*xg(2,3)+xg(1,3)*xg(2,2)-
     &        2.d0*xg(1,2)*xg(2,3)/3.d0)
     &        -2.d0*(xs(1,2)*xg(2,3)+xg(1,2)*xs(2,3))/3.d0
     &        +f0*2.d0*(xx(1,2)*xg(2,3)+xg(1,2)*xx(2,3))/3.d0
     &        -d1*xn(1,2)*xn(2,3)-d2*(xn(1,2)*xd(2,3)+
     &        xd(1,2)*xn(2,3))/2.d0+xk*vj2*xg(1,2)*xg(2,3)
     &        -xk*(vj2-1.d0)*(xg(1,2)*xg(2,3)+xg(1,3)*xg(2,2))/2.d0
         elas(20)=(umb-f0*umbb)*(xg(1,2)*xg(3,3)+xg(1,3)*xg(3,2)-
     &        2.d0*xg(1,3)*xg(2,3)/3.d0)
     &        -2.d0*(xs(1,3)*xg(2,3)+xg(1,3)*xs(2,3))/3.d0
     &        +f0*2.d0*(xx(1,3)*xg(2,3)+xg(1,3)*xx(2,3))/3.d0
     &        -d1*xn(1,3)*xn(2,3)-d2*(xn(1,3)*xd(2,3)+
     &        xd(1,3)*xn(2,3))/2.d0+xk*vj2*xg(1,3)*xg(2,3)
     &        -xk*(vj2-1.d0)*(xg(1,2)*xg(3,3)+xg(1,3)*xg(3,2))/2.d0
         elas(21)=(umb-f0*umbb)*(xg(2,2)*xg(3,3)+xg(2,3)*xg(3,2)-
     &        2.d0*xg(2,3)*xg(2,3)/3.d0)
     &        -2.d0*(xs(2,3)*xg(2,3)+xg(2,3)*xs(2,3))/3.d0
     &        +f0*2.d0*(xx(2,3)*xg(2,3)+xg(2,3)*xx(2,3))/3.d0
     &        -d1*xn(2,3)*xn(2,3)-d2*(xn(2,3)*xd(2,3)+
     &        xd(2,3)*xn(2,3))/2.d0+xk*vj2*xg(2,3)*xg(2,3)
     &        -xk*(vj2-1.d0)*(xg(2,2)*xg(3,3)+xg(2,3)*xg(3,2))/2.d0
!
      endif
!
!        updating the plastic fields
!
      do i=1,3
         cpl(i)=cpl(i)-1.d0
      enddo
      do i=1,6
         xstate(1+i,iint,iel)=-cpl(i)/2.d0
         xstate(7+i,iint,iel)=stbl(i)
      enddo
      xstate(1,iint,iel)=ep
!
!     maximum difference in the equivalent viscoplastic strain
!     in this increment based on the viscoplastic strain rate at
!     the start and the end of the increment
!
      if((kode.eq.-52).and.(ivisco.eq.1)) then
c      if(ivisco.eq.1) then
         depvisc=max(depvisc,abs(decra(1)-dtime*xstateini(14,iint,iel)))
         xstate(14,iint,iel)=decra(1)/dtime
      endif
c      depvisc=max(depvisc,ep-epini)
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      FUNCTION torflux_deriv (x)
      USE stel_kinds
      USE vmec_input, ONLY: lRFP
      USE vmec_input, ONLY: tf => aphi
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec) :: x, torflux_deriv
      REAL(rprec), EXTERNAL :: polflux_deriv, piota
      INTEGER     :: i
C-----------------------------------------------
      IF (lRFP) THEN
!        RFP/TOKAMAK
         torflux_deriv = polflux_deriv(x)/piota(x)

      ELSE
!        TOKAMAK/STELLARATOR (default is tf(1) = 1)
         torflux_deriv = 0
         DO i = UBOUND(tf,1), LBOUND(tf,1), -1
            torflux_deriv = x*torflux_deriv + i*tf(i)
         END DO
!        torflux_deriv = 1
      END IF

      END FUNCTION torflux_deriv

      FUNCTION polflux_deriv (x)
      USE stel_kinds
      USE vmec_input, ONLY: lRFP
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec) :: x, polflux_deriv
      REAL(rprec), EXTERNAL :: torflux_deriv, piota
C-----------------------------------------------
      IF (lRFP) THEN
!        RFP/TOKAMAK
         polflux_deriv = 1

      ELSE
!        TOKAMAK/STELLARATOR: dchi/ds = iota * dphi/ds
         polflux_deriv = piota(x)*torflux_deriv(x)
      END IF

      END FUNCTION polflux_deriv

      FUNCTION torflux (x)
      USE stel_kinds
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec) :: x, torflux
      REAL(rprec) :: h, xi
      REAL(rprec), EXTERNAL :: torflux_deriv
      INTEGER     :: i
C-----------------------------------------------
      h = 1.E-2_dp*x
      torflux = 0
      DO i=1,101
         xi = (i-1)*h
         torflux = torflux + torflux_deriv(xi)
      END DO
      torflux = torflux-0.5_dp*(torflux_deriv(0._dp)+torflux_deriv(x))
      torflux = h*torflux

      END FUNCTION torflux

      FUNCTION polflux (x)
      USE stel_kinds
!     USE vmec_input, ONLY: af => achi
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      REAL(rprec) :: x, polflux
      REAL(rprec) :: h, xi
      REAL(rprec), EXTERNAL :: polflux_deriv
      INTEGER     :: i
C-----------------------------------------------
      h = 1.E-2_dp*x
      polflux = 0
      DO i=1,101
         xi = (i-1)*h
         polflux = polflux + polflux_deriv(xi)
      END DO
      polflux = polflux-0.5_dp*(polflux_deriv(0._dp)+polflux_deriv(x))
      polflux = h*polflux

      END FUNCTION polflux

!  function piota moved to a separate file, piota.f. J Hanson, 2010-03-16

ext: f


text:
!! Convenient definitions of numerical constants.  This should not
!! be considered an exhaustive list of such constants.
module bl_constants_module

  use bl_types
  use amrex_constants_module

  implicit none

end module bl_constants_module

max_stars_repo_name: ylunalin/amrex

max_issues_repo_name: ylunalin/amrex

max_forks_repo_name: ylunalin/amrex


! (C) Copyright 2017-2020 UCAR
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.

module fv3jedi_traj_mod

! fv3-jedi-lm
use fv3jedi_lm_utils_mod, only: fv3jedi_traj => fv3jedi_lm_traj, deallocate_traj

! fv3-jedi
use fv3jedi_kinds_mod,    only: kind_real
use fv3jedi_state_mod,    only: fv3jedi_state

! --------------------------------------------------------------------------------------------------

implicit none
private
public :: fv3jedi_traj
public :: set
public :: wipe

! --------------------------------------------------------------------------------------------------

contains

! --------------------------------------------------------------------------------------------------

subroutine set(self, state)

implicit none
type(fv3jedi_traj),  intent(inout) :: self
type(fv3jedi_state), intent(in)    :: state

integer :: isc,iec,jsc,jec,npz

! Pointers to rank 2 state
real(kind=kind_real), allocatable, dimension(:,:,:) :: u_tmp
real(kind=kind_real), allocatable, dimension(:,:,:) :: v_tmp

real(kind=kind_real), pointer, dimension(:,:,:) :: phis
real(kind=kind_real), pointer, dimension(:,:,:) :: frocean
real(kind=kind_real), pointer, dimension(:,:,:) :: frland
real(kind=kind_real), pointer, dimension(:,:,:) :: varflt
real(kind=kind_real), pointer, dimension(:,:,:) :: ustar
real(kind=kind_real), pointer, dimension(:,:,:) :: bstar
real(kind=kind_real), pointer, dimension(:,:,:) :: zpbl
real(kind=kind_real), pointer, dimension(:,:,:) :: cm
real(kind=kind_real), pointer, dimension(:,:,:) :: ct
real(kind=kind_real), pointer, dimension(:,:,:) :: cq
real(kind=kind_real), pointer, dimension(:,:,:) :: kcbl
real(kind=kind_real), pointer, dimension(:,:,:) :: tsm
real(kind=kind_real), pointer, dimension(:,:,:) :: khl
real(kind=kind_real), pointer, dimension(:,:,:) :: khu

isc = state%isc
iec = state%iec
jsc = state%jsc
jec = state%jec
npz = state%npz

! Allocate traj
allocate(self%u      (isc:iec, jsc:jec, npz))
allocate(self%v      (isc:iec, jsc:jec, npz))
allocate(self%ua     (isc:iec, jsc:jec, npz))
allocate(self%va     (isc:iec, jsc:jec, npz))
allocate(self%t      (isc:iec, jsc:jec, npz))
allocate(self%delp   (isc:iec, jsc:jec, npz))
allocate(self%qv     (isc:iec, jsc:jec, npz))
allocate(self%qi     (isc:iec, jsc:jec, npz))
allocate(self%ql     (isc:iec, jsc:jec, npz))
allocate(self%o3     (isc:iec, jsc:jec, npz))
allocate(self%w      (isc:iec, jsc:jec, npz))
allocate(self%delz   (isc:iec, jsc:jec, npz))
allocate(self%qls    (isc:iec, jsc:jec, npz))
allocate(self%qcn    (isc:iec, jsc:jec, npz))
allocate(self%cfcn   (isc:iec, jsc:jec, npz))
allocate(self%phis   (isc:iec, jsc:jec))
allocate(self%frocean(isc:iec, jsc:jec))
allocate(self%frland (isc:iec, jsc:jec))
allocate(self%varflt (isc:iec, jsc:jec))
allocate(self%ustar  (isc:iec, jsc:jec))
allocate(self%bstar  (isc:iec, jsc:jec))
allocate(self%zpbl   (isc:iec, jsc:jec))
allocate(self%cm     (isc:iec, jsc:jec))
allocate(self%ct     (isc:iec, jsc:jec))
allocate(self%cq     (isc:iec, jsc:jec))
allocate(self%kcbl   (isc:iec, jsc:jec))
allocate(self%ts     (isc:iec, jsc:jec))
allocate(self%khl    (isc:iec, jsc:jec))
allocate(self%khu    (isc:iec, jsc:jec))

!Initialize all to zero incase not in state
self%u       = 0.0_kind_real
self%v       = 0.0_kind_real
self%ua      = 0.0_kind_real
self%va      = 0.0_kind_real
self%t       = 0.0_kind_real
self%delp    = 0.0_kind_real
self%qv      = 0.0_kind_real
self%qi      = 0.0_kind_real
self%ql      = 0.0_kind_real
self%o3      = 0.0_kind_real
self%w       = 0.0_kind_real
self%delz    = 0.0_kind_real
self%qls     = 0.0_kind_real
self%qcn     = 0.0_kind_real
self%cfcn    = 0.0_kind_real
self%phis    = 0.0_kind_real
self%frocean = 0.0_kind_real
self%frland  = 0.0_kind_real
self%varflt  = 0.0_kind_real
self%ustar   = 0.0_kind_real
self%bstar   = 0.0_kind_real
self%zpbl    = 0.0_kind_real
self%cm      = 0.0_kind_real
self%ct      = 0.0_kind_real
self%cq      = 0.0_kind_real
self%kcbl    = 0.0_kind_real
self%ts      = 0.0_kind_real
self%khl     = 0.0_kind_real
self%khu     = 0.0_kind_real

! Copy mandatory parts of the trajecotry
allocate(u_tmp(isc:iec  , jsc:jec+1, npz))
allocate(v_tmp(isc:iec+1, jsc:jec  , npz))

call state%get_field('ud'  , u_tmp     )
call state%get_field('vd'  , v_tmp     )
call state%get_field('t'   , self%t    )
call state%get_field('delp', self%delp )
call state%get_field('sphum'   , self%qv   )

self%u = u_tmp(isc:iec, jsc:jec, :)
self%v = v_tmp(isc:iec, jsc:jec, :)

deallocate(u_tmp, v_tmp)

! Copy optional parts of the trajecotry (Rank 3)
if (state%has_field('ua'     )) call state%get_field('ua'     , self%ua  )
if (state%has_field('va'     )) call state%get_field('va'     , self%va  )
if (state%has_field('ice_wat')) call state%get_field('ice_wat', self%qi  )
if (state%has_field('liq_wat')) call state%get_field('liq_wat', self%ql  )
if (state%has_field('o3mr'   )) call state%get_field('o3mr'   , self%o3  )
if (state%has_field('w'      )) call state%get_field('w'      , self%w   )
if (state%has_field('delz'   )) call state%get_field('delz'   , self%delz)
if (state%has_field('qls'    )) call state%get_field('qls'    , self%qls )
if (state%has_field('qcn'    )) call state%get_field('qcn'    , self%qcn )
if (state%has_field('cfcn'   )) call state%get_field('cfcn'   , self%cfcn)

! Copy optional parts of the trajecotry (Rank 2)
if (state%has_field('phis')) then
  call state%get_field('phis', phis)
  self%phis = phis(:,:,1)
endif
if (state%has_field('frocean')) then
  call state%get_field('frocean', frocean)
  self%frocean = frocean(:,:,1)
endif
if (state%has_field('frland')) then
  call state%get_field('frland', frland)
  self%frland = frland(:,:,1)
endif
if (state%has_field('varflt')) then
  call state%get_field('varflt', varflt)
  self%varflt = varflt(:,:,1)
endif
if (state%has_field('ustar')) then
  call state%get_field('ustar', ustar)
  self%ustar = ustar(:,:,1)
endif
if (state%has_field('bstar')) then
  call state%get_field('bstar', bstar)
  self%bstar = bstar(:,:,1)
endif
if (state%has_field('zpbl')) then
  call state%get_field('zpbl', zpbl)
  self%zpbl = zpbl(:,:,1)
endif
if (state%has_field('cm')) then
  call state%get_field('cm', cm)
  self%cm = cm(:,:,1)
endif
if (state%has_field('ct')) then
  call state%get_field('ct', ct)
  self%ct = ct(:,:,1)
endif
if (state%has_field('cq')) then
  call state%get_field('cq', cq)
  self%cq = cq(:,:,1)
endif
if (state%has_field('kcbl')) then
  call state%get_field('kcbl', kcbl)
  self%kcbl = kcbl(:,:,1)
endif
if (state%has_field('tsm')) then
  call state%get_field('tsm', tsm)
  self%ts = tsm(:,:,1)
endif
if (state%has_field('khl')) then
  call state%get_field('khl', khl)
  self%khl = khl(:,:,1)
endif
if (state%has_field('khu')) then
  call state%get_field('khu', khu)
  self%khu = khu(:,:,1)
endif

end subroutine set

! --------------------------------------------------------------------------------------------------

subroutine wipe(self)

implicit none
type(fv3jedi_traj), pointer :: self

call deallocate_traj(self)

end subroutine wipe

! --------------------------------------------------------------------------------------------------

end module fv3jedi_traj_mod


text:
SUBROUTINE idomats(g)

!-------------------------------------------------------------
!
!  Integral discrete ordinates matrices
!
!   Control the construction of jmat, kmat, jpsi, kpsi
!   Matrices nned to be made only once, then idot.f90 will
!    use them to solve the system of equations
!
!-------------------------------------------------------------

USE invar
USE solvar
IMPLICIT NONE
INTEGER, INTENT(IN) :: g
INTEGER :: ieq, i, j, k, t, u, v, info, m, tt
REAL*8 :: symm, one, zero
REAL*8, DIMENSION(:), ALLOCATABLE :: q, wrk, symmvec
REAL*8, DIMENSION(:,:), ALLOCATABLE :: tmat

! Initialize the matrices
kmat = 0.0
jpsi = 0.0
kpsi = 0.0

! Set one and zero
one = 1.0
zero = 0.0

! Allocate the iteration Jacobian
ALLOCATE(jmat(neq,neq))
jmat = 0.0

! Allocate the temp source array
ALLOCATE(q(neq))

!-------------------------------------------------------------------
! Can now start the single sweep for constructing all matrices
CALL matsweep(g)

! Begin constructing the RHS
DO t = 0, lambda
   DO u = 0, lambda
      DO v = 0, lambda
         DO k = 1, nz
            DO j = 1, ny
               DO i = 1, nx
                  ieq = ((i-1) + (j-1)*nx + (k-1)*nx*ny)*ordcb + ordsq*t + order*u + v + 1
                  src(ieq) = s(i,j,k,t,u,v,g)/sigs((mat(i,j,k)),g,g)
               END DO
            END DO
         END DO
      END DO
   END DO
END DO

! Multiply the source vector by the Jacobi iteration matrix
IF (tpose == 1) THEN
   DO j = 1, neq
      q(j) = DOT_PRODUCT(jmat(:,j),src)
   END DO
ELSE
   q = MATMUL(jmat,src)
END IF

! Form the LHS matrix
jmat = -jmat
DO tt = 1, neq
   jmat(tt,tt) = 1.0 + jmat(tt,tt)
END DO

!------------------------------------------------------------------
IF (sym == 1 .OR. idos == 1 .OR. idos == 2) THEN
   ! Symmetrize jmat
   ALLOCATE(symmvec(neq))
   IF (tpose == 1) THEN
      DO k = 1, nz
         DO j = 1, ny
            DO i = 1, nx
               m = mat(i,j,k)
               symm = sigs(m,g,g)*dx(i)*dy(j)*dz(k)
               DO t = 0, lambda
                  DO u = 0, lambda
                     DO v = 0, lambda
                        ieq = ((k-1)*nx*ny + (j-1)*nx + (i-1))*ordcb + ordsq*t + order*u + v + 1
                        IF (order == 1) symmvec(ieq) = symm
                        IF (order /= 1) symmvec(ieq) = symm*(2.0*t+1.0)*(2.0*u+1.0)*(2.0*v+1.0)
                        DO tt = 1, neq
                           jmat(tt,ieq) = symmvec(ieq)*jmat(tt,ieq)
                        END DO
                     END DO
                  END DO
               END DO
            END DO
         END DO
      END DO
   ELSE
      DO k = 1, nz
         DO j = 1, ny
            DO i = 1, nx
               m = mat(i,j,k)
               symm = sigs(m,g,g)*dx(i)*dy(j)*dz(k)
               DO t = 0, lambda
                  DO u = 0, lambda
                     DO v = 0, lambda
                        ieq = ((k-1)*nx*ny + (j-1)*nx + (i-1))*ordcb + ordsq*t + order*u + v + 1
                        IF (order == 1) symmvec(ieq) = symm
                        IF (order /= 1) symmvec(ieq) = symm*(2.0*t+1.0)*(2.0*u+1.0)*(2.0*v+1.0)
                        DO tt = 1, neq
                           jmat(ieq,tt) = symmvec(ieq)*jmat(ieq,tt)
                        END DO
                     END DO
                  END DO
               END DO
            END DO
         END DO
      END DO
   END IF
END IF

!------------------------------------------------------------------------
! Other operations depend on the IDO Solver value: 0-Direct, 1-CG, 2-SOR
! If idos=0, need to factor jmat
IF (idos == 0) THEN

   ! Call LAPACK routines according to sym flag
   IF (sym == 0) THEN

      ! Allocate the piv vector as it's now needed
      ALLOCATE(piv(neq))

      ! Always factor the jmat matrix
      CALL dgetrf(neq,neq,jmat,neq,piv,info)
      IF (info /= 0) THEN
         WRITE (8,'(//,1X,A)') "ERROR: jmat either has illegal value or is singular, cannot factor."
         STOP
      END IF

      ! Now check the invf to see what to do next based on invf
      IF (invf == 0) THEN
         ! Set sv with q. Have needed values src, sv, jmat, kmat, jpsi, kpsi, piv
         sv = q
         DEALLOCATE(q)

      ELSE IF (invf == 1) THEN
         ALLOCATE(wrk(neq))
         CALL dgetri(neq,jmat,neq,piv,wrk,neq,info)
         IF (info /= 0) THEN
            WRITE (8,'(//,1X,A)') "ERROR: jmat either has illegal value or is singular, cannot invert."
            STOP
         END IF

         ! Save the matrix vector product of the inverted matrix and the modified source
         ! Then save the product of the inverted matrix and kmat
         IF (tpose == 1) THEN
            DO j = 1, neq
               sv(j) = DOT_PRODUCT(jmat(:,j),q)
            END DO
            ! kmat*jmat
            ALLOCATE(tmat(bcs,neq))
            DO t = 1, 8
               CALL DGEMM('n','n',bcs,neq,neq,one,kmat(:,:,t),bcs,jmat,neq,zero,tmat,bcs)
               kmat(:,:,t) = tmat
            END DO
         ELSE
            sv = MATMUL(jmat,q)
            ALLOCATE(tmat(neq,bcs))
            DO t = 1, 8
               CALL DGEMM('n','n',neq,bcs,neq,one,jmat,neq,kmat(:,:,t),neq,zero,tmat,neq)
               kmat(:,:,t) = tmat
            END DO
         END IF

         ! Now have the five needed values: src, sv, kmat, jpsi, kpsi
         DEALLOCATE(jmat,tmat)
         DEALLOCATE(piv,q,wrk)
      END IF

!---------------------------------------------------------------
! Symmetric direct
   ELSE

      ! Always factor the new LHS matrix with symmetric routine
      CALL dpotrf('U',neq,jmat,neq,info)
      IF (info /= 0) THEN
         WRITE (8,'(//,1X,A)') "ERROR: jmat either has illegal value or is singular, cannot factor."
         STOP
      END IF

      ! And always apply symmvec to q and to kmat
      q = symmvec*q
      IF (tpose == 1) THEN
         DO t = 1, 8
            DO j = 1, neq
               kmat(:,j,t) = symmvec(j)*kmat(:,j,t)
            END DO
         END DO
      ELSE
         DO t = 1, 8
            DO i = 1, neq
               kmat(i,:,t) = symmvec(i)*kmat(i,:,t)
            END DO
         END DO
      END IF

      ! Now check to see what to do next based on invf
      IF (invf == 0) THEN
         ! Set sv as q. Have needed values src, sv, jmat, kmat, jpsi, kpsi
         sv = q
         DEALLOCATE(q,symmvec)

      ELSE IF (invf == 1) THEN
         CALL dpotri('U',neq,jmat,neq,info)
         IF (info /= 0) THEN
            WRITE (8,'(//,1X,A)') "ERROR: jmat either has illegal value or is singular, cannot invert."
            STOP
         END IF

         ! Symmetric inverter only returns triangular portion of invert, copy rest
         DO j = 1, neq
            DO i = 1, j
               jmat(j,i) = jmat(i,j)
            END DO
         END DO

         ! Save sv as product of jmat, q. Reset kmat as jmat*kmat.
         IF (tpose == 1) THEN
            DO j = 1, neq
               sv(j) = DOT_PRODUCT(jmat(:,j),q)
            END DO
            ! Aided a bit by the symmetry of jmat, use that feature.
            ! tmat comes out as the normal kmat shape instead of transposed
            ALLOCATE(tmat(neq,bcs))
            DO t = 1, 8
               CALL DGEMM('n','t',neq,bcs,neq,one,jmat,neq,kmat(:,:,t),bcs,zero,tmat,neq)
               kmat(:,:,t) = TRANSPOSE(tmat)
            END DO
         ELSE
            sv = MATMUL(jmat,q)
            ALLOCATE(tmat(neq,bcs))
            DO t = 1, 8
               CALL DGEMM('n','n',neq,bcs,neq,one,jmat,neq,kmat(:,:,t),neq,zero,tmat,neq)
               kmat(:,:,t) = tmat
            END DO
         END IF

         ! Now have the five needed values: src, sv, kmat, jpsi, kpsi
         DEALLOCATE(jmat,tmat,q,symmvec)
      END IF
   END IF

!-----------------------------------------------------------------------
! CG or SOR
ELSE
   ! Reset kmat with the symmvec
   IF (tpose == 1) THEN
      DO t = 1, 8
         DO j = 1, neq
            kmat(:,j,t) = symmvec(j)*kmat(:,j,t)
         END DO
      END DO
   ELSE
      DO t = 1, 8
         DO i = 1, neq
            kmat(i,:,t) = symmvec(i)*kmat(i,:,t)
         END DO
      END DO
   END IF

   ! Set the sv vector
   sv = symmvec*q

   ! If a preconditioner is to be used, the inverse diagonal will be needed
   IF (idos == 1 .AND. pcf == 1) THEN
      ALLOCATE(dv(neq))
      DO i = 1, neq
         dv(i) = 1.0/jmat(i,i)
      END DO
      ! Check if the type is SSOR
      IF (pcty == 2) dv = sorw*dv
   END IF

   ! Now have the six needed values: src, sv, jmat, kmat, jpsi, kpsi (and dv)
   DEALLOCATE(q,symmvec)
END IF

RETURN
END SUBROUTINE idomats

max_stars_repo_name: lanl/pidots

max_issues_repo_name: lanl/pidots

max_forks_repo_name: lanl/pidots


text:
** Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
** See https://llvm.org/LICENSE.txt for license information.
** SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

*   Intrinsics and generics - AINT, DINT, ..., CABS.

	parameter(n = 23)
	integer rslts(n), expect(n)
	real  rrslts(n), rexpect(n)
	equivalence (rslts, rrslts), (expect, rexpect)

	double precision d3, d349, d36
	complex c43

	data i234 / 234 /
	data x349, xn39 / 3.49, -3.9 /
	data d3, d349, d36 / 3d0, 3.49d0, 3.6d0 /
	data c43  / (4.0, 3.0) /

c ---- tests 1 - 5:  AINT, DINT:

	data (rexpect(i), i = 1, 5) /2.0, -3.0, -2.0, 3.0, 2.0 /
	rrslts(1) = aint(2.1)
	rrslts(2) = aint(xn39)
	rrslts(3) = sngl( aint(-2.9D0) )
	rrslts(4) = aint(d36)
	rrslts(5) = sngl( dint(2.0d0) )

c ---- tests 6 - 13: ANINT, DNINT, NINT, IDNINT:

	data (rexpect(i), i = 6, 9) / 3.0, -4.0, 249653.0, 3.0 /
	rrslts(6) = anint( 2.5 )
	rrslts(7) = anint( xn39 )
	rrslts(8) = sngl( anint(249653.1d0) )
	rrslts(9) = sngl( dnint(d349) )

	data (expect(i), i = 10, 13) / 23, -4, 4, 4 /
	rslts(10) = nint(22.5)
	rslts(11) = nint(xn39)
	rslts(12) = nint(d36 + .41)
	rslts(13) = idnint(.39 + d36)

c ---- tests 14 - 23: IABS, ABS, DABS, CABS:

	data (expect(i), i = 14, 16) / 5, 234, 234 /
	rslts(14) = iabs(-5)
	rslts(15) = iabs(i234)
	rslts(16) = abs( - i234 )

	data (rexpect(i), i = 17, 21) /
     +        23.47, 3.49, 2.341d-3, 2.0, +2.3d0 /
	rrslts(17) = abs(-23.47)
	rrslts(18) = abs(x349)
	rrslts(19) = sngl( abs(2.341d-3) )
	rrslts(20) = sngl( abs(1.0 - d3) )
	rrslts(21) = dabs(-2.3d0)

	data expect(22), expect(23) / 130, 50 /
	rslts(22) = nint(10 * cabs((-12.0, 5.0)))
	rslts(23) = nint(10 * abs( c43 ))

c ---- check results:

	call check(rslts, expect, n)

	end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang