{"id":253,"date":"2026-06-29T21:09:04","date_gmt":"2026-06-29T13:09:04","guid":{"rendered":"https:\/\/numsimlab.com\/?p=253"},"modified":"2026-06-29T21:09:43","modified_gmt":"2026-06-29T13:09:43","slug":"%e6%99%b6%e4%bd%93%e5%a1%91%e6%80%a7%e5%8a%9b%e5%ad%a6-huang-vumat","status":"publish","type":"post","link":"https:\/\/numsimlab.com\/?p=253","title":{"rendered":"\u6676\u4f53\u5851\u6027\u529b\u5b66-huang-vumat-\u6e90\u7801"},"content":{"rendered":"\n
      SUBROUTINE VUMAT(NBLOCK, NDIR, NSHR, NSTATEV, NFIELDV,\n     1  NPROPS,LANNEAL,STEPTIME, TOTALTIME, DT, CMNAME, COORDMP,\n     2  CHARLENGTH,PROPS, DENSITY, STRAININC, RELSPININC,\n     3  TEMPOLD, STRETCHOLD, DEFGRADOLD, FIELDOLD,\n     4  STRESSOLD, STATEOLD, ENERINTERNOLD, ENERINELASOLD,\n     5  TEMPNEW, STRETCHNEW, DEFGRADNEW, FIELDNEW,\n     6  STRESSNEW, STATENEW, ENERINTERNNEW,ENERINELASNEW)\nC     -------------------VUMAT Interface Variable Description-----------------------\nC     NBLOCK\u2014\u2014Number of material integration points processed in this VUMAT call\nC     NDIR\u2014\u2014Number of diagonal tensor components for stress\/strain; 3 for plane\/3D problems\nC     NSHR\u2014\u2014Number of off-diagonal shear tensor components; 1 for plane, 3 for 3D\nC     NSTATEV\u2014\u2014Number of user-defined state dependent variables (SDVs)\nC     NFIELDV\u2014\u2014Number of external user field variables\nC     NPROPS\u2014\u2014Number of material constitutive parameters\nC     LANNEAL\u2014\u2014Annealing flag to reinitialize internal state variables\nC     STEPTIME\u2014\u2014Time value at the start of current increment\nC     TOTALTIME\u2014\u2014Total accumulated analysis time\nC     DT\u2014\u2014Time increment size of current step\nC     CMNAME\u2014\u2014Material name character string for distinguishing multiple materials\nC     COORDMP\u2014\u2014Spatial coordinates of material integration points\nC     CHARLENGTH\u2014\u2014Characteristic element length at material point\nC     PROPS\u2014\u2014Array storing user-defined material constants (same as UMAT material table)\nC     DENSITY\u2014\u2014Material mass density defined in material input block\nC     STRAININC\u2014\u2014Incremental strain tensor at each material point\nC     RELSPININC\u2014\u2014Incremental relative rotation tensor under reference rotating coordinate system\nC     TEMPOLD\u2014\u2014Temperature at the beginning of current increment\nC     STRETCHOLD\u2014\u2014Left stretch tensor U at increment start; suffix OLD = start of increment, NEW = end of increment\nC     DEFGRADOLD\u2014\u2014Deformation gradient tensor at increment start\nC     FIELDOLD\u2014\u2014User external field variables at increment start\nC     STRESSOLD\u2014\u2014Cauchy stress tensor at increment start\nC     STATEOLD\u2014\u2014State dependent variables at increment start\nC     ENERINTERNOLD\u2014\u2014Internal energy density at increment start\nC     ENERINELASOLD\u2014\u2014Elastic strain energy density at increment start\nC     TEMPNEW\u2014\u2014Temperature at the end of current increment\nC     STRETCHNEW\u2014\u2014Left stretch tensor U at increment end\nC     DEFGRADNEW\u2014\u2014Deformation gradient tensor at increment end\nC     FIELDNEW\u2014\u2014User external field variables at increment end\nC     STRESSNEW\u2014\u2014Updated Cauchy stress tensor at increment end (output variable)\nC     STATENEW\u2014\u2014Updated state dependent variables at increment end (output variable)\nC     ENERINTERNNEW\u2014\u2014Internal energy density at increment end\nC     ENERINELASNEW\u2014\u2014Elastic strain energy density at increment end\nC     -------------------------------------------------------------------------------\n      INCLUDE 'VABA_PARAM.INC'\n\nC VARIABLE DECLARATION\n      DIMENSION PROPS(NPROPS), DENSITY(NBLOCK), COORDMP(NBLOCK,*),\n     1  CHARLENGTH(NBLOCK), STRAININC(NBLOCK, NDIR+NSHR),\n     2  RELSPININC(NBLOCK,NSHR), TEMPOLD(NBLOCK),\n     3  STRETCHOLD(NBLOCK,NDIR+NSHR),\n     4  DEFGRADOLD(NBLOCK, NDIR+NSHR+NSHR),\n     5  FIELDOLD(NBLOCK,NFIELDV), STRESSOLD(NBLOCK,NDIR+NSHR),\n     6  STATEOLD(NBLOCK,NSTATEV), ENERINTERNOLD(NBLOCK),\n     7  ENERINELASOLD(NBLOCK), TEMPNEW(NBLOCK),\n     8  STRETCHNEW(NBLOCK,NDIR+NSHR),\n     9  DEFGRADNEW(NBLOCK,NDIR+NSHR+NSHR),\n     1  FIELDNEW(NBLOCK,NFIELDV),\n     2  STRESSNEW(NBLOCK,NDIR+NSHR), STATENEW(NBLOCK,NSTATEV),\n     3  ENERINTERNNEW(NBLOCK), ENERINELASNEW(NBLOCK)\n\n      CHARACTER*80 CMNAME\n\n!---------------------------------------------------------------------\n! LOCAL VARIABLE DECLARATION\n!---------------------------------------------------------------------\n\n      INTEGER ZERO, ONE, TWO, NTENS,\n     &        NDI, NSTATV,\n     &        I, NSHRUMAT, NPROPSUMAT\n\n      DOUBLE PRECISION     DTIME\n\n      DOUBLE PRECISION STRESS(NDIR+NSHR), STATEV(NSTATEV),\n     &        STRAN(NDIR+NSHR), DSTRAN(NDIR+NSHR),\n     &        TIME(2),\n     &        DFGRD0(3,3), DFGRD1(3,3)\n\n      DOUBLE PRECISION PROPSUMAT(NPROPS)\n\n      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0)\n\n!*************************************************************************\nC Initialize specific SDVs at zero total time (initial analysis step)\n      IF (TOTALTIME .EQ. ZERO) THEN\n        DO KM = 1, NBLOCK\n           STATENEW(KM,2) = ZERO\n           STATENEW(KM,7) = STATEOLD(KM,7)\n        ENDDO\n      ENDIF\n!*************************************************************************\n\nC Loop over all material integration points in current block\n      DO 100 KM = 1,NBLOCK\n\nC Map old stress and strain increment to local UMAT-style array\n        DO I = 1, NDIR\n          STRESS(I) = STRESSOLD(KM,I)\n          DSTRAN(I) = STRAININC(KM,I)\n        ENDDO\n        STRESS(4) = STRESSOLD(KM,4)\n        DSTRAN(4) = TWO * STRAININC(KM,4)\n        IF (NSHR .GT. 1) THEN\n          STRESS(5) = STRESSOLD(KM,6)\n          DSTRAN(5) = TWO * STRAININC(KM,6)\n          STRESS(6)  = STRESSOLD(KM,5)\n          DSTRAN(6) = TWO * STRAININC(KM,5)\n        ENDIF\n\nC Assign deformation gradient at increment end DFGRD1 (3x3 full tensor)\n      DFGRD1(1,1) = DEFGRADNEW(KM,1)\n      DFGRD1(2,2) = DEFGRADNEW(KM,2)\n      DFGRD1(3,3) = DEFGRADNEW(KM,3)\n      DFGRD1(1,2) = DEFGRADNEW(KM,4)\n      DFGRD1(2,3) = DEFGRADNEW(KM,5)\n      DFGRD1(3,1) = DEFGRADNEW(KM,6)\n      DFGRD1(2,1) = DEFGRADNEW(KM,7)\n      DFGRD1(3,2) = DEFGRADNEW(KM,8)\n      DFGRD1(1,3) = DEFGRADNEW(KM,9)\n\nC Assign deformation gradient at increment start DFGRD0 (3x3 full tensor)\n      DFGRD0(1,1) = DEFGRADOLD(KM,1)\n      DFGRD0(2,2) = DEFGRADOLD(KM,2)\n      DFGRD0(3,3) = DEFGRADOLD(KM,3)\n      DFGRD0(1,2) = DEFGRADOLD(KM,4)\n      DFGRD0(2,3) = DEFGRADOLD(KM,5)\n      DFGRD0(3,1) = DEFGRADOLD(KM,6)\n      DFGRD0(2,1) = DEFGRADOLD(KM,7)\n      DFGRD0(3,2) = DEFGRADOLD(KM,8)\n      DFGRD0(1,3) = DEFGRADOLD(KM,9)\n\nC Pass global VUMAT variables to local UMAT-style variables for crystal plasticity subroutine\n        NSTATV = NSTATEV\n        STATEV(:) = STATEOLD(KM,:)\n        DTIME  = DT\n        NDI    = NDIR\n        NSHRUMAT   = NSHR\n        NTENS  = NDIR + NSHR\n        PROPSUMAT  = PROPS\n        NPROPSUMAT = NPROPS\n        TIME(1)=TOTALTIME\n        TIME(2)=STEPTIME\n\nC Call Huang Yonggang single crystal plasticity constitutive subroutine to update stress and state variables\n      CALL CRYSTALPLASTICITY(STRESS,STATEV,STRAN,DSTRAN,\n     &  TIME,DTIME,CMNAME,NDI,NSHRUMAT,NTENS,NSTATV,PROPSUMAT,\n     &  NPROPSUMAT,DFGRD0,DFGRD1)\n\nC Write updated stress back to VUMAT output array STRESSNEW\n        DO I = 1, NDIR\n          STRESSNEW(KM,I) = STRESS(I)\n        ENDDO\n\n      STRESSNEW(KM,4) = STRESS(4)\n\n      IF( NSHRUMAT .GT. 1 ) THEN\n        STRESSNEW(KM,5)  = STRESS(6)\n        STRESSNEW(KM,6)  = STRESS(5)\n      ENDIF\n\nC Write updated state dependent variables back to VUMAT output array STATENEW\n      STATENEW(KM,:) = STATEV(:)\n\n  100 CONTINUE ! End of material point loop\n\n      RETURN\n      END\n\nC     Huang Yonggang Single Crystal Plasticity Constitutive Subroutine\nC     Trimmed version retaining only stress and state variable update logic\n      SUBROUTINE CRYSTALPLASTICITY(STRESS,STATEV,STRAN,DSTRAN,\n     2 TIME,DTIME,CMNAME,NDI,NSHR,NTENS,NSTATV,PROPS,NPROPS,\n     3 DFGRD0,DFGRD1)\n\nC-----  Single precision compilation note for Cray machines:\nC     (1) Delete statement \"IMPLICIT*8 (A-H,O-Z)\";\nC     (2) Change \"REAL*8 FUNCTION\" to \"FUNCTION\";\nC     (3) Replace double precision intrinsic DSIGN with SIGN.\nC\nC-----  Internal Auxiliary Subroutines:\nC\nC       ROTATION     -- Construct crystal orientation rotation matrix;\nC                       computes direction cosines of cubic crystal [100], [010], [001]\nC                       axes under global coordinate system at initial state\nC\nC       SLIPSYS      -- Generate independent slip systems, unit slip direction vectors\nC                       and unit slip plane normal vectors for cubic crystals at initial state\nC\nC       GSLPINIT     -- Assign initial slip system critical shear strength values\nC\nC       STRAINRATE   -- Evaluate slip shear strain rate based on resolved shear stress\nC                       and current slip system hardening strength following power-law viscoplasticity\nC\nC       LATENTHARDEN -- Assemble self-hardening and latent hardening interaction matrix\nC\nC       ITERATION    -- Construct Jacobian arrays for Newton-Rhapson implicit iteration\nC\nC       LUDCMP       -- Perform LU matrix decomposition for linear system solving\nC\nC       LUBKSB       -- Solve linear system via precomputed LU decomposition\nC\nC\nC-----  Internal Auxiliary Function:\nC\nC       F -- Slip system shear strain rate function (power-law viscoplasticity)\nC\nC-----  Subroutine Input\/Output Variables:\nC\nC       STRESS -- Cauchy stress tensor (INPUT & OUTPUT)\nC                 Finite deformation framework adopts true Cauchy stress\nC       STATEV -- State dependent solution variables (INPUT & OUTPUT)\nC\nC-----  Passed auxiliary variables for reference:\nC\nC       STRAN  -- Integral logarithmic strain tensor for finite deformation\nC                 Equivalent to time integral of symmetric velocity gradient\nC       DSTRAN -- Incremental strain tensor\nC       CMNAME -- Material name string defined in *MATERIAL keyword block\nC       NDI    -- Count of direct normal stress tensor components\nC       NSHR   -- Count of engineering shear stress tensor components\nC       NTENS  -- Total tensor component count = NDI + NSHR\nC       NSTATV -- Total number of state dependent variables defined via *DEPVAR\nC       PROPS  -- User material constants defined under *USER MATERIAL keyword\nC       NPROPS -- Total number of user material constants\nC\nC-----  Constitutive Theory Overview:\nC     This subroutine implements finite deformation single crystal plasticity for ABAQUS.\nC     Crystal plastic slip follows Schmid's resolved shear stress criterion.\nC     Total strain increment decomposes additively into elastic lattice stretch strain and plastic slip strain.\nC     Elastic strain increment corresponds to pure lattice stretching; plastic strain is the superposition\nC     of shear slip over all activated slip systems.\nC     Slip shear strain increment is a power-law function of resolved shear stress normalized by slip system strength.\nC     Slip system hardening strength increment couples with accumulated shear strain via self\/latent hardening interaction.\nC\nC-----  Time Integration Scheme:\nC     Implicit backward integration algorithm proposed by Peirce, Shih & Needleman (1984) is adopted.\nC     Optional nested Newton-Rhapson iteration is available to converge stress and internal state variables per increment.\nC\nC-----  Crystal System Restriction:\nC     Core implementation for single cubic crystals (FCC\/BCC). Extensions to HCP, tetragonal, orthotropic\nC     lattices only require modifications to ROTATION and SLIPSYS subroutines to incorporate lattice aspect ratios.\nC\nC-----  Critical User Setup Requirements:\nC\nC     (1) Minimum required state variable count NSTATV:\nC         NSTATV >= 10 * NSLPTL + 5\nC         NSLPTL = total independent slip systems across all slip families\nC         Slip systems (s,-m), (-s,m), (-s,-m) are treated as dependent duplicates of (s,m) and excluded\nC         Cubic slip family examples: {110}<111> contains 12 independent slip systems\nC         If additional constitutive state parameters are required (e.g. Zarka model), extend to:\nC         NSTATV >= NPARMT + 10 * NSLPTL + 5\nC\nC     (2) Tangent stiffness matrix asymmetry:\nC         Latent hardening introduces non-symmetric consistent tangent stiffness.\nC         Must declare keyword \"UNSYMM\" under *USER MATERIAL in ABAQUS input deck.\nC\n      PARAMETER (ND=150)\nC-----  ND defines maximum array dimension for slip system storage\nC     Default value 150 supports up to three cubic slip families fully activated.\nC     Reduce ND to NSLPTL if fewer slip families are used (e.g. ND=12 for single {110}<111> family).\nC\n      include 'aba_param.inc'\nC\n      CHARACTER*8 CMNAME\n      EXTERNAL F\n\n      DIMENSION STRESS(NTENS),STATEV(NSTATV),\n     2 STRAN(NTENS),DSTRAN(NTENS),TIME(2),\n     3 PROPS(NPROPS),DROT(3,3),DFGRD0(3,3),DFGRD1(3,3)\n\n      DIMENSION ISPDIR(3), ISPNOR(3), NSLIP(3),\n     2          SLPDIR(3,ND), SLPNOR(3,ND), SLPDEF(6,ND),\n     3          SLPSPN(3,ND), DSPDIR(3,ND), DSPNOR(3,ND),\n     4          DLOCAL(6,6), D(6,6), ROTD(6,6), ROTATE(3,3),\n     5          FSLIP(ND), DFDXSP(ND), DDEMSD(6,ND),\n     6          H(ND,ND), DDGDDE(ND,6),\n     7          DSTRES(6), DELATS(6), DSPIN(3), DVGRAD(3,3),\n     8          DGAMMA(ND), DTAUSP(ND), DGSLIP(ND),\n     9          WORKST(ND,ND), INDX(ND), TERM(3,3), TRM0(3,3), ITRM(3)\n\n      DIMENSION FSLIP1(ND), STRES1(6), GAMMA1(ND), TAUSP1(ND),\n     2          GSLP1(ND), SPNOR1(3,ND), SPDIR1(3,ND), DDSDE1(6,6),\n     3          DSOLD(6), DGAMOD(ND), DTAUOD(ND), DGSPOD(ND),\n     4          DSPNRO(3,ND), DSPDRO(3,ND),\n     5          DHDGDG(ND,ND)\n      DOUBLE PRECISION LV(3,3),\n     4  LVT(3,3),WV(3,3),ROTA(3,3),ROTAT(3,3),DETADG(3,3),\n     5  DGINV(3,3),TERMIDW(3,3),TERMIPW(3,3),TERMIDWINV(3,3)\n\nC-----  NSLIP  -- Independent slip system count per slip family\nC-----  SLPDIR -- Unit slip direction vectors in initial crystal local coordinate system\nC-----  SLPNOR -- Unit slip plane normal vectors in initial crystal local coordinate system\nC-----  SLPDEF -- Slip deformation tensor (Schmid factor matrix, Voigt 6-component format)\nC                 SLPDEF(1,i) = SLPDIR(1,i)*SLPNOR(1,i)\nC                 SLPDEF(2,i) = SLPDIR(2,i)*SLPNOR(2,i)\nC                 SLPDEF(3,i) = SLPDIR(3,i)*SLPNOR(3,i)\nC                 SLPDEF(4,i) = SLPDIR(1,i)*SLPNOR(2,i)+SLPDIR(2,i)*SLPNOR(1,i)\nC                 SLPDEF(5,i) = SLPDIR(1,i)*SLPNOR(3,i)+SLPDIR(3,i)*SLPNOR(1,i)\nC                 SLPDEF(6,i) = SLPDIR(2,i)*SLPNOR(3,i)+SLPDIR(3,i)*SLPNOR(2,i)\nC                 Index i denotes ith independent slip system\nC-----  SLPSPN -- Slip spin tensor components (only required for finite rotation formulation)\nC                 SLPSPN(1,i) = 0.5*(SLPDIR(1,i)*SLPNOR(2,i)-SLPDIR(2,i)*SLPNOR(1,i))\nC                 SLPSPN(2,i) = 0.5*(SLPDIR(3,i)*SLPNOR(1,i)-SLPDIR(1,i)*SLPNOR(3,i))\nC                 SLPSPN(3,i) = 0.5*(SLPDIR(2,i)*SLPNOR(3,i)-SLPDIR(3,i)*SLPNOR(2,i))\nC-----  DSPDIR -- Incremental update of slip direction vectors under finite rotation\nC-----  DSPNOR -- Incremental update of slip plane normal vectors under finite rotation\nC\nC-----  DLOCAL -- Anisotropic elastic stiffness matrix defined in crystal local coordinate system\nC-----  D      -- Anisotropic elastic stiffness matrix rotated to global Cartesian coordinate system\nC-----  ROTD   -- Voigt rotation transformation matrix mapping DLOCAL to global D\nC\nC-----  ROTATE -- Crystal orientation rotation matrix; stores direction cosines of crystal [100], [010], [001] axes\nC                 relative to global X\/Y\/Z axes at initial state\nC\nC-----  FSLIP  -- Current shear strain rate magnitude for each slip system\nC-----  DFDXSP -- Derivative dF\/dX where X = resolved shear stress \/ slip system hardening strength\nC\nC-----  DDEMSD -- Double dot product of elastic stiffness tensor with Schmid tensor, plus\nC                 spin-stress coupling term exclusively for finite rotation kinematics\nC\nC-----  H      -- Slip system hardening interaction matrix\nC                 H(i,i) = Self-hardening modulus of ith slip system\nC                 H(i,j) = Latent hardening modulus on system i induced by slip on system j (i\u2260j)\nC\nC-----  DDGDDE -- Derivative of slip shear strain increment with respect to macroscopic strain increment\nC\nC-----  DSTRES -- Jaumann corotational stress increment tensor co-rotated with material spin\nC-----  DELATS -- Lattice elastic stretch strain increment (macro strain minus plastic slip strain)\nC                 DELATS(1-3) = Normal elastic strain increments\nC                 DELATS(4-6) = Engineering elastic shear strain increments\nC-----  DSPIN  -- Material element spin increment tensor components\nC                 DSPIN(1) = 12 component of spin tensor\nC                 DSPIN(2) = 31 component of spin tensor\nC                 DSPIN(3) = 23 component of spin tensor\nC\nC-----  DVGRAD -- Incremental velocity gradient tensor (velocity gradient \u00d7 time increment)\nC\nC-----  DGAMMA -- Incremental plastic shear strain magnitude for each slip system\nC-----  DTAUSP -- Incremental resolved shear stress on each slip system\nC-----  DGSLIP -- Incremental hardening strength increase for each slip system\nC\nC-----  Iteration temporary storage arrays:\nC            FSLIP1, STRES1, GAMMA1, TAUSP1, GSLP1 , SPNOR1, SPDIR1,\nC            DDSDE1, DSOLD , DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO,\nC            DHDGDG\nC\nC-----  STATEV State Variable Storage Layout Definition:\nC            NSLPTL = total independent slip systems across all slip families\nC\nC       STATEV(1          : NSLPTL)    : Current slip system hardening strength g_\u03b1\nC       STATEV(NSLPTL+1   : 2*NSLPTL)  : Accumulated shear strain \u03b3_\u03b1 on each slip system\nC       STATEV(2*NSLPTL+1 : 3*NSLPTL)  : Current resolved shear stress \u03c4_\u03b1 on each slip system\nC\nC       STATEV(3*NSLPTL+1 : 6*NSLPTL)  : Current updated slip plane normal vectors m_\u03b1 (3 components per system)\nC       STATEV(6*NSLPTL+1 : 9*NSLPTL)  : Current updated slip direction vectors s_\u03b1 (3 components per system)\nC\nC       STATEV(9*NSLPTL+1 : 10*NSLPTL) : Cumulative absolute shear strain |\u03b3_\u03b1| for each individual slip system\nC\nC       STATEV(10*NSLPTL+1)             : Global total cumulative absolute shear strain sum(|\u03b3_\u03b1|) over all slip systems\nC\nC       STATEV(10*NSLPTL+2 : NSTATV-4)  : User-extended auxiliary constitutive state parameters (optional)\nC\nC       STATEV(NSTATV-3)               : Slip system count of first slip family\nC       STATEV(NSTATV-2)               : Slip system count of second slip family\nC       STATEV(NSTATV-1)               : Slip system count of third slip family\nC       STATEV(NSTATV)                 : Total independent slip system count NSLPTL\nC\nC-----  PROPS Material Constant Array Layout Definition:\nC\nC       PROPS(1) - PROPS(21) -- Anisotropic elastic stiffness constants\nC\nC            Isotropic elastic model: PROPS(i)=0 for i>2\nC                          PROPS(1) = Young's Modulus E\nC                          PROPS(2) = Poisson's Ratio \u03bd\nC\nC            Cubic elastic model: PROPS(i)=0 for i>3\nC                          PROPS(1) = C11\nC                          PROPS(2) = C12\nC                          PROPS(3) = C44\nC\nC            Orthotropic elastic model: PROPS(1)-PROPS(9) match ABAQUS orthotropic elastic input order\nC                          D1111, D1122, D2222, D1133, D2233, D3333, D1212, D1313, D2323\nC\nC            General fully anisotropic elastic model: PROPS(1)-PROPS(21) full 4th-order stiffness tensor Voigt components\nC\nC\nC       PROPS(25) - PROPS(56) -- Slip family definition parameters for cubic crystal\nC\nC            PROPS(25) -- Number of distinct slip families (maximum 3, input as floating point value e.g. 3.0)\nC\nC            PROPS(33) - PROPS(35) -- Miller indices of reference slip plane normal for slip family 1 e.g. (1,1,0)\nC            PROPS(36) - PROPS(38) -- Miller indices of reference slip direction for slip family 1 e.g. [1,-1,1]\nC\nC            PROPS(41) - PROPS(43) -- Reference slip plane normal Miller indices for slip family 2\nC            PROPS(44) - PROPS(46) -- Reference slip direction Miller indices for slip family 2\nC\nC            PROPS(49) - PROPS(51) -- Reference slip plane normal Miller indices for slip family 3\nC            PROPS(52) - PROPS(54) -- Reference slip direction Miller indices for slip family 3\nC\nC\nC       PROPS(57) - PROPS(72) -- Crystal orientation definition parameters\nC            Two non-parallel vectors required to construct orientation rotation matrix\nC\nC            PROPS(57) - PROPS(59) -- Local crystal coordinate vector 1 Miller indices e.g. [1,1,0]\nC            PROPS(60) - PROPS(62) -- Global Cartesian coordinate vector 1 components (non-unit vector allowed)\nC\nC            PROPS(65) - PROPS(67) -- Local crystal coordinate vector 2 Miller indices\nC            PROPS(68) - PROPS(70) -- Global Cartesian coordinate vector 2 components\nC\nC\nC       PROPS(73) - PROPS(96) -- Viscoplastic power-law slip rate parameters for each slip family\nC\nC            PROPS(73) - PROPS(80) -- Power-law parameters for slip family 1\nC            PROPS(81) - PROPS(88) -- Power-law parameters for slip family 2\nC            PROPS(89) - PROPS(96) -- Power-law parameters for slip family 3\nC\nC\nC       PROPS(97) - PROPS(144) -- Self & latent hardening law parameters for each slip family\nC\nC            PROPS(97) - PROPS(104)-- Self-hardening parameters slip family 1\nC            PROPS(105)- PROPS(112)-- Latent hardening interaction parameters slip family 1\nC\nC            PROPS(113)- PROPS(120)-- Self-hardening parameters slip family 2\nC            PROPS(121)- PROPS(128)-- Latent hardening interaction parameters slip family 2\nC\nC            PROPS(129)- PROPS(136)-- Self-hardening parameters slip family 3\nC            PROPS(137)- PROPS(144)-- Latent hardening interaction parameters slip family 3\nC\nC\nC       PROPS(145)- PROPS(152)-- Time integration and finite deformation switch parameters\nC\nC            PROPS(145) -- Implicit integration weighting factor \u03b8 (0 \u2264 \u03b8 \u2264 1)\nC                          \u03b8=0: Explicit forward Euler integration\nC                          \u03b8=0.5: Recommended midpoint integration\nC                          \u03b8=1.0: Fully implicit backward Euler integration\nC\nC            PROPS(146) -- Finite geometry activation flag NLGEOM\nC                          0.0 = Small deformation infinitesimal theory\nC                          Non-zero = Finite rotation & finite strain kinematics (requires *NLGEOM step keyword)\nC\nC\nC       PROPS(153)- PROPS(160)-- Newton iteration control parameters\nC\nC            PROPS(153) -- Iteration enable flag ITRATN\nC                          0.0 = No nested iteration, single step solve\nC                          Non-zero = Activate Newton-Rhapson iteration\nC\nC            PROPS(154) -- Maximum allowed iteration count ITRMAX\nC\nC            PROPS(155) -- Convergence tolerance GAMERR for slip shear strain residual\nC\nC     ***************************************\nC        Supplementary routine to compute incremental rotation tensor DROT\nC     ***************************************\nC     Compute deformation gradient increment DFGRAD1 - DFGRD0\n      DO I = 1,3\n          DO J = 1,3\n              DETADG(I,J) = DFGRD1(I,J)-DFGRD0(I,J)\n          END DO\n      END DO\n\nC     Copy end-of-step deformation gradient for matrix inversion\n      DO I  = 1,3\n          DO J = 1,3\n              DGINV(I,J) = DFGRD1(I,J)\n          END DO\n      END DO\n\nC     Call subroutine to calculate inverse of deformation gradient\n      CALL GET_INV_DET(DGINV,GARB,0)\n\nC     Compute incremental velocity gradient LV = \u0394F \u00b7 F\u207b\u00b9\n      DO I = 1,3\n          DO J = 1,3\n              LV(I,J) = 0.0D0\n              DO K = 1,3\n                  LV(I,J) = LV(I,J)+DETADG(I,K)*DGINV(K,J)\n              END DO\n          END DO\n      END DO\n\nC     Extract skew-symmetric spin tensor component WV = \u03a9\u00b7dt\n      DO I =1,3\n          WV(I,I) = 0.0D0\n      END DO\n\n      WV(1,2) =0.5D0*( LV(1,2)-LV(2,1))\n      WV(2,1) = -WV(1,2)\n      WV(1,3) =0.5D0*( LV(1,3)-LV(3,1))\n      WV(3,1) = -WV(1,3)\n      WV(2,3) =0.5D0*( LV(2,3)-LV(3,2))\n      WV(3,2) = -WV(2,3)\n\nC     Zero temporary matrix storage\n      CALL CLEAR(TERMIDW,9)\n      CALL CLEAR(TERMIPW,9)\n\nC     Construct I \u00b1 0.5\u03a9dt matrix for midpoint rotation integration\n      DO I = 1,3\n          TERMIDW(I,I) = 1.0D0\n          TERMIPW(I,I) = 1.0D0\n          DO J = 1,3\n              TERMIDW(I,J) = TERMIDW(I,J)-0.5D0*WV(I,J)\n              TERMIPW(I,J) = TERMIPW(I,J)+0.5D0*WV(I,J)\n          END DO\n      END DO\n\nC     Copy matrix for inversion\n      DO I=1,3\n          DO J = 1,3\n              TERMIDWINV(I,J) = TERMIDW(I,J)\n          END DO\n      END DO\n\nC     Compute inverse of (I + 0.5\u03a9dt)\n      CALL GET_INV_DET(TERMIDWINV,GARB,0)\n\nC     Calculate incremental rotation tensor DROT = (I+0.5\u03a9dt)\u207b\u00b9 \u00b7 (I-0.5\u03a9dt)\n      DO I = 1,3\n          DO J = 1,3\n              ROTA(I,J) = 0.0D0\n              DO K =1,3\n                  ROTA(I,J) = ROTA(I,J)+TERMIDWINV(I,K)*TERMIPW(K,J)\n              END DO\n          END DO\n      END DO\n\nC     Assign final incremental rotation matrix to DROT\n      DO I = 1,3\n          DO J = 1,3\n              DROT(I,J) = ROTA(I,J)\n          END DO\n      END DO\n\nC-----  Construct anisotropic elastic stiffness matrix DLOCAL in crystal local coordinate system\n      DO J=1,6\n         DO I=1,6\n            DLOCAL(I,J)=0.\n         END DO\n      END DO\n\nC     Check if fully anisotropic elastic constants are provided\n      CHECK=0.\n      DO J=10,21\n         CHECK=CHECK+ABS(PROPS(J))\n      END DO\n\n      IF (CHECK.EQ.0.) THEN\nC     Check orthotropic elastic input\n         DO J=4,9\n            CHECK=CHECK+ABS(PROPS(J))\n         END DO\n\n         IF (CHECK.EQ.0.) THEN\nC     Isotropic or cubic elastic model\n            IF (PROPS(3).EQ.0.) THEN\nC-----  Isotropic elastic constitutive model\n               GSHEAR=PROPS(1)\/2.\/(1.+PROPS(2))\n               E11=2.*GSHEAR*(1.-PROPS(2))\/(1.-2.*PROPS(2))\n               E12=2.*GSHEAR*PROPS(2)\/(1.-2.*PROPS(2))\n\n               DO J=1,3\n                  DLOCAL(J,J)=E11\n                  DO I=1,3\n                     IF (I.NE.J) DLOCAL(I,J)=E12\n                  END DO\n                  DLOCAL(J+3,J+3)=GSHEAR\n               END DO\n\n            ELSE\nC-----  Cubic elastic stiffness matrix C11, C12, C44\n               DO J=1,3\n                  DLOCAL(J,J)=PROPS(1)\n                  DO I=1,3\n                     IF (I.NE.J) DLOCAL(I,J)=PROPS(2)\n                  END DO\n                  DLOCAL(J+3,J+3)=PROPS(3)\n               END DO\n            END IF\n\n         ELSE\nC-----  Orthotropic elastic stiffness matrix\n            DLOCAL(1,1)=PROPS(1)\n            DLOCAL(1,2)=PROPS(2)\n            DLOCAL(2,1)=PROPS(2)\n            DLOCAL(2,2)=PROPS(3)\n\n            DLOCAL(1,3)=PROPS(4)\n            DLOCAL(3,1)=PROPS(4)\n            DLOCAL(2,3)=PROPS(5)\n            DLOCAL(3,2)=PROPS(5)\n            DLOCAL(3,3)=PROPS(6)\n\n            DLOCAL(4,4)=PROPS(7)\n            DLOCAL(5,5)=PROPS(8)\n            DLOCAL(6,6)=PROPS(9)\n\n         END IF\n\n      ELSE\nC-----  General fully anisotropic elastic stiffness matrix (Voigt symmetric storage)\n         ID=0\n         DO J=1,6\n            DO I=1,J\n               ID=ID+1\n               DLOCAL(I,J)=PROPS(ID)\n               DLOCAL(J,I)=DLOCAL(I,J)\n            END DO\n         END DO\n      END IF\n\nC-----  Call orientation subroutine to build crystal rotation matrix ROTATE\n      CALL ROTATION (PROPS(57), ROTATE)\n\nC-----  Construct Voigt rotation transformation matrix ROTD for 6-component stiffness rotation\n      DO J=1,3\n         J1=1+J\/3\n         J2=2+J\/2\n         DO I=1,3\n            I1=1+I\/3\n            I2=2+I\/2\n            ROTD(I,J)=ROTATE(I,J)**2\n            ROTD(I,J+3)=2.*ROTATE(I,J1)*ROTATE(I,J2)\n            ROTD(I+3,J)=ROTATE(I1,J)*ROTATE(I2,J)\n            ROTD(I+3,J+3)=ROTATE(I1,J1)*ROTATE(I2,J2)+\n     2                    ROTATE(I1,J2)*ROTATE(I2,J1)\n         END DO\n      END DO\n\nC-----  Rotate local crystal stiffness DLOCAL to global Cartesian stiffness D\nC     D = ROTD \u00b7 DLOCAL \u00b7 ROTD\u1d40\n      DO J=1,6\n         DO I=1,6\n            D(I,J)=0.\n         END DO\n      END DO\n\n      DO J=1,6\n         DO I=1,J\n            DO K=1,6\n               DO L=1,6\n                  D(I,J)=D(I,J)+DLOCAL(K,L)*ROTD(I,K)*ROTD(J,L)\n               END DO\n            END DO\n            D(J,I)=D(I,J)\n         END DO\n      END DO\n\nC-----  Read number of independent slip families NSET\n      NSET=NINT(PROPS(25))\n      IF (NSET.LT.1) THEN\n         WRITE (6,*) '***ERROR - zero slip families defined'\n         STOP\n      ELSE IF (NSET.GT.3) THEN\n         WRITE (6,*)\n     2     '***ERROR - maximum supported slip family count is 3'\n         STOP\n      END IF\n\nC-----  Implicit integration weighting factor \u03b8\n      THETA=PROPS(145)\n\nC-----  Finite geometry switch NLGEOM\n      IF (PROPS(146).EQ.0.) THEN\n         NLGEOM=0\n      ELSE\n         NLGEOM=1\n      END IF\n\nC-----  Newton iteration enable flag ITRATN\n      IF (PROPS(153).EQ.0.) THEN\n         ITRATN=0\n      ELSE\n         ITRATN=1\n      END IF\n\n      ITRMAX=NINT(PROPS(154))\n      GAMERR=PROPS(155)\n\nC     Initialize iteration residual storage\n      NITRTN=-1\n      DO I=1,NTENS\n         DSOLD(I)=0.\n      END DO\n      DO J=1,ND\n         DGAMOD(J)=0.\n         DTAUOD(J)=0.\n         DGSPOD(J)=0.\n         DO I=1,3\n            DSPNRO(I,J)=0.\n            DSPDRO(I,J)=0.\n         END DO\n      END DO\n\nC-----  Compute material spin increment DSPIN from incremental rotation matrix DROT (finite deformation only)\n      IF (NLGEOM.NE.0) THEN\n         DO J=1,3\n            DO I=1,3\n               TERM(I,J)=DROT(J,I)\n               TRM0(I,J)=DROT(J,I)\n            END DO\n            TERM(J,J)=TERM(J,J)+1.D0\n            TRM0(J,J)=TRM0(J,J)-1.D0\n         END DO\n         CALL LUDCMP (TERM, 3, 3, ITRM, DDCMP)\n         DO J=1,3\n            CALL LUBKSB (TERM, 3, 3, ITRM, TRM0(1,J))\n         END DO\n         DSPIN(1)=TRM0(2,1)-TRM0(1,2)\n         DSPIN(2)=TRM0(1,3)-TRM0(3,1)\n         DSPIN(3)=TRM0(3,2)-TRM0(2,3)\n      END IF\n\nC-----  Volumetric strain increment trace\n      DEV=0.D0\n      DO I=1,NDI\n         DEV=DEV+DSTRAN(I)\n      END DO\n\nC-----  Newton-Rhapson iteration loop entry label\n1000  CONTINUE\n      NITRTN=NITRTN+1\n\nC-----  Initialize slip system data at first analysis increment (TOTALTIME=0)\n      IF (STATEV(1).EQ.0.) THEN\n         NSLPTL=0\n         DO I=1,NSET\n            ISPNOR(1)=NINT(PROPS(25+8*I))\n            ISPNOR(2)=NINT(PROPS(26+8*I))\n            ISPNOR(3)=NINT(PROPS(27+8*I))\n            ISPDIR(1)=NINT(PROPS(28+8*I))\n            ISPDIR(2)=NINT(PROPS(29+8*I))\n            ISPDIR(3)=NINT(PROPS(30+8*I))\nC Generate all independent slip systems for current slip family\n            CALL SLIPSYS (ISPDIR, ISPNOR, NSLIP(I), SLPDIR(1,NSLPTL+1),\n     2                    SLPNOR(1,NSLPTL+1), ROTATE)\n            NSLPTL=NSLPTL+NSLIP(I)\n         END DO\nC Check array dimension ND is sufficient for total slip systems\n         IF (ND.LT.NSLPTL) THEN\n            WRITE (6,*)\n     2 '***ERROR - Parameter ND smaller than total independent slip systems NSLPTL'\n            STOP\n         END IF\n\nC-----  Assemble Schmid slip deformation tensor SLPDEF for all slip systems\n         DO J=1,NSLPTL\n            SLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)\n            SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)\n            SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)\n            SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)\n            SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)\n            SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J)\n         END DO\n\nC-----  Write slip system metadata to state variables\n         STATEV(NSTATV)=FLOAT(NSLPTL)\n         DO I=1,NSET\n            STATEV(NSTATV-4+I)=FLOAT(NSLIP(I))\n         END DO\n\nC Store initial slip plane normals and slip directions to SDVs\n         IDNOR=3*NSLPTL\n         IDDIR=6*NSLPTL\n         DO J=1,NSLPTL\n            DO I=1,3\n               IDNOR=IDNOR+1\n               STATEV(IDNOR)=SLPNOR(I,J)\n               IDDIR=IDDIR+1\n               STATEV(IDDIR)=SLPDIR(I,J)\n            END DO\n         END DO\n\nC Assign initial critical slip strength g0 via GSLPINIT\n         CALL GSLPINIT (STATEV(1), NSLIP, NSLPTL, NSET, PROPS(97))\n\nC Initialize accumulated shear strain and cumulative absolute slip\n         DO I=1,NSLPTL\n            STATEV(NSLPTL+I)=0.\n            STATEV(9*NSLPTL+I)=0.\n         END DO\n         STATEV(10*NSLPTL+1)=0.\n\nC Calculate initial resolved shear stress \u03c4_\u03b1 = \u03c3 : S_\u03b1\n         DO I=1,NSLPTL\n            TERM1=0.\n            DO J=1,NTENS\n               IF (J.LE.NDI) THEN\n                  TERM1=TERM1+SLPDEF(J,I)*STRESS(J)\n               ELSE\n                  TERM1=TERM1+SLPDEF(J-NDI+3,I)*STRESS(J)\n               END IF\n            END DO\n            STATEV(2*NSLPTL+I)=TERM1\n         END DO\n\n      ELSE\nC-----  Post-initialization: Read existing slip system data from state variables\n         NSLPTL=NINT(STATEV(NSTATV))\n         DO I=1,NSET\n            NSLIP(I)=NINT(STATEV(NSTATV-4+I))\n         END DO\n\nC Recover current slip plane normals and slip directions from SDVs\n         IDNOR=3*NSLPTL\n         IDDIR=6*NSLPTL\n         DO J=1,NSLPTL\n            DO I=1,3\n               IDNOR=IDNOR+1\n               SLPNOR(I,J)=STATEV(IDNOR)\n               IDDIR=IDDIR+1\n               SLPDIR(I,J)=STATEV(IDDIR)\n            END DO\n         END DO\n\nC Rebuild Schmid tensor SLPDEF from updated slip vectors\n         DO J=1,NSLPTL\n            SLPDEF(1,J)=SLPDIR(1,J)*SLPNOR(1,J)\n            SLPDEF(2,J)=SLPDIR(2,J)*SLPNOR(2,J)\n            SLPDEF(3,J)=SLPDIR(3,J)*SLPNOR(3,J)\n            SLPDEF(4,J)=SLPDIR(1,J)*SLPNOR(2,J)+SLPDIR(2,J)*SLPNOR(1,J)\n            SLPDEF(5,J)=SLPDIR(1,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(1,J)\n            SLPDEF(6,J)=SLPDIR(2,J)*SLPNOR(3,J)+SLPDIR(3,J)*SLPNOR(2,J)\n         END DO\n      END IF\n\nC-----  Construct slip spin tensor SLPSPN (finite rotation kinematics only)\n      IF (NLGEOM.NE.0) THEN\n         DO J=1,NSLPTL\n            SLPSPN(1,J)=0.5*(SLPDIR(1,J)*SLPNOR(2,J)-\n     2                       SLPDIR(2,J)*SLPNOR(1,J))\n            SLPSPN(2,J)=0.5*(SLPDIR(3,J)*SLPNOR(1,J)-\n     2                       SLPDIR(1,J)*SLPNOR(3,J))\n            SLPSPN(3,J)=0.5*(SLPDIR(2,J)*SLPNOR(3,J)-\n     2                       SLPDIR(3,J)*SLPNOR(2,J))\n         END DO\n      END IF\n\nC-----  Compute DDEMSD coupling matrix D:S_\u03b1 + \u03c3\u00d7W_\u03b1 (finite rotation extra term)\n      DO J=1,NSLPTL\n         DO I=1,6\n            DDEMSD(I,J)=0.\n            DO K=1,6\n               DDEMSD(I,J)=DDEMSD(I,J)+D(K,I)*SLPDEF(K,J)\n            END DO\n         END DO\n      END IF\n\n      IF (NLGEOM.NE.0) THEN\n         DO J=1,NSLPTL\n            DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(1,J)*STRESS(1)\n            DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(2,J)*STRESS(1)\n            IF (NDI.GT.1) THEN\n               DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(1,J)*STRESS(2)\n               DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(3,J)*STRESS(2)\n            END IF\n            IF (NDI.GT.2) THEN\n               DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(2,J)*STRESS(3)\n               DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(3,J)*STRESS(3)\n            END IF\n            IF (NSHR.GE.1) THEN\n               DDEMSD(1,J)=DDEMSD(1,J)+SLPSPN(1,J)*STRESS(NDI+1)\n               DDEMSD(2,J)=DDEMSD(2,J)-SLPSPN(1,J)*STRESS(NDI+1)\n               DDEMSD(5,J)=DDEMSD(5,J)-SLPSPN(3,J)*STRESS(NDI+1)\n               DDEMSD(6,J)=DDEMSD(6,J)+SLPSPN(2,J)*STRESS(NDI+1)\n            END IF\n            IF (NSHR.GE.2) THEN\n               DDEMSD(1,J)=DDEMSD(1,J)-SLPSPN(2,J)*STRESS(NDI+2)\n               DDEMSD(3,J)=DDEMSD(3,J)+SLPSPN(2,J)*STRESS(NDI+2)\n               DDEMSD(4,J)=DDEMSD(4,J)+SLPSPN(3,J)*STRESS(NDI+2)\n               DDEMSD(6,J)=DDEMSD(6,J)-SLPSPN(1,J)*STRESS(NDI+2)\n            END IF\n            IF (NSHR.EQ.3) THEN\n               DDEMSD(2,J)=DDEMSD(2,J)+SLPSPN(3,J)*STRESS(NDI+3)\n               DDEMSD(3,J)=DDEMSD(3,J)-SLPSPN(3,J)*STRESS(NDI+3)\n               DDEMSD(4,J)=DDEMSD(4,J)-SLPSPN(2,J)*STRESS(NDI+3)\n               DDEMSD(5,J)=DDEMSD(5,J)+SLPSPN(1,J)*STRESS(NDI+3)\n            END IF\n         END DO\n      END IF\n\nC-----  Evaluate slip shear strain rate FSLIP and derivative dF\/dX via STRAINRATE\n      ID=1\n      DO I=1,NSET\n         IF (I.GT.1) ID=ID+NSLIP(I-1)\n         CALL STRAINRATE (STATEV(NSLPTL+ID), STATEV(2*NSLPTL+ID),\n     2                    STATEV(ID), NSLIP(I), FSLIP(ID), DFDXSP(ID),\n     3                    PROPS(65+8*I))\n      END DO\n\nC-----  Assemble self\/latent hardening interaction matrix H\n       CALL LATENTHARDEN (STATEV(NSLPTL+1), STATEV(2*NSLPTL+1),\n     2                   STATEV(1), STATEV(9*NSLPTL+1),\n     3                   STATEV(10*NSLPTL+1), NSLIP, NSLPTL,\n     4                   NSET, H(1,1), PROPS(97), ND)\n\nC-----  Build Jacobian matrix for slip shear strain increment solve\n      TERM1=THETA*DTIME\n      DO I=1,NSLPTL\n         TAUSLP=STATEV(2*NSLPTL+I)\n         GSLIP=STATEV(I)\n         X=TAUSLP\/GSLIP\n         TERM2=TERM1*DFDXSP(I)\/GSLIP\n         TERM3=TERM1*X*DFDXSP(I)\/GSLIP\n         DO J=1,NSLPTL\n            TERM4=0.\n            DO K=1,6\n               TERM4=TERM4+DDEMSD(K,I)*SLPDEF(K,J)\n            END DO\n            WORKST(I,J)=TERM2*TERM4+H(I,J)*TERM3*DSIGN(1.D0,FSLIP(J))\n            IF (NITRTN.GT.0) WORKST(I,J)=WORKST(I,J)+TERM3*DHDGDG(I,J)\n         END DO\n         WORKST(I,I)=WORKST(I,I)+1.\n      END DO\n\nC LU decomposition of slip system Jacobian matrix\n      CALL LUDCMP (WORKST, NSLPTL, ND, INDX, DDCMP)\n\n\nC-----  Increment of shear strain in a slip system: DGAMMA\n      TERM1=THETA*DTIME\n      DO I=1,NSLPTL\n\n         IF (NITRTN.EQ.0) THEN\n            TAUSLP=STATEV(2*NSLPTL+I)\n            GSLIP=STATEV(I)\n            X=TAUSLP\/GSLIP\n            TERM2=TERM1*DFDXSP(I)\/GSLIP\n\n            DGAMMA(I)=0.\n            DO J=1,NDI\n               DGAMMA(I)=DGAMMA(I)+DDEMSD(J,I)*DSTRAN(J)\n            END DO\n\n            IF (NSHR.GT.0) THEN\n               DO J=1,NSHR\n                  DGAMMA(I)=DGAMMA(I)+DDEMSD(J+3,I)*DSTRAN(J+NDI)\n               END DO\n            END IF\n\n            DGAMMA(I)=DGAMMA(I)*TERM2+FSLIP(I)*DTIME\n\n         ELSE\n            DGAMMA(I)=TERM1*(FSLIP(I)-FSLIP1(I))+FSLIP1(I)*DTIME\n     2                -DGAMOD(I)\n\n         END IF\n\n      END DO\n\n      CALL LUBKSB (WORKST, NSLPTL, ND, INDX, DGAMMA)\n\n      DO I=1,NSLPTL\n         DGAMMA(I)=DGAMMA(I)+DGAMOD(I)\n      END DO\n\nC-----  Update the shear strain in a slip system: STATEV(NSLPTL+1) - \nC     STATEV(2*NSLPTL)\nC\n      DO I=1,NSLPTL\n         STATEV(NSLPTL+I)=STATEV(NSLPTL+I)+DGAMMA(I)-DGAMOD(I)  ! \u66f4\u65b0\u5207\u5e94\u53d8\n      END DO\n\nC-----  Increment of current strength in a slip system: DGSLIP\n      DO I=1,NSLPTL\n         DGSLIP(I)=0.\n         DO J=1,NSLPTL\n            DGSLIP(I)=DGSLIP(I)+H(I,J)*ABS(DGAMMA(J))  !\u6c42\u6ed1\u79fb\u7cfb\u5f3a\u5ea6\u589e\u91cf\n         END DO\n      END DO\n\nC-----  Update the current strength in a slip system: STATEV(1) - \nC     STATEV(NSLPTL)\nC\n      DO I=1,NSLPTL\n         STATEV(I)=STATEV(I)+DGSLIP(I)-DGSPOD(I)  !\u5e76\u66f4\u65b0\u5f3a\u5ea6\n      END DO\n\nC-----  Increment of strain associated with lattice stretching: DELATS\n      DO J=1,6\n         DELATS(J)=0.\n      END DO\n\n      DO J=1,3\n         IF (J.LE.NDI) DELATS(J)=DSTRAN(J)\n         DO I=1,NSLPTL\n            DELATS(J)=DELATS(J)-SLPDEF(J,I)*DGAMMA(I)\n         END DO\n      END DO\n\n      DO J=1,3\n         IF (J.LE.NSHR) DELATS(J+3)=DSTRAN(J+NDI)\n         DO I=1,NSLPTL\n            DELATS(J+3)=DELATS(J+3)-SLPDEF(J+3,I)*DGAMMA(I)\n         END DO\n      END DO\n\nC-----  Increment of deformation gradient associated with lattice \nC     stretching in the current state, i.e. the velocity gradient \nC     (associated with lattice stretching) times the increment of time:\nC     DVGRAD (only needed for finite rotation)\nC\n      IF (NLGEOM.NE.0) THEN\n         DO J=1,3\n            DO I=1,3\n               IF (I.EQ.J) THEN\n                  DVGRAD(I,J)=DELATS(I)\n               ELSE\n                  DVGRAD(I,J)=DELATS(I+J+1)\n               END IF\n            END DO\n         END DO\n\n         DO J=1,3\n            DO I=1,J\n               IF (J.GT.I) THEN\n                  IJ2=I+J-2\n                  IF (MOD(IJ2,2).EQ.1) THEN\n                     TERM1=1.\n                  ELSE\n                     TERM1=-1.\n                  END IF\n\n                  DVGRAD(I,J)=DVGRAD(I,J)+TERM1*DSPIN(IJ2)\n                  DVGRAD(J,I)=DVGRAD(J,I)-TERM1*DSPIN(IJ2)\n\n                  DO K=1,NSLPTL\n                     DVGRAD(I,J)=DVGRAD(I,J)-TERM1*DGAMMA(K)*\n     2                                       SLPSPN(IJ2,K)\n                     DVGRAD(J,I)=DVGRAD(J,I)+TERM1*DGAMMA(K)*\n     2                                       SLPSPN(IJ2,K)\n                  END DO\n               END IF\n\n            END DO\n         END DO\n\n      END IF\n\nC-----  Increment of resolved shear stress in a slip system: DTAUSP\n      DO I=1,NSLPTL\n         DTAUSP(I)=0.\n         DO J=1,6\n            DTAUSP(I)=DTAUSP(I)+DDEMSD(J,I)*DELATS(J)\n         END DO\n      END DO\n\nC-----  Update the resolved shear stress in a slip system: \nC     STATEV(2*NSLPTL+1) - STATEV(3*NSLPTL)\nC\n      DO I=1,NSLPTL\n         STATEV(2*NSLPTL+I)=STATEV(2*NSLPTL+I)+DTAUSP(I)-DTAUOD(I)\n      END DO\n\nC-----  Increment of stress: DSTRES\n      IF (NLGEOM.EQ.0) THEN\n         DO I=1,NTENS\n            DSTRES(I)=0.\n         END DO\n      ELSE\n         DO I=1,NTENS\n            DSTRES(I)=-STRESS(I)*DEV\n         END DO\n      END IF\n\n      DO I=1,NDI\n         DO J=1,NDI\n            DSTRES(I)=DSTRES(I)+D(I,J)*DSTRAN(J)\n         END DO\n\n         IF (NSHR.GT.0) THEN\n            DO J=1,NSHR\n               DSTRES(I)=DSTRES(I)+D(I,J+3)*DSTRAN(J+NDI)\n            END DO\n         END IF\n\n         DO J=1,NSLPTL\n            DSTRES(I)=DSTRES(I)-DDEMSD(I,J)*DGAMMA(J)\n         END DO\n      END DO\n\n      IF (NSHR.GT.0) THEN\n         DO I=1,NSHR\n\n            DO J=1,NDI\n               DSTRES(I+NDI)=DSTRES(I+NDI)+D(I+3,J)*DSTRAN(J)\n            END DO\n\n            DO J=1,NSHR\n               DSTRES(I+NDI)=DSTRES(I+NDI)+D(I+3,J+3)*DSTRAN(J+NDI)\n            END DO\n\n            DO J=1,NSLPTL\n               DSTRES(I+NDI)=DSTRES(I+NDI)-DDEMSD(I+3,J)*DGAMMA(J)\n            END DO\n\n         END DO\n      END IF\n\nC-----  Update the stress: STRESS\n      DO I=1,NTENS\n         STRESS(I)=STRESS(I)+DSTRES(I)-DSOLD(I)\n      END DO\n\nC-----  Increment of normal to a slip plane and a slip direction (only \nC     needed for finite rotation)\nC\n      IF (NLGEOM.NE.0) THEN\n         DO J=1,NSLPTL\n            DO I=1,3\n               DSPNOR(I,J)=0.\n               DSPDIR(I,J)=0.\n\n               DO K=1,3\n                  DSPNOR(I,J)=DSPNOR(I,J)-SLPNOR(K,J)*DVGRAD(K,I)\n                  DSPDIR(I,J)=DSPDIR(I,J)+SLPDIR(K,J)*DVGRAD(I,K)\n               END DO\n\n            END DO\n         END DO\n\nC-----  Update the normal to a slip plane and a slip direction (only \nC     needed for finite rotation)\nC\n         IDNOR=3*NSLPTL\n         IDDIR=6*NSLPTL\n         DO J=1,NSLPTL\n            DO I=1,3\n               IDNOR=IDNOR+1\n               STATEV(IDNOR)=STATEV(IDNOR)+DSPNOR(I,J)-DSPNRO(I,J)\n\n               IDDIR=IDDIR+1\n               STATEV(IDDIR)=STATEV(IDDIR)+DSPDIR(I,J)-DSPDRO(I,J)\n            END DO\n         END DO\n\n      END IF\n\nC-----  Iteration ?\n      IF (ITRATN.NE.0) THEN\n\nC-----  Save solutions (without iteration):\nC            Shear strain-rate in a slip system FSLIP1\nC            Current strength in a slip system GSLP1\nC            Shear strain in a slip system GAMMA1\nC            Resolved shear stress in a slip system TAUSP1\nC            Normal to a slip plane SPNOR1\nC            Slip direction SPDIR1\nC            Stress STRES1\nC            Jacobian matrix DDSDE1\nC\n         IF (NITRTN.EQ.0) THEN\n\n            IDNOR=3*NSLPTL\n            IDDIR=6*NSLPTL\n            DO J=1,NSLPTL\n               FSLIP1(J)=FSLIP(J)\n               GSLP1(J)=STATEV(J)\n               GAMMA1(J)=STATEV(NSLPTL+J)\n               TAUSP1(J)=STATEV(2*NSLPTL+J)\n               DO I=1,3\n                  IDNOR=IDNOR+1\n                  SPNOR1(I,J)=STATEV(IDNOR)\n\n                  IDDIR=IDDIR+1\n                  SPDIR1(I,J)=STATEV(IDDIR)\n               END DO\n            END DO\n         END IF\n\nC-----  Increments of stress DSOLD, and solution dependent state \nC     variables DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO (for the next \nC     iteration)\nC\n         DO I=1,NTENS\n            DSOLD(I)=DSTRES(I)\n         END DO\n\n         DO J=1,NSLPTL\n            DGAMOD(J)=DGAMMA(J)\n            DTAUOD(J)=DTAUSP(J)\n            DGSPOD(J)=DGSLIP(J)\n            DO I=1,3\n               DSPNRO(I,J)=DSPNOR(I,J)\n               DSPDRO(I,J)=DSPDIR(I,J)\n            END DO\n         END DO\n\nC-----  Check if the iteration solution converges\n         IDBACK=0\n         ID=0\n         DO I=1,NSET\n            DO J=1,NSLIP(I)\n               ID=ID+1\n               X=STATEV(2*NSLPTL+ID)\/STATEV(ID)\n               RESIDU=THETA*DTIME*F(X,PROPS(65+8*I))+DTIME*(1.0-THETA)*\n     2                FSLIP1(ID)-DGAMMA(ID)\n               IF (ABS(RESIDU).GT.GAMERR) IDBACK=1\n            END DO\n         END DO\n\n         IF (IDBACK.NE.0.AND.NITRTN.LT.ITRMAX) THEN\nC-----  Iteration: arrays for iteration\nCFIXA\n            CALL ITERATION (STATEV(NSLPTL+1), STATEV(2*NSLPTL+1), \n     2                      STATEV(1), STATEV(9*NSLPTL+1), \n     3                      STATEV(10*NSLPTL+1), NSLPTL, \n     4                      NSET, NSLIP, ND, PROPS(97), DGAMOD,\n     5                      DHDGDG)\nCFIXB\n\n            GO TO 1000\n\n         ELSE IF (NITRTN.GE.ITRMAX) THEN\nC-----  Solution not converge within maximum number of iteration (the \nC     solution without iteration will be used)\n\n            IDNOR=3*NSLPTL\n            IDDIR=6*NSLPTL\n            DO J=1,NSLPTL\n               STATEV(J)=GSLP1(J)\n               STATEV(NSLPTL+J)=GAMMA1(J)\n               STATEV(2*NSLPTL+J)=TAUSP1(J)\n\n               DO I=1,3\n                  IDNOR=IDNOR+1\n                  STATEV(IDNOR)=SPNOR1(I,J)\n\n                  IDDIR=IDDIR+1\n                  STATEV(IDDIR)=SPDIR1(I,J)\n               END DO\n            END DO\n\n         END IF\n\n      END IF\n\nC-----  Total cumulative shear strains on all slip systems (sum of the \nC       absolute values of shear strains in all slip systems)\nCFIX--  Total cumulative shear strains on each slip system (sum of the \nCFIX    absolute values of shear strains in each individual slip system)\nC\n      DO I=1,NSLPTL\nCFIXA\n         STATEV(10*NSLPTL+1)=STATEV(10*NSLPTL+1)+ABS(DGAMMA(I))\n         STATEV(9*NSLPTL+I)=STATEV(9*NSLPTL+I)+ABS(DGAMMA(I))\nCFIXB\n      END DO\n\n      RETURN\n      END\n\n\nC---------------------\u4ee5\u4e0b\u662fVUMAT\u8c03\u7528\u7684\u51fd\u6570--------------------------------\nC     \n\n      SUBROUTINE ROTATION (PROP, ROTATE)\n\nC-----  This subroutine calculates the rotation matrix, i.e. the \nC     direction cosines of cubic crystal [100], [010] and [001] \nC     directions in global system\n\nC-----  The rotation matrix is stored in the array ROTATE.\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      DIMENSION PROP(16), ROTATE(3,3), TERM1(3,3), TERM2(3,3), INDX(3) \n\nC-----  Subroutines:\nC\nC       CROSS  -- cross product of two vectors\nC\nC       LUDCMP -- LU decomposition\nC\nC       LUBKSB -- linear equation solver based on LU decomposition \nC                 method (must call LUDCMP first)\n\n\nC-----  PROP -- constants characterizing the crystal orientation \nC               (INPUT)\nC\nC            PROP(1) - PROP(3) -- direction of the first vector in \nC                                 local cubic crystal system\nC            PROP(4) - PROP(6) -- direction of the first vector in \nC                                 global system\nC\nC            PROP(9) - PROP(11)-- direction of the second vector in \nC                                 local cubic crystal system\nC            PROP(12)- PROP(14)-- direction of the second vector in \nC                                 global system\nC\nC-----  ROTATE -- rotation matrix (OUTPUT):\nC\nC            ROTATE(i,1) -- direction cosines of direction [1 0 0] in \nC                           local cubic crystal system\nC            ROTATE(i,2) -- direction cosines of direction [0 1 0] in \nC                           local cubic crystal system\nC            ROTATE(i,3) -- direction cosines of direction [0 0 1] in \nC                           local cubic crystal system\n\nC-----  local matrix: TERM1\n      CALL CROSS (PROP(1), PROP(9), TERM1, ANGLE1)\n\nC-----  LU decomposition of TERM1\n      CALL LUDCMP (TERM1, 3, 3, INDX, DCMP)\n\nC-----  inverse matrix of TERM1: TERM2\n      DO J=1,3\n         DO I=1,3\n            IF (I.EQ.J) THEN\n               TERM2(I,J)=1.\n            ELSE\n               TERM2(I,J)=0.\n            END IF\n         END DO\n      END DO\n\n      DO J=1,3\n         CALL LUBKSB (TERM1, 3, 3, INDX, TERM2(1,J))\n      END DO\n\nC-----  global matrix: TERM1\n      CALL CROSS (PROP(4), PROP(12), TERM1, ANGLE2)\n\nC-----  Check: the angle between first and second vector in local and \nC     global systems must be the same.  The relative difference must be\nC     less than 0.1%.\nC\n      IF (ABS(ANGLE1\/ANGLE2-1.).GT.0.001) THEN \n         WRITE (6,*) \n     2      '***ERROR - angles between two vectors are not the same'\n         STOP\n      END IF\n\nC-----  rotation matrix: ROTATE\n      DO J=1,3\n         DO I=1,3\n            ROTATE(I,J)=0.\n            DO K=1,3\n               ROTATE(I,J)=ROTATE(I,J)+TERM1(I,K)*TERM2(K,J)\n            END DO\n         END DO\n      END DO\n\n      RETURN\n      END\n\n\nC-----------------------------------\n\n\n           SUBROUTINE CROSS (A, B, C, ANGLE)\n\nC-----  (1) normalize vectors A and B to unit vectors\nC       (2) store A, B and A*B (cross product) in C\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\n           DIMENSION A(3), B(3), C(3,3)\n\n           SUM1=SQRT(A(1)**2+A(2)**2+A(3)**2)\n           SUM2=SQRT(B(1)**2+B(2)**2+B(3)**2)\n\n           IF (SUM1.EQ.0.) THEN\n              WRITE (6,*) '***ERROR - first vector is zero'\n              STOP\n           ELSE\n              DO I=1,3\n                 C(I,1)=A(I)\/SUM1\n              END DO\n           END IF\n\n           IF (SUM2.EQ.0.) THEN\n              WRITE (6,*) '***ERROR - second vector is zero'\n              STOP\n           ELSE\n              DO I=1,3\n                 C(I,2)=B(I)\/SUM2\n              END DO\n           END IF\n\n           ANGLE=0.\n           DO I=1,3\n              ANGLE=ANGLE+C(I,1)*C(I,2)\n           END DO\n           ANGLE=ACOS(ANGLE)\n\n           C(1,3)=C(2,1)*C(3,2)-C(3,1)*C(2,2)\n           C(2,3)=C(3,1)*C(1,2)-C(1,1)*C(3,2)\n           C(3,3)=C(1,1)*C(2,2)-C(2,1)*C(1,2)\n           SUM3=SQRT(C(1,3)**2+C(2,3)**2+C(3,3)**2)\n           IF (SUM3.LT.1.E-8) THEN\n              WRITE (6,*) \n     2           '***ERROR - first and second vectors are parallel'\n               STOP\n            END IF\n\n           RETURN\n           END\n\n\nC----------------------------------------------------------------------\n\n\n      SUBROUTINE SLIPSYS (ISPDIR, ISPNOR, NSLIP, SLPDIR, SLPNOR, \n     2                    ROTATE)\n\nC-----  This subroutine generates all slip systems in the same set for \nC     a CUBIC crystal.  For other crystals (e.g., HCP, Tetragonal, \nC     Orthotropic, ...), it has to be modified to include the effect of\nC     crystal aspect ratio.\n\nC-----  Denote s as a slip direction and m as normal to a slip plane.  \nC     In a cubic crystal, (s,-m), (-s,m) and (-s,-m) are NOT considered\nC     independent of (s,m).\n\nC-----  Subroutines:  LINE1 and LINE\n\nC-----  Variables:\nC\nC     ISPDIR -- a typical slip direction in this set of slip systems \nC               (integer)  (INPUT)\nC     ISPNOR -- a typical normal to slip plane in this set of slip \nC               systems (integer)  (INPUT)\nC     NSLIP  -- number of independent slip systems in this set \nC               (OUTPUT)\nC     SLPDIR -- unit vectors of all slip directions  (OUTPUT)\nC     SLPNOR -- unit normals to all slip planes  (OUTPUT)\nC     ROTATE -- rotation matrix (INPUT)\nC          ROTATE(i,1) -- direction cosines of [100] in global system\nC          ROTATE(i,2) -- direction cosines of [010] in global system\nC          ROTATE(i,3) -- direction cosines of [001] in global system\nC\nC     NSPDIR -- number of all possible slip directions in this set\nC     NSPNOR -- number of all possible slip planes in this set\nC     IWKDIR -- all possible slip directions (integer)\nC     IWKNOR -- all possible slip planes (integer)\n\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      DIMENSION ISPDIR(3), ISPNOR(3), SLPDIR(3,50), SLPNOR(3,50), \n     *          ROTATE(3,3), IWKDIR(3,24), IWKNOR(3,24), TERM(3)\n\n      NSLIP=0\n      NSPDIR=0\n      NSPNOR=0\n\nC-----  Generating all possible slip directions in this set\nC\nC       Denote the slip direction by [lmn].  I1 is the minimum of the \nC     absolute value of l, m and n, I3 is the maximum and I2 is the \nC     mode, e.g. (1 -3 2), I1=1, I2=2 and I3=3.  I1<=I2<=I3.\n\n      I1=MIN(IABS(ISPDIR(1)),IABS(ISPDIR(2)),IABS(ISPDIR(3)))\n      I3=MAX(IABS(ISPDIR(1)),IABS(ISPDIR(2)),IABS(ISPDIR(3)))\n      I2=IABS(ISPDIR(1))+IABS(ISPDIR(2))+IABS(ISPDIR(3))-I1-I3\n\n      RMODIR=SQRT(FLOAT(I1*I1+I2*I2+I3*I3))\n\nC     I1=I2=I3=0\n      IF (I3.EQ.0) THEN \n         WRITE (6,*) '***ERROR - slip direction is [000]'\n         STOP\n\nC     I1=I2=0, I3>0   ---   [001] type\n      ELSE IF (I2.EQ.0) THEN\n         NSPDIR=3\n         DO J=1,3\n            DO I=1,3\n               IWKDIR(I,J)=0\n               IF (I.EQ.J) IWKDIR(I,J)=I3\n            END DO\n         END DO\n\nC     I1=0, I3>=I2>0\n      ELSE IF (I1.EQ.0) THEN\n\nC        I1=0, I3=I2>0   ---   [011] type\n         IF (I2.EQ.I3) THEN\n            NSPDIR=6\n            DO J=1,6\n               DO I=1,3\n                  IWKDIR(I,J)=I2\n                  IF (I.EQ.J.OR.J-I.EQ.3) IWKDIR(I,J)=0\n                  IWKDIR(1,6)=-I2\n                  IWKDIR(2,4)=-I2\n                  IWKDIR(3,5)=-I2\n               END DO\n            END DO\n\nC        I1=0, I3>I2>0   ---   [012] type\n         ELSE\n            NSPDIR=12\n            CALL LINE1 (I2, I3, IWKDIR(1,1), 1)\n            CALL LINE1 (I3, I2, IWKDIR(1,3), 1)\n            CALL LINE1 (I2, I3, IWKDIR(1,5), 2)\n            CALL LINE1 (I3, I2, IWKDIR(1,7), 2)\n            CALL LINE1 (I2, I3, IWKDIR(1,9), 3)\n            CALL LINE1 (I3, I2, IWKDIR(1,11), 3)\n\n         END IF\n\nC     I1=I2=I3>0   ---   [111] type\n      ELSE IF (I1.EQ.I3) THEN\n         NSPDIR=4\n         CALL LINE (I1, I1, I1, IWKDIR)\n\nC     I3>I2=I1>0   ---   [112] type\n      ELSE IF (I1.EQ.I2) THEN\n         NSPDIR=12\n         CALL LINE (I1, I1, I3, IWKDIR(1,1))\n         CALL LINE (I1, I3, I1, IWKDIR(1,5))\n         CALL LINE (I3, I1, I1, IWKDIR(1,9))\n\nC     I3=I2>I1>0   ---   [122] type\n      ELSE IF (I2.EQ.I3) THEN\n         NSPDIR=12\n         CALL LINE (I1, I2, I2, IWKDIR(1,1))\n         CALL LINE (I2, I1, I2, IWKDIR(1,5))\n         CALL LINE (I2, I2, I1, IWKDIR(1,9))\n\nC     I3>I2>I1>0   ---   [123] type\n      ELSE\n         NSPDIR=24\n         CALL LINE (I1, I2, I3, IWKDIR(1,1))\n         CALL LINE (I3, I1, I2, IWKDIR(1,5))\n         CALL LINE (I2, I3, I1, IWKDIR(1,9))\n         CALL LINE (I1, I3, I2, IWKDIR(1,13))\n         CALL LINE (I2, I1, I3, IWKDIR(1,17))\n         CALL LINE (I3, I2, I1, IWKDIR(1,21))\n\n      END IF\n\nC-----  Generating all possible slip planes in this set\nC\nC       Denote the normal to slip plane by (pqr).  J1 is the minimum of\nC     the absolute value of p, q and r, J3 is the maximum and J2 is the\nC     mode, e.g. (1 -2 1), J1=1, J2=1 and J3=2.  J1<=J2<=J3.\n\n      J1=MIN(IABS(ISPNOR(1)),IABS(ISPNOR(2)),IABS(ISPNOR(3)))\n      J3=MAX(IABS(ISPNOR(1)),IABS(ISPNOR(2)),IABS(ISPNOR(3)))\n      J2=IABS(ISPNOR(1))+IABS(ISPNOR(2))+IABS(ISPNOR(3))-J1-J3\n\n      RMONOR=SQRT(FLOAT(J1*J1+J2*J2+J3*J3))\n\n      IF (J3.EQ.0) THEN \n         WRITE (6,*) '***ERROR - slip plane is [000]'\n         STOP\n\nC     (001) type\n      ELSE IF (J2.EQ.0) THEN\n         NSPNOR=3\n         DO J=1,3\n            DO I=1,3\n               IWKNOR(I,J)=0\n               IF (I.EQ.J) IWKNOR(I,J)=J3\n            END DO\n         END DO\n\n      ELSE IF (J1.EQ.0) THEN\n\nC     (011) type\n         IF (J2.EQ.J3) THEN\n            NSPNOR=6\n            DO J=1,6\n               DO I=1,3\n                  IWKNOR(I,J)=J2\n                  IF (I.EQ.J.OR.J-I.EQ.3) IWKNOR(I,J)=0\n                  IWKNOR(1,6)=-J2\n                  IWKNOR(2,4)=-J2\n                  IWKNOR(3,5)=-J2\n               END DO\n            END DO\n\nC     (012) type\n         ELSE\n            NSPNOR=12\n            CALL LINE1 (J2, J3, IWKNOR(1,1), 1)\n            CALL LINE1 (J3, J2, IWKNOR(1,3), 1)\n            CALL LINE1 (J2, J3, IWKNOR(1,5), 2)\n            CALL LINE1 (J3, J2, IWKNOR(1,7), 2)\n            CALL LINE1 (J2, J3, IWKNOR(1,9), 3)\n            CALL LINE1 (J3, J2, IWKNOR(1,11), 3)\n\n         END IF\n\nC     (111) type\n      ELSE IF (J1.EQ.J3) THEN\n         NSPNOR=4\n         CALL LINE (J1, J1, J1, IWKNOR)\n\nC     (112) type\n      ELSE IF (J1.EQ.J2) THEN\n         NSPNOR=12\n         CALL LINE (J1, J1, J3, IWKNOR(1,1))\n         CALL LINE (J1, J3, J1, IWKNOR(1,5))\n         CALL LINE (J3, J1, J1, IWKNOR(1,9))\n\nC     (122) type\n      ELSE IF (J2.EQ.J3) THEN\n         NSPNOR=12\n         CALL LINE (J1, J2, J2, IWKNOR(1,1))\n         CALL LINE (J2, J1, J2, IWKNOR(1,5))\n         CALL LINE (J2, J2, J1, IWKNOR(1,9))\n\nC     (123) type\n      ELSE\n         NSPNOR=24\n         CALL LINE (J1, J2, J3, IWKNOR(1,1))\n         CALL LINE (J3, J1, J2, IWKNOR(1,5))\n         CALL LINE (J2, J3, J1, IWKNOR(1,9))\n         CALL LINE (J1, J3, J2, IWKNOR(1,13))\n         CALL LINE (J2, J1, J3, IWKNOR(1,17))\n         CALL LINE (J3, J2, J1, IWKNOR(1,21))\n\n      END IF\n\nC-----  Generating all slip systems in this set\nC\nC-----  Unit vectors in slip directions: SLPDIR, and unit normals to \nC     slip planes: SLPNOR in local cubic crystal system\nC\n      WRITE (6,*) '          '\n      WRITE (6,*) ' #          Slip plane          Slip direction'\n\n      DO J=1,NSPNOR\n         DO I=1,NSPDIR\n\n            IDOT=0\n            DO K=1,3\n               IDOT=IDOT+IWKDIR(K,I)*IWKNOR(K,J)\n            END DO\n\n            IF (IDOT.EQ.0) THEN\n               NSLIP=NSLIP+1\n               DO K=1,3\n                  SLPDIR(K,NSLIP)=IWKDIR(K,I)\/RMODIR\n                  SLPNOR(K,NSLIP)=IWKNOR(K,J)\/RMONOR\n               END DO\n\n               WRITE (6,10) NSLIP, \n     2                      (IWKNOR(K,J),K=1,3), (IWKDIR(K,I),K=1,3)\n\n            END IF\n\n         END DO\n      END DO\n10    FORMAT(1X,I2,9X,'(',3(1X,I2),1X,')',10X,'[',3(1X,I2),1X,']')\n\n      WRITE (6,*) 'Number of slip systems in this set = ',NSLIP\n      WRITE (6,*) '          '\n\n      IF (NSLIP.EQ.0) THEN\n         WRITE (6,*) \n     *      'There is no slip direction normal to the slip planes!'\n         STOP\n\n      ELSE\n\nC-----  Unit vectors in slip directions: SLPDIR, and unit normals to \nC     slip planes: SLPNOR in global system\nC\n         DO J=1,NSLIP\n            DO I=1,3\n               TERM(I)=0.\n               DO K=1,3\n                  TERM(I)=TERM(I)+ROTATE(I,K)*SLPDIR(K,J)\n               END DO\n            END DO\n            DO I=1,3\n               SLPDIR(I,J)=TERM(I)\n            END DO\n\n            DO I=1,3\n               TERM(I)=0.\n               DO K=1,3\n                  TERM(I)=TERM(I)+ROTATE(I,K)*SLPNOR(K,J)\n               END DO\n            END DO\n            DO I=1,3\n               SLPNOR(I,J)=TERM(I)\n            END DO\n         END DO\n\n      END IF\n\n      RETURN\n      END\n\n\nC----------------------------------\n\n\n           SUBROUTINE LINE (I1, I2, I3, IARRAY)\n\nC-----  Generating all possible slip directions <lmn> (or slip planes \nC     {lmn}) for a cubic crystal, where l,m,n are not zeros.\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\n           DIMENSION IARRAY(3,4)\n\n           DO J=1,4\n              IARRAY(1,J)=I1\n              IARRAY(2,J)=I2\n              IARRAY(3,J)=I3\n           END DO\n\n           DO I=1,3\n              DO J=1,4\n                 IF (J.EQ.I+1) IARRAY(I,J)=-IARRAY(I,J)\n              END DO\n           END DO\n\n           RETURN\n           END\n\n\nC-----------------------------------\n\n\n           SUBROUTINE LINE1 (J1, J2, IARRAY, ID)\n\nC-----  Generating all possible slip directions <0mn> (or slip planes \nC     {0mn}) for a cubic crystal, where m,n are not zeros and m does \nC     not equal n.\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\n           DIMENSION IARRAY(3,2)\n\n           IARRAY(ID,1)=0\n           IARRAY(ID,2)=0\n\n           ID1=ID+1\n           IF (ID1.GT.3) ID1=ID1-3\n           IARRAY(ID1,1)=J1\n           IARRAY(ID1,2)=J1\n\n           ID2=ID+2\n           IF (ID2.GT.3) ID2=ID2-3\n           IARRAY(ID2,1)=J2\n           IARRAY(ID2,2)=-J2\n  \n           RETURN\n           END\n\n\nC----------------------------------------------------------------------\n\n\n      SUBROUTINE GSLPINIT (GSLIP0, NSLIP, NSLPTL, NSET, PROP)\n\nC-----  This subroutine calculates the initial value of current \nC     strength for each slip system in a rate-dependent single crystal.\nC     Two sets of initial values, proposed by Asaro, Pierce et al, and \nC     by Bassani, respectively, are used here.  Both sets assume that \nC     the initial values for all slip systems are the same (initially \nC     isotropic).\n\nC-----  These initial values are assumed the same for all slip systems \nC     in each set, though they could be different from set to set, e.g.\nC     <110>{111} and <110>{100}.\n\nC-----  Users who want to use their own initial values may change the \nC     function subprogram GSLP0.  The parameters characterizing these \nC     initial values are passed into GSLP0 through array PROP.\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      EXTERNAL GSLP0\n      DIMENSION GSLIP0(NSLPTL), NSLIP(NSET), PROP(16,NSET)\n\nC-----  Function subprograms:\nC\nC       GSLP0 -- User-supplied function subprogram given the initial \nC                value of current strength at initial state\n\nC-----  Variables:\nC\nC     GSLIP0 -- initial value of current strength (OUTPUT)\nC\nC     NSLIP  -- number of slip systems in each set (INPUT)\nC     NSLPTL -- total number of slip systems in all the sets (INPUT)\nC     NSET   -- number of sets of slip systems (INPUT)\nC\nC     PROP   -- material constants characterizing the initial value of \nC               current strength (INPUT)\nC\nC               For Asaro, Pierce et al's law \nC               PROP(1,i) -- initial hardening modulus H0 in the ith \nC                            set of slip systems\nC               PROP(2,i) -- saturation stress TAUs in the ith set of  \nC                            slip systems\nC               PROP(3,i) -- initial critical resolved shear stress \nC                            TAU0 in the ith set of slip systems\nC\nC               For Bassani's law \nC               PROP(1,i) -- initial hardening modulus H0 in the ith \nC                            set of slip systems\nC               PROP(2,i) -- stage I stress TAUI in the ith set of  \nC                            slip systems (or the breakthrough stress \nC                            where large plastic flow initiates)\nC               PROP(3,i) -- initial critical resolved shear stress \nC                            TAU0 in the ith set of slip systems\nC\n\n      ID=0\n      DO I=1,NSET\n         ISET=I\n         DO J=1,NSLIP(I)\n            ID=ID+1\n            GSLIP0(ID)=GSLP0(NSLPTL,NSET,NSLIP,PROP(1,I),ID,ISET)\n         END DO\n      END DO\n\n      RETURN\n      END\n\n\nC----------------------------------\n\n\nC-----  Use single precision on cray\nC\n           REAL*8 FUNCTION GSLP0(NSLPTL,NSET,NSLIP,PROP,ISLIP,ISET)\n\nC-----     User-supplied function subprogram given the initial value of\nC        current strength at initial state\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\n           DIMENSION NSLIP(NSET), PROP(16)\n\n           GSLP0=PROP(3)\n\n           RETURN\n           END\n\n\nC----------------------------------------------------------------------\n\n\n      SUBROUTINE STRAINRATE (GAMMAR, TAUSLP, GSLIP, NSLIP, FSLIP, \n     2                       DFDXSP, PROP)\n\nC-----  This subroutine calculates the shear strain-rate in each slip \nC     system for a rate-dependent single crystal.  The POWER LAW \nC     relation between shear strain-rate and resolved shear stress \nC     proposed by Hutchinson, Pan and Rice, is used here.\n\nC-----  The power law exponents are assumed the same for all slip \nC     systems in each set, though they could be different from set to \nC     set, e.g. <110>{111} and <110>{100}.  The strain-rate coefficient\nC     in front of the power law form are also assumed the same for all \nC     slip systems in each set. \n\nC-----  Users who want to use their own constitutive relation may \nC     change the function subprograms F and its derivative DFDX, \nC     where F is the strain hardening law, dGAMMA\/dt = F(X), \nC     X=TAUSLP\/GSLIP.  The parameters characterizing F are passed into \nC     F and DFDX through array PROP.\n\nC-----  Function subprograms:\nC\nC       F    -- User-supplied function subprogram which gives shear \nC               strain-rate for each slip system based on current \nC               values of resolved shear stress and current strength\nC\nC       DFDX -- User-supplied function subprogram dF\/dX, where x is the\nC               ratio of resolved shear stress over current strength\n\nC-----  Variables:\nC\nC     GAMMAR  -- shear strain in each slip system at the start of time \nC               step  (INPUT)\nC     TAUSLP -- resolved shear stress in each slip system (INPUT)\nC     GSLIP  -- current strength (INPUT)\nC     NSLIP  -- number of slip systems in this set (INPUT)\nC\nC     FSLIP  -- current value of F for each slip system (OUTPUT)\nC     DFDXSP -- current value of DFDX for each slip system (OUTPUT)\nC\nC     PROP   -- material constants characterizing the strain hardening \nC               law (INPUT)\nC\nC               For the current power law strain hardening law \nC               PROP(1) -- power law hardening exponent\nC               PROP(1) = infinity corresponds to a rate-independent \nC               material\nC               PROP(2) -- coefficient in front of power law hardening\n\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      EXTERNAL F, DFDX\n      DIMENSION GAMMAR(NSLIP), TAUSLP(NSLIP), GSLIP(NSLIP), \n     2          FSLIP(NSLIP), DFDXSP(NSLIP), PROP(8)\n\n      DO I=1,NSLIP\n         X=TAUSLP(I)\/GSLIP(I)\n         FSLIP(I)=F(X,PROP)\n         DFDXSP(I)=DFDX(X,PROP)\n      END DO\n\n      RETURN\n      END\n\n\nC-----------------------------------\n\n\nC-----  Use single precision on cray\nC\n           REAL*8 FUNCTION F(X,PROP)\n\nC-----     User-supplied function subprogram which gives shear \nC        strain-rate for each slip system based on current values of \nC        resolved shear stress and current strength\nC\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\n           DIMENSION PROP(8)\n\n           F=PROP(2)*(ABS(X))**PROP(1)*DSIGN(1.D0,X)\n\n           RETURN\n           END\n\n\nC-----------------------------------\n\n\nC-----  Use single precision on cray\nC\n           REAL*8 FUNCTION DFDX(X,PROP)\n\nC-----     User-supplied function subprogram dF\/dX, where x is the \nC        ratio of resolved shear stress over current strength\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\n           DIMENSION PROP(8)\n\n           DFDX=PROP(1)*PROP(2)*(ABS(X))**(PROP(1)-1.)\n\n           RETURN\n           END\n\n\nC----------------------------------------------------------------------\n\nCFIXA\n      SUBROUTINE LATENTHARDEN (GAMMAR, TAUSLP, GSLIP, GMSLTL, GAMTOL, \n     2                         NSLIP, NSLPTL, NSET, H, PROP, ND)\nCFIXB\n\nC-----  This subroutine calculates the current self- and latent-\nC     hardening moduli for all slip systems in a rate-dependent single \nC     crystal.  Two kinds of hardening law are used here.  The first \nC     law, proposed by Asaro, and Pierce et al, assumes a HYPER SECANT \nC     relation between self- and latent-hardening moduli and overall \nC     shear strain.  The Bauschinger effect has been neglected.  The \nC     second is Bassani's hardening law, which gives an explicit \nC     expression of slip interactions between slip systems.  The \nC     classical three stage hardening for FCC single crystal could be \nC     simulated.\n\nC-----  The hardening coefficients are assumed the same for all slip \nC     systems in each set, though they could be different from set to \nC     set, e.g. <110>{111} and <110>{100}.\n\nC-----  Users who want to use their own self- and latent-hardening law \nC     may change the function subprograms HSELF (self hardening) and \nC     HLATNT (latent hardening).  The parameters characterizing these \nC     hardening laws are passed into HSELF and HLATNT through array \nC     PROP.\n\n\nC-----  Function subprograms:\nC\nC       HSELF  -- User-supplied self-hardening function in a slip \nC                 system\nC\nC       HLATNT -- User-supplied latent-hardening function\n\nC-----  Variables:\nC\nC     GAMMAR  -- shear strain in all slip systems at the start of time \nC               step  (INPUT)\nC     TAUSLP -- resolved shear stress in all slip systems (INPUT)\nC     GSLIP  -- current strength (INPUT)\nCFIX  GMSLTL -- total cumulative shear strains on each individual slip system \nCFIX            (INPUT)\nC     GAMTOL -- total cumulative shear strains over all slip systems \nC               (INPUT)\nC     NSLIP  -- number of slip systems in each set (INPUT)\nC     NSLPTL -- total number of slip systems in all the sets (INPUT)\nC     NSET   -- number of sets of slip systems (INPUT)\nC\nC     H      -- current value of self- and latent-hardening moduli \nC               (OUTPUT)\nC               H(i,i) -- self-hardening modulus of the ith slip system\nC                         (no sum over i)\nC               H(i,j) -- latent-hardening molulus of the ith slip \nC                         system due to a slip in the jth slip system \nC                         (i not equal j)\nC\nC     PROP   -- material constants characterizing the self- and latent-\nC               hardening law (INPUT)\nC\nC               For the HYPER SECANT hardening law \nC               PROP(1,i) -- initial hardening modulus H0 in the ith \nC                            set of slip systems\nC               PROP(2,i) -- saturation stress TAUs in the ith set of  \nC                            slip systems\nC               PROP(3,i) -- initial critical resolved shear stress \nC                            TAU0 in the ith set of slip systems\nC               PROP(9,i) -- ratio of latent to self-hardening Q in the\nC                            ith set of slip systems\nC               PROP(10,i)-- ratio of latent-hardening from other sets \nC                            of slip systems to self-hardening in the \nC                            ith set of slip systems Q1\nC\nC               For Bassani's hardening law \nC               PROP(1,i) -- initial hardening modulus H0 in the ith \nC                            set of slip systems\nC               PROP(2,i) -- stage I stress TAUI in the ith set of  \nC                            slip systems (or the breakthrough stress \nC                            where large plastic flow initiates)\nC               PROP(3,i) -- initial critical resolved shear stress \nC                            TAU0 in the ith set of slip systems\nC               PROP(4,i) -- hardening modulus during easy glide Hs in \nC                            the ith set of slip systems\nC               PROP(5,i) -- amount of slip Gamma0 after which a given \nC                            interaction between slip systems in the \nC                            ith set reaches peak strength\nC               PROP(6,i) -- amount of slip Gamma0 after which a given \nC                            interaction between slip systems in the \nC                            ith set and jth set (i not equal j) \nC                            reaches peak strength\nC               PROP(7,i) -- representing the magnitude of the strength\nC                            of interaction in the ith set of slip \nC                            system\nC               PROP(8,i) -- representing the magnitude of the strength\nC                            of interaction between the ith set and jth\nC                            set of system\nC               PROP(9,i) -- ratio of latent to self-hardening Q in the\nC                            ith set of slip systems\nC               PROP(10,i)-- ratio of latent-hardening from other sets \nC                            of slip systems to self-hardening in the \nC                            ith set of slip systems Q1\nC\nC     ND     -- leading dimension of arrays defined in subroutine UMAT \nC               (INPUT) \n\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      EXTERNAL HSELF, HLATNT\nCFIXA\n      DIMENSION GAMMAR(NSLPTL), TAUSLP(NSLPTL), GMSLTL(NSLPTL),\n     2          GSLIP(NSLPTL), NSLIP(NSET), PROP(16,NSET), \n     3          H(ND,NSLPTL)\nCFIXB\n\n      CHECK=0.\n      DO I=1,NSET\n         DO J=4,8\n            CHECK=CHECK+ABS(PROP(J,I))\n         END DO\n      END DO\n\nC-----  CHECK=0   --  HYPER SECANT hardening law\nC       otherwise --  Bassani's hardening law\n\n      ISELF=0\n      DO I=1,NSET\n         ISET=I\n         DO J=1,NSLIP(I)\n            ISELF=ISELF+1\n\n            DO LATENT=1,NSLPTL\n               IF (LATENT.EQ.ISELF) THEN\nCFIXA\n                  H(LATENT,ISELF)=HSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,\n     2                                  NSET,NSLIP,PROP(1,I),CHECK,\n     3                                  ISELF,ISET)\nCFIXB\n               ELSE\nCFIXA\n                  H(LATENT,ISELF)=HLATNT(GAMMAR,GMSLTL,GAMTOL,NSLPTL,\n     2                                   NSET,NSLIP,PROP(1,I),CHECK,\n     3                                   ISELF,ISET,LATENT)\nCFIXB\n\n               END IF\n            END DO\n\n         END DO\n      END DO\n\n      RETURN\n      END\n\n\nC-----------------------------------\n\n\nC-----  Use single precision on cray\nCFIXA\n           REAL*8 FUNCTION HSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,\n     2                           NSLIP,PROP,CHECK,ISELF,ISET)\nCFIXB\n\nC-----     User-supplied self-hardening function in a slip system\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\nCFIXA\n           DIMENSION GAMMAR(NSLPTL), NSLIP(NSET), PROP(16),\n     2               GMSLTL(NSLPTL)\nCFIXB\n\n           IF (CHECK.EQ.0.) THEN\n\nC-----  HYPER SECANT hardening law by Asaro, Pierce et al\n              TERM1=PROP(1)*GAMTOL\/(PROP(2)-PROP(3))\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              HSELF=PROP(1)*TERM2**2\n\n           ELSE\n\nC-----  Bassani's hardening law\nCFIXA\n              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)\/(PROP(2)-PROP(3))\nCFIXB\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)\n\n              ID=0\n              G=1.\n              DO I=1,NSET\n                 IF (I.EQ.ISET) THEN\n                    GAMMA0=PROP(5)\n                    FAB=PROP(7)\n                 ELSE\n                    GAMMA0=PROP(6)\n                    FAB=PROP(8)\n                 END IF\n\n                 DO J=1,NSLIP(I)\n                    ID=ID+1\n                    IF (ID.NE.ISELF) THEN\nCFIXA\n\t\t       G=G+FAB*TANH(GMSLTL(ID)\/GAMMA0)\nCFIXB\n\t\t    END IF\n\n                 END DO\n              END DO\n\n              HSELF=F*G\n\n           END IF\n\n           RETURN\n           END\n\n\nC-----------------------------------\n\n\nC-----  Use single precision on cray\nCFIXA\n           REAL*8 FUNCTION HLATNT(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,\n     2                            NSLIP,PROP,CHECK,ISELF,ISET,LATENT)\nCFIXB\n\nC-----     User-supplied latent-hardening function\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\nCFIXA\n           DIMENSION GAMMAR(NSLPTL), NSLIP(NSET), PROP(16),\n     2               GMSLTL(NSLPTL)\nCFIXB\n\n           ILOWER=0\n           IUPPER=NSLIP(1)\n           IF (ISET.GT.1) THEN\n              DO K=2,ISET\n                 ILOWER=ILOWER+NSLIP(K-1)\n                 IUPPER=IUPPER+NSLIP(K)\n              END DO\n           END IF\n\n           IF (LATENT.GT.ILOWER.AND.LATENT.LE.IUPPER) THEN\n              Q=PROP(9)\n           ELSE\n              Q=PROP(10)\n           END IF\n\n           IF (CHECK.EQ.0.) THEN\n\nC-----  HYPER SECANT hardening law by Asaro, Pierce et al\n              TERM1=PROP(1)*GAMTOL\/(PROP(2)-PROP(3))\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              HLATNT=PROP(1)*TERM2**2*Q\n\n           ELSE\n\nC-----  Bassani's hardening law\nCFIXA\n              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)\/(PROP(2)-PROP(3))\nCFIXB\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)\n\n              ID=0\n              G=1.\n              DO I=1,NSET\n                 IF (I.EQ.ISET) THEN\n                    GAMMA0=PROP(5)\n                    FAB=PROP(7)\n                 ELSE\n                    GAMMA0=PROP(6)\n                    FAB=PROP(8)\n                 END IF\n\n                 DO J=1,NSLIP(I)\n                    ID=ID+1\n                    IF (ID.NE.ISELF) THEN\nCFIXA\n\t\t       G=G+FAB*TANH(GMSLTL(ID)\/GAMMA0)\nCFIXB\n\t\t    END IF\n\n                 END DO\n              END DO\n\n              HLATNT=F*G*Q\n\n           END IF\n\n           RETURN\n           END\n\n\nC----------------------------------------------------------------------\n\nCFIXA\n      SUBROUTINE ITERATION (GAMMAR, TAUSLP, GSLIP, GMSLTL, GAMTOL, \n     2                      NSLPTL, NSET, NSLIP, ND, PROP, DGAMOD, \n     3                      DHDGDG)\nCFIXB\n\nC-----  This subroutine generates arrays for the Newton-Rhapson \nC     iteration method.\n\nC-----  Users who want to use their own self- and latent-hardening law \nC     may change the function subprograms DHSELF (self hardening) and \nC     DHLATN (latent hardening).  The parameters characterizing these \nC     hardening laws are passed into DHSELF and DHLATN through array \nC     PROP.\n\n\nC-----  Function subprograms:\nC\nC       DHSELF -- User-supplied function of the derivative of self-\nC                 hardening moduli\nC\nC       DHLATN -- User-supplied function of the derivative of latent-\nC                 hardening moduli\n\nC-----  Variables:\nC\nC     GAMMAR  -- shear strain in all slip systems at the start of time \nC               step  (INPUT)\nC     TAUSLP -- resolved shear stress in all slip systems (INPUT)\nC     GSLIP  -- current strength (INPUT)\nCFIX  GMSLTL -- total cumulative shear strains on each individual slip system \nCFIX            (INPUT)\nC     GAMTOL -- total cumulative shear strains over all slip systems \nC               (INPUT)\nC     NSLPTL -- total number of slip systems in all the sets (INPUT)\nC     NSET   -- number of sets of slip systems (INPUT)\nC     NSLIP  -- number of slip systems in each set (INPUT)\nC     ND     -- leading dimension of arrays defined in subroutine UMAT \nC               (INPUT) \nC\nC     PROP   -- material constants characterizing the self- and latent-\nC               hardening law (INPUT)\nC\nC               For the HYPER SECANT hardening law \nC               PROP(1,i) -- initial hardening modulus H0 in the ith \nC                            set of slip systems\nC               PROP(2,i) -- saturation stress TAUs in the ith set of  \nC                            slip systems\nC               PROP(3,i) -- initial critical resolved shear stress \nC                            TAU0 in the ith set of slip systems\nC               PROP(9,i) -- ratio of latent to self-hardening Q in the\nC                            ith set of slip systems\nC               PROP(10,i)-- ratio of latent-hardening from other sets \nC                            of slip systems to self-hardening in the \nC                            ith set of slip systems Q1\nC\nC               For Bassani's hardening law \nC               PROP(1,i) -- initial hardening modulus H0 in the ith \nC                            set of slip systems\nC               PROP(2,i) -- stage I stress TAUI in the ith set of  \nC                            slip systems (or the breakthrough stress \nC                            where large plastic flow initiates)\nC               PROP(3,i) -- initial critical resolved shear stress \nC                            TAU0 in the ith set of slip systems\nC               PROP(4,i) -- hardening modulus during easy glide Hs in \nC                            the ith set of slip systems\nC               PROP(5,i) -- amount of slip Gamma0 after which a given \nC                            interaction between slip systems in the \nC                            ith set reaches peak strength\nC               PROP(6,i) -- amount of slip Gamma0 after which a given \nC                            interaction between slip systems in the \nC                            ith set and jth set (i not equal j) \nC                            reaches peak strength\nC               PROP(7,i) -- representing the magnitude of the strength\nC                            of interaction in the ith set of slip \nC                            system\nC               PROP(8,i) -- representing the magnitude of the strength\nC                            of interaction between the ith set and jth\nC                            set of system\nC               PROP(9,i) -- ratio of latent to self-hardening Q in the\nC                            ith set of slip systems\nC               PROP(10,i)-- ratio of latent-hardening from other sets \nC                            of slip systems to self-hardening in the \nC                            ith set of slip systems Q1\nC\nC-----  Arrays for iteration:\nC\nC       DGAMOD (INPUT)\nC\nC       DHDGDG (OUTPUT)\nC\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      EXTERNAL DHSELF, DHLATN\nCFIXA\n      DIMENSION GAMMAR(NSLPTL), TAUSLP(NSLPTL), GMSLTL(NSLPTL),\n     2          GSLIP(NSLPTL), NSLIP(NSET), PROP(16,NSET), \n     3          DGAMOD(NSLPTL), DHDGDG(ND,NSLPTL)\nCFIXB\n\n      CHECK=0.\n      DO I=1,NSET\n         DO J=4,8\n            CHECK=CHECK+ABS(PROP(J,I))\n         END DO\n      END DO\n\nC-----  CHECK=0   --  HYPER SECANT hardening law\nC       otherwise --  Bassani's hardening law\n\n      ISELF=0\n      DO I=1,NSET\n         ISET=I\n         DO J=1,NSLIP(I)\n            ISELF=ISELF+1\n\n            DO KDERIV=1,NSLPTL\n               DHDGDG(ISELF,KDERIV)=0.\n\n               DO LATENT=1,NSLPTL\n                  IF (LATENT.EQ.ISELF) THEN\nCFIXA\n                     DHDG=DHSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,\n     2                           NSLIP,PROP(1,I),CHECK,ISELF,ISET,\n     3                           KDERIV)\nCFIXB\n                  ELSE\nCFIXA\n                     DHDG=DHLATN(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,\n     2                           NSLIP,PROP(1,I),CHECK,ISELF,ISET,\n     3                           LATENT,KDERIV)\nCFIXB\n                  END IF\n\n                  DHDGDG(ISELF,KDERIV)=DHDGDG(ISELF,KDERIV)+\n     2                                 DHDG*ABS(DGAMOD(LATENT))\n               END DO\n\n            END DO\n         END DO\n      END DO\n\n      RETURN\n      END\n\n\nC-----------------------------------\n\n\nC-----  Use single precision on cray\nCFIXA\n           REAL*8 FUNCTION DHSELF(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,\n     2                            NSLIP,PROP,CHECK,ISELF,ISET,\n     3                            KDERIV)\nCFIXB\n\nC-----  User-supplied function of the derivative of self-hardening\nC     moduli\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\nCFIXA\n           DIMENSION GAMMAR(NSLPTL), GMSLTL(NSLPTL), \n     2               NSLIP(NSET), PROP(16)\nCFIXB\n\n           IF (CHECK.EQ.0.) THEN\n\nC-----  HYPER SECANT hardening law by Asaro, Pierce et al\n              TERM1=PROP(1)*GAMTOL\/(PROP(2)-PROP(3))\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              TERM3=PROP(1)\/(PROP(2)-PROP(3))*DSIGN(1.D0,GAMMAR(KDERIV))\n              DHSELF=-2.*PROP(1)*TERM2**2*TANH(TERM1)*TERM3\n\n           ELSE\n\nC-----  Bassani's hardening law\nCFIXA\n              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)\/(PROP(2)-PROP(3))\nCFIXB\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              TERM3=(PROP(1)-PROP(4))\/(PROP(2)-PROP(3))\n\n              IF (KDERIV.EQ.ISELF) THEN\n                 F=-2.*(PROP(1)-PROP(4))*TERM2**2*TANH(TERM1)*TERM3\n                 ID=0\n                 G=1.\n                 DO I=1,NSET\n                    IF (I.EQ.ISET) THEN\n                       GAMMA0=PROP(5)\n                       FAB=PROP(7)\n                    ELSE\n                       GAMMA0=PROP(6)\n                       FAB=PROP(8)\n                    END IF\n\n                    DO J=1,NSLIP(I)\n                       ID=ID+1\nCFIXA\n                       IF (ID.NE.ISELF) G=G+FAB*TANH(GMSLTL(ID)\/GAMMA0)\nCFIXB\n                    END DO\n                 END DO\n\n              ELSE\n                 F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)\n                 ILOWER=0\n                 IUPPER=NSLIP(1)\n                 IF (ISET.GT.1) THEN\n                    DO K=2,ISET\n                       ILOWER=ILOWER+NSLIP(K-1)\n                       IUPPER=IUPPER+NSLIP(K)\n                    END DO\n                 END IF\n\n                 IF (KDERIV.GT.ILOWER.AND.KDERIV.LE.IUPPER) THEN\n                    GAMMA0=PROP(5)\n                    FAB=PROP(7)\n                 ELSE\n                    GAMMA0=PROP(6)\n                    FAB=PROP(8)\n                 END IF\n\nCFIXA\n                 TERM4=GMSLTL(KDERIV)\/GAMMA0\nCFIXB\n                 TERM5=2.*EXP(-TERM4)\/(1.+EXP(-2.*TERM4))\n                 G=FAB\/GAMMA0*TERM5**2\n\n              END IF\n\n              DHSELF=F*G\n\n           END IF\n\n           RETURN\n           END\n\n\nC-----------------------------------\n\n\nC-----  Use single precision on cray\nCFIXA\n           REAL*8 FUNCTION DHLATN(GAMMAR,GMSLTL,GAMTOL,NSLPTL,NSET,\n     2                            NSLIP,PROP,CHECK,ISELF,ISET,LATENT,\n     3                            KDERIV)\nCFIXB\n\nC-----  User-supplied function of the derivative of latent-hardening \nC     moduli\n\nC-----  Use single precision on cray\nC\n           IMPLICIT REAL*8 (A-H,O-Z)\nCFIXA\n           DIMENSION GAMMAR(NSLPTL), GMSLTL(NSLPTL), NSLIP(NSET), \n     2               PROP(16)\nCFIXB\n\n           ILOWER=0\n           IUPPER=NSLIP(1)\n           IF (ISET.GT.1) THEN\n              DO K=2,ISET\n                 ILOWER=ILOWER+NSLIP(K-1)\n                 IUPPER=IUPPER+NSLIP(K)\n              END DO\n           END IF\n\n           IF (LATENT.GT.ILOWER.AND.LATENT.LE.IUPPER) THEN\n              Q=PROP(9)\n           ELSE\n              Q=PROP(10)\n           END IF\n\n           IF (CHECK.EQ.0.) THEN\n\nC-----  HYPER SECANT hardening law by Asaro, Pierce et al\n              TERM1=PROP(1)*GAMTOL\/(PROP(2)-PROP(3))\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              TERM3=PROP(1)\/(PROP(2)-PROP(3))*DSIGN(1.D0,GAMMAR(KDERIV))\n              DHLATN=-2.*PROP(1)*TERM2**2*TANH(TERM1)*TERM3*Q\n\n           ELSE\n\nC-----  Bassani's hardening law\nCFIXA\n              TERM1=(PROP(1)-PROP(4))*GMSLTL(ISELF)\/(PROP(2)-PROP(3))\nCFIXB\n              TERM2=2.*EXP(-TERM1)\/(1.+EXP(-2.*TERM1))\n              TERM3=(PROP(1)-PROP(4))\/(PROP(2)-PROP(3))\n\n              IF (KDERIV.EQ.ISELF) THEN\n                 F=-2.*(PROP(1)-PROP(4))*TERM2**2*TANH(TERM1)*TERM3\n                 ID=0\n                 G=1.\n                 DO I=1,NSET\n                    IF (I.EQ.ISET) THEN\n                       GAMMA0=PROP(5)\n                       FAB=PROP(7)\n                    ELSE\n                       GAMMA0=PROP(6)\n                       FAB=PROP(8)\n                    END IF\n\n                    DO J=1,NSLIP(I)\n                       ID=ID+1\nCFIXA\n                       IF (ID.NE.ISELF) G=G+FAB*TANH(GMSLTL(ID)\/GAMMA0)\nCFIXB\n                    END DO\n                 END DO\n\n              ELSE\n                 F=(PROP(1)-PROP(4))*TERM2**2+PROP(4)\n                 ILOWER=0\n                 IUPPER=NSLIP(1)\n                 IF (ISET.GT.1) THEN\n                    DO K=2,ISET\n                       ILOWER=ILOWER+NSLIP(K-1)\n                       IUPPER=IUPPER+NSLIP(K)\n                    END DO\n                 END IF\n\n                 IF (KDERIV.GT.ILOWER.AND.KDERIV.LE.IUPPER) THEN\n                    GAMMA0=PROP(5)\n                    FAB=PROP(7)\n                 ELSE\n                    GAMMA0=PROP(6)\n                    FAB=PROP(8)\n                 END IF\nCFIXA\n                 TERM4=GMSLTL(KDERIV)\/GAMMA0\nCFIXB\n                 TERM5=2.*EXP(-TERM4)\/(1.+EXP(-2.*TERM4))\n                 G=FAB\/GAMMA0*TERM5**2\n\n              END IF\n\n              DHLATN=F*G*Q\n\n           END IF\n\n           RETURN\n           END\n\n\nC----------------------------------------------------------------------\n\n\n      SUBROUTINE LUDCMP (A, N, NP, INDX, D)\n\nC-----  LU decomposition\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      PARAMETER (NMAX=200, TINY=1.0E-20)\n      DIMENSION A(NP,NP), INDX(N), VV(NMAX)\n\n      D=1.\n      DO I=1,N\n         AAMAX=0.\n\n         DO J=1,N\n            IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))\n         END DO\n\n         IF (AAMAX.EQ.0.) PAUSE 'Singular matrix.'\n         VV(I)=1.\/AAMAX\n      END DO\n\n      DO J=1,N\n         DO I=1,J-1\n            SUM=A(I,J)\n\n            DO K=1,I-1\n               SUM=SUM-A(I,K)*A(K,J)\n            END DO\n\n            A(I,J)=SUM\n         END DO\n         AAMAX=0.\n\n         DO I=J,N\n            SUM=A(I,J)\n\n            DO K=1,J-1\n               SUM=SUM-A(I,K)*A(K,J)\n            END DO\n\n            A(I,J)=SUM\n            DUM=VV(I)*ABS(SUM)\n            IF (DUM.GE.AAMAX) THEN\n               IMAX=I\n               AAMAX=DUM\n            END IF\n         END DO\n\n         IF (J.NE.IMAX) THEN\n            DO K=1,N\n               DUM=A(IMAX,K)\n               A(IMAX,K)=A(J,K)\n               A(J,K)=DUM\n            END DO\n\n            D=-D\n            VV(IMAX)=VV(J)\n         END IF\n\n         INDX(J)=IMAX\n         IF (A(J,J).EQ.0.) A(J,J)=TINY\n         IF (J.NE.N) THEN\n            DUM=1.\/A(J,J)\n            DO I=J+1,N\n               A(I,J)=A(I,J)*DUM\n            END DO\n         END IF\n\n      END DO\n\n      RETURN\n      END\n\n\nC----------------------------------------------------------------------\n\n\n      SUBROUTINE LUBKSB (A, N, NP, INDX, B)\n\nC-----  Linear equation solver based on LU decomposition\n\nC-----  Use single precision on cray\nC\n      IMPLICIT REAL*8 (A-H,O-Z)\n      DIMENSION A(NP,NP), INDX(N), B(N)\n\n      II=0\n      DO I=1,N\n         LL=INDX(I)\n         SUM=B(LL)\n         B(LL)=B(I)\n\n         IF (II.NE.0) THEN\n            DO J=II,I-1\n               SUM=SUM-A(I,J)*B(J)\n            END DO\n         ELSE IF (SUM.NE.0.) THEN\n            II=I\n         END IF\n\n         B(I)=SUM\n      END DO\n\n      DO I=N,1,-1\n         SUM=B(I)\n\n         IF (I.LT.N) THEN\n            DO J=I+1,N\n               SUM=SUM-A(I,J)*B(J)\n            END DO\n         END IF\n\n         B(I)=SUM\/A(I,I)\n      END DO\n\n      RETURN\n      END\n      \nC     ********************* 25. GET_INV_DET   *********************\n\nc   this subroutine inverts a given matrix and return it  \n      \n      SUBROUTINE GET_INV_DET(XS,DET,NFLAG)\n                                      \n      IMPLICIT REAL*8 (A-H,O-Z)                                         \n                 \n      DIMENSION XS(3,3),A(3,3) \n      \n      DO 10 I=1,3                                                       \n      I1=I+1                                                            \n      I2=I+2                                                            \n      IF(I1.GT.3) I1=I1-3                                               \n      IF(I2.GT.3) I2=I2-3                                               \n      DO 10 J=1,3                                                       \n          J1=J+1                                                  \n          J2=J+2                                                  \n          IF(J1.GT.3) J1=J1-3                                     \n          IF(J2.GT.3) J2=J2-3                                     \n   10         A(I,J)=XS(I1,J1)*XS(I2,J2)-XS(I1,J2)*XS(I2,J1)\n                      \n      DET=0.D0 \n                                                             \n      DO 20 I=1,3                                                       \n   20     DET=DET+A(1,I)*XS(1,I) \n   \n      NE=1  \n\n      IF (NFLAG.EQ.0.AND.DET.LE.0.0D0) THEN\n                DET=0.00001\n      !WRITE(*,*) 'ERROR IN THE JACOBAIN, DETERMINANT', DET\nC      PAUSE'ERROR IN THE JACOBIAN' \n      ENDIF         \n                                           \n      DO 30  I=1,3                                                      \n      DO 30 J=1,3                                                       \n   30         XS(I,J)=A(J,I)\/DET \n                                                  \n      RETURN                                                            \n      END   \nC     ********************* 29. CLEAR    *********************\n\nc     this subroutine initializes a real matrix to zero\n\n      SUBROUTINE CLEAR(A,N)\n      \n      IMPLICIT REAL*8 (A-H,O-Z)\n      \n      DIMENSION A(N)\n      \n      DO 10 I=1,N\n 10       A(I)=0.D0 \n \n      RETURN\n      END \n<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[7,8],"tags":[],"class_list":["post-253","post","type-post","status-publish","format-standard","hentry","category-abaqus","category-abaqus-subroutine"],"blocksy_meta":[],"_links":{"self":[{"href":"https:\/\/numsimlab.com\/index.php?rest_route=\/wp\/v2\/posts\/253","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/numsimlab.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/numsimlab.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/numsimlab.com\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/numsimlab.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=253"}],"version-history":[{"count":2,"href":"https:\/\/numsimlab.com\/index.php?rest_route=\/wp\/v2\/posts\/253\/revisions"}],"predecessor-version":[{"id":257,"href":"https:\/\/numsimlab.com\/index.php?rest_route=\/wp\/v2\/posts\/253\/revisions\/257"}],"wp:attachment":[{"href":"https:\/\/numsimlab.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=253"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/numsimlab.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=253"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/numsimlab.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=253"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}