using improved dEq0/dNeq instead of dEq/dNeq with 0.0 argument

code/crystallite.f90

@@ -508,7 +508,7 @@ end subroutine crystallite_init
 subroutine crystallite_stressAndItsTangent(updateJaco)
  use prec, only: &
    tol_math_check, &
-   dNeq
+   dNeq0
  use numerics, only: &
    subStepMinCryst, &
    subStepSizeCryst, &
@@ -776,7 +776,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
          endif
        else
          subFracIntermediate = maxval(crystallite_subFrac, mask=.not.crystallite_localPlasticity)
-         if (dNeq(subFracIntermediate,0.0_pReal)) then
+         if (dNeq0(subFracIntermediate)) then
            crystallite_neighborEnforcedCutback = .false.  ! look for ips that require a cutback because of a nonconverged neighbor
            !$OMP PARALLEL
            !$OMP DO PRIVATE(neighboring_e,neighboring_i)
@@ -817,7 +817,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
            !$OMP DO PRIVATE(neighboring_e,neighboring_i)
            do e = FEsolving_execElem(1),FEsolving_execElem(2)
              do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
-               if (.not. crystallite_localPlasticity(1,i,e) .and. dNeq(crystallite_subFrac(1,i,e),0.0_pReal)) then
+               if (.not. crystallite_localPlasticity(1,i,e) .and. dNeq0(crystallite_subFrac(1,i,e))) then
                  do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
                    neighboring_e = mesh_ipNeighborhood(1,n,i,e)
                    neighboring_i = mesh_ipNeighborhood(2,n,i,e)
@@ -3070,7 +3070,7 @@ end subroutine crystallite_integrateStateFPI
 logical function crystallite_stateJump(ipc,ip,el)
  use prec, only: &
    prec_isNaN, &
-   dNeq
+   dNeq0
  use debug, only: &
    debug_level, &
    debug_crystallite, &
@@ -3122,7 +3122,7 @@ logical function crystallite_stateJump(ipc,ip,el)
  enddo
 
 #ifndef _OPENMP
