c     ****************************************************************          
c     *                                                              *          
c     *                      subroutine chkcrack                     *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 12/6/20 rhd                *          
c     *                                                              *          
c     *        This routine drives the check on conditions to tell   *          
c     *        if crack growth will occur.  This calls routines for  *          
c     *        both crack growth and discrete crack growth.          *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
c                                                                               
c                                                                               
      subroutine chkcrack( step, iter )   
c                                      
      use global_data, only : out
      use damage_data, only : crack_growth_type
c                                                          
      implicit none                                                    
c                                                                               
      integer, intent(in) :: step, iter
      logical :: debug, ldebug                                                            
c                                                                               
      debug  = .false.    
      ldebug = .false.                                                       
c                                                                               
      if( ldebug ) write(out,*) '>>> entered chkcrack ...'                        
c                                                                               
c              Return if no crack growth specified.                             
c                                                                               
c              crack_growth_type:                                               
c                                                                               
c                 (1)  element extinction by critical porosity                  
c                      in Gurson-Tvergaard material model                       
c                      or the extended Gurson-Tvergaard model                   
c                 (2)  node release on attainment of critical                   
c                      ctoa                                                     
c                 (3)  element extinction by stress-modified                    
c                      critical strain criterion. any material                  
c                      model.                                                   
c                 (4)  crack growth using cohesive elements              
c                 (5)  user requested by input command. element list
c                      provided in module damage_data       
c                                                                               
c              Lower level routines are separated into other .f files by        
c              the type of crack growth.                                        
c                                                                               
c                                                                               
      select case( crack_growth_type )                                           
      case(0)                                                                  
         if(debug) write (out,*) '>>>> no crack growth specified.'             
      case(1)                                                                  
         if(debug) write (out,*) '>>>> use gurson crack growth.'               
         call chk_elem_kill( debug, step, iter )                                
      case(2)                                                                  
         if(debug) write (out,*) '>>>> use node release crack growth.'         
         call chk_node_release( debug, step, iter )                             
      case(3)                                                                  
         if(ldebug) write (out,*) '>>>> use smcs crack growth.'                 
         call chk_elem_kill( debug, step, iter )                                
      case(4)    !  cohesive                                                               
         call chk_elem_kill( debug, step, iter )                                
         return                                                                 
      case(5)                                                                  
         if(debug) write (out,*) '>>>> user specified kill element.'      
         call chk_elem_kill( debug, step, iter )                                
      end select                                                                
c                                                                               
      if( debug ) call dam_debug    
      if( ldebug ) then 
        write(out,*) '>>> leaving chkcrack ...'                       
        write(out,*) ' '
      end if
c                                                                               
      return                                                                    
      end 

c     ****************************************************************          
c     *                                                              *          
c     *                      subroutine chk_elem_kill                *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 9/15/21 rhd                *          
c     *                                                              *          
c     *        checks conditions to see if crack growth will occur   *          
c     *        by element extinction. if so start the deletion.      *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
c                                                                               
      subroutine chk_elem_kill( debug, step, iter )
c                             
      use global_data, only : out, noelem, current_load_time_step, 
     &                        iprops, props, nelblk, elblks, num_error,
     &                        mxgp
      use dam_param_code, only : dam_param, growth_set_dbar
      use elem_extinct_data, only : dam_state, smcs_start_kill_step,                 
     &                              kill_order_list, dam_print_list,
     &                              smcs_d_values,
     &                              smcs_eps_plas_at_death,
     &                              smcs_stress_at_death,
     &                              Oddy_metrics
      use main_data, only : output_packets, packet_file_no
      use damage_data, only : no_killed_elems, dam_ptr, max_dam_state,
     &     crack_growth_type, num_user_kill_elems, release_type, 
     &     user_kill_list_now, all_elems_killed, killed_element_limit,
     &     num_elements_killed, kill_order, smcs_deleted_list_file_flag, 
     &     smcs_deleted_list_file_name, num_kill_order_list, smcs_type,
     &     stop_killed_elist_length, deleted_elist_to_stop, 
     &     distortion_plastic_limit, smcs_type_5_tp_critical,
     &     use_mesh_regularization,
     &     tol_regular =>
     &     tolerance_mesh_regularization, use_distortion_metric, 
     &     Oddy_critical_ratio    
      use constants                   
c                                                         
      implicit none
c
      integer, intent(in) :: step, iter
      logical, intent(in) :: debug                                                                                                                                   
c                                                                               
c               local declarations                                               
c
      integer, parameter :: max_local_list=500 ! max elements that
c                                                can be killed this call    
c                                       
      integer :: local_packet(max_local_list), elem, elem_ptr,
     &           local_count, blk,  cohes_type, count, 
     &           element, felem, i, span
      integer, parameter :: local_pkt_type(5) = (/ 20, 0, 21, 10, 0 / )                   
      logical :: blk_killed, killed_this_time, distortion_killed,           
     &           kill_the_elem, ldummy1, ldummy2,
     &           local_debug, stop_solution, standard_kill, 
     &           elements_have_been_killed,  
     &           fgm_cohes, ext_gurson, option_exponential,                        
     &           option_ppr, local_write_packets, do_kill_in_order,                
     &           found_exponential, found_ppr, found_cavit,                        
     &           option_cavit, elem_deleted                                                     
      logical, external :: scan_entry_in_list, chkcrack_is_elem_deleted                                                                                                                       
      double precision ::                                                        
     &   porosity, eps_plas, eps_crit, d_old, d_new,           
     &   values(20), local_status(max_local_list,3), orig_poros,                
     &   ext_shape, ext_spacing, avg_triaxiality, avg_zeta,
     &   avg_bar_theta, tear_param, sig_mises, mises_at_death,
     &   sig_1, sig_1_at_death, sig_mean
c        
      local_debug = .false.
      if( local_debug ) then      
         write(out,*) '... entered chk_elem_kill ...'                                                    
         write(out,*) '    ... no_killed_elems:',no_killed_elems
      end if                                                                    
c                                                                               
c              1) Loop for all elements. process killable elements.
c                 start kill as required.           
c
c                 standard kill: immediately delete element. no
c                                [Ke], stresses, properties, histories.
c                                element forces save for prescribed
c                                reduction to zero in future steps.
c                 use_mesh_regularization: element stiffness, stresses
c                                gradually released to zero following
c                                a user-defined damage function.      
c                                                                               
      killed_this_time  = .false.                                                
      local_count       = 0                                                           
      found_exponential = .false.                                               
      found_ppr         = .false.                                               
      found_cavit       = .false.    
      stop_solution     = .false.   
      standard_kill     = .not. use_mesh_regularization                                       
c     
      do elem = 1, noelem
c          
         elem_ptr = dam_ptr(elem)                                               
         if( elem_ptr == 0 ) cycle     ! element not killable                                      
         if( local_debug ) write(out,*) 'killable element, dam state:',
     &                                  elem, dam_state(elem_ptr)
c                                                                               
c                If we previously started killing element, update it's             
c                unloading state vector or damage parameter, 
c                and then skip the rest of the          
c                loop. we do this for traction-separation law as well           
c                since other routines look at dam_state.                        
c              
         if( dam_state(elem_ptr) > 0 ) then ! already being killed
            if( use_mesh_regularization ) then
               call chk_elem_kill_upmr
               cycle 
            end if
            if( standard_kill ) then
               if( dam_state(elem_ptr) <=  max_dam_state )                        
     &             dam_state(elem_ptr) = dam_state(elem_ptr) + 1                  
               cycle   
            end if
            write(out,9100) elem
            call die_abort
         end if                                                                 
c                                                                               
c                 calculate damage parameter(s) for element,                    
c                 then check it. if exceeded, set-up data structures for        
c                 internal force release using fixed no.                        
c                 of load steps or the traction-separation law.                 
c                                                                               
         select case( crack_growth_type )
         case( 1 )  ! GTN
            call dam_param(
     &         elem=elem, kill_now=kill_the_elem, debug=debug,
     &         porosity=porosity, eps_plas=eps_plas, eps_crit=eps_crit,
     &         sig_mean=sig_mean, sig_mises=sig_mises, 
     &         ext_gurson=ext_gurson, ext_shape=ext_shape,
     &         ext_spacing=ext_spacing )
            distortion_killed = .false.
            if( .not. kill_the_elem ) then ! distortion may kill it
               if( use_distortion_metric ) call chk_gt_kill_distortion
            end if      
         case( 3 )  ! SMCS   
            call dam_param(
     &         elem=elem, kill_now=kill_the_elem, debug=debug,              
     &         eps_plas=eps_plas, eps_crit=eps_crit, 
     &         sig_mises=sig_mises,          
     &         avg_triaxiality=avg_triaxiality, avg_zeta=avg_zeta, 
     &         avg_bar_theta=avg_bar_theta, tear_param=tear_param,
     &         max_princ_stress=sig_1 )
         case( 4 )
            call dam_param_cohes( elem, kill_the_elem, debug,                  
     &                            values, 2 )  
         case( 5 )
            kill_the_elem = .false.
            do i = 1, num_user_kill_elems
             if( elem  == user_kill_list_now(i) ) then
                 kill_the_elem = .true.
                 exit
             end if
            end do
         end select 
c                              
         killed_this_time = killed_this_time .or. kill_the_elem                 
         if( .not. kill_the_elem ) cycle                                                                
c                                                                               
c                  for standard element deletion, the element is 
c                  immediately removed from the mesh by zeroing 
c                  strains, stresses, history, stiffness. forces 
c                  exerted by element on model nodes are gradually
c                  relaxed to zero using one of several methods.
c
c                  for mesh regularization, the element stiffness and
c                  internal forces are degraded towards zero with 
c                  increasing deformation based on a damage 
c                  parameter (d). When d > (1.0-tol), the element
c                  stiffness, strains, stresses, .. are zeroed and no
c                  further element computations occur
c
         local_count = local_count + 1                                          
         if( crack_growth_type == 1 ) call chk_output_kill_messages_1
         if( crack_growth_type == 3 ) call chk_output_kill_messages_3
         if( crack_growth_type == 4 ) call chk_output_kill_messages_4
         if( crack_growth_type == 5 ) call chk_output_kill_messages_5
c
         if( output_packets ) call chk_elem_kill_make_packets
c                                                                               
c                 initialize global variables for 1st killed element.
c                 for standard kill process (no damage tracking),  
c                 null props and history for this newly killed element          
c                 so subsequent stiffness will be zero. zero stresses,          
c                 start the element force reduction process, constrain          
c                 completely uncoupled nodes if required.                       
c     
         if( no_killed_elems ) then     
            call allocate_damage( 4 ) ! setup to find newly free nodes                                         
            if( release_type .eq. 2 ) call allocate_damage( 5 )                 
            call chk_dam_init2( debug )                                             
         end if                                                                 
         num_elements_killed = num_elements_killed + 1                       
         if( standard_kill ) then
           call chk_update_killed_energy( elem )                                   
           call chk_kill_element_now( elem, debug )                                    
           call chk_store_ifv( elem, elem_ptr, debug )                             
           call chk_update_node_elecnt( elem, debug )                              
           if( release_type .eq. 2 ) call growth_set_dbar( elem,              
     &             elem_ptr, debug, -1 ) 
         end if
         if( use_mesh_regularization ) then
            mises_at_death = sig_mises
            sig_1_at_death = sig_1
            smcs_eps_plas_at_death(elem_ptr) = eps_plas
            smcs_stress_at_death(elem_ptr) = sig_1_at_death
            smcs_start_kill_step(elem_ptr) = current_load_time_step
            dam_state(elem_ptr) = 1
            call chk_get_d( elem, zero, sig_1_at_death, d_new, out )
            d_old = smcs_d_values(elem_ptr)
            if( d_new < d_old .and. iter > 1) then
              write(out,9020) elem, d_old, d_new
              num_error = num_error + 1 ! stop at next compute cmd
            end if
            smcs_d_values(elem_ptr) = d_new
         end if
c
c                 the user may have provided a list if elements.
c                 if this newly killed element is is the list,
c                 stop solution at cleanup in this routine.
c
         if( stop_killed_elist_length > 0 ) then
           if( scan_entry_in_list( elem, deleted_elist_to_stop, 
     &         stop_killed_elist_length ) ) stop_solution = .true.
         end if

