mat104c_nodam_nice.F

Go to the documentation of this file.

Copyright>        OpenRadioss
Copyright>        Copyright (C) 1986-2025 Altair Engineering Inc.
Copyright>
Copyright>        This program is free software: you can redistribute it and/or modify
Copyright>        it under the terms of the GNU Affero General Public License as published by
Copyright>        the Free Software Foundation, either version 3 of the License, or
Copyright>        (at your option) any later version.
Copyright>
Copyright>        This program is distributed in the hope that it will be useful,
Copyright>        but WITHOUT ANY WARRANTY; without even the implied warranty of
Copyright>        MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
Copyright>        GNU Affero General Public License for more details.
Copyright>
Copyright>        You should have received a copy of the GNU Affero General Public License
Copyright>        along with this program.  If not, see <https://www.gnu.org/licenses/>.
Copyright>
Copyright>
Copyright>        Commercial Alternative: Altair Radioss Software
Copyright>
Copyright>        As an alternative to this open-source version, Altair also offers Altair Radioss
Copyright>        software under a commercial license.  Contact Altair to discuss further if the
Copyright>        commercial version may interest you: https://www.altair.com/radioss/.
!||====================================================================
!||    mat104c_nodam_nice   ../engine/source/materials/mat/mat104/mat104c_nodam_nice.F
!||--- called by ------------------------------------------------------
!||    sigeps104c           ../engine/source/materials/mat/mat104/sigeps104c.f
!||====================================================================
      SUBROUTINE mat104c_nodam_nice(
     1     NEL     ,NGL     ,NUPARAM ,NUVAR   , 
     2     TIME    ,TIMESTEP,UPARAM  ,UVAR    ,JTHE    ,OFF     ,
     3     GS      ,RHO     ,PLA     ,DPLA    ,EPSD    ,SOUNDSP ,
     4     DEPSXX  ,DEPSYY  ,DEPSXY  ,DEPSYZ  ,DEPSZX  ,
     5     SIGOXX  ,SIGOYY  ,SIGOXY  ,SIGOYZ  ,SIGOZX  ,
     6     SIGNXX  ,SIGNYY  ,SIGNXY  ,SIGNYZ  ,SIGNZX  ,THKLY   ,
     7     THK     ,SIGY    ,ET      ,TEMPEL  ,DPLANL  ,TEMP    ,
     8     SEQ     ,INLOC   )
      !=======================================================================
      !      Implicit types
      !=======================================================================
#include      "implicit_f.inc"
      !=======================================================================
      !      Common
      !=======================================================================
      !=======================================================================
      !      Dummy arguments
      !=======================================================================
      INTEGER NEL,NUPARAM,NUVAR,JTHE,INLOC
      INTEGER ,DIMENSION(NEL), INTENT(IN)    :: NGL
      my_real 
     .   TIME,TIMESTEP
      my_real,DIMENSION(NUPARAM), INTENT(IN) :: 
     .   UPARAM
      my_real,DIMENSION(NEL), INTENT(IN)     :: 
     .   RHO,TEMPEL,DPLANL,
     .   depsxx,depsyy,depsxy,depsyz,depszx,
     .   sigoxx,sigoyy,sigoxy,sigoyz,sigozx,
     .   gs , thkly
c
      my_real ,DIMENSION(NEL), INTENT(OUT)   :: 
     .   soundsp,sigy,et,
     .   signxx,signyy,signxy,signyz,signzx
c
      my_real ,DIMENSION(NEL), INTENT(INOUT)       :: 
     .   pla,dpla,epsd,off,thk,temp,seq
      my_real ,DIMENSION(NEL,NUVAR), INTENT(INOUT) :: 
     .   uvar
c      
      !=======================================================================
      !      Local Variables
      !=======================================================================
      INTEGER I,K,II,NINDX,INDEX(NEL),NICE