- if (any(dNeq(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c),0.0_pReal)) &
+ if (any(dNeq0(plasticState(p)%deltaState(1:mySizePlasticDeltaState,c))) &
      .and. iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
      .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
              .or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
@@ -3178,7 +3178,7 @@ logical function crystallite_integrateStress(&
  use prec, only:         pLongInt, &
                          tol_math_check, &
                          prec_isNaN, &
-                         dEq
+                         dEq0
  use numerics, only:     nStress, &
                          aTol_crystalliteStress, &
                          rTol_crystalliteStress, &
@@ -3325,7 +3325,7 @@ logical function crystallite_integrateStress(&
  !* inversion of Fp_current...
 
  invFp_current = math_inv33(Fp_current)
- failedInversionFp: if (all(dEq(invFp_current,0.0_pReal))) then
+ failedInversionFp: if (all(dEq0(invFp_current))) then
 #ifndef _OPENMP
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
      write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
@@ -3341,7 +3341,7 @@ logical function crystallite_integrateStress(&
  !* inversion of Fi_current...
 
  invFi_current = math_inv33(Fi_current)
- failedInversionFi: if (all(dEq(invFi_current,0.0_pReal))) then
+ failedInversionFi: if (all(dEq0(invFi_current))) then
 #ifndef _OPENMP
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
      write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip ipc ',&
@@ -3600,7 +3600,7 @@ logical function crystallite_integrateStress(&
  invFp_new = math_mul33x33(invFp_current,B)
  invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal)                               ! regularize by det
  Fp_new = math_inv33(invFp_new)
- failedInversionInvFp: if (all(dEq(Fp_new,0.0_pReal))) then
+ failedInversionInvFp: if (all(dEq0(Fp_new))) then
 #ifndef _OPENMP
    if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
      write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip ipc ',&

code/homogenization_RGC.f90

@@ -387,7 +387,7 @@ end subroutine homogenization_RGC_partitionDeformation
 !--------------------------------------------------------------------------------------------------
 function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  use prec, only: &
-   dEq
+   dEq0
  use debug, only: &
    debug_level, &
    debug_homogenization,&
@@ -443,7 +443,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
  real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
  real(pReal), dimension(:), allocatable   :: resid,relax,p_relax,p_resid,drelax
  
- zeroTimeStep: if(dEq(dt,0.0_pReal)) then
+ zeroTimeStep: if(dEq0(dt)) then
    homogenization_RGC_updateState = .true.                                                          ! pretend everything is fine and return
    return                                                                    
  endif zeroTimeStep

code/plastic_disloUCLA.f90

@@ -974,7 +974,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent
 subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
    tol_math_check, &
-   dEq
+   dEq0
  use math, only: &
    pi
  use material, only: &
@@ -1112,7 +1112,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
      !* Dipole formation
      EdgeDipMinDistance = &
        plastic_disloUCLA_CEdgeDipMinDistance(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)
-     if (dEq(tau_slip_pos,0.0_pReal)) then
+     if (dEq0(tau_slip_pos)) then
        DotRhoDipFormation = 0.0_pReal
      else
        EdgeDipDistance = &
@@ -1140,7 +1140,7 @@ subroutine plastic_disloUCLA_dotState(Tstar_v,Temperature,ipc,ip,el)
         plastic_disloUCLA_CAtomicVolume(instance)*plastic_disloUCLA_burgersPerSlipSystem(j,instance)**(3.0_pReal)
      VacancyDiffusion = &
         plastic_disloUCLA_D0(instance)*exp(-plastic_disloUCLA_Qsd(instance)/(kB*Temperature))
-     if (dEq(tau_slip_pos,0.0_pReal)) then
+     if (dEq0(tau_slip_pos)) then
        DotRhoEdgeDipClimb = 0.0_pReal
      else
        ClimbVelocity = &
@@ -1174,7 +1174,7 @@ end subroutine plastic_disloUCLA_dotState
 function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
    tol_math_check, &
-   dEq
+   dEq0
  use math, only: &
    pi
  use material, only: &
@@ -1402,7 +1402,7 @@ function plastic_disloUCLA_postResults(Tstar_v,Temperature,ipc,ip,el)
           c = c + ns
         elseif(plastic_disloUCLA_outputID(o,instance) == stress_exponent_ID) then
           do j = 1_pInt, ns
-            if (dEq(gdot_slip_pos(j)+gdot_slip_neg(j),0.0_pReal)) then
+            if (dEq0(gdot_slip_pos(j)+gdot_slip_neg(j))) then
               plastic_disloUCLA_postResults(c+j) = 0.0_pReal
             else
               plastic_disloUCLA_postResults(c+j) = (tau_slip_pos(j)+tau_slip_neg(j))/&

code/plastic_dislotwin.f90

@@ -200,7 +200,8 @@ contains
 subroutine plastic_dislotwin_init(fileUnit)
  use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
  use prec, only: &
-   dEq, &
+   dEq0, &
+   dNeq0, &
    dNeq
  use debug, only: &
    debug_level,&
@@ -748,8 +749,8 @@ subroutine plastic_dislotwin_init(fileUnit)
       if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) &
         call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
       if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then
-        if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
+        if (dEq0(plastic_dislotwin_SFE_0K(instance)) .and. &
-            dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
+            dEq0(plastic_dislotwin_dSFE_dT(instance)) .and. &
                             lattice_structure(phase) == LATTICE_fcc_ID) &
           call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
         if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) &
@@ -758,8 +759,8 @@ subroutine plastic_dislotwin_init(fileUnit)
           call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOTWIN_label//')')
       endif
       if (sum(plastic_dislotwin_Ntrans(:,instance)) > 0_pInt) then
-        if (dEq(plastic_dislotwin_SFE_0K(instance), 0.0_pReal) .and. &
+        if (dEq0(plastic_dislotwin_SFE_0K(instance)) .and. &
-            dEq(plastic_dislotwin_dSFE_dT(instance),0.0_pReal) .and. &
+            dEq0(plastic_dislotwin_dSFE_dT(instance)) .and. &
                             lattice_structure(phase) == LATTICE_fcc_ID) &
           call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
         if (plastic_dislotwin_aTolTransFrac(instance) <= 0.0_pReal) &
@@ -772,7 +773,7 @@ subroutine plastic_dislotwin_init(fileUnit)
       if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
           plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) &
         call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
-      if (dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 0.0_pReal) .and. &
+      if (dNeq0(plastic_dislotwin_dipoleFormationFactor(instance)) .and. &
           dNeq(plastic_dislotwin_dipoleFormationFactor(instance), 1.0_pReal)) &
         call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')')
       if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
@@ -1624,7 +1625,7 @@ end subroutine plastic_dislotwin_microstructure
 subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
    tol_math_check, &
-   dNeq
+   dNeq0
  use math, only: &
    math_Plain3333to99, &
    math_Mandel6to33, &
@@ -1771,8 +1772,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
  
 !--------------------------------------------------------------------------------------------------
 ! Shear banding (shearband) part
- if(dNeq(plastic_dislotwin_sbVelocity(instance),  0.0_pReal) .and. &
+ if(dNeq0(plastic_dislotwin_sbVelocity(instance)) .and. dNeq0(plastic_dislotwin_sbResistance(instance))) then
-    dNeq(plastic_dislotwin_sbResistance(instance),0.0_pReal)) then
    gdot_sb = 0.0_pReal
    dgdot_dtausb = 0.0_pReal
    call math_eigenValuesVectorsSym(math_Mandel6to33(Tstar_v),eigValues,eigVectors,error)
@@ -1939,7 +1939,7 @@ end subroutine plastic_dislotwin_LpAndItsTangent
 subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
    tol_math_check, &
-   dEq
+   dEq0
  use math, only: &
    pi
  use material, only: &
@@ -2040,7 +2040,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
       !* Dipole formation
       EdgeDipMinDistance = &
         plastic_dislotwin_CEdgeDipMinDistance(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)
-      if (dEq(tau_slip(j),0.0_pReal)) then
+      if (dEq0(tau_slip(j))) then
         DotRhoDipFormation = 0.0_pReal
       else
         EdgeDipDistance = &
@@ -2068,10 +2068,10 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
         plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal)
       VacancyDiffusion = &
         plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature))