c                         
      end do  ! over all model elements                                         
c                                                                               
c                                                                               
c         (2) packet output for newly killed elements at                        
c             for gurson, extended-gurson, smcs, cohesive                       
c                                                                               
      local_write_packets = output_packets .and. local_count .gt. 0             
      if( local_write_packets ) call chk_elem_kill_output_packets
c                                                                               
c         (3) are all elements in model now killed ?
c             std kill means all killable elements have dam state >0
c             regularization kill means damage > one
c          
      all_elems_killed = .true.
      do elem = 1, noelem
         elem_deleted = chkcrack_is_elem_deleted( elem )
         if( elem_deleted ) cycle
         all_elems_killed = .false.
         exit
      end do
c                                                                               
c         (4) if a killing order has been specified then make sure that         
c             no element "holes" have been left.                                
c 
      if( local_debug )write(out,*) '.. kill_order, killed this time:',
     &      kill_order, killed_this_time                                                     
      do_kill_in_order = kill_order .and. killed_this_time                      
      if( do_kill_in_order ) call chk_elem_kill_in_order 
c                                                                               
c                                                                               
c         (5) deprecated code
c                                                                               
c         (6) re-check for free nodes if any elements have ever                 
c             been killed. Just makes triple sure that there are                
c             no free nodes.                                                    
c              
      elements_have_been_killed =  .not. no_killed_elems                                                                 
      if( elements_have_been_killed ) call chk_free_nodes( debug )                 
c                                                                               
c         (7) if MPI, then all processors need to know what blocks              
c             and elements have been killed.  send dam_state,...          
c             Also have processors kill any elements            
c             which they own which should be killed.                            
c             this is a dummy for non-MPI                                      
c                                                                               
      call wmpi_send_growth( killed_this_time )                                
c                                                                               
c         (8) if limit on number of killed elements is exceeded,
c             stop the solution.
c
      if( num_elements_killed >= killed_element_limit ) then
          write(out,9000) killed_element_limit
          call store( ' ','kill_limit_restart.db', ldummy1, ldummy2 )
          call warp3d_normal_stop
      end if
c                                                                               
c         (9) the user may have provided a list of elements
c             to stop the analysis if <any> of them are killed
c             stop the solution.
c
      if( stop_solution ) then
          write(out,9010) 
          call store( ' ','kill_limit_restart.db', ldummy1, ldummy2 )
          call warp3d_normal_stop
      end if
c      
      return    