c
      my_real 
     .   young,bulk,lam,g,g2,nu,cdr,kdr,hard,yld0,qvoce,bvoce,jcc,
     .   epsp0,mtemp,tini,tref,eta,cp,dpis,dpad,asrate,afiltr,a11,a12,
     .   nnu, normsig
      my_real 
     .   h,ldav,trdfds,sigvm,sigdr2,yld2i,omega,
     .   dtemp,fcosh,fsinh,dpdt,dlam,ddep,depxx,depyy
      my_real 
     .   dpxx,dpyy,dpxy,dpzz,dsdrdj2,dsdrdj3,
     .   dj3dsxx,dj3dsyy,dj3dsxy,dj3dszz,
     .   dj2dsxx,dj2dsyy,dj2dsxy,
     .   dfdsxx,dfdsyy,dfdsxy,
     .   normxx,normyy,normxy,normzz,
     .   denom,dphi,dphi_dtrsig,sig_dfdsig,dfdsig2,
     .   dphi_dsig,dphi_dyld,dphi_dfdr,dpdt_nl,
     .   dyld_dpla,dyld_dtemp,dtemp_dlam,dlam_nl
c
      my_real, DIMENSION(NEL) ::
     .   dsigxx,dsigyy,dsigxy,trsig,trdep,
     .   sxx,syy,sxy,szz,sigm,j2,j3,sigdr,yld,weitemp,
     .   hardp,fhard,frate,ftherm,fdr,phi0,phi,dpla_dlam,
     .   dezz
      !=======================================================================
      !             DRUCKER - VOCE - JOHNSON-COOK MATERIAL LAW
      !=======================================================================
      !UVAR(1)   PHI    YIELD FUNCTION VALUE
      !=======================================================================
c      
      !=======================================================================
      ! - INITIALISATION OF COMPUTATION ON TIME STEP
      !=======================================================================
      ! Recovering model parameters
      ! Elastic parameters      
      young   = uparam(1)  ! Young modulus
      bulk    = uparam(2)  ! Bulk modulus 
      g       = uparam(3)  ! Shear modulus 
      g2      = uparam(4)  ! 2*Shear modulus 
      lam     = uparam(5)  ! Lambda Hooke parameter 
      nu      = uparam(6)  ! Poisson ration 
      nnu     = uparam(7)  ! NU/(1-NU)
      a11     = uparam(9)  ! Diagonal term, elastic matrix in plane stress
      a12     = uparam(10) ! Non-diagonal term, elastic matrix in plane stress   
      ! Plastic criterion and hardening parameters [Drucker, 1948]
      nice    = nint(uparam(11)) ! Choice of the Nice method
      cdr     = uparam(12) ! Drucker coefficient
      kdr     = uparam(13) ! Drucker 1/K coefficient
      tini    = uparam(14) ! Initial temperature
      hard    = uparam(15) ! Linear hardening
      yld0    = uparam(16) ! Initial yield stress
      qvoce   = uparam(17) ! 1st Voce parameter
      bvoce   = uparam(18) ! 2nd Voce parameter
      ! Strain-rate dependence parameters
      jcc     = uparam(20) ! Johnson-Cook strain rate coefficient
      epsp0   = uparam(21) ! Johnson-Cook reference strain rate
      ! Thermal softening and self-heating parameters
      mtemp   = uparam(22) ! Thermal softening slope
      tref    = uparam(23) ! Reference temperature
      eta     = uparam(24) ! Taylor-Quinney coefficient
      cp      = uparam(25) ! Thermal mass capacity
      dpis    = uparam(26) ! Isothermal plastic strain rate
      dpad    = uparam(27) ! Adiabatic plastic strain rate
      ! Plastic strain-rate filtering parameters
      asrate  = uparam(28) ! Plastic strain rate filtering frequency
      afiltr  = asrate*timestep/(asrate*timestep + one)
c      
      ! Recovering internal variables
      DO i=1,nel
        IF (off(i) < em01) off(i) = zero
        IF (off(i) < one)  off(i) = off(i)*four_over_5
        ! User inputs
        phi0(i)  = uvar(i,1) ! Previous yield function value
        ! Standard inputs
        dpla(i)  = zero      ! Initialization of the plastic strain increment
        et(i)    = one       ! Initialization of hourglass coefficient
        hardp(i) = zero      ! Initialization of hardening modulus
      ENDDO