-      if (dEq(tau_slip(j),0.0_pReal)) then
+      if (dEq0(tau_slip(j))) then
         DotRhoEdgeDipClimb = 0.0_pReal
       else
-        if (dEq(EdgeDipDistance-EdgeDipMinDistance,0.0_pReal)) then
+        if (dEq0(EdgeDipDistance-EdgeDipMinDistance)) then
           DotRhoEdgeDipClimb = 0.0_pReal
         else
           ClimbVelocity = 3.0_pReal*lattice_mu(ph)*VacancyDiffusion*AtomicVolume/ &
@@ -2187,7 +2187,7 @@ end subroutine plastic_dislotwin_dotState
 function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
  use prec, only: &
    tol_math_check, &
-   dEq
+   dEq0
  use math, only: &
    pi, &
    math_Mandel6to33, &
@@ -2503,7 +2503,7 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
 
              !* Stress exponent
              plastic_dislotwin_postResults(c+j) = &
-               merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq(gdot_slip(j),0.0_pReal))
+               merge(0.0_pReal,(tau/gdot_slip(j))*dgdot_dtauslip,dEq0(gdot_slip(j)))
          enddo ; enddo
          c = c + ns
       case (sb_eigenvalues_ID)

code/plastic_isotropic.f90

@@ -517,7 +517,7 @@ end subroutine plastic_isotropic_LiAndItsTangent
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
  use prec, only: &
-   dEq
+   dEq0
  use math, only: &
    math_mul6x6
  use material, only: &
@@ -564,7 +564,7 @@ subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
 !--------------------------------------------------------------------------------------------------
 ! hardening coefficient
  if (abs(gamma_dot) > 1e-12_pReal) then
-   if (dEq(param(instance)%tausat_SinhFitA,0.0_pReal)) then
+   if (dEq0(param(instance)%tausat_SinhFitA)) then
      saturation = param(instance)%tausat
    else
      saturation = (  param(instance)%tausat &

code/plastic_nonlocal.f90

@@ -1542,7 +1542,7 @@ end subroutine plastic_nonlocal_aTolState
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_nonlocal_microstructure(Fe, Fp, ip, el)
 use prec, only: &
-  dEq
+  dEq0
 use IO, only: &
   IO_error
 use math, only: &
@@ -1786,7 +1786,7 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
                                   - neighbor_rhoExcess(c,s,neighbors(2))
       enddo
       invConnections = math_inv33(connections)
-      if (all(dEq(invConnections,0.0_pReal))) &
+      if (all(dEq0(invConnections))) &
         call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
       rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
                                          math_mul33x3(invConnections,rhoExcessDifferences))
@@ -2195,7 +2195,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el)
 use prec, only: &