c
 9000 format(/,1x,">>>>> Number of currently killed (deleted)",
     & " elements has reached the ",
     &       /,7x,"user-specified limit of: ",i6,".",
     &       /,7x,"Restart file: kill_limit_restart.db written.",
     &       /,7x,"Job ended normally."// )
 9010 format(/,1x,">>>>> An element just killed appears in the user",
     & " defined list to stop the solution.",
     &       /,7x,"Restart file: kill_limit_restart.db written.",
     &       /,7x,"Job ended normally."// )
 9020 format(/,1x,">>>>> Inconsistent condition for element: ",i8,
     &       /,1x,"      new damage paramter < old damage parameter",
     &       /,1x,"      values: ", 2f8.4,
     &       /,1x,"      in routine chk_elem_kill", 
     &       /,1x,"      Execution will stop at next Compute cmd",
     &       /)
 9100 format('>> FATAL ERROR: routine ',
     & 'chk_elem_kill. @ 1',
     & 10x,'inconsistent condition. for element:',i8,
     & ' job terminated.'//)
c                                                                               
c
      contains
c     ========
c
c
c     ****************************************************************
c     *                                                              *
c     *          internal subroutine chk_gt_kill_distortion          *
c     *                                                              *
c     *                   last modified : 8/29/21 rhd                *
c     *                                                              *
c     *     Gurson element did not fail by critical porosity.        *
c     *     Check if it fails by distortion metrics                  *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine chk_gt_kill_distortion
c
      implicit none
c
      logical :: local_debug, kill_by_Oddy, kill_immediately,
     &           kill_by_plastic_limit
      double precision :: total_plastic_strain, max_Oddy
c
      local_debug = elem == -1
      if( local_debug ) write(out,9000) elem
c
c              current plastic strain
c
      total_plastic_strain = eps_plas
      if( eps_plas < zero ) total_plastic_strain = zero ! roundoff
c
c               could also delete by distortion
c               limit or plastic strain limit.
c
      max_Oddy = Oddy_metrics(elem_ptr,2) / Oddy_metrics(elem_ptr,1)
      kill_by_Oddy = max_Oddy >= Oddy_critical_ratio 
      kill_by_plastic_limit = total_plastic_strain >= 
     &                        distortion_plastic_limit
      kill_immediately = kill_by_Oddy  .or. kill_by_plastic_limit
      if( local_debug ) then
        write(out,9020) total_plastic_strain, max_Oddy
      end if  
      if( .not. kill_immediately ) return
c      
      kill_the_elem = .true.
      distortion_killed = .true. 
      return
c
 9000 format(5x,"... entered chk_xxx: ",i8)
 9010 format(10x,'eps_plas ',f12.8)
 9020 format(10x,'eps_plas, gt_max_Oddy: ',f8.5,f7.2)
 9200 format(/,' @2  >> element removed by distortion metrics: ',i7,
     & '  max Oddy ratio: ',f7.2,2x,'plastic strain: ',
     & f4.2,2x,a)
c
      end subroutine chk_gt_kill_distortion
cc     ****************************************************************
c     *                                                              *
c     *          internal subroutine chk_elem_kill_upmr              *
c     *                                                              *
c     *                   last modified : 8/28/21 rhd                *
c     *                                                              *
c     *     element in deletion via SMCS using mesh regularization.  *
c     *     update damage parameter. if fully deleted, updated       *
c     *     other data structures                                    *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine chk_elem_kill_upmr
c
      implicit none
c
c              element is in progress of being deleted. mesh
c              regularization is being used. Update the damage
c              parameter "d" for the element. It was zeroed
c              upon triggering deletion. When d >=1, the
c              element is to be fully deleted as in the standard
c              killing process.
c
      logical :: local_debug, kill_by_Oddy, kill_by_d_value, 
     &           kill_immediately, kill_by_plastic_limit
      double precision ::  avg_eps_plas, avg_sig_mises, d_now, d_old,
     &                     deps_plas, stress_at_death, avg_sig_1,
     &                     max_Oddy, total_plastic_strain, max_eps_plas
      character(len=50) :: distortion_msg
c
      local_debug = elem == -1
      if( local_debug ) write(out,9000) elem
c
      if( smcs_d_values(elem_ptr) > tol_regular ) return ! already fully killed
c
c              current plastic strain and mises stress averages
c              for element
c
      call chk_get_ele_vals( elem, avg_eps_plas, avg_sig_mises, 
     &                       avg_sig_1, max_eps_plas )
      deps_plas = avg_eps_plas - smcs_eps_plas_at_death(elem_ptr)
      if( local_debug ) then
        write(out,9010) avg_eps_plas, deps_plas, 
     &                  smcs_eps_plas_at_death(elem_ptr)
      end if
      if( deps_plas < zero ) deps_plas = zero ! round-off @ just killed
c
c               get current damage parameter value "d". decide to
c               fully delete element. could also delete by distortion
c               limit or plastic strain limit.
c
      stress_at_death = smcs_stress_at_death(elem_ptr)
      call chk_get_d( elem, deps_plas, stress_at_death, d_now, out )
      d_old = smcs_d_values(elem_ptr)
      if( d_now < d_old .and. iter > 1 ) then
         write(out,9040) elem, d_old, d_now
         num_error = num_error + 1 ! stop at next compute cmd
      end if
      smcs_d_values(elem_ptr) = d_now
      if( local_debug ) write(out,9020) d_old, d_now, stress_at_death
c
      if( use_distortion_metric ) then
         max_Oddy = Oddy_metrics(elem_ptr,2) / Oddy_metrics(elem_ptr,1)
         total_plastic_strain = max_eps_plas
      else
         max_Oddy = zero
         total_plastic_strain = -one ! actual value ust be larger
      end if 
c
      kill_by_Oddy = max_Oddy >= Oddy_critical_ratio 
      kill_by_d_value = d_now > tol_regular
      kill_by_plastic_limit = total_plastic_strain >= 
     &                        distortion_plastic_limit
c      
      kill_immediately = kill_by_Oddy  .or. kill_by_d_value .or. 
     &                   kill_by_plastic_limit
      if( .not. kill_immediately ) return
      smcs_d_values(elem_ptr) = hundred
c
c               fully kill the element. 
c               data structures like a standard (immediate) deletion
c               process.
c                   - energy ??
c                   - zero properties, stresses, history
c                   - decreases by 1 the count of elements attached to
c                     the now fully deleted element
c                   - if a node now has no elements attached, add
c                     constraints to prevent rbm
c
      if( local_debug )  write(out,9030)
      call chk_update_killed_energy( elem )     
      call chk_kill_element_now( elem, debug )                                    
      call chk_update_node_elecnt( elem, debug )  
      call chk_free_nodes( .false. ) ! debug flag  
      if( use_distortion_metric ) then
        distortion_msg = ' '
        if( kill_by_Oddy ) distortion_msg = '** Oddy exceeds limit'
        if( kill_by_plastic_limit ) distortion_msg =
     &                '** eps plastic exceeds limit'
        if( kill_by_Oddy .and. kill_by_plastic_limit )
     &     distortion_msg = '** Oddy & eps plastic exceed limits'
        write(out,9200) elem, d_now, max_Oddy, total_plastic_strain, 
     &                  distortion_msg(1:len(distortion_msg))
      else
        write(out,9210) elem, d_now
      end if  
c               
      return
c
 9000 format(5x"... entered chk_elem_kill_upmr. elem: ",i8)
 9010 format(10x,'avg eps pls, deps_pls, eps_pls @ death: ',3f12.8)
 
 9020 format(10x,'d_old: ',f6.3,'d_now: ',f6.3,' sig_1 @ death:',f8.3)
 9030 format(10x,'ready to leave. fully deleting element')
 9040 format(/,1x,">>>>> Inconsistent condition for element: ",i8,
     &       /,1x,"      new damage paramter < old damage parameter",
     &       /,1x,"      values: ", 2f8.4,
     &       /,1x,"      in routine chk_elem_kill_upmr", 
     &       /,1x,"      Execution will stop at next Compute cmd",
     &       /)
 9100 format('>> FATAL ERROR: routine ',
     & 'chk_elem_kill_update_mesh_regularization.',
     & 10x,'inconsistent condition. for element:',i8,
     & ' job terminated.'//)
 9200 format(/,' >> element removed by regularization: ',i7,
     & '  d: ', f5.2,'  max Oddy ratio: ',f7.2,2x,'plastic strain: ',
     & f4.2,2x,a)
 9210 format(/,' >> element removed by regularization: ',i7,
     & '  d: ', f5.2 )
c
      end subroutine chk_elem_kill_upmr
c     ****************************************************************
c     *                                                              *
c     *          internal subroutine chk_get_ele_vals                *
c     *                                                              *
c     *                   last modified : 8/26/21 rhd                *
c     *                                                              *
c     *         service routine to get strain-stress values for a    *
c     *         single element                                       *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine chk_get_ele_vals( elem, avg_eps_plas, avg_sig_mises,
     &                             avg_sig_1, max_eps_plas  )
c
      implicit none
c
c              return current values of average plastic strain and 
c              average mises stress for element
c
      integer, intent(in) :: elem
      double precision, intent(out) :: avg_eps_plas, avg_sig_mises,
     &                                 avg_sig_1, max_eps_plas
c
      integer :: elem_type, ngp
      logical :: threed_solid_elem, l1, l2, l3, l4, l5, l6, l7, l8,
     &           l9, l10
      logical :: ldebug 
      double precision :: workvec(27)
c
      ldebug = .false. ! elem .eq. 4 .and. step > 57
      if( ldebug ) 
     &    write(out,9000) elem
c
      elem_type = iprops(1,elem)
      ngp       = iprops(6,elem)    
c
      if( step < 2 ) then  ! sanity check. step passed to chkcrack
        call set_element_type( elem_type, threed_solid_elem, l1, l2, l3,
     &                       l4, l5, l6, l7, l8, l9, l10 )  
        if( .not. threed_solid_elem ) then
          write(out,9100) elem
          call die_abort
        end if
      end if
c 
      call mm_return_values( "plastic_strain", elem, workvec, ngp )
      avg_eps_plas = sum( workvec(1:ngp) ) * ( one / dble( ngp ) )
      max_eps_plas = maxval( workvec(1:ngp) )
c
      call mm_return_values( "avg_mises", elem, workvec, ngp )
      avg_sig_mises = workvec(1)
      call mm_return_values( "avg_princ_stress", elem, workvec, ngp )
      avg_sig_1 = workvec(1)
      if( ldebug ) then
        write(out,9010) avg_eps_plas, avg_sig_mises
      end if
      return
c
 9000 format(5x"... enter chk_get_ele_vals. elem: ",i8)
 9010 format(10x,'... leave chk_get_ele_vals with:',
     &   /,10x,'avg eps pls:',f12.8,' avg_sig_mises:',f8.3)
 9100 format('>> FATAL ERROR: routine ',
     & 'chk_get_ele_vals.',
     & 10x,'inconsistent condition. for element:',i8,
     & ' job terminated.'//)
c
      end subroutine chk_get_ele_vals
c
c
c     ****************************************************************
c     *                                                              *
c     *        internal subroutine chk_output_kill_messages_1        *
c     *                                                              *
c     *                   last modified : 9/19/21 rhd                *
c     *                                                              *
c     *                     Gurson crack growth                      *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_output_kill_messages_1
c
      implicit none
c
      integer :: device
      logical :: fexists, kill_by_Oddy, kill_by_plastic_limit
      double precision ::  total_plastic_strain, max_Oddy
      character(len=50) :: distortion_msg      
c
      if( ext_gurson ) then
         write(out,9002) elem, porosity, ext_spacing, ext_shape    
         return
      endif
c     
      if( distortion_killed ) then
        total_plastic_strain = eps_plas
        if( eps_plas < zero ) total_plastic_strain = zero ! roundoff
        max_Oddy = Oddy_metrics(elem_ptr,2) / Oddy_metrics(elem_ptr,1)
        kill_by_Oddy = max_Oddy >= Oddy_critical_ratio 
        kill_by_plastic_limit = total_plastic_strain >= 
     &                          distortion_plastic_limit
        distortion_msg = ' '
        if( kill_by_Oddy ) distortion_msg = '** Oddy exceeds limit'
        if( kill_by_plastic_limit ) distortion_msg =
     &                '** eps plastic exceeds limit'
        if( kill_by_Oddy .and. kill_by_plastic_limit )
     &     distortion_msg = '** Oddy & eps plastic exceed limits'
        write(out,9200) elem, porosity, max_Oddy, total_plastic_strain, 
     &                  distortion_msg(1:len(distortion_msg))
      else  !  killed by porosity
      if( .not. use_distortion_metric ) write(out,9000) elem, porosity
      if( use_distortion_metric )  write(out,9010) elem, porosity,
     &              Oddy_metrics(elem_ptr,2) / Oddy_metrics(elem_ptr,1)   
      end if 
c   
c              use same file name for list of deleted elements
c              as used for scs since both cannot be used in
c              same model 
      
      if( .not. smcs_deleted_list_file_flag ) return ! delete elements file
      inquire( file=smcs_deleted_list_file_name, exist=fexists )
      open( newunit=device,file=smcs_deleted_list_file_name,
     &      form='formatted', status='unknown', 
     &      access='sequential', position='append' )   
      if( .not. fexists ) then
        if( .not. use_distortion_metric ) write(device,9018)
        if( use_distortion_metric ) write(device,9019)
      end if   
      if( .not. use_distortion_metric )
     &  write(device,9020) current_load_time_step,
     &                   elem, porosity, sig_mean, sig_mises, eps_plas
      if( use_distortion_metric ) then
         max_Oddy = Oddy_metrics(elem_ptr,2) / Oddy_metrics(elem_ptr,1)
         write(device,9020) current_load_time_step,
     &                   elem, porosity, sig_mean, sig_mises, eps_plas,
     &                   max_Oddy
      end if         
      close(unit=device)    
                 

 9000 format(/,' @1  >> element death invoked for element: ',i7,                  
     & '.   f: ',f6.3)                                                          
 9002 format(/,'   >> element death invoked for element: ',i7,                  
     & '.   f, W, X: ',3f8.3)                                                   
 9010 format(/,' @1  >> element death invoked for element: ',i7,                  
     & '.   f: ',f6.3,2x,'Oddy metric: ',f8.3)                                                          
 9018 format('!',/,'!   Gurson quantities when elements deleted',/,
     & '!   (element deletion begins at listed step number)',
     & /,'!',/,
     & 6x,'step     elem   porosity     mean stress    mises stress',
     & '  plastic strain')
 9019 format('!',/,'!  Gurson quantities when elements deleted',/,
     & '!   (element deletion begins at listed step number)',
     & /,'!',/,
     & 6x,'step     elem   porosity     mean stress    mises stress',
     & '  plastic strain  Oddy metric')
 9020 format(2x,i8,i9,3x,f8.5,2x,e14.6,2x,e14.6,2x,e14.6,2x,f7.2)    
 9200 format(/,' @2  >> element removed by distortion metrics: ',i7,
     & '  f: ', f5.3,' max Oddy ratio: ',f7.2,2x,'plastic strain: ',
     & f4.2,2x,a)
c
      end subroutine chk_output_kill_messages_1
c
c
c     ****************************************************************
c     *                                                              *
c     *        internal subroutine chk_output_kill_messages_3        *
c     *                                                              *
c     *                   last modified : 8/29/21 rhd                *
c     *                                                              *
c     *                     SMCS crack growth                        *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_output_kill_messages_3
c
      implicit none
c
      integer :: device
      logical :: fexists
c
      if( smcs_type <= 4 ) then
        write(out,9010) elem, eps_plas, eps_crit, avg_triaxiality,
     &                avg_zeta, avg_bar_theta
      end if
      if( smcs_type == 5 ) then
        write(out,9025) elem, tear_param, smcs_type_5_tp_critical
      end if
      if( .not. smcs_deleted_list_file_flag ) return ! delete elements file
      inquire( file=smcs_deleted_list_file_name, exist=fexists )
      open( newunit=device,file=smcs_deleted_list_file_name,
     &      form='formatted', status='unknown', 
     &      access='sequential', position='append' )   
      if( .not. fexists ) write(device,9018)
      write(device,9020) current_load_time_step,
     &                   elem, eps_plas, eps_crit,
     &                   avg_triaxiality, avg_zeta, avg_bar_theta,
     &                   tear_param 
      close(unit=device)    
c
      return
c
 9010 format(/,' >> element death option invoked before next step',             
     &       /,'    element: ',i7,' is now killed.',                            
     &       /,'    plastic strain: ',f8.5,' is > limit of: ',f8.5,
     &       /,'    plastic strain weighted values: bar_T, ',
     &         'bar_xi, bar_theta:',f6.3,2x,f6.3,2x,f6.3)    
 9018 format('!',/,'!  SMCS quantities when critical plastic ',
     & 'strain reached for elements', /,
     & '!   (element deletion begins at listed step number)',
     & /,'!',/,
     & 3x,'step     elem   ebarp  ebarp_crit bar_T  bar_xi',
     &  '   bar_theta   tear_param')  
 9020 format(2x,2i8,f8.5,1x,f8.5,1x,f6.3,2x,f6.3,2x,f6.3,4x,f8.3)      
 9025 format(/,' >> element death option invoked before next step',             
     &       /,'    element: ',i7,' is now killed.',                            
     &       /,'    tearing parameter: ',f9.3,' is > limit of: ',
     &        f9.3 )
c 
      end subroutine chk_output_kill_messages_3
c
c
c     ****************************************************************
c     *                                                              *
c     *        internal subroutine chk_output_kill_messages_4        *
c     *                                                              *
c     *                   last modified : 8/29/21 rhd                *
c     *                                                              *
c     *                    cohesive element/material growth          *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_output_kill_messages_4
c
      implicit none
c
c                 start killing this element right now. record status               
c                 for subsequent packet output if neeeded. write usual              
c                 output about element death.                                   
c
      cohes_type         = iprops(27,elem)                               
      option_exponential = cohes_type .eq. 4                             
      option_ppr         = cohes_type .eq. 6                             
      option_cavit       = cohes_type .eq. 7                             
      found_exponential  = found_exponential .or.                         
     &                     option_exponential                             
      found_ppr   = found_ppr .or. option_ppr                            
      found_cavit = found_cavit .or. option_cavit                        
      count = 0                                                          
      if( found_exponential ) count = count + 1                          
      if( found_ppr ) count = count + 1                                  
      if( found_cavit ) count = count + 1                                
      if( count .gt. 1 ) then                                            
        write(out,9400) elem                                            
        call die_gracefully                                              
      end if                                                             
      if( option_exponential )                                           
     &   write(out,9200) elem, values(6)/values(7), values(8)            
      if( option_ppr )   write(out,9220) elem                            
      if( option_cavit ) write(out,9230) elem  
c
      return
c
 9200 format(/,'   >> element death invoked for element: ',i7,                  
     & '.   Deff/Dpeak: ',f5.2,' Teff/Tpeak: ',f5.2)                            
 9220 format(/,'   >> element death invoked for element: ',i7,                  
     & '.  PPR cohesive option')                                                
 9230 format(/,'   >> element death invoked for element: ',i7,                  
     & '.  cavit cohesive option')                                              
 9400 format(/,'FATAL ERROR: mixed cohesive options not allowed',               
     & /,      '             at present in WARP3D with crack growth.'           
     & /,      '             found while processing element: ',i7,              
     & /,      '             job aborted.' )  
c 
      end subroutine chk_output_kill_messages_4
c     ****************************************************************
c     *                                                              *
c     *        internal subroutine chk_output_kill_messages_5        *
c     *                                                              *
c     *                   last modified : 8/29/21 rhd                *
c     *                                                              *
c     *                 user directed element deletion               *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_output_kill_messages_5
c
      implicit none
c
      write(out,9240) elem                         
c
      return
c
 9240 format(/,'   >> element death invoked for element: ',i7,                  
     & '.  user-directed removal')                                              
c
      end subroutine chk_output_kill_messages_5
c     ****************************************************************
c     *                                                              *
c     *        internal subroutine chk_elem_kill_output_packets      *
c     *                                                              *
c     *                   last modified : 3/7/21 rhd                 *
c     *                                                              *
c     *     packet output for new element starting deletion          *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_elem_kill_output_packets
c      implicit none
c
      write(packet_file_no) local_pkt_type(crack_growth_type),                  
     &                      local_count, step, iter                             
      do elem = 1, local_count                                                  
        element = local_packet(elem)                                            
        select case( crack_growth_type )
        case( 1 )
          orig_poros = props(26,element)                                        
          write(packet_file_no) element,orig_poros,                             
     &                  local_status(elem,1),                                   
     &                  local_status(elem,2),local_status(elem,3)               
        case( 3 )
          write(packet_file_no) element, local_status(elem,1),                  
     &                                  local_status(elem,2)                    
        case( 4 ) 
          write(packet_file_no) element, local_status(elem,1),                  
     &          local_status(elem,2), int(local_status(elem,3))                 
        end select
      end do   
c                                                                 
      return
      end subroutine chk_elem_kill_output_packets
c
c     ****************************************************************
c     *                                                              *
c     *        internal subroutine chk_elem_kill_in_order            *
c     *                                                              *
c     *                   last modified : 4/13/21 rhd                *
c     *                                                              *
c     *     after starting deletion of new element, see if other     *
c     *     elements must also be killed to satisfy the user input   *
c     *     command to kill in order                                 *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_elem_kill_in_order
c
      use constants
c
      implicit none
c
      integer :: i, index_killed, elem, elem_ptr
      logical :: killed_found_in_order_list
      logical, parameter :: ldebug = .false.
      double precision :: eps_plas_at_death, mises_at_death_local,
     &                    sig_1_at_death_local, max_eps_plas
c
      if( ldebug ) write (out,*) '>>> chk_elem_kill_in_order....'       
c
c              have started deaths for 1 or more elements just now.
c              user has provided a kill-in-order list.
c              start at end of user supplied order list. work backwards.
c              find first killed element if any. kill all elements from 
c              that point to beginning of list unless already 
c              being killed.
c            
      killed_found_in_order_list = .false.
c
      do i = num_kill_order_list, 1, -1                                   
        elem     = kill_order_list(i)                                    
        elem_ptr = dam_ptr(elem)                                                
        if( dam_state(elem_ptr) == 0 ) cycle ! not yet killed
        killed_found_in_order_list = .true.
        index_killed = i
        exit
      end do ! i
c
      if( .not. killed_found_in_order_list ) return 
c
c              start killing elements from beginning of in-order list up to
c              first element in list already being killed.
c              
      do i = 1, index_killed-1
        elem     = kill_order_list(i)                                    
        elem_ptr = dam_ptr(elem)                                                
        if( dam_state(elem_ptr) .ne. 0 ) cycle ! already being killed
        write(out,9100) elem                                              
        num_elements_killed = num_elements_killed + 1  
        if( standard_kill ) then                    
          call chk_update_killed_energy( elem )                                  
          call chk_kill_element_now( elem, ldebug )                                   
          call chk_store_ifv( elem, elem_ptr, ldebug )                            
          call chk_update_node_elecnt( elem, ldebug )                             
          if ( release_type .eq. 2 )                                         
     &       call growth_set_dbar( elem, elem_ptr, ldebug, -1 ) 
        elseif( use_mesh_regularization ) then
          call chk_get_ele_vals( elem, eps_plas_at_death,
     &                           mises_at_death_local, 
     &                           sig_1_at_death_local,
     &                           max_eps_plas )
          smcs_eps_plas_at_death(elem_ptr) = eps_plas_at_death
          smcs_stress_at_death(elem_ptr) = sig_1_at_death_local
          dam_state(elem_ptr) = 1
          call chk_get_d( elem, zero, sig_1_at_death_local, 
     &                    smcs_d_values(elem_ptr), out )
        else
          write(out,9200) elem
          call die_abort
        end if                                                                  
      end do  ! i 
c   
      return
c
 9100 format(/,' >> element death option invoked before next step',             
     &       /,'    element: ',i7,' is now killed per user request',               
     &       /,'    via kill in order input command')    
 9200 format('>>>> FATAL ERROR: chk_elem_kill_in_order. elem',i8,
     &     /,10x,'Job terminated...',//)             
                                         
c
      end subroutine   chk_elem_kill_in_order 
c                                                             
c     ****************************************************************
c     *                                                              *
c     *        internal subroutine chk_elem_kill_make_packets        *
c     *                                                              *
c     *                   last modified : 3/7/21 rhd                 *
c     *                                                              *
c     *     new element starting deletion. output packets            *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_elem_kill_make_packets
c
      implicit none
c
      if( local_count .gt. max_local_list ) then                           
         write(out,9300)                                                   
         call die_gracefully                                               
      end if                                                               
      select case( crack_growth_type )
      case( 1 )
           local_packet(local_count) = elem                                
           local_status(local_count,1) = porosity                          
           local_status(local_count,2) = ext_spacing                       
           local_status(local_count,3) = ext_shape                         
      case( 3 ) 
           local_packet(local_count) = elem                                
           local_status(local_count,1) = eps_plas                          
           local_status(local_count,2) = eps_crit                          
      case( 4 )
           if( option_exponential ) then                                   
             local_packet(local_count) = elem                              
             local_status(local_count,1) = values(6)/values(7)             
             local_status(local_count,2) = values(8)                       
             local_status(local_count,3) = 4.0                             
           end if                                                          
           if( option_ppr ) then                                           
             local_packet(local_count) = elem                              
             local_status(local_count,1) = values(5)                       
             local_status(local_count,2) = values(6)                       
             local_status(local_count,3) = 6.0                             
           end if                                                          
           if( option_cavit) then                                          
             write(out,*) ' cavit death packets not implemented'             
             call die_abort                                                
           end if  
      end select 
c                                                       
      return
c
 9300 format(/,'FATAL ERROR: list length exceeded in chk_elem_kill',            
     & /,      '             job aborted.' )                                    
c
      end subroutine chk_elem_kill_make_packets
      end subroutine chk_elem_kill                                                                      
c
c     ****************************************************************          
c     *                                                              *          
c     *                      subroutine chk_dam_init2                *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 12/6/20 rhd                *          
c     *                                                              *          
c     *        called when 1st element is deleted to setup.          *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
      subroutine chk_dam_init2( debug ) 
c                                            
      use global_data, only : out, nonode, noelem, iprops, c ! coords
c                                                                               
      use elem_extinct_data, only : dam_node_elecnt, dam_face_nodes             
      use main_data, only : incmap, incid, elems_to_blocks,                     
     &                      inverse_incidences                                  
      use elem_block_data, only : cdest_blocks                                  
      use damage_data, only : no_killed_elems, release_type, dam_ptr, 
     &                        crk_pln_normal_idx, crack_plane_coord,
     &                        gurson_cell_size                                                           
c                                                                               
      implicit none                                                    
c                                                                               
      logical, intent(in) :: debug
c
      integer :: i, elem, elem_ptr, incptr, num_enodes, k1, k2,
     &           face_count, blk, rel_elem, enode, snode
      logical, parameter :: local_debug = .false.                                            
      double precision :: node_coord(3), coord
      double precision, parameter :: plane_tol = 0.01d0                                 
      integer, dimension(:,:), pointer :: cdest                                 
c                                                                               
      if( debug ) write(out,*) '>>>> in chk_dam_init2'                             
c                                                                               
      no_killed_elems = .false.                                                 
      do i = 1, nonode                                                          
         dam_node_elecnt(i) = inverse_incidences(i)%element_count               
      end do                                                                    
      if( debug ) then                                                         
         do i = 1, nonode                                                       
            write(out,9030) i, dam_node_elecnt(i)                              
         end do                                                                 
      end if                                                                    
c                                                                               
c               set up the force release by traction-displacement               
c               separation law if that option is being used.                    
c               for each killable element on the crack plane,                   
c               find the 4 nodes on the face opposite the crack                 
c               plane face used to define the element deformation.              
c               we examine only the first 8 nodes of the element                
c               (corner nodes).                                                 
c                                                                               
      if( release_type .ne. 2 ) then
        if( debug ) write(out,*) '<<<< leaving chkdam_init2'                        
        return                                                                    
      end if
c
      do elem  = 1, noelem                                                      
        elem_ptr = dam_ptr(elem)                                                
        if( elem_ptr .eq. 0 ) cycle                                            
        incptr     = incmap(elem)-1                                             
        num_enodes = iprops(2,elem)                                             
        k1         = num_enodes                                                 
        k2         = 2*num_enodes                                               
        face_count = 0                                                          
        blk        = elems_to_blocks(elem,1)                                    
        rel_elem   = elems_to_blocks(elem,2)                                    
        cdest      => cdest_blocks(blk)%ptr                                     
        do enode = 1, 8                                                         
          snode         = incid(incptr+enode)                                   
          node_coord(1) = c( cdest(enode,rel_elem) )                            
          node_coord(2) = c( cdest(k1+enode,rel_elem) )                         
          node_coord(3) = c( cdest(k2+enode,rel_elem) )                         
          coord         = node_coord(crk_pln_normal_idx)                        
          if( abs(coord-crack_plane_coord) .le. plane_tol*                     
     &        gurson_cell_size ) cycle                                         
          face_count = face_count + 1                                           
          dam_face_nodes(face_count,elem_ptr) = snode                           
        end do   ! enode                                                               
        if( face_count .ne. 4 ) then                                           
             write(out,*) 'FATAL error 1 in dam_init2'                          
             write(out,*) 'invalid plane definition for growth'                 
             write(out,*)  elem, elem_ptr, num_enodes,face_count                
             call die_gracefully                                                
             stop                                                               
        end if                                                                  
      end do ! elem                                                                  
c                                                                              
      if( local_debug ) then                                                  
       write(out,*) '> element release type 2. face node table'                
       do elem  = 1, noelem                                                    
        elem_ptr = dam_ptr(elem)                                               
        if( elem_ptr .eq. 0 ) cycle                                           
        write(out,9050) elem, dam_face_nodes(1,elem_ptr),                      
     &       dam_face_nodes(2,elem_ptr), dam_face_nodes(3,elem_ptr),         
     &       dam_face_nodes(4,elem_ptr)                                      
       end do                                                                  
      end if                                                                   
c                                                                                                                                             
      if( debug ) write(out,*) '<<<< leaving chk_dam_init2'                        
      return                                                                    
 9030 format(' dam_node_elecnt(',i4,')=',i2)                                    
 9050 format(' elem: ',i7, '  face nodes: ',8i8)                                
      end                                                                       
c     ****************************************************************          
c     *                                                              *          
c     *                      subroutine chk_free_nodes               *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 12/6/20 rhd                *          
c     *                                                              *          
c     *        check each node in the structure to                   *          
c     *        see if they are no longer connected to any element.   *          
c     *        if so, it constrains them.                            *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
      subroutine chk_free_nodes( debug )                                        
      use global_data ! old common.main
c                                                                               
      use elem_extinct_data, only : dam_node_elecnt                             
      use main_data,         only : cnstrn, cnstrn_in                          
      use damage_data, only : csttail  
      use constants                                                    
c                                                                               
      implicit none                                                    
c              
      logical, intent(in) :: debug
c
      integer :: node, num_dof, j, glbal_dof
c
      if( debug ) write(out,*) '>>> in chk_free_nodes <<<'                     
c                                                                               
c              Loop over each node -- if a node is now attatched to             
c              no elements, constrain it.                                       
c              Note: this could be more efficient -- we are looping             
c                over all the nodes each time, and we could just                
c                constrain the newly zeroed nodes.  This, however,              
c                will also take care of the problem of a user                   
c                putting in a new set of constraints.                           
c                                                                               
      do node = 1, nonode                                                       
       if( dam_node_elecnt(node) .ne. 0 ) cycle                             
       if( debug ) then                                                   
         write(out,*) 'free node ',node                                   
         write (out,9010) node,cnstrn(dstmap(node)),                      
     &          cnstrn(dstmap(node)+1),cnstrn(dstmap(node)+2)               
       end if                                                              
       do j = 1, 3    ! num_dof                                                   
         glbal_dof            = dstmap(node) + j-1                        
         cnstrn(glbal_dof)    = zero                                      
         cnstrn_in(glbal_dof) = zero                                      
         if( cstmap(glbal_dof) .eq. 0 ) then                             
           cstmap(csttail) = glbal_dof                                   
           csttail         = glbal_dof                                   
         end if                                                           
       end do                                                              
       cstmap(csttail) = -1                                                
      end do                                                                    
c                                                                               
      if ( debug ) write(out,*) '>>> in chk_free_nodes <<<'                     
      return                                                                    
 9010 format('  old constraints for node ',i7,' :',3(1x,e14.6))                 
      end                                                                       
c     ****************************************************************          
c     *                                                              *          
c     *                      subroutine dam_print                    *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 12/6/20 rhd                *          
c     *                                                              *          
c     *     print status of killable elements or released nodes at   *
c     *     the beginning of a load step.                            *       
c     *     optional smcs states output for smcs model
c     *                                                              *          
c     ****************************************************************          
c                                                                               
c                                                                               
      subroutine dam_print( step, iter ) 
c                                       
      use damage_data, only : crack_growth_type, print_status, 
     &                        print_top_list, const_front, smcs_states                                                           
c
      implicit none    
c
      integer, intent(in) :: step, iter  
c
      logical :: doprint                                              
c                                                                               
c           check to make sure crack growth is on and some type of
c           output is requested
c                                                                               
      select case( crack_growth_type )                                           
      case(0)       ! no action                                                             
      case(1)                                                                  
         doprint = print_status .or. print_top_list
         if( .not. doprint ) return
         call dam_print_elem1( step, iter )                                  
      case(2)
         doprint = print_status .or. print_top_list
         if( .not. doprint ) return
         if( const_front ) then                                                
           call dam_print_front( step, iter)                                
         else                                                                 
           call dam_print_node( step, iter )                                
         end if    
      case(3)                                                            
         doprint = print_status .or. print_top_list .or. smcs_states
         if( .not. doprint ) return
         call dam_print_elem3( step, iter )                                  
      case(4)
         doprint = print_status .or. print_top_list
         if( .not. doprint ) return
         call dam_print_elem4( step, iter )                                  
      end select                                                                    
c                                                                               
      return                                                                    
      end                                                                       
c     ****************************************************************          
c     *                                                              *          
c     *                      subroutine dam_debug                    *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 06/23/94                   *          
c     *                                                              *          
c     *     This routine prints out all of the crack growth          *          
c     *     parameters.                                              *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               

      subroutine dam_debug                                                      
      use global_data ! old common.main
      use elem_extinct_data                                                     
      use node_release_data                                                     
      use damage_data                                                           
c                                                                               
      implicit integer (a-z)                                                    
      double precision                                                          
     &     zero                                                                 
      data zero / 0.0 /                                                         
                                                                                
c                                                                               
      write (out,*)                                                             
      write (out,*) '>>>>>>>> dumping all relevent crack growth info:'          
      write (out,*)                                                             
      write (out,*) '> crack_growth_type:',crack_growth_type                    
c                                                                               
c              -- Gurson crack growth                                           
c                                                                               
      if( crack_growth_type.eq.1 ) then                                         
         write (out,*) '> max_dam_state:', max_dam_state                        
         write (out,*) '> num_kill_elem:', num_kill_elem                        
         write (out,*) '> dam_ptr:'                                             
         do i=1, noelem                                                         
            write (out,'(" dam_ptr(",i4,")=",i4)') i,dam_ptr(i)                 
         end do                                                                 
         write (out,*) '> dam_state:'                                           
         do i = 1, num_kill_elem                                                
            write (out,'("dam_state(",i3,")=",i3)')i,dam_state(i)               
         end do                                                                 
         write (out,*) '> dam_ifv:'                                             
         do i=1, num_kill_elem                                                  
            write (out,*) 'ifv entries for ptr:',i                              
            write(out,9000) (dam_ifv(j,i), j=1, mxedof)                         
         end do                                                                 
         write (out,'(" > porosity_limit:",e14.6)') porosity_limit              
         if( kill_order ) then                                                  
            write (out,*) '> kill_order_list:'                                  
            do i = 1, num_kill_order_list                                       
               write (out,'("  entry:",i7," elem:",i6)')i,                      
     &              kill_order_list(i)                                          
            end do                                                              
         end if                                                                 
         if( .not.no_killed_elems ) then                                        
            write(out,*) '> elems had been killed'                              
            write(out,*) '> dam_node_elecnt:'                                   
            do i = 1, nonode                                                    
               write (out,'(" dam_node_elecnt(",i5,")=",i3)')i,                 
     &              dam_node_elecnt(i)                                          
            end do                                                              
            write (out,*) '> dam_dbar_elems(current height, fraction)'          
            if( release_type .eq. 2 ) then                                      
               do i = 1, num_kill_elem                                          
                  write (out,'(e13.6,2x,e13.6)') dam_dbar_elems(1,i),           
     &                 dam_dbar_elems(2,i)                                      
               end do                                                           
            end if                                                              
         else                                                                   
            write (out,*) '> no elems have been killed.'                        
         end if                                                                 
c                                                                               
c              -- Node release crack growth                                     
c                                                                               
      else if (crack_growth_type.eq.2) then                                     
c                                                                               
         write (out,*) '> Number of crack plane nodes:',                        
     &        num_crack_plane_nodes                                             
         write (out,'(" > half of CTOA for release:",e13.6)')                   
     &        critical_angle                                                    
         write (out,*) '> entry, node, state, ifv:'                             
         do i=1,num_crack_plane_nodes                                           
            write (out,'(2x,i7,i6,i6,3x,e13.6)')i,crack_plane_nodes(i),         
     &           crkpln_nodes_state(i), crkpln_nodes_react(i)                   
         enddo                                                                  
         write (out,*) '> global nodes => crack plane nodes'                    
         do i=1, nonode                                                         
            write (out,*) ' ',i,':',inv_crkpln_nodes(i)                         
         enddo                                                                  
         write (out,*) '> crack front list:'                                    
         call write_list                                                        
         write (out,*) '> neighbor_nodes:'                                      
         do i=1, num_crack_plane_nodes                                          
            write (out,*) ' ',i,':  for node:',crack_plane_nodes(i),            
     &           ' # neighbors=', num_neighbors(i),' which are:'                
            write (out,'(7x,12i5)') (neighbor_nodes(j,i),                       
     &           j=1,num_neighbors(i))                                          
         enddo                                                                  
c                                                                               
         if (release_type .eq. 2) then                                          
            write(out,'(" > charlen,rel.frac,rel.height:",3e13.6)')             
     &           char_length, release_fraction, release_height                  
            write (out,*) '> node_release_frac:'                                
            do i = 1, num_crack_plane_nodes                                     
               write (out,'(2x,i7,2x,e13.6)')i,node_release_frac(i)             
            enddo                                                               
         endif                                                                  
c                                                                               
         if ( overshoot_control_crk_grth ) then                                 
            write (out,*) '> overshoot control is on:'                          
            write (out,'("   -control load fact:",e13.6)')                      
     &           control_load_fact                                              
            write (out,'("   -old load fact:",e13.6)') old_load_fact            
            write (out,'("   -min load fact:",e13.6)') min_load_fact            
            write (out,'("   -overshoot_limit:",e13.6)')                        
     &           overshoot_limit                                                
            write (out,'("   -CTOA range:",e13.6)') CTOA_range                  
            if (overshoot_allocated) then                                       
               write (out,*) '   -overshoot stuff is allocated'                 
               write (out,'("   -old angles:")')                                
               write (out,'("     entry  angles")')                             
               do i=1, num_crack_plane_nodes                                    
                  write (out,'(6x,i4,1x,12e13.6)') i,                           
     &                 (old_angles_at_front(i,j),j=1,12)                        
               enddo                                                            
            else                                                                
               write (out,*) '   -overshoot stuff is not allocated'             
            endif                                                               
         else                                                                   
            write (out,*) '> no overshoot control'                              
         endif                                                                  
c                                                                               
         if ( load_size_control_crk_grth ) then                                 
            write (out,*) '> load size control is on'                           
            write (out,'("   -perm load fact:",e13.6)') perm_load_fact          
            write (out,'("   -min steps:",i3)') min_steps_for_release           
         else                                                                   
            write (out,*) '> no load size control'                              
         endif                                                                  
c                                                                               
      if (const_front) then                                                     
       write (out,*) '> const front is set.'                                    
       write (out,*) '  -master nodes:'                                         
       do i=1, num_crack_fronts                                                 
         write (out,'(5x,i7)') master_nodes(i)                                  
       enddo                                                                    
       write (out,*) '  -crack_front_list:'                                     
       do i=1, num_crack_fronts * num_nodes_grwinc                              
        write (out,'(5x,i3,":",40i7)') i,                                       
     &             (crack_front_list(i,j),j=1,num_nodes_thick)                  
       enddo                                                                    
c                                                                               
       write (out,*) '> measured CTOA'                                          
       write (out,'("   -init_ctoa_dist:",e13.6)') init_ctoa_dist               
       write (out,'("   -ctoa_dist:",e13.6)') ctoa_dist                         
       write (out,*) '  -num_nodes_back =',num_nodes_back                       
       write (out,*) '  -master list (first 20 nodes):'                         
       do i=1, num_crack_fronts                                                 
         write (out,'(5x,i7,":",20i6)') i, (master_lines(i,j),                  
     &              j = 1, 20)                                                  
       enddo                                                                    
      endif                                                                     
c                                                                               
         write (out,*)                                                          
c                                                                               
c              -- No crack growth                                               
c                                                                               
      else if (crack_growth_type.eq.0) then                                     
         write (out,*) '> no crack growth specified'                            
c                                                                               
c              -- Error                                                         
c                                                                               
      else                                                                      
         write (out,*) ' Incorrect crack_growth_type:',                         
     &        crack_growth_type                                                 
      endif                                                                     
c                                                                               
      write (out,*)                                                             
      write (out,*) '<<<<<<<< finished dumping crack growth info.'              
      write (out,*)                                                             
 9999 return                                                                    
c                                                                               
 9000 format(5(3x,4e14.6,/))                                                    
      end                                                                       
c                                                                           
c     ****************************************************************          
c     *                                                              *          
c     *               service function chk_killed                    *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 7/7/21 rhd                 *          
c     *                                                              *          
c     *     given an element, returns true if the                    *          
c     *     element has been or is in-progress of being killed       *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
      logical function chk_killed( elem ) result( t_chk_killed )                                      
c                                                                               
      use elem_extinct_data, only : dam_state                                   
      use damage_data, only : dam_ptr                                                          
c                                                                               
      implicit none
c
      integer, intent(in) :: elem
c
      integer :: elem_ptr                                         
c                                                                               
      t_chk_killed = .false.  
      if( .not. allocated( dam_ptr ) ) return                                                    
      elem_ptr =  dam_ptr(elem)    ! assume checking elsewhere                                            
      if( elem_ptr == 0 ) return ! element is not killable                                             
      t_chk_killed = dam_state( elem_ptr) > 0
c                                                                               
      return                                                                    
      end                                                                       
c
c     ****************************************************************
c     *                                                              *
c     *                    subroutine chk_killed_blk                 *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 3/28/21 rhd                *
c     *                                                              *
c     *     for the specified block of elements, return a logical    *
c     *     vector indicating if each element is killed or active.   *
c     *     "killed" can have two meanings based on type             *
c     *     also a flag if the whole block is killed                 *
c     *                                                              *
c     *     service routine for tanstf and do_nleps_block            *
c     *                                                              *
c     ****************************************************************
c
      subroutine chk_killed_blk( span, felem, step, block_no, 
     &                           status_vec, block_killed, type )
c
      use global_data, only : out, mxvl 
      use elem_extinct_data, only : dam_state, smcs_d_values
      use damage_data, only : growth_by_kill, dam_ptr, 
     &                        use_mesh_regularization,
     &                        tol_regular =>
     &                        tolerance_mesh_regularization     
      use constants
c
      implicit none
c
c              parameters
c
      integer, intent(in)  :: span, block_no, felem, type, step
      logical, intent(out) :: status_vec(mxvl), block_killed
c      
c              locals
c
      integer :: elem, elem_ptr, i
      logical :: standard_kill_method
      logical, parameter :: ldebug = .false.
      double precision :: chk_value
c
      block_killed = .false.
      status_vec   = .false. ! all entries

      if( .not. growth_by_kill ) return
c 
      if( step <= 2 ) then
        if( .not. allocated(dam_state) ) then         ! sanity check                                
          write(out,9100) 1                                                         
          call die_abort                                                          
        end if 
      end if
c
      if( dam_ptr(felem) == 0 ) return ! no killable elems in blk
c
c               standard_kill_method = no mesh regularization
c               for standard, element is killed immediately. state > 0
c
c               for mesh regularization, have to check current value
c               of the damage parameter 'd'. if > 0, element is 
c               being killed (and d<1.0) or has been 
c               fully killed (d>=1.0)
c
c               type = 1 used in tanstf. mark element killed if
c                      d > 0.0. tanstf will use either the [Ke] saved
c                      at death degraded by 1.0 - d or set [Ke] =0
c                      if d >= 1.0
c               type = 2 used in drive_eps_sig_internal_forces.
c                      mark element killed only if d > 1.0.
c                      for d < 1.0, the element strains, stresses
c                      and internal forces are computed as usual.
c                      routine addifv will degrade them by 1 - d before
c                      assembly into structure. if d >= 1.0 we don't
c                      want to process strains, stresses, internal 
c                      forces for the element.
c 
c
c               run sanity check. every element in blk must be 
c               killable. 
c
      if( step <= 2 ) then  
         do  i = 1, span  
           elem = felem + i - 1
           elem_ptr = dam_ptr( elem )
           if( elem_ptr == 0 ) then ! bad data structures
             write(out,9100) 3                                                         
             call die_abort                                                          
           end if 
         end do 
      end if 
c
      standard_kill_method = .not. use_mesh_regularization
c
      if( standard_kill_method ) then
        do i = 1, span  
          elem = felem + i -1
          elem_ptr = dam_ptr( elem )
          status_vec(i) = dam_state(elem_ptr) > 0
        end do
        if( all( status_vec(1:span) ) ) block_killed = .true.
        return
      end if
c
      if( .not. use_mesh_regularization ) then
          write(out,9100) 4
          call die_abort
      end if
c
      if( .not. allocated( smcs_d_values ) ) then ! really bad
          write(out,9100) 5
          call die_abort
      end if
c
      if( type < 1 .or. type > 2 ) then
          write(out,9100) 6
          call die_abort
      end if
c 
      chk_value = zero
      if( type == 2 ) chk_value = tol_regular
c
      do i = 1, span  
        elem = felem + i - 1
        elem_ptr = dam_ptr( elem )
        status_vec(i) = smcs_d_values(elem_ptr) > chk_value
      end do
c
      if( all( status_vec(1:span) ) ) block_killed = .true.
c
      return       
c      
 9100 format(/,'FATAL ERROR: chk_killed_blk. Contact WARP3D group',             
     &       /,'             Job terminated at ',i1,//)                         
c
      end
c     ****************************************************************          
c     *                                                              *          
c     *                     subroutine damage_update                 *          
c     *                                                              *          
c     *                       written by : rhd                       *          
c     *                                                              *          
c     *                   last modified : 12/18/2020 rhd             *          
c     *                                                              *          
c     *        certain damage models need variables updated outside  *
c     *        element stress/history data (e.g. smcs )              *
c     *                                                              *          
c     ****************************************************************          
c                                                                               
c                                                                               
      subroutine damage_update( now_step, now_iter )
c
      use global_data, only : out, noelem
      use elem_extinct_data, only : dam_state
      use damage_data, only :  crack_growth_type, dam_ptr        
      use dam_param_code, only : dam_param_3_get_values 
c                                             
      implicit none
c
      integer :: now_step, now_iter
c                                       
      integer :: elem, elem_ptr, dowhat
      logical :: debug, local_debug, process, ld1                                                   
      double precision :: d1, d2, d3, d4, d5, d6, d7, d8, d9, d10
c        
      debug = .false.
      local_debug = .false. ! debug
      if( local_debug ) then      
         write(out,*) '... entered damage_update ...'                                                    
      end if                                                                    
c 
c              only damage model using SMCS requires this update
c              using data outside the stress/history data.
c
c              the Gurson model, for example, only needs the instantaneous
c              values stored in the element stress/history data to
c              determine if element needs to be killed, a step
c              size reduction is needed, values for printing, etc.
c      
      process = crack_growth_type .eq. 3
      if( .not. process ) return 
c
c              Loop for all elements. process killable elements.
c                                                                               
c$OMP PARALLEL DO PRIVATE( elem, elem_ptr, dowhat, d1, d2, d3, d4, 
c$OMP&                     d5, d6, d7, d8, d9, ld1, d10 )
      do elem = 1, noelem
         elem_ptr = dam_ptr(elem)                                               
         if( elem_ptr .eq. 0 ) cycle     ! element not killable                                      
         if( dam_state(elem_ptr) .gt. 0 ) cycle ! element killed already    
         if( local_debug ) write(out,*) 'killable element is:',elem                           
c                                                                               
c                 calculate/update damage parameter(s) for element
c                 if this damage model requires it. 
c
c                 only SMCS needs this now to support plastic strain
c                 weighted averages.                   
c                                                                               
         select case( crack_growth_type )
         case( 3 ) !  SMCS   
          dowhat = 4
          call dam_param_3_get_values( elem, debug, d1, d2, d3, d4,            
     &                                 d5, d6, d7, d8, dowhat, d9,
     &                                 ld1, d10 )                
         end select 
c                         
      end do  ! over all model elements                                         
c$OMP END PARALLEL DO
c      
      return    
      end
c     ****************************************************************          
c     *                                                              *          
c     *                      subroutine chk_store_ifv                *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 03/9/21 rhd                *          
c     *                                                              *          
c     *        setup to gradually reduce newly killed element        *
c     *        internal forces -- for standard death process         *
c     *        use_mesh_regularization -- nothing to do.             *
c     *                                                              *          
c     ****************************************************************          
c                                                                               
      subroutine chk_store_ifv( elem, elem_ptr, gdebug )  
c                           
      use global_data, only : out, iprops, load                                                                               
      use elem_extinct_data, only : dam_ifv, dam_state                          
      use elem_block_data,   only : edest_blocks                                
      use main_data,         only : elems_to_blocks                             
      use damage_data, only : release_type, use_mesh_regularization  
      use constants                                    
c                                                                               
      implicit none
c
      integer, intent(in) :: elem, elem_ptr
      logical, intent(in) :: gdebug
c
      integer :: num_edof, i, blk, rel_elem                                                     
      integer, dimension(:,:), pointer :: edest
      logical :: ldebug                                 
c     
      ldebug = .false.
c                                                                          
      if( ldebug ) write (*,*) '>>>>> Entering chk_store_ifv'      
c
      if( use_mesh_regularization ) return                 
c                                                                               
c            initialize state of release for element                            
c                                                                               
      dam_state(elem_ptr) = 1                                                   
c                                                                               
c            subtract the element internal forces from the                      
c            structure load vector. for all killable elements                   
c            which are not in the process of releasing, we                      
c            grab their internal forces at start of each load                   
c            step in another routine so they are available                      
c            here.                                                              
c                                                                               
      num_edof = iprops(2,elem) * iprops(4,elem)                                
c                                                                               
c            for the traction separation release, provide                       
c            an option to immediately                                           
c            drop the internal forces by a fraction, then release the           
c            remaining forces linearly with opening displacement.               
c            current reduction is zero.                                         
c                                                                               
      if( release_type .eq. 2 ) then                                           
        do i = 1, num_edof                                                      
            dam_ifv(i,elem_ptr) = one * dam_ifv(i,elem_ptr)                     
        end do                                                                  
      end if                                                                    
c                                                                               
      blk         = elems_to_blocks(elem,1)                                     
      rel_elem    = elems_to_blocks(elem,2)                                     
      edest       => edest_blocks(blk)%ptr                                      
      do i = 1, num_edof                                                        
        load(edest(i,rel_elem)) = load(edest(i,rel_elem)) -                     
     &                            dam_ifv(i,elem_ptr)                           
      end do                                                                    
c                                                                               
      if( ldebug ) then                                                         
         write (*,*) '>>>>> load in storeifv for elem:',elem                    
         do i = 1, num_edof                                                     
            write (*,'("new load:",i5,1x,e14.6)') edest(i,rel_elem),            
     &        load(edest(i,rel_elem))                                           
         end do                                                                 
         write (*,*) '>>>>>> Leaving chk_store_ifv'                                  
      end if                                                                    
c                                                                               
      return                                                                    
      end                                                                       
c     ****************************************************************          
c     *                                                              *          
c     *                  subroutine chk_kill_element_now             *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 9/19/21 rhd                *          
c     *                                                              *          
c     *        This routine actually kills an element in the model.  *          
c     *        Material properties and history data for the element  *          
c     *        are modified so that subsequent element [k]s will be  *          
c     *        zero. the existing stresses are zeroed. the props and *          
c     *        history changes are such that stresses will remain    *          
c     *        zero for all subsequent steps.                        *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
c                                                                               
c                                                                               
      subroutine chk_kill_element_now( ielem, debug ) 
c                                  
      use global_data, only : out, iprops, props, nstrs, nstr,
     &                             current_load_time_step
      use main_data, only : elems_to_blocks, cohesive_ele_types                                                              
      use elem_block_data, only : history_blocks, urcs_n_blocks,                
     &                            urcs_n1_blocks, history_blk_list,
     &                            eps_n_blocks, eps_n1_blocks 
      use damage_data, only : dam_ptr, crack_growth_type,
     &                        use_mesh_regularization,
     &                        use_distortion_metric,
     &                        smcs_removed_list_file_flag,
     &                        smcs_removed_list_file_name
      use elem_extinct_data, only : smcs_start_kill_step, Oddy_metrics
      use constants           
c
      implicit none   
c   
      integer, intent(in) :: ielem
      logical, intent(in) :: debug
c 
c              locals
c
      integer :: elem_type, material_type, ngp, blk, rel_elem, device,
     &           hist_size, hist_offset, gp, hisloc, start, end, i,
     &           se, ss, elem_ptr, npts
      double precision :: porosity
      double precision, dimension(:), pointer :: history, urcs_n1,
     &                                           urcs_n, eps_n, eps_n1
      logical :: gurson, regular_material, cohesive_elem, 
     &           output_fully_killed_msg, fexists
      character(len=80) :: file_name
c                                                                               
      if( debug ) write (out,*) '>>>> in chk_kill_element'                           
c                                                                               
c              for regular elements (hex, tet, et.c),                           
c              set young's modulus and poisson's ratio to zero for              
c              element so that any future request for a linear stiffness        
c              computes zero. for cohesive elements, we set a row of            
c              element props table = 1 so the element and cohesive model        
c              routines know element is killed.                                 
c                                                                               
      elem_type     = iprops(1,ielem)   
      material_type = iprops(25,ielem)      
      gurson        = material_type .eq. 3  ! also could just be mises                                      
      regular_material = .not. gurson
      cohesive_elem = cohesive_ele_types(elem_type)    
c                         
      if( .not. cohesive_elem ) then                                           
         props(7,ielem) = 0.0    !   E and nu                                                
         props(8,ielem) = 0.0                                                  
      else  ! cohesive elements                                                                    
         props(7:9,ielem) = 0.0
         props(20,ielem)  = 0.0
         props(21,ielem)  = 0.0
         iprops(32,ielem) = 1                                                   
      end if                                                                    
c                                                                              
c              for each gauss point of the element:                             
c               (a)  reset history data to a null state characteristic          
c                    of an unstressed element                                   
c               (b)  set stresses, streains at n and n+1 to zero.                        
c                                                                               
      ngp           = iprops(6,ielem)                                           
      blk           = elems_to_blocks(ielem,1)                                  
      rel_elem      = elems_to_blocks(ielem,2)                                  
      hist_size     = history_blk_list(blk)                                     
      hist_offset   = (rel_elem-1)*hist_size*ngp                                
      history       => history_blocks(blk)%ptr                                  
      urcs_n1       => urcs_n1_blocks(blk)%ptr                                  
      urcs_n        => urcs_n_blocks(blk)%ptr                                   
      urcs_n1       => urcs_n1_blocks(blk)%ptr                                  
      eps_n         => eps_n_blocks(blk)%ptr                                   
      eps_n1        => eps_n1_blocks(blk)%ptr                                   
c    
      if( gurson ) then ! don't zero porosity
         do gp = 1, ngp                                                            
           hisloc = hist_offset + (gp-1)*hist_size 
           porosity = history(hisloc+5)                          
           history(hisloc+1:hisloc+hist_size) = zero 
           history(hisloc+5) = porosity  
         end do                                         
      end if
c
      if( cohesive_elem ) then
         do gp = 1, ngp                                                            
            hisloc = hist_offset + (gp-1)*hist_size                                
            history(hisloc+1) = zero                                            
            history(hisloc+2) = zero                                            
            history(hisloc+3) = zero                                            
         end do                                                                    
      end if
c
      if( regular_material ) then
         do gp = 1, ngp                                                            
           hisloc = hist_offset + (gp-1)*hist_size 
           history(hisloc+1:hisloc+hist_size) = zero 
         end do
      end if
c                                                                               
      start = (rel_elem-1)*nstrs*ngp + 1                                    
      end   = start + ngp*nstrs - 1                                     
      do i = start, end                                                
        urcs_n1(i) = zero                                                      
        urcs_n(i)  = zero                                                      
      end do  
c                                                                  
      start = (rel_elem-1)*nstr*ngp + 1                                    
      end   = start + ngp*nstr - 1                                     
      do i = start, end                                                
        eps_n1(i) = zero                                                      
        eps_n(i)  = zero                                                      
      end do                                                                    
c 
      output_fully_killed_msg = .true.
      if( crack_growth_type .ne. 3 ) return
      if( .not. use_mesh_regularization ) return
      if( .not. output_fully_killed_msg ) return    
      elem_ptr = dam_ptr(ielem)                                               
      if( elem_ptr == 0 ) then
         write(out,9000) ielem
         call die_abort
      end if
c      
      if( .not. smcs_removed_list_file_flag ) return ! remove elements file
      inquire( file=smcs_removed_list_file_name, exist=fexists )
      open( newunit=device,file=smcs_removed_list_file_name,
     &      form='formatted', status='unknown', 
     &      access='sequential', position='append' )   
      if( .not. fexists ) write(device,9018)
      se = current_load_time_step
      ss = smcs_start_kill_step(elem_ptr)
      if( use_distortion_metric ) then
         write(device,9020) ielem, ss, se, se-ss, 
     &    Oddy_metrics(elem_ptr,2) / Oddy_metrics(elem_ptr,1)
      else
         write(device,9020) ielem, ss, se, se-ss
      end if   
      close(unit=device)    
c                                                                          
      if( debug ) write(out,*) '>>>> leaving chk_kill_element_now'                      
c                                                                               
      return
c
 9000 format('>> FATAL ERROR: routine ',
     & 'chk_kill_element_now',
     & 10x,'inconsistent condition. for element:',i8,
     & ' job terminated.'//)
 9018 format('!',/,
     & '!  SMCS steps numbers for element removal', /,
     & 6x,'elem  start step   end step # steps release',
     & '  max Oddy ratio')  
 9020 format(2x,i8,4x,i6,5x,i6,8x,i6,4x,f9.2)
 9100 format(10x,8f5.2)
c
      end                                                                       
c     ****************************************************************          
c     *                                                              *          
c     *                subroutine chk_update_killed_energy           *          
c     *                                                              *          
c     *                       written by : rhd                       *          
c     *                                                              *          
c     *                   last modified : 03/9/21 rhd                *          
c     *                                                              *          
c     *        An element is being killed. We need to save the       *          
c     *        internal and plastic work right now so it can be      *          
c     *        included in future work totals. this is because all   *          
c     *        such data is about to be zeroed.                      *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
c                                                                               
c                                                                               
      subroutine chk_update_killed_energy( ielem )  
c                                
      use global_data, only : out, iprops, nstrs, killed_ele_int_work,
     &                        killed_ele_pls_work
      use main_data, only : elems_to_blocks, cohesive_ele_types
      use elem_block_data, only : urcs_n_blocks, element_vol_blocks
      use constants
c            
      implicit none
c
      integer :: ielem
c
      integer :: elem_type, ngp, blk, rel_elem, sig_start, gp
      double precision :: sum_internal_energy, sum_plastic_energy,                       
     &                    element_volume, rgp                                                  
      double precision, pointer :: urcs_n(:) 
      logical :: cohesive_elem   
      logical, parameter :: debug = .false.                                           
c                                                                               
      if( debug ) write (out,*) '>>>> in update_killed_energy'                   
c                                                                               
c              for regular elements (hex, tet, et.c),                           
c              compute the average internal energy and plastic energy           
c              at the gauss points. these values are multipled by element       
c              volume and added to accumluated totals for killed                
c              elements.                                                        
c                                                                               
      elem_type     = iprops(1,ielem)                                           
      cohesive_elem = cohesive_ele_types(elem_type)                             
      if( cohesive_elem ) return                                               
      ngp           = iprops(6,ielem)                                           
      blk           = elems_to_blocks(ielem,1)                                  
      rel_elem      = elems_to_blocks(ielem,2)                                  
      urcs_n        => urcs_n_blocks(blk)%ptr                                   
      sig_start     = (rel_elem-1) * nstrs * ngp                                
c                                                                               
      sum_internal_energy = zero                                                
      sum_plastic_energy  = zero                                                
c                                                                               
      do gp = 1, ngp                                                            
        sum_internal_energy = sum_internal_energy + urcs_n(sig_start+7)         
        sum_plastic_energy  = sum_plastic_energy  + urcs_n(sig_start+8)         
        sig_start           = sig_start + nstrs                                 
      end do                                                                    
c                                                                               
      rgp                 = dble(ngp)                                           
      sum_internal_energy = sum_internal_energy / rgp                           
      sum_plastic_energy  = sum_plastic_energy  / rgp                           
c                                                                               
      element_volume      = element_vol_blocks(blk)%ptr(rel_elem)               
      killed_ele_int_work = killed_ele_int_work +                               
     &                      element_volume * sum_internal_energy                
      killed_ele_pls_work = killed_ele_pls_work +                               
     &                      element_volume * sum_plastic_energy                 
c                                                                               
      if ( debug ) write (out,*) '>>>> leaving update_killed_energy'              
c                                                                               
      return                                                                    
      end                                                                       
                                                                                
c     ****************************************************************          
c     *                                                              *          
c     *                   subroutine chk_update_node_elecnt          *          
c     *                                                              *          
c     *                       written by : ag                        *          
c     *                                                              *          
c     *                   last modified : 03/9/21 rhd                *          
c     *                                                              *          
c     *        Updates the vector that holds the number of           *          
c     *        elements connected to each node after a node is       *          
c     *        killed.                                               *          
c     *                                                              *          
c     ****************************************************************          
c                                                                               
c                                                                               
      subroutine chk_update_node_elecnt( ielem, debug ) 
c                            
      use global_data, only : out, iprops
      use elem_extinct_data, only : dam_node_elecnt                             
      use main_data, only : incid, incmap, inverse_incidences                   
c                                                                               
      implicit none
c
      integer, intent(in) :: ielem
      logical, intent(in) :: debug                                                             
c 
      integer :: nnode, i, node, elem, ndof                                                                              
c
      if( debug ) write (out,*) '>>>> in update_node_elecnt '                    
c                                                                               
c                loop over nodes connected to killed element.  Get node         
c                number.  Then decrease the count of the 
c                elements connected to that node by one.                                           
c                
      nnode = iprops(2,ielem)                                                   
      do i = 1, nnode                                                           
         node = incid(incmap(ielem)+(i-1))                                      
         elem = inverse_incidences(node)%element_list(1)                        
         ndof = iprops(4,elem)                                                  
         dam_node_elecnt(node)= dam_node_elecnt(node) - 1                       
      end do                                                                    
      if( debug ) write (out,*) '>>>> leaving update_node_elecnt'                  
c                                                                               
      return                                                                    
 9000 format ('>> node ',i2,' of elem ',i7,' is ',i6)                           
 9010 format ('>>   elems connected to node is ',i2)                            
 9020 format('old constraints for node ',i7,' :',3(1x,e14.6)) 
c                  
      end                                                                       
c
c     ****************************************************************
c     *                                                              *
c     *                    subroutine chk_get_d                      *
c     *                                                              *
c     *                   last modified : 4/30/21 rhd                *
c     *                                                              *
c     *     get current value of damage parameter "d" for mesh       *
c     *     regularization                                           *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine chk_get_d( elem, deps_plas, stress_at_death, d_now,
     &                      out )
      use damage_data, only : tol_regular =>
     &                        tolerance_mesh_regularization,
     &                        regular_npoints, regular_length,
     &                        regular_up_max, regular_points, ! 10x2
     &                        regular_type, regular_alpha,
     &                        regular_Gf, regular_m_power
      use constants
c
      implicit none
c
      integer, intent(in) :: out, elem
      double precision, intent(in) :: deps_plas, stress_at_death
      double precision, intent(out) :: d_now
      double precision, external :: chklint
c
      double precision :: u_plastic, u_normed, Gf, m, up, uf, T0,
     &                    T1, T2, Gtilde, Gratio, alpha
      logical :: local_debug
c
      local_debug = elem == -1
      u_plastic = regular_length * deps_plas
c
      select case( regular_type )
      case( 1 )
        u_normed =  u_plastic / regular_up_max
        d_now = chklint( u_normed, regular_npoints, 
     &                   regular_points(1,1), regular_points(1,2) )
        if( u_plastic >= regular_up_max ) d_now = one 
      case( 2 )
        if( stress_at_death < zero ) then
          write(out,9300) elem, stress_at_death
          call die_abort
        end if 
        Gf = regular_Gf
        m  = regular_m_power
        T0 = stress_at_death
        up = u_plastic
        uf = Gf * (m+one) / m / T0
        T1 = one - (up/uf)**m / (m + one)
        Gtilde = T0 * up * T1
        d_now = Gtilde / Gf
        if( up >= uf ) d_now = one ! happens w/ power law + large jumps
        if( local_debug ) write(out,9100) elem, u_plastic, regular_Gf,
     &                    m, uf, Gtilde, stress_at_death, d_now
      case( 3 )
        if( stress_at_death < zero ) then
          write(out,9300) elem, stress_at_death
          call die_abort
        end if 
        Gf     = regular_Gf
        m      = regular_m_power
        alpha  = regular_alpha
        T0     = stress_at_death
        up     = u_plastic
        uf     = Gf * (m+one) / m / T0
        T1     = one - (up/uf)**m / (m + one)
        Gtilde = T0 * up * T1
        Gratio = Gtilde / Gf
        T2     = one - exp( -alpha * Gratio )
        d_now  = T2 / ( one - exp(-alpha))
        if( up >= uf ) d_now = one ! happens w/ power law + large jumps
        if( local_debug ) write(out,9200) elem, u_plastic, regular_Gf,
     &                    m, alpha, uf, Gtilde, Gratio, stress_at_death,
     &                    d_now
      case default
         write(out,9000)
         call die_abort
      end select
c
      if( d_now < zero ) d_now = zero
      if( d_now > tol_regular ) d_now = one 
      if( local_debug ) write(out,9020) d_now
c
      return
 9000 format('>> FATAL ERROR: routine ',
     & 'chk_get_d.',
     & 10x,'inconsistent condition.',
     & ' job terminated.'//)
 9020 format(10x,'leaving chk_get_d w/ d_now: ',f6.3)
 9100 format(/,10x,'... compute d for element: ',i8,
     & /,15x,'u plastic:         ',f10.6,
     & /,15x,'Gf:                ',f10.6,
     & /,15x,'m-power:           ',f10.6,
     & /,15x,'computed uf:       ',f10.6,
     & /,15x,'computed G-tilde:  ',f10.5,
     & /,15x,'stress @ death:    ',f10.3,
     & /,15x,'d now:             ',f10.3 ) 
 9200 format(/,10x,'... compute d for element: ',i8,
     & /,15x,'u plastic:         ',f10.6,
     & /,15x,'Gf:                ',f10.6,
     & /,15x,'m-power:           ',f10.6,
     & /,15x,'alpha:             ',f10.6,
     & /,15x,'computed uf:       ',f10.6,
     & /,15x,'computed G-tilde:  ',f10.5,
     & /,15x,'Gratio:            ',f10.5,
     & /,15x,'stress @ death:    ',f10.3,
     & /,15x,'d now:             ',f10.3 ) 
 9300 format('>> FATAL ERROR: routine ',
     & 'chk_get_d. the stress @ element death < 0. this is an',
     & 10x,'inconsistent condition.',
     & ' job terminated.'//)
c
      end
c     ****************************************************************
c     *                                                              *
c     *         double precision function chklint                    *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified: 4/1/21 rhd                  *
c     *                                                              *
c     *    execute linear interpolation on a tabular function        *
c     *    of a single variable where the x values are sorted in     *
c     *    increasing value but are not req'd to be uniformly spaced *
c     *                                                              *
c     ****************************************************************
c
      function chklint( xvalue, n, x, y ) result( value )
      implicit none
c
      integer :: n
      double precision :: xvalue, x(n), y(n), value
c
      integer :: i, point
      double precision :: x1, x2, y1, y2
c
      if( xvalue .le. x(1) ) then
        value = y(1)
        return
      end if
c
      if( xvalue .ge. x(n) ) then
        value = y(n)
        return
      end if
c
      do point = 2, n
        if( xvalue .gt. x(point) ) cycle
        x1 = x(point-1)
        x2 = x(point)
        y1 = y(point-1)
        y2 = y(point)
        value = y1 + (xvalue-x1)*(y2-y1)/(x2-x1)
        return
      end do
c
      end
c
c     ****************************************************************
c     *                                                              *
c     *          logical function chkcrack_is_elem_deleted           *
c     *                                                              *
c     *                   last modified : 3/28/21 rhd                *
c     *                                                              *
c     *     has element been (fully) deleted from the solution?      *
c     *                                                              *
c     ****************************************************************
c
c
      logical function chkcrack_is_elem_deleted( elem ) result( test )
c
      use global_data, only : out
      use dam_param_code, only : dam_param
      use elem_extinct_data, only : dam_state, smcs_d_values
      use damage_data, only : dam_ptr,
     &                        use_mesh_regularization,
     &                        tol_regular =>
     &                        tolerance_mesh_regularization     
c
      implicit none
c
      integer, intent(in) :: elem
c
      integer :: elem_ptr
      logical :: std_kill  
c
      test = .false.
      elem_ptr = dam_ptr(elem)   ! assumes much checking before here
c                                            
      if( elem_ptr == 0 ) return ! element not killable        
      if( dam_state(elem_ptr) == 0 ) return ! not yet started killing
c
      std_kill = .not. use_mesh_regularization
      if( std_kill ) then
        test = .true.
        return ! dam_state /= 0 element actually deleted
      end if
c
      if( use_mesh_regularization ) then
        test = smcs_d_values(elem_ptr) > tol_regular ! fully deleted 
        return
      end if
c
      write(out,9000) elem
      call die_abort
c
 9000 format('>> FATAL ERROR: routine ',
     & 'chkcrack_is_elem_deleted.',
     & 10x,'inconsistent condition. for element:',i8,
     & ' job terminated.'//)
c
      end function chkcrack_is_elem_deleted
c
c     ****************************************************************
c     *                                                              *
c     *          subroutine chkcrack_Oddy_needed                     *
c     *                                                              *
c     *            last modified : 8/29/21 rhd                       *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine chkcrack_Oddy_needed( felem, lflag )
c
      use global_data, only : out
      use elem_extinct_data, only : dam_state
      use damage_data, only : dam_ptr, crack_growth_type, 
     &                        use_distortion_metric
      use constants
c
      implicit none
c
      integer, intent(in) :: felem
      logical, intent(out) :: lflag
c
      integer :: elem, elem_ptr
      logical :: process
c
      lflag = .false.
c
      if( .not. use_distortion_metric ) return
      process = crack_growth_type .eq. 1  .or. crack_growth_type .eq. 3
      if( .not. process ) return  ! no element deletion growth
      elem = felem
      elem_ptr = dam_ptr(elem)                                               
      if( elem_ptr == 0 ) return ! no elements in block killable
      lflag = .true. ! killable element in block
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *          subroutine chkcrack_Oddy_initial_print              *
c     *                                                              *
c     *            last modified : 9/10/21 rhd                       *
c     *                                                              *
c     *     after step = 1, iter = 1, print a file of Oddy metrics   *
c     *     if requested by the user. file name is fixed             *
c     *                                                              *
c     ****************************************************************
c
      subroutine chkcrack_Oddy_initial_print
c
      use global_data, only : out, noelem
      use damage_data, only :  dam_ptr, Oddy_print_initial
      use elem_extinct_data, only : oddy_metrics_initial
      use constants
c
      integer :: elem, elem_ptr, fileno, element_min, element_max,
     &           element_max_ratio
      logical, save :: already_printed = .false.
      real :: Oddy_min, Oddy_max, global_min, global_max, 
     &        global_max_ratio
      character(len=80) :: fname
c
      if( .not. Oddy_print_initial ) return
      if( already_printed ) return ! only print file once
      already_printed = .true.
c
      if( .not. allocated( dam_ptr ) ) then
         write(out,9900) 1
         call die_abort
      end if
      if( .not. allocated( Oddy_metrics_initial ) ) then
         write(out,9900) 2
         call die_abort
      end if
c          
      fname = "Oddy_initial_metrics_killable_elements"
      open(newunit=fileno,file=fname,status="unknown")
      write(fileno,9010) 
c
      global_min = 1000000.0
      global_max = -100000.0
      global_max_ratio = -100000.0
c
      do elem = 1, noelem
       elem_ptr = dam_ptr(elem)                                               
       if( elem_ptr == 0 ) cycle 
       Oddy_min = Oddy_metrics_initial(elem_ptr,1)
       Oddy_max = Oddy_metrics_initial(elem_ptr,2)
       if( Oddy_min < global_min ) then
         global_min = Oddy_min
         element_min = elem
       end if
       if( Oddy_max > global_max ) then
         global_max = Oddy_max
         element_max = elem
       end if
       ele_ratio = Oddy_max / oddy_min
       if( ele_ratio > global_max_ratio ) then
         global_max_ratio = ele_ratio
         element_max_ratio = elem
       end if
       write(fileno,9000) elem, Oddy_min, Oddy_max, ele_ratio
      end do
c
      write(fileno,9020) global_min, element_min, global_max,
     &                   element_max, global_max_ratio, 
     &                   element_max_ratio
      write(fileno,*) " "
c
      close(unit=fileno,status='keep')
      deallocate( Oddy_metrics_initial ) ! not saved in restart
c
      return
c
 9000 format(2x,i8,2x,f12.3,2x,f12.3,2x,f12.3)
 9010 format(/," Oddy metrics for killable elements @ t=0",
     & /,/,'    element      Oddy min      Oddy max      max/min')
 9020 format(/,5x,'minimum:       ',f12.3,' @ element: ',i8,/,
     &         5x,'maximum:       ',f12.3,' @ element: ',i8,/,
     &         5x,'maximum ratio: ',f12.3,' @ element: ',i8,/)
 9900 format('>> FATAL ERROR: @',i2,' in chkcrack_Oddy_initial_print',
     & /,    '                Job terminated',//)           
c
      end
c
c     ****************************************************************
c     *                                                              *
c     *          subroutine chkcrack_Oddy                            *
c     *                                                              *
c     *            last modified : 9/10/21 rhd                       *
c     *                                                              *
c     *     compute/store the Oddy distortion paramter for block     *
c     *     of killable elements. only called if calcs needed        *
c     *                                                              *
c     ****************************************************************
c
      subroutine chkcrack_Oddy( step, itype, span, felem, gpn, detJ,
     &                          jac )
c
      use global_data, only : out, mxvl, iprops, props
      use dam_param_code, only : dam_param
      use elem_extinct_data, only : dam_state, Oddy_metrics, 
     &    oddy_metrics_initial
      use damage_data, only : dam_ptr, Oddy_print_initial
      use constants
c
c              compute/store the current Oddy paramter for a block 
c              of killable elements at gpn.
c              allocate space during step 1 as needed. 
c              called by gtmat while processing blocks in parallel
c              to compute various required deformation gradients F.
c              gtmat1 comoutes coordinate jacobian (and det) for
c              current element (deformed) shape.
c              this routines runs on each thread.
c
      implicit none
c
      integer, intent(in) :: step, itype, span, felem, gpn
      double precision, intent(in) :: detJ(mxvl), jac(mxvl,3,3)
c
      integer :: i, elem, elem_ptr
      logical :: process
      logical, parameter :: ldebug1 = .false., ldebug2 = .false.,
     &                      ldebug3 = .false., print = .false.
      logical, save :: header = .true.
      double precision :: Oddy_value, Oddy_prior_value, Oddy_new
      real :: ymod, old_value, new_value, old_min, new_min, old_max,
     &        new_max
c
      if( ldebug1 ) write(out,9000) step, itype, felem, gpn, span
      if( .not. allocated( dam_ptr ) ) then
         write(out,8900)
         call die_abort
      end if
c
c              use special loop for step 1. undeformed Jac passed
c
      if( itype == 1 ) then
        do i = 1, span
          elem = felem + i - 1
          elem_ptr = dam_ptr(elem)                                               
          if( elem_ptr == 0 ) cycle ! really fatal error
          call chkcrack_Oddy_compute( i, Oddy_value )
          old_value = Oddy_metrics(elem_ptr,1)
          new_value = min( sngl(Oddy_value), old_value )
          if( new_value <= 0.0) new_value = 1.0 ! cubical elems
          Oddy_metrics(elem_ptr,1) = new_value
          Oddy_metrics(elem_ptr,2) = 0.0
          if( ldebug1 .and. elem == 1 ) then
              write(out,9010) elem, detJ(i), Oddy_value, new_value
              write(out,9012) jac(i,1:3,1:3)
          end if
          if( Oddy_print_initial ) then
            old_min = Oddy_metrics_initial(elem_ptr,1)
            new_min = min( old_min, sngl(Oddy_value) )
            old_max = Oddy_metrics_initial(elem_ptr,2)
            new_max = max( old_max, sngl(Oddy_value) )
            Oddy_metrics_initial(elem_ptr,1) = new_min
            Oddy_metrics_initial(elem_ptr,2) = new_max
          end if
          if( .not. print ) cycle
          if( header ) write(out,9310)
          header = .false.
          write(out,9320) elem, gpn, new_min, new_max
        end do
        return
      end if
c
c              regular loop after step 1. all data allocated. put 
c              current Oddy value in row 2 of table. only need to 
c              process elements being killed and not yet fully 
c              deleted from model
c
      if( itype /= 2 ) then
        write(out,9300) 
        call die_abort
      end if
c        
      do i = 1, span
        elem = felem + i - 1
        ymod = props(7,elem)  ! has been zeroed for fully killed element
        if( ymod <= 0.0 ) cycle  ! props is single precision
        elem_ptr = dam_ptr(elem) 
        call chkcrack_Oddy_compute( i, Oddy_value )
        Oddy_prior_value = dble( Oddy_metrics(elem_ptr,2) )
        if( Oddy_value <= zero ) Oddy_value = one
        Oddy_new = max( Oddy_prior_value, Oddy_value )
        Oddy_metrics(elem_ptr,2) = sngl(Oddy_new) 
        if( step <4 .and. ldebug3 .and. elem == 1 ) then
            write(out,9200) elem, gpn, Oddy_prior_value, Oddy_value, 
     &                      Oddy_new
            write(out,9012) jac(i,1:3,1:3)
        end if    
      end do
c
      return 
c
 8900 format(/1x,'>>>>> FATAL ERROR: routine chkcrack_Oddy.'
     &  /,'                    should not be called',
     &  /,'        Execution terminated.',/)
 9000 format(/,'.... Enter chkcrack_Oddy:',/,
     &  10x,'step, itype, felem, gpn, span: ',i7,i3,i7,i3,i4) 
 9010 format(10x,'elem, detJ, Oddy value, Oddy new:',i8,3f10.4)
 9012 format(15x,3f10.5)
 9100 format(10x, i8, 3f10.4 )
 9200 format('.... Oddy trace. elem, gp, Oddy prior, Oddy now, ',
     &   ' Oddy new: ', i8,i4,3f10.4) 
 9300 format(/,">>>>>  FATAL ERROR: routine chkcrack_Oddy",/,
     &         "         Job terminated.")
 9310 format(5x,'--------   Oddy metric values undeformed element  ',
     &   '--------',/,
     &       15x,'Element     Gauss Pt     min Oddy value   ',
     &       'max Oddy value') 
 9320 format(15x, i7, 8x, i2, 9x, f8.3, 9x, f8.3)          
c
      contains
c     --------
c
      subroutine chkcrack_Oddy_compute( relem, Oddy_value )
      implicit none
      integer, intent(in) :: relem
      double precision, intent(out) :: Oddy_value
c
      double precision :: tens(3,3), trans_tens(3,3), C(3,3), 
     &                    Cvec(9), t1, t2
      equivalence ( C, Cvec )
c
      tens(1:3,1:3) = jac(relem,1:3,1:3) / detJ(relem)**third 
      trans_tens = transpose( tens )
      C = matmul( trans_tens, tens )
      t2 = ( C(1,1) + C(2,2) + C(3,3) )**2 / three
      t1 = dot_product( Cvec, Cvec ) ! tensor contraction C:C
      Oddy_value = one +  (t1 - t2)
      if( ldebug2 .and. elem == 1 ) then
         write(out,9000)
         write(out,9010) C
         write(out,9015) t1, t2
      end if
c
      return
 9000 format(10x,"... [C] ....:")
 9010 format(15x,3f10.5)
 9015 format(10x,"... t1,t2: ", 2f10.5)
c
      end subroutine chkcrack_Oddy_compute
      end subroutine chkcrack_Oddy
        