c      
      ! Initialization of temperature and self-heating weight factor
      IF (time == zero) THEN   
        temp(1:nel) = tini 
      ENDIF 
      IF (cp > zero) THEN     
        IF (jthe == 0) THEN     
          IF (inloc == 0) THEN
            DO i=1,nel
              ! Computation of the weight factor
              IF (epsd(i) < dpis) THEN
                weitemp(i) = zero
              ELSEIF (epsd(i) > dpad) THEN
                weitemp(i) = one
              ELSE
                weitemp(i) = ((epsd(i)-dpis)**2 )
     .                  * (three*dpad - two*epsd(i) - dpis)
     .                  / ((dpad-dpis)**3)
              ENDIF
            ENDDO
          ENDIF
        ELSE
          temp(1:nel) = tempel(1:nel)
        ENDIF
      ENDIF
c      
      ! Computation of the initial yield stress
      DO i = 1,nel
        ! a) - Hardening law
        fhard(i) = yld0 + hard*pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
        ! b) - Correction factor for strain-rate dependence (Johnson-Cook)
        frate(i)  = one
        IF (epsd(i) > epsp0) frate(i) = frate(i) + jcc*log(epsd(i)/epsp0)   
        ! c) - Correction factor for thermal softening      
        ftherm(i) = one
        IF (cp > zero) ftherm(i) = one - mtemp * (temp(i) - tref)
        ! d) - Computation of the yield function and value check
        yld(i) = fhard(i)*frate(i)*ftherm(i)
        ! e) - Checking values
        yld(i) = max(em10, yld(i))
      ENDDO
c
      !========================================================================
      ! - COMPUTATION OF TRIAL VALUES
      !========================================================================
      DO i=1,nel
c
        ! Computation of the trial stress tensor
        signxx(i) = sigoxx(i) + a11*depsxx(i) + a12*depsyy(i)
        signyy(i) = sigoyy(i) + a11*depsyy(i) + a12*depsxx(i)
        signxy(i) = sigoxy(i) + depsxy(i)*g
        signyz(i) = sigoyz(i) + depsyz(i)*gs(i)
        signzx(i) = sigozx(i) + depszx(i)*gs(i)
        ! Computation of the trace of the trial stress tensor
        trsig(i)  = signxx(i) + signyy(i) 
        sigm(i)   = -trsig(i) * third
        ! Computation of the deviatoric trial stress tensor
        sxx(i)    = signxx(i) + sigm(i)
        syy(i)    = signyy(i) + sigm(i)
        szz(i)    = sigm(i)
        sxy(i)    = signxy(i)
        dezz(i)   = -nnu*depsxx(i) - nnu*depsyy(i)
        ! Second deviatoric invariant
        j2(i) = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2 
        ! Third deviatoric invariant
        j3(i)  = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2 
        ! Drucker equivalent stress
        fdr(i) = j2(i)**3 - cdr*(j3(i)**2)
        ! Checking equivalent stress values
        IF (fdr(i) > zero) THEN
          sigdr(i) = kdr * exp((one/six)*log(fdr(i)))  ! FDR(I)**(1/6)
        ELSE
          sigdr(i) = zero
        ENDIF
      ENDDO
c
      !========================================================================
      ! - COMPUTATION OF YIELD FONCTION
      !========================================================================
      phi(1:nel) = (sigdr(1:nel) / yld(1:nel))**2 - one
c
      ! Checking plastic behavior for all elements
      nindx = 0
      DO i=1,nel         
        IF (phi(i) >= zero .AND. off(i) == one) THEN
          nindx=nindx+1
          index(nindx)=i
        ENDIF
      ENDDO
c      
      !====================================================================
      ! - PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
      !====================================================================
#include "vectorize.inc" 
      DO ii=1,nindx   
c      
        ! Number of the element with plastic behaviour                                               
        i = index(ii)
c        
        ! Computation of the trial stress increment
        dsigxx(i) = a11*depsxx(i) + a12*depsyy(i)
        dsigyy(i) = a11*depsyy(i) + a12*depsxx(i)
        dsigxy(i) = depsxy(i)*g  
        dlam      = zero 