-  dNeq
+  dNeq0
 use debug,    only: debug_level, &
                     debug_constitutive, &
                     debug_levelBasic, &
@@ -2318,8 +2318,8 @@ dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) / (4.0_pReal * pi * abs
 
 
 forall (c = 1_pInt:2_pInt)
-  where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
+  where(dNeq0(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+abs(rhoSgl(1:ns,2*c+3))&
-             abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
+                                     +abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)))) &
     dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & 
                        + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
                        dUpper(1:ns,c))
@@ -2331,7 +2331,7 @@ deltaDUpper = dUpper - dUpperOld
 !*** dissociation by stress increase
 deltaRhoDipole2SingleStress = 0.0_pReal
 forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
-                                        dNeq(dUpperOld(s,c) - dLower(s,c),0.0_pReal)) &
+                                        dNeq0(dUpperOld(s,c) - dLower(s,c))) &
   deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
                                            / (dUpperOld(s,c) - dLower(s,c))
 
@@ -2379,8 +2379,9 @@ subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
                                      timestep,subfrac, ip,el)
 
 use prec,     only: DAMASK_NaN, &
+                    dNeq0, &
                     dNeq, &
-                    dEq
+                    dEq0
 use numerics, only: numerics_integrationMode, &
                     numerics_timeSyncing
 use IO,       only: IO_error
@@ -2614,8 +2615,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
 dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
                / (4.0_pReal * pi * abs(tau))
 forall (c = 1_pInt:2_pInt)
-  where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
+  where(dNeq0(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+abs(rhoSgl(1:ns,2*c+3))&
-             abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
+                                     +abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)))) &
     dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & 
                        + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
                        dUpper(1:ns,c))
@@ -2827,11 +2828,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
       my_rhoSgl = rhoSgl
       my_v = v
       if(numerics_timeSyncing) then
-        if (dEq(subfrac(1_pInt,ip,el),0.0_pReal)) then
+        if (dEq0(subfrac(1_pInt,ip,el))) then
           my_rhoSgl = rhoSgl0
           my_v = v0
         elseif (neighbor_n > 0_pInt) then
-          if (dEq(subfrac(1_pInt,neighbor_ip,neighbor_el),0.0_pReal)) then
+          if (dEq0(subfrac(1_pInt,neighbor_ip,neighbor_el))) then
             my_rhoSgl = rhoSgl0
             my_v = v0
           endif
@@ -3170,7 +3171,7 @@ end subroutine plastic_nonlocal_updateCompatibility
 !*********************************************************************
 function plastic_nonlocal_dislocationstress(Fe, ip, el)
 use prec, only: &
-  dEq
+  dEq0
 use math,     only: math_mul33x33, &
                     math_mul33x3, &
                     math_inv33, &
@@ -3383,7 +3384,7 @@ if (.not. phase_localPlasticity(ph)) then
                     Rsquare = R * R
                     Rcube = Rsquare * R 
                     denominator = R * (R + flipSign * lambda)
-                    if (dEq(denominator,0.0_pReal)) exit ipLoop
+                    if (dEq0(denominator)) exit ipLoop
                       
                     sigma(1,1) = sigma(1,1) - real(side,pReal) &
                                             * flipSign * z / denominator &
@@ -3428,7 +3429,7 @@ if (.not. phase_localPlasticity(ph)) then
                     Rsquare = R * R
                     Rcube = Rsquare * R 
                     denominator = R * (R + flipSign * lambda)
-                    if (dEq(denominator,0.0_pReal)) exit ipLoop
+                    if (dEq0(denominator)) exit ipLoop
                     
                     sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z &
                                                                * (1.0_pReal - lattice_nu(ph)) / denominator &
@@ -3518,7 +3519,7 @@ end function plastic_nonlocal_dislocationstress
 !--------------------------------------------------------------------------------------------------
 function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
  use prec, only: &
-   dNeq
+   dNeq0
  use math, only: &
    math_mul6x6, &
    math_mul33x3, &
@@ -3635,8 +3636,8 @@ dUpper(1:ns,1) = lattice_mu(ph) * burgers(1:ns,instance) &
 dUpper(1:ns,2) = lattice_mu(ph) * burgers(1:ns,instance) &
                / (4.0_pReal * pi * abs(tau))
 forall (c = 1_pInt:2_pInt)
