C DDSIRANK  SOURCE    AM        12/02/28    21:15:09     7291           
      SUBROUTINE ddsigRank(sigma,d2gt,i1,i2,i3,i4,i5,i6,lcp,
     .                     parahot3,idimpara3)  
                                              
c     This subroutine calculates the second derivative of the Rankine
c     plastic function with respect to sigma                                          
c     gt = ft = S1 = largest principal stress

      IMPLICIT REAL*8 (A-B,D-H,O-Z)
      implicit integer (I-K,M,N)
      implicit logical (L)
      implicit character*10 (C)        

****  dimension sigma(i6),dgt(i6),d2gt(i6,i6)
      dimension sigma(i6),dgt(6),d2gt(i6,i6)
****  dimension rcossigmapr(i3,i3),P1(3,3),pi1(3)
      dimension rcossigmapr(3,3),P1(3,3),pi1(3)
      dimension sigmacp(3),vloc(3),vloc2(1,3),vloc1(1)
      dimension vloc3(3,3),d2gtcp(3,3),parahot3(idimpara3)
      dimension sigmaim(6),sigmail(6),rhi(6),dgtm(6),dgtl(6)



*****
*     MESSAGES D'ERREUR ( SUPPRESSION DES AUTOMATIC OBJECTS)
      IF(I3.GT.3) PRINT *, ' DDSIGRANK - ERREUR  I2 = ', I2, ' > 2 '
      IF(I6.GT.6) PRINT *, ' DDSIGRANK - ERREUR  I6 = ', I6, ' > 6 '
*****


      if (lcp) then
c       In plane stress, the second derivative can be calculated analytically
        do jloc=1,6
          do iloc=1,6
            d2gt(iloc,jloc) = 0.d0
          enddo
        enddo
        P1(1,1) = 0.5d0
        P1(2,2) = 0.5d0
        P1(1,2) = -0.5d0
        P1(2,1) = -0.5d0
        do iloc=1,2
          P1(iloc,3) = 0.d0
          P1(3,iloc) = 0.d0
        enddo
        P1(3,3) = 2.d0
        pi1(1) = 1.d0
        pi1(2) = 1.d0
        pi1(3) = 0.d0
        sigmacp(1) = sigma(1)
        sigmacp(2) = sigma(2)
        sigmacp(3) = sigma(4)
        call mulAB(P1,sigmacp,vloc,3,3,3,1)
        call mulATB(sigmacp,vloc,vloc1,3,1,3,1)
        rloc = 0.5d0 * vloc1(1)
        psi1 = dsqrt(rloc)
        call mulATB(sigmacp,P1,vloc2,3,1,3,3)
        call mulAB(vloc,vloc2,vloc3,3,1,1,3)
               
        test = 0.5d0*(sigmacp(1)+sigmacp(2))
        psi1lim=max(1.E-7*test,1.d-3)
        if (dabs(psi1).le.psi1lim) then
          continue
c         d2gt = 0
        else
          do jloc=1,3
            do iloc=1,3
              d2gtcp(iloc,jloc)=(P1(iloc,jloc)/(2.0d0*psi1)) 
     .                 - (vloc3(iloc,jloc)/(4.0d0*psi1*psi1*psi1))
            enddo
          enddo
          do iloc=1,2
            d2gt(1,iloc) = d2gtcp(1,iloc)
            d2gt(2,iloc) = d2gtcp(2,iloc)
            d2gt(4,iloc) = d2gtcp(3,iloc)
            d2gt(iloc,4) = d2gtcp(iloc,3)
          end do
          d2gt(4,4) = d2gtcp(3,3)
        endif
      else
c       In full 3D, numerical differentiation
        do jloc=1,6
          do iloc=1,6
            d2gt(iloc,jloc) = 0.d0
          enddo
        enddo
        rh = 1.0d-5
        do iloc=1,6
c         ith component of sigma => ith column of d2gt
c         calculation of   sigma + hi * ei   and   sigma - hi * ei
          do jloc=1,6
            rhi(jloc) = 0.d0
          enddo
          rhi(iloc) = 1.0d0 * rh
          do jloc=1,6
            sigmaim(jloc) = sigma(jloc) + rhi(jloc)
            sigmail(jloc) = sigma(jloc) - rhi(jloc)
          enddo
c         Calculation of q(x+hi*ei) = dft(x+hi*ei)/dsig
          call dsigRank(sigmaim,dgtm,3,6,parahot3,idimpara3,
     .                  rkappait,lcp)
          call dsigRank(sigmail,dgtl,3,6,parahot3,idimpara3,
     .                  rkappait,lcp)
          do jloc=1,6
            d2gt(jloc,iloc) = (dgtm(jloc) - dgtl(jloc))/(2.0d0*rh)
          enddo
        enddo
c       Calculation of f(x+hi*ei)        
      endif
     
      RETURN
      END