c        
        ! Computation of Drucker equivalent stress of the previous stress tensor
        trsig(i)  = sigoxx(i) + sigoyy(i)
        sigm(i)   = -trsig(i) * third
        sxx(i)    = sigoxx(i) + sigm(i)
        syy(i)    = sigoyy(i) + sigm(i)
        szz(i)    = sigm(i)
        sxy(i)    = sigoxy(i)
        j2(i)     = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2) + sxy(i)**2 
        j3(i)     = sxx(i)*syy(i)*szz(i) - szz(i)*sxy(i)**2          
        fdr(i)    = j2(i)**3 - cdr*(j3(i)**2)
        sigdr(i)  = kdr * exp((one/six)*log(fdr(i))) 
c        
        ! Note     : in this part, the purpose is to compute in one iteration
        ! a plastic multiplier allowing to update internal variables to satisfy
        ! the consistency condition. 
        ! Its expression is : DLAMBDA = - (PHI_OLD + DPHI) / DPHI_DLAMBDA
        ! -> PHI_OLD : old value of yield function (known)
        ! -> DPHI : yield function prediction (to compute)
        ! -> DPHI_DLAMBDA : derivative of PHI with respect to DLAMBDA by taking
        !                   into account of internal variables kinetic : 
        !                   plasticity, temperature and damage (to compute)
c        
        ! 1 - Computation of DPHI_DSIG the normal to the yield surface
        !-------------------------------------------------------------  
c        
        ! Derivative with respect to the equivalent stress
        yld2i       = one/(yld(i)**2)
        dphi_dsig   = two*sigdr(i)*yld2i
c        
        ! Computation of the Eulerian norm of the stress tensor
        normsig = sqrt(sigoxx(i)*sigoxx(i)
     .               + sigoyy(i)*sigoyy(i)
     .          + two*sigoxy(i)*sigoxy(i))
        normsig = max(normsig,one)
c     
        ! Derivative with respect to Fdr 
        fdr(i)     =  (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)       
        dphi_dfdr  =  dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))             
        dsdrdj2    =  dphi_dfdr*three*(j2(i)/(normsig**2))**2                             
        dsdrdj3    = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))                   
        ! dJ3/dSig                                                        
        dj3dsxx =  two_third*(syy(i)*szz(i))/(normsig**2)
     .              -  third*(sxx(i)*szz(i))/(normsig**2)
     .              -  third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)                 
        dj3dsyy =    - third*(syy(i)*szz(i))/(normsig**2)
     .           + two_third*(sxx(i)*szz(i))/(normsig**2)
     .              -  third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)             
        dj3dszz =    - third*(syy(i)*szz(i))/(normsig**2)
     .               - third*(sxx(i)*szz(i))/(normsig**2)
     .           + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2) 
        dj3dsxy =  two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)               
        ! dPhi/dSig
        normxx  = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
        normyy  = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
        normzz  = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
        normxy  = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
c        
        ! Restoring previous value of the yield function
        phi(i) = phi0(i)
c
        ! Computation of yield surface trial increment DPHI       
        dphi = normxx * dsigxx(i)  
     .       + normyy * dsigyy(i)   
     .       + normxy * dsigxy(i)      
c             
        ! 2 - Computation of DPHI_DLAMBDA
        !---------------------------------------------------------
c        
        !   a) Derivative with respect stress increments tensor DSIG
        !   --------------------------------------------------------
        dfdsig2 = normxx * (a11*normxx + a12*normyy)
     .          + normyy * (a11*normyy + a12*normxx)
     .          + normxy * normxy * g
c     
        !   b) Derivatives with respect to plastic strain P 
        !   ------------------------------------------------
c        
        !     i) Derivative of the yield stress with respect to plastic strain dYLD / dPLA
        !     ----------------------------------------------------------------------------
        hardp(i)   = hard + qvoce*bvoce*exp(-bvoce*pla(i))      
        dyld_dpla  = hardp(i)*frate(i)*ftherm(i)       
c              
        !     ii) Derivative of dPLA with respect to DLAM
        !     -------------------------------------------
        sig_dfdsig = sigoxx(i) * normxx
     .             + sigoyy(i) * normyy
     .             + sigoxy(i) * normxy        
        dpla_dlam(i) = sig_dfdsig / yld(i)  
c        
        !  c) Derivatives with respect to the temperature TEMP 
        !  ---------------------------------------------------
        IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
        !    i)  Derivative of the yield stress with respect to temperature dYLD / dTEMP
        !    ---------------------------------------------------------------------------          
          dyld_dtemp = -fhard(i)*frate(i)*mtemp
        !    ii) Derivative of the temperature TEMP with respect to DLAM
        !    -----------------------------------------------------------
          dtemp_dlam = weitemp(i)*(eta/(rho(i)*cp))*sig_dfdsig
        ELSE
          dyld_dtemp = zero
          dtemp_dlam = zero
        ENDIF
c        
        !  d) derivative with respect to the yield stress
        !  ----------------------------------------------
        sigdr2    =  sigdr(i)**2                        
        dphi_dyld = -two*sigdr2/(yld(i)**3)   
c        
        !  e) Derivative of PHI with respect to DLAM ( = -DENOM)
        denom = dfdsig2 - (dphi_dyld * dyld_dpla * dpla_dlam(i))
        IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN
          denom = denom - (dphi_dyld * dyld_dtemp * dtemp_dlam)
        ENDIF
        denom = sign(max(abs(denom),em20) ,denom)   
c        
        ! 3 - Computation of plastic multiplier and variables update
        !----------------------------------------------------------
c                              
        ! Computation of the plastic multiplier
        dlam = (dphi + phi(i)) / denom
        dlam = max(dlam, zero)
c        
        ! Plastic strains tensor update
        dpxx = dlam * normxx
        dpyy = dlam * normyy
        dpzz = dlam * normzz
        dpxy = dlam * normxy
c        
        ! Elasto-plastic stresses update   
        signxx(i) = sigoxx(i) + dsigxx(i) - (a11*dpxx + a12*dpyy)
        signyy(i) = sigoyy(i) + dsigyy(i) - (a11*dpyy + a12*dpxx)      
        signxy(i) = sigoxy(i) + dsigxy(i) - dpxy*g
        trsig(i)  = signxx(i) + signyy(i)
        sigm(i)   = -trsig(i) * third
        sxx(i)    = signxx(i) + sigm(i)
        syy(i)    = signyy(i) + sigm(i)
        szz(i)    = sigm(i)
        sxy(i)    = signxy(i)
c  
        ! Cumulated plastic strain and strain rate update
        ddep    = (dlam/yld(i))*sig_dfdsig
        dpla(i) = dpla(i) + ddep 
        pla(i)  = pla(i)  + ddep  
c        
        ! Drucker equivalent stress update
        j2(i)    = half*(sxx(i)**2 + syy(i)**2 + szz(i)**2 ) + sxy(i)**2 
        j3(i)    = sxx(i) * syy(i) * szz(i) - szz(i) * sxy(i)**2 
        fdr(i)   = j2(i)**3 - cdr*(j3(i)**2)
        sigdr(i) = kdr * exp((one/six)*log(fdr(i)))
c        
        ! Temperature update
        IF (jthe == 0 .AND. cp > zero .AND. inloc == 0) THEN 
          dtemp     = weitemp(i)*yld(i)*ddep*eta/(rho(i)*cp)
          temp(i)   = temp(i) + dtemp
          ftherm(i) = one - mtemp* (temp(i) - tref)
        ENDIF
c        
        ! Hardening law update
        fhard(i) = yld0 + hard * pla(i) + qvoce*(one-exp(-bvoce*pla(i)))
c        
        ! Yield stress update
        yld(i) = fhard(i) * frate(i) * ftherm(i)
        yld(i) = max(yld(i), em10)
c        
        ! Yield function value update
        sigdr2 = sigdr(i)**2
        yld2i  = one / yld(i)**2
        phi(i) = sigdr2 * yld2i - one  
c        
        ! Transverse strain update
        IF (inloc == 0) THEN
          dezz(i) = dezz(i) + nnu*dpxx + nnu*dpyy + dpzz  
        ENDIF                
      ENDDO 
      !===================================================================
      ! - END OF PLASTIC CORRECTION WITH NICE - EXPLICIT METHOD
      !===================================================================     
c      
      ! Storing new values
      DO i=1,nel  
        ! USR Outputs & coefficient for hourglass
        IF (dpla(i) > zero) THEN 
          uvar(i,1) = phi(i)  ! Yield function value
          et(i)     = hardp(i)*frate(i) / (hardp(i)*frate(i) + young)  
        ELSE
          uvar(i,1) = zero
          et(i)     = one
        ENDIF
        seq(i)     = sigdr(i)
        ! Plastic strain-rate (filtered)
        dpdt       = dpla(i) / max(em20,timestep)
        epsd(i)    = afiltr * dpdt + (one - afiltr) * epsd(i)
        ! Computation of the sound speed   
        soundsp(i) = sqrt(a11/rho(i))
        ! Storing the yield stress
        sigy(i)    = yld(i)  
        ! Non-local thickness variation + temperature
        IF ((inloc > 0).AND.(off(i) == one).AND.(dplanl(i)>zero)) THEN 
          yld2i  = one/(yld(i)**2) 
          dphi_dsig = two*sigdr(i)*yld2i   
          normsig = sqrt(signxx(i)*signxx(i)
     .                 + signyy(i)*signyy(i)
     .            + two*signxy(i)*signxy(i))
          normsig = max(normsig,one)
          fdr(i)  = (j2(i)/(normsig**2))**3 - cdr*((j3(i)/(normsig**3))**2)  
          IF (fdr(i) > zero) THEN  
            dphi_dfdr = dphi_dsig*kdr*(one/six)*exp(-(five/six)*log(fdr(i)))   
          ELSE
            dphi_dfdr = zero
          ENDIF
          dsdrdj2 =  dphi_dfdr*three*(j2(i)/(normsig**2))**2                             
          dsdrdj3 = -dphi_dfdr*two*cdr*(j3(i)/(normsig**3))                                                   
          dj3dsxx =  two_third*(syy(i)*szz(i))/(normsig**2)
     .                 - third*(sxx(i)*szz(i))/(normsig**2)
     .                 - third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)                 
          dj3dsyy =    - third*(syy(i)*szz(i))/(normsig**2)
     .             + two_third*(sxx(i)*szz(i))/(normsig**2)
     .                -  third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2)             
          dj3dszz =    - third*(syy(i)*szz(i))/(normsig**2)
     .                 - third*(sxx(i)*szz(i))/(normsig**2)
     .             + two_third*(sxx(i)*syy(i)-sxy(i)**2)/(normsig**2) 
          dj3dsxy =  two*(sxx(i)*sxy(i) + sxy(i)*syy(i))/(normsig**2)                   
          normxx  = dsdrdj2*sxx(i)/normsig + dsdrdj3*dj3dsxx
          normyy  = dsdrdj2*syy(i)/normsig + dsdrdj3*dj3dsyy
          normzz  = dsdrdj2*szz(i)/normsig + dsdrdj3*dj3dszz
          normxy  = two*dsdrdj2*sxy(i)/normsig + dsdrdj3*dj3dsxy
          sig_dfdsig = signxx(i) * normxx
     .               + signyy(i) * normyy
     .               + signxy(i) * normxy
          IF (sig_dfdsig > em01) THEN
            dlam_nl = (yld(i)*dplanl(i))/sig_dfdsig
            dezz(i) = dezz(i) + nnu*(dlam_nl*normxx)
     .                        + nnu*(dlam_nl*normyy)
     .                        + dlam_nl*normzz
          ENDIF
          IF (cp > zero .AND. jthe == 0) THEN     
            ! Computation of the weight factor
            dpdt_nl = dplanl(i)/max(timestep,em20)
            IF (dpdt_nl < dpis) THEN
              weitemp(i) = zero
            ELSEIF (dpdt_nl > dpad) THEN
              weitemp(i) = one
            ELSE
              weitemp(i) = ((dpdt_nl-dpis)**2 )
     .                  * (three*dpad - two*dpdt_nl - dpis)
     .                  / ((dpad-dpis)**3)
            ENDIF
            dtemp   = weitemp(i)*dplanl(i)*yld(i)*eta/(rho(i)*cp)
            temp(i) = temp(i) + dtemp 
          ENDIF
        ENDIF
        ! Computation of the thickness variation  
        thk(i) = thk(i) + dezz(i)*thkly(i)*off(i)  
      ENDDO
c
      END