-  where(dNeq(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+&
+  where(dNeq0(sqrt(rhoSgl(1:ns,2*c-1)+rhoSgl(1:ns,2*c)+abs(rhoSgl(1:ns,2*c+3))&
-             abs(rhoSgl(1:ns,2*c+3))+abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)),0.0_pReal)) &
+                                     +abs(rhoSgl(1:ns,2*c+4))+rhoDip(1:ns,c)))) &
     dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & 
                        + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
                        dUpper(1:ns,c))

code/plastic_phenoplus.f90

@@ -113,7 +113,7 @@ contains
 subroutine plastic_phenoplus_init(fileUnit)
  use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
  use prec, only: &
-   dEq
+   dEq0
  use debug, only: &
    debug_level, &
    debug_constitutive,&
@@ -479,8 +479,7 @@ subroutine plastic_phenoplus_init(fileUnit)
      if (any(plastic_phenoplus_tausat_slip(:,instance) <= 0.0_pReal .and. &
              plastic_phenoplus_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPLUS_label//')')
-     if (any(dEq(plastic_phenoplus_a_slip(instance),0.0_pReal) .and. &
+     if (any(dEq0(plastic_phenoplus_a_slip(instance)) .and. plastic_phenoplus_Nslip(:,instance) > 0)) &
-             plastic_phenoplus_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPLUS_label//')')
      if (any(plastic_phenoplus_tau0_twin(:,instance) < 0.0_pReal .and. &
              plastic_phenoplus_Ntwin(:,instance) > 0)) &
@@ -918,7 +917,7 @@ end subroutine plastic_phenoplus_microstructure
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
  use prec, only: &
-   dNeq
+   dNeq0
  use math, only: &
    math_Plain3333to99, &
    math_Mandel6to33
@@ -1034,7 +1033,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
                (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (dNeq(gdot_slip_pos,0.0_pReal)) then
+     if (dNeq0(gdot_slip_pos)) then
        dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenoplus_n_slip(instance)/tau_slip_pos
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -1042,7 +1041,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
                                                      nonSchmid_tensor(m,n,1)
      endif
 
-     if (dNeq(gdot_slip_neg,0.0_pReal)) then
+     if (dNeq0(gdot_slip_neg)) then
        dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenoplus_n_slip(instance)/tau_slip_neg
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -1069,7 +1068,7 @@ subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
      Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (dNeq(gdot_twin,0.0_pReal)) then
+     if (dNeq0(gdot_twin)) then
        dgdot_dtautwin = gdot_twin*plastic_phenoplus_n_twin(instance)/tau_twin
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &

code/plastic_phenopowerlaw.f90

@@ -125,7 +125,7 @@ contains
 subroutine plastic_phenopowerlaw_init(fileUnit)
  use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
  use prec, only: &
-   dEq
+   dEq0
  use debug, only: &
    debug_level, &
    debug_constitutive,&
@@ -485,8 +485,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit)
      if (any(plastic_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. &
              plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
-     if (any(dEq(plastic_phenopowerlaw_a_slip(instance),0.0_pReal) .and. &
+     if (any(dEq0(plastic_phenopowerlaw_a_slip(instance)) .and. plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
-             plastic_phenopowerlaw_Nslip(:,instance) > 0)) &
        call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
      if (any(plastic_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. &
              plastic_phenopowerlaw_Ntwin(:,instance) > 0)) &
@@ -769,7 +768,7 @@ end subroutine plastic_phenopowerlaw_aTolState
 !--------------------------------------------------------------------------------------------------
 subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
  use prec, only: &
-   dNeq
+   dNeq0
  use math, only: &
    math_Plain3333to99, &
    math_Mandel6to33
@@ -859,7 +858,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
                (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (dNeq(gdot_slip_pos,0.0_pReal)) then
+     if (dNeq0(gdot_slip_pos)) then
        dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenopowerlaw_n_slip(instance)/tau_slip_pos
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -867,7 +866,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
                                                      nonSchmid_tensor(m,n,1)
      endif
 
-     if (dNeq(gdot_slip_neg,0.0_pReal)) then
+     if (dNeq0(gdot_slip_neg)) then
        dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
@@ -894,7 +893,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,
      Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
 
      ! Calculation of the tangent of Lp
-     if (dNeq(gdot_twin,0.0_pReal)) then
+     if (dNeq0(gdot_twin)) then
        dgdot_dtautwin = gdot_twin*plastic_phenopowerlaw_n_twin(instance)/tau_twin
        forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
          dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &