text:
       SUBROUTINE BREAD (IG,INV,II3,NORIG,KG)
C
C      THIS ROUTINE IS USED ONLY IN BANDIT MODULE
C      IT READS THE CONNECTING ELEMENTS AND GENEL ELEM. FROM GEOM2 FILE
C      AND PREPROCESS THE MPC CARDS AND THE RIGID ELEMENTS FROM GEOM4
C
C      REVISED BY G.CHAN/UNISYS
C      12/89, TO INCLUDE NEW RIGID ELEMENTS CRROD, CRBAR, CRTRPLT,
C      CRBE1, CREB2, CRBE3 AND CRSPLINE
C      03/92, TO INCLUDE DUMMY ELEMENTS, CDUM1,...,CDUM9
C
      IMPLICIT INTEGER (A-Z)
      LOGICAL          DEBUG
      DIMENSION        SUB(2),   XXX(3),   IZ(3),    KG(7),    IG(1),
     1                 NORIG(1), INV(II3,1)
      CHARACTER        UFM*23,   UWM*25
      COMMON /XMSSG /  UFM,      UWM
      COMMON /BANDA /  IBUF1,    NOMPC
      COMMON /BANDB /  NBITIN,   KORE,     IFL,      NGRID,    IPNW(2),
     1                 KDIM
      COMMON /BANDD /  DUM6(6),  NEL,      NEQ,      NEQR
      COMMON /BANDS /  NN(10)
      COMMON /GEOMX /  GEOM1,    GEOM2,    GEOM4,    SCR1
      COMMON /NAMES /  RD,       RDREW,    WRT,      WRTREW,   REW
      COMMON /SYSTEM/  IBUF,     NOUT,     DUM43(43),KDUM(9)
      COMMON /GPTA1 /  NE,       LAST,     INCR,     KE(1)
      COMMON /ZZZZZZ/  Z(1)
      DATA             CRIGDR,   CRIGD1,   CRIGD2,   CRIGD3,   GENEL  /
     1                 8210,     5310,     5410,     8310,     4301   /
      DATA             CHBDY,    PLOTEL,   CRROD,    CRBAR,    CRTRPT /
     1                 4208,     5201,     6510,     6610,     6710   /
      DATA             CRBE1,    CRBE2,    CRBE3,    CRSPLN,   MSET   /
     1                 6810,     6910,     7010,     7110,     4HMSET /
      DATA             SUB,                MPC,      MAXMPC,   DEBUG  /
     1                 4HBREA,   4HD   ,   4901,     150,      .FALSE./
C
C
C     CHECK THE PRESENCE OF GEOM2 FILE
C
      KG(1) = GEOM2
      CALL RDTRL (KG(1))
      J = KG(2) + KG(3) + KG(4) + KG(5) + KG(6) + KG(7)
      IF (KG(1).LT.0 .OR. J.EQ.0) GO TO 370
      DO 10 I = 1,7
 10   KG(I) = 0
C
C     UPDATE /GPTA1/ IF DUMMY ELEMENTS ARE PRESENT
C
      DO 15 I = 1,9
      IF (KDUM(I) .EQ. 0) GO TO 15
      K = KDUM(I)/10000000
      L = (KDUM(I)-K*10000000)/10000
      J = (I+51)*INCR
      KE(J+ 6) = 2 + K + L
      KE(J+10) = K
 15   CONTINUE
C
C     CHECK THE PRESENCE OF MPC CARDS AND RIGID ELEMENTS.  SAVE THEIR
C     GRID DATA IN SCR1 FILE FOR TIGER AND UPDATE NEQ AND NEQR COUNTERS
C
      IF (NOMPC .EQ. 0) GO TO 200
      Z(1) = GEOM4
      CALL RDTRL (Z(1))
      J = 0
      DO 20 I = 2,7
 20   J = J + Z(I)
      IF (Z(1).LT.0 .OR. J.EQ.0) GO TO 200
C
      IBUF2 = IBUF1 - IBUF
      CALL OPEN (*290,SCR1,Z(IBUF2),WRTREW)
      IFILE  = GEOM4
      CALL PRELOC (*190,Z(IBUF1),GEOM4)
C
      IF (NOMPC .EQ. 1) GO TO 40
C
      XXX(1) = MPC
      XXX(2) = XXX(1)/100
      CALL LOCATE (*40,Z(IBUF1),XXX,J)
 25   J = 1
      CALL READ (*300,*40,GEOM4,IZ,1,0,M)
 30   J = J + 1
      CALL READ (*300,*40,GEOM4,KG(J),3,0,M)
      IF (KG(J) .NE. -1) IF (J+3-MAXMPC) 30,30,320
      J = J - 1
      KG(1) = J - 1
      CALL WRITE (SCR1,KG,J,1)
      NEQ = NEQ + 1
      GO TO 25
C
C     LOCATE ANY CRIGDR AND CRROD ELEMENTS, AND SAVE THE GRID DATA IN
C     SCR1. (DEPENDENT GRID FIRST, AND ONE INDEPENDENT GRID LAST)
C
C     FOR ALL RIGID ELEMENTS, THE FIRST WORD OF KG ARRAY CONTAINS
C     (NO. OF DEPENDENT + INDEP. GRIDS)*1000 + (NO. OF INDEP. GRIDS)
C     THE DATA IN SCR1 WILL BE PROCESSED BY TIGER
C
 40   IF (NOMPC .EQ. 3) GO TO 180
      XXX(1) = CRIGDR
 50   XXX(2) = XXX(1)/100
      CALL LOCATE (*60,Z(IBUF1),XXX,J)
 55   CALL READ (*300,*60,GEOM4,IZ,1,0,M)
      CALL READ (*300,*60,GEOM4,KG(3),3,0,M)
      KG(1) = 2*1000 + 1
      KG(2) = KG(4)
      CALL WRITE (SCR1,KG,3,1)
      NEQR  = NEQR + 1
      GO TO 55
C
 60   IF (XXX(1) .EQ. CRROD) GO TO 70
      XXX(1) = CRROD
      GO TO 50
C
C     LOCATE ANY CRIGD1, CRIGD2  AND CRBE2  ELEMENTS, AND SAVE GRID
C     DATA IN SCR1. PUT THE ONE INDEPENDENT GRID LAST
C
 70   XXX(1) = CRIGD1
 75   XXX(2) = XXX(1)/100
      CALL LOCATE (*90,Z(IBUF1),XXX,J)
 80   J = 1
      CALL READ (*300,*90,GEOM4,IZ,2,0,M)
      IZ2 = IZ(2)
 85   J = J + 1
      CALL READ (*300,*90,GEOM4,KG(J),1,0,M)
      CALL READ (*300,*90,GEOM4,   0,-6,0,M)
      IF (KG(J) .NE. -1) IF (J-MAXMPC) 85,85,320
      KG(J) = IZ2
      KG(1) = (J-1)*1000 + 1
      CALL WRITE (SCR1,KG,J,1)
      NEQR  = NEQR + 1
      GO TO 80
 90   IF (XXX(1) .EQ. CRBE2) GO TO 110
C
C     LOCATE ANY CRIGD2 ELEMENT
C
      IF (XXX(1) .EQ. CRIGD2) GO TO 100
      XXX(1) = CRIGD2
      GO TO 75
C
C     LOCATE ANY CRBE2 ELEMENT
C
 100  XXX(1) = CRBE2
      GO TO 75
C
C     LOCATE ANY CRIGD3, CRBE1, CRBAR AND CRTRPLT ELEMENTS, AND SAVE
C     GRID DATA IN SCR1 FILE. PUT THE INDEPENDENT GRID LAST
C
 110  XXX(1) = CRBAR
      ASSIGN 115 TO IRTN
      GO TO 130
 115  XXX(1) = CRTRPT
      ASSIGN 120 TO IRTN
      GO TO 130
 120  XXX(1) = CRBE1
      ASSIGN 125 TO IRTN
      GO TO 130
 125  XXX(1) = CRIGD3
      ASSIGN 150 TO IRTN
 130  XXX(2) = XXX(1)/100
      CALL LOCATE (*145,Z(IBUF1),XXX,J)
 133  J = 2
      K = 1
      CALL READ (*300,*145,GEOM4,IZ,1,0,M)
 135  CALL READ (*300,*145,GEOM4,IZ(K),1,0,M)
      IF (IZ(K) .EQ. MSET) GO TO 137
      CALL READ (*300,*145,GEOM4,0,-6,0,M)
      K = K + 1
      IF (K .GT. 999) GO TO 340
      GO TO 135
 137  CALL READ (*300,*145,GEOM4,KG(J),1,0,M)
      CALL READ (*300,*145,GEOM4,0,-6,0,M)
      IF (KG(J) .EQ. -1) GO TO 140
      J = J + 1
      IF (J .GT. MAXMPC) GO TO 320
      GO TO 137
 140  K = K - 1
      DO 142 I = 1,K
      KG(J) = IZ(I)
 142  J = J + 1
      J = J - 1
      KG(1) = (J-1)*1000 + K
      CALL WRITE (SCR1,KG,J,1)
      NEQR = NEQR + 1
      GO TO 133
C
C     LOCATE ANY CRSPLINE ELEMENTS, AND SAVE GRID DATA IN SCR1 FILE.
C     PUT THE INDEPENDENT GRIDS LAST
 145  GO TO IRTN, (115,120,125,150)
C
C     LOCATE ANY CRBE3 ELEMENTS, AND SAVE GRID DATA IN SCR1 FILE. PUT
C     THE INDEPENDENT GRID LAST
C
 150  XXX(1) = CRBE3
      XXX(2) = XXX(1)/100
      CALL LOCATE (*165,Z(IBUF1),XXX,J)
 151  CALL READ (*300,*165,GEOM4,IZ,3,0,M)
      IZ2 = IZ(2)
      J = 2
      CALL READ (*300,*165,GEOM4,0,-2,0,M)
 153  CALL READ (*300,*165,GEOM4,KG(J),1,0,M)
      K = -KG(J)
      IF (K .GT. 0) GO TO (155,157,160) K
      J = J + 1
      IF (J-MAXMPC) 153,153,320
 155  CALL READ (*300,*165,GEOM4,I,1,0,M)
      IF (I .EQ. -2) GO TO 157
      CALL READ (*300,*165,GEOM4,0,-1,0,M)
      GO TO 153
 157  CALL READ (*300,*165,GEOM4,KG(J),1,0,M)
      IF (KG(J) .LT. 0) GO TO 160
      CALL READ (*300,*165,GEOM4,0,-1,0,M)
      J = J + 1
      GO TO 157
 160  KG(J) = IZ2
      KG(1) = (J-1)*1000 + 1
      CALL WRITE (SCR1,KG,J,1)
      NEQR = NEQR + 1
      GO TO 151
C
C     LOCATE ANY CRSPLINE ELEMENTS, AND SAVE GRID DATA IN SCR1 FILE.
C     PUT THE INDEPENDENT GRIDS LAST
C
 165  XXX(1) = CRSPLN
      XXX(2) = XXX(1)/100
      CALL LOCATE (*180,Z(IBUF1),XXX,J)
 167  CALL READ (*300,*180,GEOM4,IZ,3,0,M)
      K = 1
      IZ(K) = IZ(3)
      J = 1
 170  J = J + 1
 173  CALL READ (*300,*175,GEOM4,KG(J),2,0,M)
      IF (KG(J) .EQ. -1) GO TO 175
      IF (J+2 .GT. MAXMPC) GO TO 320
      IF (KG(J+1) .NE.  0) GO TO 170
      K = K + 1
      IF (K .GT. 999) GO TO 340
      IZ(K) = KG(J)
      GO TO 173
 175  DO 177 I = 1,K
      KG(J) = IZ(I)
 177  J = J + 1
      J = J - 1
      KG(1) = (J-1)*1000 + K
      CALL WRITE (SCR1,KG,J,1)
      NEQR = NEQR + 1
      GO TO 167
C
 180  DO 185 K = 1,MAXMPC
 185  KG(K) = 0
 190  CALL CLOSE (GEOM4,REW)
      CALL CLOSE (SCR1,REW)
C
C     PROCESS ELEMENT CARDS AND FILL UP CONNECTION TABLE IG
C
 200  IFILE = GEOM2
      CALL PRELOC (*300,Z(IBUF1),GEOM2)
      IELEM = 1 - INCR
 205  IELEM = IELEM + INCR
      IF (IELEM .GT. LAST) GO TO 250
      IF (KE(IELEM+3) .EQ. CHBDY ) GO TO 205
      IF (KE(IELEM+3) .EQ. PLOTEL) GO TO 205
      SCALAR = KE(IELEM+10)
      IF (SCALAR .EQ. -1) GO TO 205
      CALL LOCATE (*205,Z(IBUF1),KE(IELEM+3),J)
      NWDS  = KE(IELEM+ 5)
      NGPTS = KE(IELEM+ 9)
      NGPT1 = KE(IELEM+12)
      NCON  = NGPTS
 210  CALL READ (*300,*205,GEOM2,KG(1),NWDS,0,M)
      IF (SCALAR .EQ. 0) GO TO 220
      IF (KG(5).EQ.0 .OR. KG(6).EQ.0) GO TO 210
C     THE ABOVE CONDITIONS HOLD TRUE FOR CDAMPI, CELASI, AND CMASSI
C     WHERE I = 1,2
 220  NEL = NEL + 1
      CALL SCAT (KG(NGPT1),NCON,INV,II3,NORIG)
      IF (NGRID .EQ. -1) GO TO 270
      IF (NCON  .LE.  1) GO TO 240
      NGPT2 = NGPT1 + NCON - 1
      K = NGPT2 - 1
      DO 230 I = NGPT1,K
      L = I + 1
      DO 230 J = L,NGPT2
 230  CALL SETIG (KG(I),KG(J),IG,NORIG)
 240  IF (IELEM-LAST) 210,210,255
C
C     SPECIAL TREATMENT FOR GENERAL ELEM.
C     (LIMITED TO KDIM*4 GRID POINTS PER GENEL)
C
 250  XXX(1) = GENEL
      XXX(2) = XXX(1)/100
      CALL LOCATE (*270,Z(IBUF1),XXX,J)
      KDIM4  = KDIM*4
 255  NTOT   = 0
      CALL READ (*300,*270,GEOM2,K,1,0,M)
      K    = 0
      KGPV = 0
      GO TO 263
 260  IF (KG(NCON) .EQ. KGPV) GO TO 265
      KGPV = KG(NCON)
 263  NTOT = NTOT + 1
      IF (NTOT .LT. KDIM4) NCON = NTOT
 265  CALL READ (*300,*270,GEOM2,KG(NCON),2,0,M)
      IF (KG(NCON) .NE. -1) IF (KG(NCON+1)) 260,  265,  260
C                                           GRD  SCALAR GRD
C                                           PT.   PT.   PT.
      K = K + 1
      XXX(K) = KG(NCON+1)
      IF (K .LT. 2) GO TO 265
      NCON = NCON - 1
      M    = XXX(1)
      NWDS = 1 + (M*M-M)/2 + M
      CALL READ (*300,*270,GEOM2,K,-NWDS,0,M)
      CALL READ (*300,*270,GEOM2,K,   1,0,M)
      NGPT1 = 1
      IF (K .EQ. 0) GO TO 220
      NWDS  = M*XXX(2)
      CALL READ (*300,*270,GEOM2,K,-NWDS,0,M)
      GO TO 220
 270  CALL CLOSE (GEOM2,REW)
      IF (NTOT .GT. KDIM4) GO TO 330
      IF (.NOT.DEBUG) RETURN
C
      M = NN(1)
      WRITE  (NOUT,280) NN
      WRITE  (NOUT,285) ((INV(I,J),J=1,2),I=1,M)
 280  FORMAT (//21H /BANDS/ FROM BREAD =,10I8)
 285  FORMAT (/12H TABLE INV =,(/10X,2I8))
      RETURN
C
 290  IFILE = SCR1
 300  CALL MESAGE (-1,IFILE,SUB)
 320  WRITE  (NOUT,325) UWM,IZ(1),MAXMPC
 325  FORMAT (A25,', MPC SET (OR CRIGID ID)',I9,
     1        ' IS TOO LONG,  ONLY THE FIRST',I4, /5X,
     2        ' GRID POINTS ARE USED IN THE BANDIT COMPUTATION')
      GO TO 180
 330  WRITE  (NOUT,335) UFM,NTOT
 335  FORMAT (A23,', GENEL ELEMENT HAS TOO MANY GRID POINTS,',I7)
      J = NTOT/400 + 1
      IF (J .LE. 9) WRITE (NOUT,336) J
 336  FORMAT (5X,'USER NEEDS TO ADD A ''NASTRAN BANDTDIM=',I1,
     1       ''' CARD AND RERUN JOB')
      GO TO 350
 340  WRITE  (NOUT,345)
 345  FORMAT ('0*** MORE THAN 1000 INDEPENDENT GRID POINTS USED IN A ',
     1        'RIGID ELEMENT')
 350  CALL MESAGE (-61,0,0)
C
 370  NGRID = 0
      RETURN
      END


      FUNCTION BESSY (N,X)
! ------------------------------------------------------------------
!     This subroutine calculates the second kind Bessel Function of
!     integer order N, for any real X. We use here the classical
!     recursive formula.
! ------------------------------------------------------------------
      IMPLICIT NONE
      REAL *8 X,BESSY,BESSY0,BESSY1,TOX,BY,BYM,BYP
      INTEGER N,J
      IF (N.EQ.0) THEN
      BESSY = BESSY0(X)
      RETURN
      ENDIF
      IF (N.EQ.1) THEN
      BESSY = BESSY1(X)
      RETURN
      ENDIF
      IF (X.EQ.0.) THEN
      BESSY = -1.E30
      RETURN
      ENDIF
      TOX = 2./X
      BY  = BESSY1(X)
      BYM = BESSY0(X)
      DO 11 J = 1,N-1
      BYP = J*TOX*BY-BYM
      BYM = BY
      BY  = BYP
   11 CONTINUE
      BESSY = BY
      RETURN
      END


! { dg-do run }
! { dg-options "-fdump-tree-original" }
!
! Test the fix for PR66079. The original problem was with the first
! allocate statement. The rest of this testcase fixes problems found
! whilst working on it!
!
! Reported by Damian Rouson  <damian@sourceryinstitute.org>
!
  type subdata
    integer, allocatable :: b
  endtype
!  block
    call newRealVec
!  end block
contains
  subroutine newRealVec
    type(subdata), allocatable :: d, e, f
    character(:), allocatable :: g, h, i
    character(8), allocatable :: j
    allocate(d,source=subdata(1)) ! memory was lost, now OK
    allocate(e,source=d) ! OK
    allocate(f,source=create (99)) ! memory was lost, now OK
    if (d%b .ne. 1) STOP 1
    if (e%b .ne. 1) STOP 2
    if (f%b .ne. 99) STOP 3
    allocate (g, source = greeting1("good day"))
    if (g .ne. "good day") STOP 4
    allocate (h, source = greeting2("hello"))
    if (h .ne. "hello") STOP 5
    allocate (i, source = greeting3("hiya!"))
    if (i .ne. "hiya!") STOP 6
    call greeting4 (j, "Goodbye ") ! Test that dummy arguments are OK
    if (j .ne. "Goodbye ") STOP 7
  end subroutine

  function create (arg) result(res)
    integer :: arg
    type(subdata), allocatable :: res, res1
    allocate(res, res1, source = subdata(arg))
  end function

  function greeting1 (arg) result(res) ! memory was lost, now OK
    character(*) :: arg
    Character(:), allocatable :: res
    allocate(res, source = arg)
  end function

  function greeting2 (arg) result(res)
    character(5) :: arg
    Character(:), allocatable :: res
    allocate(res, source = arg)
  end function

  function greeting3 (arg) result(res)
    character(5) :: arg
    Character(5), allocatable :: res, res1
    allocate(res, res1, source = arg) ! Caused an ICE
    if (res1 .ne. res) STOP 8
  end function

  subroutine greeting4 (res, arg)
    character(8), intent(in) :: arg
    Character(8), allocatable, intent(out) :: res
    allocate(res, source = arg) ! Caused an ICE
  end subroutine
end
! { dg-final { scan-tree-dump-times "builtin_malloc" 20 "original" } }
! { dg-final { scan-tree-dump-times "builtin_free" 21 "original" } }


text:
c ********************************************************************
c *                                                                  *
c *  routines to support crack growth by critical ctoa criterion     *
c *                                                                  *
c ********************************************************************
c
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine dam_print_node               *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *     This routine prints out the status of the crack front    *
c     *     nodes for a crack growth by node release model at the    *
c     *     beginning of a load step.                                *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine dam_print_node( step, iter )
c
      use global_data, only : out, dstmap
c
      use node_release_data, only : crack_front_nodes, inv_crkpln_nodes,
     &     num_neighbors, neighbor_nodes, crkpln_nodes_state
      use main_data, only : cnstrn, cnstrn_in, output_packets,
     &                      packet_file_no
      use damage_data
      use constants
c
      implicit none
c
      integer :: step, iter
c
      integer, parameter :: max_local_list=300
      integer :: i, use_init, pointer, max_node,
     &           max_neighbor, node, node_data_entry, neighbor,
     &           neighbor_node, num_local_list
      integer :: local_node_list(max_local_list,3)
      double precision :: angle, max_angle, denom_angle,
     &                    local_angles_list(max_local_list)
      character(len=1), parameter :: star = "*"
      logical :: first_entry
c
c           print the header
c
      write(out,*) ' '
      write(out,*) ' ********************************************* '
      write(out,*) ' ***         Crack Front Status            *** '
      write(out,*) ' ********************************************* '
      write(out,*)
c
      write(out,*) 'crack front node  neighbor  CTOA(degrees)',
     &     '  CTOA/critical CTOA'
      write(out,*) '----------------  --------  -------------',
     &     '  ------------------'
c
c          set up for parsing the crack front node list
c
      pointer = crack_front_start
      if( pointer .eq. -1 ) then
         write(out,*) '>>>>> No crack front nodes'
         go to 9999
      end if
      max_angle      = zero
      max_node       = 0
      max_neighbor   = 0
      num_local_list = 0
c
c          check next crack front node
c
 10   continue
      first_entry     = .true.
      node            = crack_front_nodes(pointer,1)
      node_data_entry = inv_crkpln_nodes(node)
c
      do neighbor = 1, num_neighbors(node_data_entry)
c
c                find out if current neighbor is constrained -- if so,
c                skip to next node.  if not, must compute the angle
c
         neighbor_node = neighbor_nodes(neighbor,node_data_entry)
         if( cnstrn(dstmap(neighbor_node)+crk_pln_normal_idx-1)
     &        .ne. d32460 ) cycle
c
c               neighbor node is unconstrained.  call routine to get
c               the angle between it and the original node. Print info.
c
c               note that the angle returned is between the line and
c               the crack plane -- to make it relate to the CTOA, we
c               must multiply it by two.
c
         call get_slope( neighbor_node, node, crk_pln_normal_idx,
     &                   angle )
         angle = angle * two
         if( angle .gt. max_angle ) then
            max_angle    = angle
            max_node     = node
            max_neighbor = neighbor_node
         end if
c
         num_local_list = num_local_list + 1
         if( num_local_list .gt. max_local_list ) then
           write(out,9030)
           call die_gracefully
         end if
c
         local_node_list(num_local_list,1) = node
         local_node_list(num_local_list,2) = neighbor_node
         local_node_list(num_local_list,3) = 0
         local_angles_list(num_local_list) = angle

         if( crkpln_nodes_state(inv_crkpln_nodes(neighbor_node) )
     &        .eq. 0 ) then
            if( first_entry ) then
               write(out,"(6x,i5,9x,i5,7x,f7.2,a1,11x,f4.1,3x,'<=')")
     &              node, neighbor_node,
     &              angle, star, angle/(init_crit_ang*two)
               first_entry = .false.
            else
               write(out,"(20x,i5,7x,f7.2,a1,11x,f4.1,3x,'<=')")
     &              neighbor_node,
     &              angle, star, angle/(init_crit_ang*two)
            end if
         else
            if( first_entry ) then
               write(out,"(6x,i5,9x,i5,7x,f7.2,12x,f4.1,3x,'<=')")
     &              node, neighbor_node,
     &              angle,angle/(critical_angle*two)
               first_entry = .false.
               local_node_list(num_local_list,3) = 1
            else
               write(out,"(20x,i5,7x,f7.2,12x,f4.1,3x,'<=')")
     &              neighbor_node,
     &              angle, angle/(critical_angle*two)
               local_node_list(num_local_list,3) = 1
            end if
         end if
c
      end do
c
c         now go to next entry in list
c
      pointer = crack_front_nodes(pointer,2)
      if( pointer .ne. -1 ) go to 10
c
      write(out,9000) max_angle, max_node, max_neighbor
      write(out,9050) num_ctoa_released_nodes

c
      if( output_packets .and. num_local_list .gt. 0  ) then
        write(packet_file_no) 9, num_local_list, step, iter
        do i = 1, num_local_list
          node        = local_node_list(i,1)
          neighbor    = local_node_list(i,2)
          use_init    = local_node_list(i,3)
          angle       = local_angles_list(i)
          denom_angle = critical_angle*two
          if ( use_init .eq. 0 ) denom_angle = init_crit_ang*two
          write(packet_file_no) node, neighbor, angle, use_init,
     &                          angle/(denom_angle)
        end do
      end if
c
 9999 continue
      write(out,*)
      return
c
 9000 format (/1x,'** Max CTOA: ',f6.2,' degrees',
     &        ' @ crack front node: ',i7,
     &        ' connected to node: ',i7,
     &      /,1x,'   NOTE: * indicates critical CTOA is for',
     &        ' initiation',
     &        / )
 9030 format(/,'FATAL ERROR: list length exceeded in',
     & ' dam_print_node. job aborted.' )
 9050 format(8x,'** Total nodes released in model: ',i5,' **')
c
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine dam_print_front              *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *     This routine prints out the status of the crack front    *
c     *     nodes for a crack growth by node release model at the    *
c     *     beginning of a load step if it uses a constant front.    *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine dam_print_front( step, iter )
c
      use global_data, only : out
      use node_release_data, only : master_lines
      use main_data, only : output_packets, packet_file_no
      use damage_data
      use constants
c
      implicit none
c
      integer :: step, iter
c
      integer, parameter ::  max_local_list = 200
      double precision :: angle, max_angle,
     &                    local_angles_list(max_local_list)
      integer :: i, local_master_list(max_local_list), max_node,
     &           max_neighbor, num_local_list, num_line, node, idummy
      character(len=3) :: special_char
      logical :: use_init
c
c           print the header
c
      write(out,*) ' '
      write(out,*) ' ********************************************* '
      write(out,*) ' ***         Crack Front Status            *** '
      write(out,*) ' ********************************************* '
      write(out,*)
c
      write(out,*) 'master node       CTOA(degrees)',
     &     '  CTOA/critical CTOA'
      write(out,*) '-----------       -------------',
     &     '  ------------------'
c
c          loop over the entries in master_lines
c
      max_angle      = zero
      max_node       = 0
      max_neighbor   = 0
      num_local_list = 0
c
      do num_line = 1, num_crack_fronts
c
         call get_slope_master_line( num_line, use_init, angle,
     &                               idummy )
         if( angle .eq. zero ) cycle
         node  = master_lines( num_line, 1)
         angle = angle * two
         if( angle .gt. max_angle ) then
            max_angle = angle
            max_node = node
         end if
         special_char(1:3) = ' '
         if( use_init ) special_char(1:3) = '(*)'
         if( use_init ) then
          write(out,"(6x,i5,7x,f7.2,a3,11x,f4.1,3x,'<=')")
     &        node, angle, special_char, angle/(init_crit_ang*two)
         else
          write(out,"(6x,i5,7x,f7.2,a3,11x,f4.1,3x,'<=')")
     &        node, angle, special_char, angle/(critical_angle*two)
         end if
c
c          save into vector for packet output after loop
c
         num_local_list = num_local_list + 1
         if( num_local_list .gt. max_local_list ) then
           write(out,9030)
           call die_gracefully
         end if
         local_master_list(num_local_list) = node
         local_angles_list(num_local_list) = angle
c
      end do
c
      write(out,9000) max_angle, max_node
      write(out,9050) num_ctoa_released_nodes
      write(out,*)
c
      if( output_packets .and. num_local_list .gt. 0  ) then
        write(packet_file_no) 8, num_local_list, step, iter
        do i = 1, num_local_list
          node  = local_master_list(i)
          angle = local_angles_list(i)
          if( use_init ) then
           write(packet_file_no) node, angle, 0,
     &                           angle/(init_crit_ang*two)
          else
           write(packet_file_no) node, angle, 1,
     &                           angle/(critical_angle*two)
          end if
        end do
      end if

      return
 9000 format (/1x,'** Max CTOA: ',f6.2,' degrees',
     &        ' @ crack front node: ',i7,
     &      /,1x,'   NOTE: * indicates critical CTOA is for',
     &        ' initiation',
     &        / )
 9030 format(/,'FATAL ERROR: list length exceeded in',
     & ' dam_print_front. job aborted.' )
 9050 format(8x,'** Total nodes released in model: ',i5,' **')
c
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine chk_node_release             *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *        This routine updates previously released nodes and    *
c     *        checks the current crack front for crack growth       *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine chk_node_release( debug, step, iter )
c
      use global_data, only : out, growth_k_flag
      use node_release_data, only: crkpln_nodes_state
      use damage_data
c
      implicit none
c
      integer :: step, iter
      logical :: debug
c
      integer :: i
      logical :: killed_this_time
c
c
c     This routine:
c        1) updates information for previously released nodes by
c           - incrementing the state variable for each released node
c           - eliminiating any new constraints placed on previously
c              released degrees of freedom
c
c        2) finds new nodes for releaseing and executes release through
c            two passes of the linked list for the crack front:
c               Pass 1: compute CTOA, mark node for release if
c                         CTOA > critical angle and alter data structures
c                  for the released node
c               Pass 2: update crack front list to reflect newly released
c                       nodes
c
c
      if( debug ) write(out,*) '>>>> entering chk_node_release'
c
c
c        1) Update previously released nodes
c
c             loop through all crack plane nodes -- if node
c             has been released, increment state by one. This
c             keeps a record of the load steps since release
c             of the node, used by the load size reduction
c             mechanism.
c
c             Even though the exact state value is not used for
c             traction separation, we still increment the state
c             value since it shows if a node has been released
c             or not.
c
      if( debug ) write(out,*) '>>>>>> Checking released nodes '
c
      do i = 1, num_crack_plane_nodes
         if( crkpln_nodes_state(i) .gt. 0 )
     &        crkpln_nodes_state(i) = crkpln_nodes_state(i) + 1
      end do
c
c
c        2) Check nodes on crack front for release.  If any are released,
c           then update the crack front list to reflect the changes.
c           Note that if we are using constant front growth
c           we call a different check front routine.
c
      if( debug ) write(out,*) '>>>>>> Checking crack front nodes'
      if( const_front ) then
         call chk_crack_front_const( killed_this_time, debug,
     &                               step, iter )
      else
         call chk_crack_front( killed_this_time, debug, step, iter )
      end if
c
      if( killed_this_time ) growth_k_flag = .true.
c
      enforce_node_release = .false.
c
      if( debug ) write(out,*) '<<<< leaving chk_node_release'
c
      return
      end
c     ****************************************************************
c     *                                                              *
c     *                      subroutine chk_reconstraint             *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *        This routine checks the current constraints to see    *
c     *        if a previously released node has been re-constrained *
c     *        in the direction normal to the crack plane            *
c     *        by user-input constraints.  If a released node has    *
c     *        been re-constrained, unconstrain it.                  *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine chk_reconstraint
      use global_data, only : csthed, dstmap, cstmap
c
      use node_release_data, only : inv_crkpln_nodes,
     &     crkpln_nodes_state, crkpln_nodes_react
      use main_data, only : invdst, cnstrn_in, cnstrn
      use damage_data
      use constants
c
      implicit none
c
      integer :: cst_ptr, above_ptr, node, crkpln_ptr
      character(len=1) :: dums
      real :: dumr
      double precision :: dumd
c
c             traverse the singly-linked list of constraints in the model. If
c             the user has input new constraints, then some of the previously
c             released nodes may be re-constrained.  If a released node
c             has been re-constrained, remove the constraint.
c
c                set up pointers for constraint linked list
c
      cst_ptr   = csthed
      above_ptr = -1
c
c                enter top of constraint linked list loop.
c
 10   continue
      node       = invdst(cst_ptr)
      crkpln_ptr = inv_crkpln_nodes(node)
c
c                    if node is not on crack plane, or the dof is not in
c                    the crack plane direction, or the node is not yet
c                    killed, then go to next entry
c
c      if ( (crkpln_ptr .eq. 0) .or.
c     &     (cst_ptr .ne. dstmap(node) + crk_pln_normal_idx -1) .or.
c     &     (crkpln_nodes_state(crkpln_ptr) .eq. 0) ) go to 20
      if( crkpln_ptr .eq. 0 ) go to 20
      if( cst_ptr .ne. dstmap(node) + crk_pln_normal_idx -1 ) go to 20
      if( crkpln_nodes_state(crkpln_ptr) .eq. 0 ) go to 20
c
c                    the node has been re-constrained and must be
c                    re-released.  write a warning message,
c                    set the constraint value to null (32460),
c                    and remove it from the linked list.  Removing a
c                    constraint from linked list automatically
c                    updates the traversal pointer to point to the next
c                    entry
c
      call errmsg( 248, node, dums, dumr, dumd )
      cnstrn(cst_ptr)    = d32460
      cnstrn_in(cst_ptr) = d32460
c
      if( above_ptr .eq.-1 ) then
c
c                           at top of list.  move head pointer.
c
         csthed          = cstmap(cst_ptr)
         cstmap(cst_ptr) = 0
         cst_ptr         = csthed
      else
c
c                           in middle of list. correct pointers
c
         cstmap(above_ptr) = cstmap(cst_ptr)
         cstmap(cst_ptr)   = 0
         if( csttail .eq. cst_ptr ) csttail = above_ptr
         cst_ptr = cstmap(above_ptr)
      end if
      go to 30
c
c                    the constraint was not corrected.  we need to
c                    move to the next entry in the list.
c
 20   continue
      above_ptr = cst_ptr
      cst_ptr   = cstmap(cst_ptr)
c
c                    check to see if we are at end of list. if not, go to
c                    next value.
c
 30   continue
      if( cst_ptr .ne. -1 ) go to 10
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine chk_crack_front              *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/28/21 rhd                *
c     *                                                              *
c     *        This routine traverses the nodes on the current       *
c     *        crack front, releasing any nodes that have (1) CTOA   *
c     *        larger than the given critical angle, or (2) the      *
c     *        user has specifically requested a node release in     *
c     *        this step.                                            *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine chk_crack_front( killed_this_time, debug, step, iter )
c
      use global_data, only : out, dstmap
c
      use node_release_data, only : inv_crkpln_nodes, num_neighbors,
     &     neighbor_nodes, crack_front_nodes, crkpln_nodes_state
      use main_data, only : cnstrn, output_packets, packet_file_no
      use damage_data
      use constants
c
      implicit none
c
      logical :: killed_this_time, debug
      integer :: step, iter
c
      integer, parameter :: max_local_list=200
      integer :: local_killed_list(max_local_list)
      integer :: i, dof,  neighbor_node, neighbor,  node_data_entry,
     &           node, node_ptr, prev_node, num_killed_now,
     &           prev_node_ptr, node_causing_stop
      double precision :: angle, crit_angle,
     &                    local_killed_angles(max_local_list)
      logical :: kill, last_node_released, ldummy1, ldummy2, test
      integer, external :: scan_entry_in_list
c
c          traverse the linked list of crack front nodes and
c          check if nodes need to be released.
c
c             set pointer to top of list, and the previous node
c             pointer to null.  Also, set the flag that shows if
c             any nodes have been released this step to false.
c
      killed_this_time = .false.
      num_killed_now   = 0
      prev_node_ptr    = -1
      node_ptr         = crack_front_start
      if( node_ptr .eq. -1 ) return
c
c             for current node, find out if it is to be
c             released. loop over the neighbors, calculating
c             the opening angle (from the crack plane to the connected
c             element edges).
c
 10   continue
c
      node            = crack_front_nodes(node_ptr,1)
      node_data_entry = inv_crkpln_nodes(node)
      kill            = .false.
      if( debug ) then
        write(out,*) ' >> Checking node ', node
        write(out,*) '    > enforced release flag: ',
     &                   enforce_node_release
      end if
c
      do neighbor = 1, num_neighbors(node_data_entry)
c
c                find out if current neighbor is constrained -- if so,
c                skip to next node.  if not, compute the angle.  If
c                angle is larger than the critical angle (within a
c                range of CTOA_range of the critical angle), signal
c                that the node should be killed.  Note that we have two
c                critical angles -- one for initiation and one for
c                crack growth.  We use the initiation angle on any
c                neighboring nodes that are unconstrained but have a state
c                of 0, indicating that they were initially unconstrained.
c
         neighbor_node = neighbor_nodes(neighbor,node_data_entry)
         dof           = dstmap(neighbor_node)+crk_pln_normal_idx-1
         if( cnstrn(dof) .ne. d32460 ) cycle
c
         call get_slope( neighbor_node, node, crk_pln_normal_idx,
     &        angle )
         if( debug ) write(out,'("     Angle is: ",e13.6)') angle
c
         if( crkpln_nodes_state(inv_crkpln_nodes(neighbor_node))
     &        .eq. 0 ) then
            if( debug ) write(out,*) ' using initiation angle'
            crit_angle = init_crit_ang * ( one - CTOA_range )
         else
            if( debug ) write(out,*) ' using crack growth angle'
            crit_angle = critical_angle * ( one - CTOA_range )
         end if
c
         if ( angle .ge. crit_angle .or. enforce_node_release ) then
            kill = .true.
            exit
         end if
c
      end do
c
c           if node should be released, then print message and release it
c
c           note that the angle as computed in WARP3D is between the
c           crack plane and the appropriate element edge.  When we output
c           we need to multiply by two so we get the total opening angle
c           for the specimen.
c
c           the user may have provided a list if nodes.
c           if this newly released node is is the list,
c           stop solution at cleanup in this routine.
c
      if( kill ) then
         if( debug ) write(out,*) ' -> Crit.CTOA reached. Kill node.'
         if( .not. killed_this_time ) write(out,9000)
         write(out,9010) node, angle * two
         killed_this_time = .true.
         call release_node( node, node_data_entry, debug )
         num_killed_now = num_killed_now + 1
         if( num_killed_now .gt. max_local_list ) then
           write(out,9030)
           call die_gracefully
         end if
         local_killed_list(num_killed_now) = node
         local_killed_angles(num_killed_now) = angle * two
         if( stop_released_nlist_length > 0 ) then
           test = scan_entry_in_list( node, released_nlist_to_stop,
     &                             stop_released_nlist_length )
           if( test ) last_node_released = .true.
                node_causing_stop = node
         end if
      end if
c
c            move to the next entry in the list.  If we are not at the end,
c            process the next node.
c
      prev_node_ptr = node_ptr
      node_ptr      = crack_front_nodes(node_ptr,2)
      if( node_ptr .ne. -1 ) go to 10
c
c            print critical CTOA if any elements were killed.  Note that
c            again we multiply the angle by two before we output so that
c            we are outputting the total critical CTOA.  Also update the
c            crack front to refect killed elements.
c
c
      if( killed_this_time ) then
         if( enforce_node_release ) then
           write(out,9025)
         else
           write(out,9020) critical_angle * two,
     &          init_crit_ang * two, CTOA_range * hundred
         end if
         call update_crack_front( debug )
         write(out,9050)  num_ctoa_released_nodes
      end if
c
      if( output_packets .and. killed_this_time ) then
        write(packet_file_no) 18, num_killed_now+1, step, iter
        write(packet_file_no) critical_angle*two,
     &        init_crit_ang*two, CTOA_range*hundred
        do i = 1, num_killed_now
          write(packet_file_no) local_killed_list(i),
     &         local_killed_angles(i)
        end do
      end if
c
      if( last_node_released ) then
          write(out,9060) node_causing_stop
          call store( ' ','kill_limit_restart.db', ldummy1, ldummy2 )
          call warp3d_normal_stop
      end if

      if( debug ) write(out,*) '>> leaving chk_crack_front'
c
c
      return
 9000 format(/,' >> node release option invoked for the following:',/)
 9010 format(  '        node: ',i7,'  CTOA: ',f6.2,' (degrees)')
 9020 format(/,'        * critical CTOA:',
     &       /,'            for growth    : ',f6.2,' (degrees)',
     &       /,'            for initiation: ',f6.2,' (degrees)',
     &       /,'            nodes released within ',f4.1,' % of ',
     &                        'critical angle.',
     &       /)
 9025 format(/,'        *** user requested node release ***', / )
 9030 format(/,'FATAL ERROR: list length exceeded in chk_crack_front',
     & /,      '             job aborted.' )
 9050 format(8x,'** Total nodes released in model: ',i5,' **')
 9060 format(/,1x,">>>>> Limiting node released: ",i7,
     &       /,7x,"Restart file: kill_limit_restart.db written.",
     &       /,7x,"Job ended normally."// )
      end
c
c
c     ****************************************************************
c     *                                                              *
c     *                subroutine chk_crack_front_const              *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/28/21 rhd                *
c     *                                                              *
c     *        This routine checks the crack front for the constant  *
c     *        front algorithm, where the CTOA is measured a given   *
c     *        distance behind the crack front.  The algorithm just  *
c     *        cycles through the master_line lists, evaluating the  *
c     *        CTOA for each master node.                            *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine chk_crack_front_const( killed_this_time, debug, step,
     &                                  iter )
      use global_data, only : out, c
      use node_release_data, only : master_lines
      use main_data, only :  output_packets, packet_file_no, crdmap
      use damage_data
      use constants
c
      implicit none
c
      logical :: killed_this_time, debug
      integer :: step, iter
c
      integer, parameter :: max_local_list=200
      integer :: local_killed_list(max_local_list)
      integer :: i, orig_node, loop, num_line, num_killed_now, idummy,
     &           mnode, node_causing_stop
      double precision :: angle, crit_angle, dist, targ_dist, xc, yc,
     &                    zc, local_killed_angles(max_local_list)
      logical :: kill, use_init, last_node_released, ldummy1,
     &           ldummy2, test
      integer, external :: scan_entry_in_list
c
      if( debug ) write(out,*) '>>> in chk_front_master_line'
c
      killed_this_time = .false.
      num_killed_now   = 0
      last_node_released  = .false.
c
c          loop over the entries in master_lines
c
      do num_line = 1, num_crack_fronts
c
c             get CTOA for each master node
c
         if( debug ) write(out,'("    check front:",i7)') num_line
         call get_slope_master_line( num_line, use_init, angle,
     &        idummy )
         if( debug ) write(out,'("     Angle is: ",e13.6)') angle
c
         if( use_init ) then
            if( debug ) write(out,*) ' using initiation angle'
            crit_angle = init_crit_ang * ( one - CTOA_range )
         else
            if( debug ) write(out,*) ' using crack growth angle'
            crit_angle = critical_angle * ( one - CTOA_range )
         end if
c
c           if node should be released, then print message and release it
c
c           note that the angle as computed in WARP3D is between the
c           crack plane and the appropriate location on the crack face.
c           When we output we need to multiply by two so we get the total
c           opening angle for the specimen.
c
c           To eliminate mesh dependence of the CTOA results, the CTOA is
c           measured at a given distance behind the crack front.  This
c           allows an arbitrary number of elements between the crack tip and
c           the point at which the CTOA is measured.  However, the large
c           amount of deformation present at the crack tip may cause a number
c           of nodes to release consecutively after the first.  To avoid this
c           "spurt growth", we automatically release all the fronts ahead of
c           the first release until the crack has grown the distance used to
c           evaluate the CTOA.  This significantly improves results.
c
c
         if( angle .ge. crit_angle .or. enforce_node_release ) then
c
c              loop over nodes within critical distance
c
            if( debug ) write(out,*) ' ->Crit.CTOA reached.Kill node.'
            loop = 0
            orig_node = master_lines(num_line,1)
c
            do while (.true.)
c
c                  set the killed_this_time flag to indicate nodes have been
c                  released.
c
               loop = loop + 1
               if( .not. killed_this_time ) then
                  write(out,9000)
                  killed_this_time = .true.
               end if
c
c                  print the releasing master nodes
c
               num_killed_now = num_killed_now + 1
               if( num_killed_now .gt. max_local_list ) then
                  write(out,9030)
                  call die_gracefully
               end if
               local_killed_angles(num_killed_now) = angle * two
c
               mnode = master_lines(num_line,1)
c
c                 the user may have provided a list if nodes.
c                 if this newly released node is is the list,
c                 stop solution at cleanup in this routine.
c
               if( stop_released_nlist_length > 0 ) then
                test = scan_entry_in_list( mnode,
     &                  released_nlist_to_stop,
     &                  stop_released_nlist_length )
                if( test ) last_node_released = .true.
                node_causing_stop = mnode
               end if
               xc = c(crdmap(mnode)+0)
               yc = c(crdmap(mnode)+1)
               zc = c(crdmap(mnode)+2)
               if( loop .eq. 1 ) then
                  write(out,9010) mnode, angle*two, step-1, xc, yc, zc
                  local_killed_list(num_killed_now) = mnode
               else
                  write(out,9015) mnode, angle*two, step-1, xc, yc, zc
                  local_killed_list(num_killed_now) = mnode
               end if
c
c                  release the front related to the master node, then
c                  update the crack front.
c
               call release_front( master_lines(num_line,1), debug )
c
               call update_crack_front( debug )
c
c                  if the master_lines entry is now zero, then the
c                  crack has coalesed or grown to a free surface:
c                exit the release loop.
c
           if( master_lines(num_line,1) .eq. 0 ) exit
c
c                  now get the distance between current master node and
c                  the old master node which initiated the growth increment.
c                  If releasing the next increment would get us closer to
c                  the target distance, then continue releasing.  Otherwise
c                  stop releasing.
c
               call get_dist_to_node ( num_line, orig_node, dist )
c
               if( use_init ) then
                  targ_dist = init_ctoa_dist
               else
                  targ_dist = ctoa_dist
               end if
               if( abs(dist - targ_dist) .lt.
     &             abs(dist + char_length - targ_dist)) exit
c
            end do
c
         end if
c
      end do
c
c            print critical CTOA if any nodes were released.  Note that
c            again we multiply the angle by two before we output so that
c            we are outputting the total critical CTOA.  Also update the
c            crack front to refect killed elements.
c
      if( killed_this_time ) then
         if( enforce_node_release ) then
           write(out,9025)
         else
           write(out,9020) critical_angle * two,
     &        init_crit_ang * two, CTOA_range * hundred
         end if
         write(out,9050)  num_ctoa_released_nodes
      end if
c
      if( output_packets .and. killed_this_time ) then
        write(packet_file_no) 19, num_killed_now+1, step, iter
        write(packet_file_no) critical_angle*two,
     &        init_crit_ang*two, CTOA_range*hundred
        do i = 1, num_killed_now
          write(packet_file_no) local_killed_list(i),
     &         local_killed_angles(i)
        end do
      end if
c
      if( last_node_released ) then
          write(out,9060) node_causing_stop
          call store( ' ','kill_limit_restart.db', ldummy1, ldummy2 )
          call warp3d_normal_stop
      end if
c
      if ( debug ) write(out,*) '>> leaving chk_front_master_line'
c
      return
 9000 format(/,' >> node release option invoked for the following:',/)
 9010 format('  master node: ',i7,' CTOA: ',f6.2,' (degrees)',
     & 2x, 'step:',i7,2x,'(X,Y,Z):',e10.3,1x,e10.3,1x,e10.3)
 9015 format('  interim node:',i7,' CTOA: ',f6.2,' (degrees)',
     & 2x, 'step:',i7,2x,'(X,Y,Z):',e10.3,1x,e10.3,1x,e10.3)
 9020 format(/,'        * critical CTOA:',
     &       /,'            for growth    : ',f6.2,' (degrees)',
     &       /,'            for initiation: ',f6.2,' (degrees)',
     &       /,'            nodes released within ',f4.1,' % of ',
     &                        'critical angle.',
     &       /)
 9025 format(/,'        *** user requested node release ***', / )
 9030 format(/,'FATAL ERROR: list length exceeded in',
     &       /,'             chk_crack_front_const. job aborted.' )
 9050 format(8x,'** Total nodes released in model: ',i5,' **')
 9060 format(/,1x,">>>>> Limiting node released: ",i7,
     &       /,7x,"Restart file: kill_limit_restart.db written.",
     &       /,7x,"Job ended normally."// )

c
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine release_node                 *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *    This routine releases a node -- removes the constraint    *
c     *    and sets up reaction force release information.           *
c     *                                                              *
c     ****************************************************************
c
      subroutine release_node( node, node_data_entry, debug )
c
      use global_data, only : cstmap, load, dstmap, ifv, csthed, out
      use node_release_data, only : crkpln_nodes_state,
     &                              crkpln_nodes_react
      use main_data, only : cnstrn_in, cnstrn
      use damage_data
      use constants
c
      implicit none
c
      integer :: node, node_data_entry
      logical :: debug
c
      integer :: dof, cst_ptr
c
      if( debug ) write(out,*) '>> releasing node..'
c
c           To release the node:
c                   a) set global flag that says "no elements
c                      have been released" to false.
c
      no_released_nodes = .false.
c
c
c                   b) remove constraint from node.  Set the
c                      constraint value to null (32460) and remove it
c                      from the constraint linked list (cstmap).
c
      dof            = dstmap(node)+crk_pln_normal_idx-1
      cnstrn(dof)    = d32460
      cnstrn_in(dof) = d32460
      cst_ptr        = csthed
c
      if( cst_ptr .eq. dof ) then
c
c                           we are killing the first entry in the list.
c                           move the list head pointer.
c
         csthed      = cstmap(dof)
         cstmap(dof) = 0
      else
c
c                           cycle through the list to find the entry.  When
c                           found, make the entry that points to the
c                           killed entry point to the next valid entry.
c
         do while (.true.)
            if( cstmap(cst_ptr) .eq. dof ) then
               cstmap(cst_ptr) = cstmap(dof)
               if( cstmap(cst_ptr) .eq. -1 ) csttail = cst_ptr
               cstmap(dof) = 0
               exit
            end if
            cst_ptr = cstmap(cst_ptr)
            if( cst_ptr .eq. -1 ) then
               write(out,*) '>>>> Fatal Error: failed to find',
     &              'constraint in constraint linked list.'
               call die_gracefully
               stop
            end if
         end do
      end if
c
c                  c) Add newly released node to the node release algortihm.
c                     Store the internal force at the node, and initialize
c                     the node's state to 1.
c
      crkpln_nodes_state(node_data_entry) = 1
      if( release_type .eq. 1 ) then
        crkpln_nodes_react(node_data_entry) = ifv(dof)
      end if
      if( release_type .eq. 2 ) then
        load(dof) = load(dof) - point_two * ifv(dof)
        crkpln_nodes_react ( node_data_entry ) = point_eight * ifv(dof)
      end if
c
      num_ctoa_released_nodes = num_ctoa_released_nodes + 1
c
 9999 continue
      if (debug) write(out,*) '<< finished releasing node..'
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine release_front                *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *    Used in the const_front algorithm. This routine, given a  *
c     *    master node, finds the attatched crack front, stores it   *
c     *    for use in the force relaxation algorithm, and releases   *
c     *    the crack front constraints on each crack front node.     *
c     *                                                              *
c     ****************************************************************
c
c
      subroutine release_front( master_node, debug )
c
      use global_data, only : out, dstmap
      use node_release_data, only : crkpln_nodes_state,
     &                              crkpln_nodes_react,
     &                              crack_front_list, inv_crkpln_nodes,
     &                              num_neighbors, neighbor_nodes
      use main_data, only : cnstrn
      use damage_data
      use constants
c
      implicit none
c
      integer :: master_node
      logical :: debug
c
      integer :: i, j, k, dof, entry, idum, list_entry, ptr, ptr_last,
     &           node, neighbor_node, neigh2, dof2,
     &           neighbor_data_entry, node_data_entry
      character(len=1) :: dums
      double precision :: dumd
      real :: dumr
c
      logical :: crack_node
      logical, external :: same_front
c
c
      if( debug ) write(out,*) '>> releasing front...'
c
c           check if master node has already been released. This can happen
c           if two cracks are just about to coalese, such that there is
c           a single line of constrained nodes surrounded by free nodes.
c           In this case, return.
c
      dof = dstmap(master_node)+crk_pln_normal_idx-1
      if( cnstrn(dof) .eq. d32460 ) goto 9999
c
c           initialize const_front_list for holding the crack front as
c           we find it.  This will be used by the reaction force
c           relaxation algorithm.
c
      entry = 0
      do while (.true.)
         entry = entry + 1
c
c               if no free space in the const_front_list
c               is found, then stop execution
c
         if( entry .gt. num_crack_fronts * num_nodes_grwinc ) then
            call errmsg(285, idum, dums, dumr, dumd)
            call die_gracefully
            stop
         endif
c
c               free space found; start the list
c
         if( crack_front_list(entry,1) .eq. 0 ) then
            list_entry = entry
            crack_front_list(list_entry,1) = master_node
            exit
         endif
      enddo
c
c           loop over the crack front node list.  For each node, release
c           it, then check each of its neighbors to find out if they are
c           nodes on the crack front.  If a neighbor is a crackfront node,
c           and it is not on the crack_front_node list, then add it to the
c           list.  In this way, the crack_front_node list will be both
c           processed such that each node on it is released, and the list
c           will be expanded as new crack front nodes are found.
c
      if( debug ) write(out,*) '> check crack_front_list'
      ptr = 0
      ptr_last = 1
      do while (.true.)
         ptr = ptr + 1
         if( debug ) write(out,*) '  - ptr = ', ptr
c
c              check if we have processed all the nodes on the list; exit
c              if we have.
c
         if( ptr .gt. num_nodes_thick ) exit
         node = crack_front_list ( entry, ptr )
         if( ptr .gt. ptr_last ) exit
         if( node .eq. 0 ) exit
         node_data_entry = inv_crkpln_nodes ( node )
c
c              node has not yet been processed.  release it.
c
         call release_node( node, node_data_entry, debug )
         if( node .ne. master_node ) write(out,9000) node,master_node
c
c              now search the node's neighbors, looking for crack front
c              nodes.  A neighbor node is a crack front node if it is
c              a constrained node and one of it's neighbors is an
c              unconstrained node.  Once neighbor nodes have been
c              located, add them to crack_front_list if they are not
c              already in the list.
c
         do i = 1, num_neighbors(node_data_entry)
c
            neighbor_node = neighbor_nodes (i,node_data_entry)
            dof = dstmap(neighbor_node)+crk_pln_normal_idx-1
            if( cnstrn(dof) .eq. d32460 ) cycle
c
c                       check if node is already assigned in
c                       const_front_list; skip it if it is
c
            do j = 1, ptr_last
               if( crack_front_list(list_entry, j )
     &             .eq. neighbor_node ) go to 100
            end do
c
c                       check the neighbors of the neighbors
c
            crack_node = .false.
            neighbor_data_entry = inv_crkpln_nodes(neighbor_node)
c
            do j = 1, num_neighbors(neighbor_data_entry)
               neigh2 = neighbor_nodes (j,neighbor_data_entry)
               dof2   = dstmap(neigh2)+crk_pln_normal_idx-1
               if( cnstrn(dof2) .ne. d32460 ) cycle
c
c                          check if node was released just this step
c
               do k = 1, ptr_last
                  if( crack_front_list(list_entry,k)
     &                .eq. neighbor_node ) cycle
               end do
c
c                          skip node if it is the node we are searching
c                          from, or if node found is not on same crack
c                          front.
c
               if( neigh2 .eq. node ) cycle
           if( .not. same_front(neighbor_node, neighbor_data_entry,
     &               node) ) cycle
c
               crack_node = .true.
               exit
            end do
c
c                       if node is a crack node, then put it in
c                       crack_front_list
c
            if( crack_node ) then
               ptr_last = ptr_last + 1
               if( debug ) write(out,*) '     cfn neighbor:',
     &                neighbor_node,' into ',ptr_last
               if( ptr_last .gt. num_nodes_thick ) then
                  call errmsg(285, idum, dums, dumr, dumd)
                  call die_gracefully
                  stop
               endif
               crack_front_list(list_entry,ptr_last) = neighbor_node
            endif
c
 100        continue
         end do
c
c              loop over to next entry in crack_front_list
c
      end do
c
c              if we are using release by steps, then we don't a
c              need the nodal list anymore.  Clear it out.
c
      if( release_type .eq. 1 ) then
         do i = 1, num_nodes_thick
            crack_front_list(list_entry,i) = 0
         end do
      endif
c
c
 9999 continue
      if( debug ) write(out,*) '<< finished releasing front..'
c
      return
 9000 format('         node: ',i7,' connected to master node: ',i7)
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine update_crack_front           *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *      This routine updates the crack front list, removing     *
c     *      entries of newly released nodes and adding to the       *
c     *      to the list their neighbors that are still constrained  *
c     *      in the direction normal to the crack plane.             *
c     *                                                              *
c     ****************************************************************
c
      subroutine update_crack_front( debug )
c
      use global_data, only : out, dstmap
      use node_release_data, only : inv_crkpln_nodes, num_neighbors,
     &                              neighbor_nodes, crack_front_nodes,
     &                              crkpln_nodes_state, master_nodes
      use main_data, only : cnstrn
      use damage_data
      use constants
c
      implicit none
c
      logical :: debug
c
      integer :: pointer, above_node_ptr, node, node_data_entry,
     &           temp_ptr, neighbor, neighbor_node, dof, new_node,
     &           master_entry, dumi
      character(len=1) :: dums
      real :: dumr
      double precision :: dumd
      logical :: next_node_found
      integer, external :: incrack_master, find_in_list
c
c
c            Traverse linked list of the crack front nodes.
c            If we find a node that has been newly killed:
c                   a) delete the node from the crack front list
c                   b) add any neighbors of this node that are still
c                      constrained to the crack front list
c
      if( debug ) write(out,*) '>>>>>> Updating linked list.'
c
c            initialize the pointers for traversing the list.
c
      pointer = crack_front_start
      above_node_ptr = -1
c
c            top of list traversal loop
c
 10   continue
c
c            a node is newly killed if its state is 1.  get the
c            state and skip node if the state is not 1.
c
      node = crack_front_nodes(pointer,1)
      node_data_entry = inv_crkpln_nodes(node)
      if( crkpln_nodes_state(node_data_entry) .ne. 1 ) goto 20
c
c                  node is newly killed.
c
c                    a) remove it from the list
c
c                       Note: by doing this, we won't need to move
c                             the pointers for the next list entry
c
      temp_ptr = crack_front_nodes(pointer,2)
      call rm_from_list(pointer, above_node_ptr)
      pointer = temp_ptr
c
c                    b) add any still constrained neighbors to
c                       crack_front_nodes
c
      next_node_found = .false.
      do neighbor = 1, num_neighbors(node_data_entry)
c
         neighbor_node = neighbor_nodes(neighbor,node_data_entry)
         dof = dstmap(neighbor_node)+crk_pln_normal_idx-1
         if( cnstrn(dof) .ne. d32460 ) then
c
c                         neighbor is constrained. check if it is
c                         already on the list.  If it isn't, add to
c                         the list.
c
c                         if we are doing const_growth, then replace the
c                         old master node in the list of master nodes with
c                         the new node.
c
           new_node = neighbor_node
           next_node_found = .true.
           if( find_in_list(neighbor_node,dumi).eq.-1 ) then
               call add_to_list( neighbor_node )
               master_entry = incrack_master(node)
               if( master_entry .gt. 0 ) then
                  master_nodes(master_entry) = neighbor_node
               endif
            endif
         end if
      end do
c
c           If we are using the constant front algorithm:
c            Since we have removed a node, we need to update the master
c            line list.  If we removed a node, but didn't add a new
c            one, then a crack has coalesed, so we zero out the
c            corresponding entry.
c
      if( const_front ) call update_master_line( node, new_node,
     &     next_node_found, debug )
c
      goto 30
c
c            node was not released in this step, so we need to
c            move the pointers to the next entry in the list.
c
 20   continue
      above_node_ptr = pointer
      pointer = crack_front_nodes(pointer,2)
c
c            if we are not at end of list, check the next node.
c
 30   continue
      if( pointer .ne. -1 ) goto 10
c
      if( debug ) write(out,*) '>>>>>> Linked list is now updated.'
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine update_master_line           *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *        This routine updates the master line list after       *
c     *        all the crack front releases have been calculated.    *
c     *        If the node in front of the old master node is        *
c     *        unconstrained, then a crack has coalesed. In this     *
c     *        case we zero out the entry in the master_line list    *
c     *        to indicate the crack front no longer exists.         *
c     *                                                              *
c     ****************************************************************
c
      subroutine update_master_line( old_node, new_node,
     &                               next_node_found, debug )
c
      use global_data, only : out
      use node_release_data, only : master_lines, old_angles_at_front
      use damage_data
      use constants
c
      implicit none
c
      integer :: old_node, new_node, num_line
      logical :: next_node_found, debug
c
      integer :: i, j
c
c          if no replacement node was found, then the crack is coalesing
c          or ending.  Zero out the corresponding entry in the master_lines
c          structure.
c
      if( .not. next_node_found ) then
c
         do i = 1, num_crack_fronts
            if( master_lines(i,1) .eq. old_node ) then
               if( debug ) write(out,*)
     &           '-> dead node, zero master_line:', old_node
               do j = 1, num_nodes_back + 1
                  master_lines(i,j) = 0
               enddo
            endif
         enddo
c
c           if the crack is growing normally, then new_node holds the
c           new crack front node.  move all of the entries in the
c           master_line list over one to make room for the new crack
c           tip node.
c
      else
c
         do num_line = 1, num_crack_fronts
            if( master_lines(num_line,1) .eq. old_node ) then
               do i = num_nodes_back,1,-1
                  master_lines(num_line,i+1) =
     &                 master_lines(num_line,i)
               enddo
               exit
            endif
         enddo
         master_lines(num_line,1) = new_node
c
      endif
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine get_slope                    *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *        This routine gets the angle between a crack front     *
c     *        node and an unconstrained neighboring crack face      *
c     *        node.  This is used in the node_release algorithm.    *
c     *                                                              *
c     ****************************************************************
c
      subroutine get_slope( node1, node2, normal, angle )
c
      use global_data, only : c, dstmap, u
      use main_data, only : crdmap
      use constants
c
      implicit none
c
      integer :: node1, node2, normal
      double precision :: angle
c
      double precision :: dist_x, dist_y, dist_z, plane_dist
c
c               calculate the distances in the three coordinate
c               directions between the two nodes.
c
      dist_x = c(crdmap(node1)) + u(dstmap(node1))
     &      - ( c(crdmap(node2)) + u(dstmap(node2)))
      dist_y = c(crdmap(node1)+1) + u(dstmap(node1)+1)
     &      - ( c(crdmap(node2)+1) + u(dstmap(node2)+1))
      dist_z = c(crdmap(node1)+2) + u(dstmap(node1)+2)
     &      - ( c(crdmap(node2)+2) + u(dstmap(node2)+2))
c
      if( normal .eq. 1 ) then
c
c               x is direction of crack plane normal
c
         plane_dist = sqrt (dist_y ** 2 + dist_z ** 2)
         angle = atan( dist_x / plane_dist ) * oneeighty / pi
c
      else if( normal .eq. 2 ) then
c
c               y is direction of crack plane normal
c
         plane_dist = sqrt (dist_x ** 2 + dist_z ** 2)
         angle = atan( dist_y / plane_dist ) * oneeighty / pi
c
      else if( normal .eq. 3 ) then
c
c               z is direction of crack plane normal
c
         plane_dist = sqrt (dist_x ** 2 + dist_y ** 2)
         angle = atan( dist_z / plane_dist ) * oneeighty / pi
c
      end if
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine get_slope_master_line        *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *        This routine gets the CTOA at a master node at a      *
c     *        defined distance behind the crack front, as defined   *
c     *        by the master_line for the node.                      *
c     *                                                              *
c     ****************************************************************
c
      subroutine get_slope_master_line( num_line, use_init, angle,
     &                                  num_elems_back)
c
      use global_data, only : out
      use node_release_data, only : master_lines, crkpln_nodes_state,
     &                              inv_crkpln_nodes
      use damage_data
      use constants
c
      implicit none
c
      integer :: num_line, num_elems_back
      double precision :: angle
      logical :: use_init
c
      integer :: i, base_node, next_node, node_idx, istate
      double precision :: dist, d_dist, height, factor, d_height,
     &                    ctoa_dist_back, dumd, old_height
      real :: dumr
      character(len=1) :: dums
      logical :: debug
c
      debug = .false.
      num_elems_back = 0
c
c           if master_lines entry is zeroed, then set angle to zero
c           and leave.
c
      if( master_lines(num_line,1) .eq. 0 ) then
         angle = zero
         return
      endif
c             determine if we should use the initial angle or continuation
c             angle. If the second node in the master_line list has never
c             been released, then the crack has not yet moved, so we can
c             use the initiation angle and distance.
c
      i = master_lines(num_line,2)
      istate = crkpln_nodes_state(inv_crkpln_nodes(i))
      if( istate .eq. 0 ) then
         use_init = .true.
         ctoa_dist_back = init_ctoa_dist
         if( debug ) write(out,*) 'Using initiation angle'
      else
         use_init = .false.
         ctoa_dist_back = ctoa_dist
         if( debug ) write(out,*) 'Using release angle'
      endif
c
c           calculate arc distance along the master line until we
c           find the nodes which bound our distance, and interpolate.
c
      dist = zero
      height = zero
      base_node = master_lines(num_line,1)
      next_node = master_lines(num_line,2)
      node_idx = 2
      old_height = zero
c
c               travel along the master line
c
      do while (.true.)
c
         num_elems_back = num_elems_back + 1
c
c                  get distance
c
         call get_dist( next_node, base_node, crk_pln_normal_idx,
     &                  d_dist, d_height )
c
c                  check distance -- if further than we want to go,
c                  calculate the angle
c
         if( dist + d_dist .ge. ctoa_dist_back ) then
c
            factor = ( ctoa_dist_back - dist ) / d_dist
            dist = dist + d_dist * factor
            height = height + d_height * factor
            angle = atan( height / dist ) * oneeighty / pi
            exit
c
         else
c
c                  we haven't passed our distance yet.  Continue travelling
c                  along the master line.
c
            dist = dist + d_dist
            height = height + d_height
            node_idx = node_idx + 1
            if( node_idx .gt. num_nodes_back + 1 ) then
               call errmsg( 290, num_nodes_back, dums, dumr, dumd )
               call die_gracefully
               stop
            endif
            base_node = next_node
            next_node = master_lines(num_line,node_idx)
            if( next_node .eq. 0 ) then
               call errmsg( 291, master_lines(num_line,1), dums, dumr,
     &              dumd )
               call die_gracefully
            endif
c
         endif
c
      end do
c
      if( debug ) then
         write(out,*) ' nodes;',base_node, next_node
         write(out,*) ' number of elements back;',num_elems_back
         write(out,*) ' dist:',dist, '  height:',height
      endif
c
      return
      end
c
c

c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine get_dist_to_node             *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *        This routine gets the distance from a given old       *
c     *        master node to the current master node.               *
c     *                                                              *
c     ****************************************************************
c
c
c
      subroutine get_dist_to_node( num_line, old_node, dist )
c
      use global_data, only : out
      use node_release_data, only : master_lines
      use damage_data
      use constants
c
      implicit none
c
      integer :: num_line, old_node
      double precision :: dist
c
      integer :: base_node, next_node, node_idx
      double precision :: d_dist, height, factor, d_height,
     &                    ctoa_dist_back, dumd
      real :: dumr
      character(len=1) :: dums
      logical :: debug, use_init
c
      debug = .false.
c
c           if master_lines entry is zeroed, then set distance to zero
c           and leave.
c
      if( master_lines(num_line,1) .eq. 0 ) then
         dist = zero
         return
      endif
c
c           calculate arc distance along the master line until we
c           find the nodes which bound our distance, and interpolate.
c
      dist = zero
      base_node = master_lines(num_line,1)
      next_node = master_lines(num_line,2)
      node_idx = 2
c
c               travel along the master line
c
      do while (.true.)
c
c                  get distance
c
         call get_dist( next_node, base_node, crk_pln_normal_idx,
     &                  d_dist, d_height )
c
c                  check distance -- if further than we want to go,
c                  calculate the angle
c
         dist = dist + d_dist
         node_idx = node_idx + 1
         if( node_idx .gt. num_nodes_back + 1 ) then
            call errmsg( 290, num_nodes_back, dums, dumr, dumd )
            call die_gracefully
         endif
         base_node = next_node
         next_node = master_lines(num_line,node_idx)
         if( base_node .eq. old_node ) then
            exit
         else if( next_node .eq. 0 ) then
            call errmsg( 291, master_lines(num_line,1), dums, dumr,
     &           dumd )
            call die_gracefully
            stop
         endif
c
      end do
c
      if( debug ) then
         write(out,*)
     &       '>>>>> in get_dist_to_node: dist to old node:',dist
      endif
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine get_dist                     *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *        This routine calculates the distance between two      *
c     *        nodes originally on the crack plane.  Distance is     *
c     *        given in terms of planar distance and height.         *
c     *                                                              *
c     ****************************************************************
c
      subroutine get_dist( node1, node2, normal, plane_dist, height )
c
      use global_data, only : c, u, dstmap
      use main_data, only : crdmap
c
      implicit none
c
      integer :: node1, node2, normal
      double precision :: height, plane_dist
c
      double precision :: dist_x, dist_y, dist_z
c
c               calculate the distances in the three coordinate
c               directions between the two nodes.
c
      dist_x = c(crdmap(node1)) + u(dstmap(node1))
     &      - ( c(crdmap(node2)) + u(dstmap(node2)))
      dist_y = c(crdmap(node1)+1) + u(dstmap(node1)+1)
     &      - ( c(crdmap(node2)+1) + u(dstmap(node2)+1))
      dist_z = c(crdmap(node1)+2) + u(dstmap(node1)+2)
     &      - ( c(crdmap(node2)+2) + u(dstmap(node2)+2))
c
      if( normal .eq. 1 ) then
c
c               x is direction of crack plane normal
c
         plane_dist = sqrt (dist_y ** 2 + dist_z ** 2)
         height = dist_x
c
      else if( normal .eq. 2 ) then
c
c               y is direction of crack plane normal
c
         plane_dist = sqrt (dist_x ** 2 + dist_z ** 2)
         height = dist_y
c
      else if( normal .eq. 3 ) then
c
c               z is direction of crack plane normal
c
         plane_dist = sqrt (dist_x ** 2 + dist_y ** 2)
         height = dist_z
c
      end if
c
      return
      end
c
c
c     ****************************************************************
c     *                                                              *
c     *                      function same_front                     *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 6/22/21 rhd (cleanup)      *
c     *                                                              *
c     *   This function, given two nodes on crack fronts, returns    *
c     *   as true if the nodes are neighbors on the same crack front.*
c     *   This function is used in the const_front algorithm.        *
c     *                                                              *
c     ****************************************************************
c
c
      logical function same_front( new_node, new_data_entry,
     &                             check_node )
c
      use global_data, only : dstmap
      use node_release_data, only : inv_crkpln_nodes, num_neighbors,
     &                              neighbor_nodes
      use main_data, only : cnstrn, cnstrn_in
      use damage_data
      use constants
c
      implicit none
c
      integer :: new_node, new_data_entry, check_node
c
      integer :: i, j, k, node1, node1_data_entry, dof, node2,
     &           node2_data_entry, node3
c
c           This function checks a newly found crack front node to see if
c           it is on the same front as the node from which the new node
c           was found.  This becomes an issue when two crack fronts are
c           back to back; in this case, the crack front location algorithm
c           can break down, finding both crack fronts.
c
c           The basic idea behind this routine is to loop back over the
c           unconstrained nodes adjoining the crack front to make sure
c           that there is a path from the new_node to the check_node
c           consisting of precisely two unconstrained nodes.
c
c           This routine operates as follows:
c
c           - find an unconstrained neighbor to the new node.
c           - loop over this neighbor's neighbors.  for each that is
c             unconstrained:
c                - check the constrained neighbors of this new node.  If
c                  one is the check_node, then new_node and check_node
c                  are on the same crack front.
c
      same_front = .false.
      do i = 1, num_neighbors(new_data_entry)
         node1 = neighbor_nodes(i,new_data_entry)
         if( node1 .eq. check_node ) cycle
         node1_data_entry = inv_crkpln_nodes(node1)
         dof = dstmap(node1) + crk_pln_normal_idx-1
         if( cnstrn(dof) .ne. d32460 ) cycle
         do j = 1, num_neighbors(node1_data_entry)
            node2 = neighbor_nodes(j,node1_data_entry)
            node2_data_entry = inv_crkpln_nodes(node2)
            dof = dstmap(node2)+crk_pln_normal_idx-1
            if( cnstrn(dof) .ne. d32460 ) cycle
            do k = 1, num_neighbors(node2_data_entry)
               node3 = neighbor_nodes(k,node2_data_entry)
               if( node3 .eq. check_node ) then
                  same_front = .true.
                  go to 9999
               end if
            end do
         end do
      end do
c
 9999 continue
      return
      end

ext: f


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     #############################################################
c     ##                                                         ##
c     ##  function random  --  portable random number generator  ##
c     ##                                                         ##
c     #############################################################
c
c
c     "random" generates a random number on [0,1] via a long
c     period generator due to L'Ecuyer with Bays-Durham shuffle
c
c     literature references:
c
c     P. L'Ecuyer, Communications of the ACM, 31, 742-774 (1988)
c
c     W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P.
c     Flannery, Numerical Recipes (Fortran), 2nd Ed., Cambridge
c     University Press, 1992, Section 7.1
c
c
      function random ()
      implicit none
      include 'sizes.i'
      include 'inform.i'
      include 'iounit.i'
      include 'keys.i'
      integer im1,ia1,iq1,ir1
      integer im2,ia2,iq2,ir2
      integer big,nshuffle
      integer imm1,ndiv
      real*8 factor
      parameter (im1=2147483563)
      parameter (ia1=40014)
      parameter (iq1=53668)
      parameter (ir1=12211)
      parameter (im2=2147483399)
      parameter (ia2=40692)
      parameter (iq2=52774)
      parameter (ir2=3791)
      parameter (big=141803398)
      parameter (nshuffle=32)
      parameter (imm1=im1-1)
      parameter (ndiv=1+imm1/nshuffle)
      parameter (factor=1.0d0/im1)
      integer i,k,iy,next
      integer seed,seed2
      integer year,month,day
      integer hour,minute,second
      integer ishuffle(nshuffle)
      real*8 random
      logical first
      character*20 keyword
      character*120 record
      character*120 string
      save first
      save seed,seed2
      save iy,ishuffle
      data first  / .true. /
c
c
c     random number seed is first set to a big number,
c     then incremented by the seconds elapsed this decade
c
      if (first) then
         first = .false.
         seed = big
         call calendar (year,month,day,hour,minute,second)
         year = mod(year,10)
         seed = seed + 32140800*year + 2678400*(month-1)
         seed = seed + 86400*(day-1) + 3600*hour
         seed = seed + 60*minute + second
c
c     search the keywords for a random number seed
c
         do i = 1, nkey
            next = 1
            record = keyline(i)
            call gettext (record,keyword,next)
            call upcase (keyword)
            if (keyword(1:11) .eq. 'RANDOMSEED ') then
               string = record(next:120)
               read (string,*,err=10)  seed
               seed = max(1,seed)
            end if
   10       continue
         end do
c
c     print the value used for the random number seed
c
         if (verbose) then
            write (iout,20)  seed
   20       format (/,' Random Number Generator Initialized',
     &                 ' with SEED :',3x,i12)
         end if
c
c     warm up and then load the shuffling table
c
         seed2 = seed
         do i = nshuffle+8, 1, -1
            k = seed / iq1
            seed = ia1 * (seed-k*iq1) - k*ir1
            if (seed .lt. 0)  seed = seed + im1
            if (i .le. nshuffle)  ishuffle(i) = seed
         end do
         iy = ishuffle(1)
      end if
c
c     get a new random number value each call
c
      k = seed / iq1
      seed = ia1*(seed-k*iq1) - k*ir1
      if (seed .lt. 0)  seed = seed + im1
      k = seed2 / iq2
      seed2 = ia2*(seed2-k*iq2) - k*ir2
      if (seed2 .lt. 0)  seed2 = seed2 + im2
      i = 1 + iy/ndiv
      iy = ishuffle(i) - seed2
      ishuffle(i) = seed
      if (iy .lt. 1)  iy = iy + imm1
      random = factor * iy
c
c     print the value of the current random number
c
c     if (debug) then
c        write (iout,30)  random
c  30    format (' RANDOM  --  The Random Number Value is',f12.8)
c     end if
      return
      end
c
c
c     ############################################################
c     ##                                                        ##
c     ##  function normal  --  random number from normal curve  ##
c     ##                                                        ##
c     ############################################################
c
c
c     "normal" generates a random number from a normal Gaussian
c     distribution with a mean of zero and a variance of one
c
c
      function normal ()
      implicit none
      include 'inform.i'
      include 'iounit.i'
      real*8 random,v1,v2,rsq
      real*8 factor,store,normal
      logical compute
      save compute,store
      data compute  / .true. /
c
c
c     get a pair of random values from the distribution
c
      if (compute) then
   10    continue
         v1 = 2.0d0 * random () - 1.0d0
         v2 = 2.0d0 * random () - 1.0d0
         rsq = v1**2 + v2**2
         if (rsq .ge. 1.0d0)  goto 10
         factor = sqrt(-2.0d0*log(rsq)/rsq)
         store = v1 * factor
         normal = v2 * factor
         compute = .false.
c
c     use the second random value computed at the last call
c
      else
         normal = store
         compute = .true.
      end if
c
c     print the value of the current random number
c
c     if (debug) then
c        write (iout,20)  normal
c  20    format (' NORMAL  --  The Random Number Value is',f12.8)
c     end if
      return
      end
c
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine ranvec  --  unit vector in random direction  ##
c     ##                                                          ##
c     ##############################################################
c
c
c     "ranvec" generates a unit vector in 3-dimensional
c     space with uniformly distributed random orientation
c
c     literature references:
c
c     G. Marsaglia, Ann. Math. Stat., 43, 645 (1972)
c
c     R. C. Rapaport, The Art of Molecular Dynamics Simulation,
c     2nd Edition, Cambridge University Press, 2004, Section 18.4
c
c
      subroutine ranvec (vector)
      implicit none
      include 'inform.i'
      include 'iounit.i'
      real*8 x,y,s
      real*8 random
      real*8 vector(3)
c
c
c     get a pair of appropriate components in the plane
c
      s = 2.0d0
      do while (s .ge. 1.0d0)
         x = 2.0d0 * random () - 1.0d0
         y = 2.0d0 * random () - 1.0d0
         s = x**2 + y**2
      end do
c
c     construct the 3-dimensional random unit vector
c
      vector(3) = 1.0d0 - 2.0d0*s
      s = 2.0d0 * sqrt(1.0d0 - s)
      vector(2) = s * y
      vector(1) = s * x
c
c     print the components of the random unit vector
c
c     if (debug) then
c        write (iout,10)  vector(1),vector(2),vector(3)
c  10    format (' RANVEC  --  The Random Vector is',3f10.4)
c     end if
      return
      end
c
c
c     ##############################################################
c     ##                                                          ##
c     ##  subroutine sphere  --  uniform set of points on sphere  ##
c     ##                                                          ##
c     ##############################################################
c
c
c     "sphere" finds a specified number of uniformly distributed
c     points on a sphere of unit radius centered at the origin
c
c     literature reference:
c
c     E. B. Saff and A. B. J. Kuijlaars, "Distributing Many
c     Points on a Sphere", The Mathematical Intelligencer,
c     19, 5-11 (1997)
c
c
      subroutine sphere (ndot,dot)
      implicit none
      include 'math.i'
      integer i,ndot
      real*8 theta,phi
      real*8 h,phiold
      real*8 tot,tot1
      real*8 dot(3,*)
c
c
c     find spherical coordinates then convert to Cartesian
c
      tot = dble(ndot)
      tot1 = dble(ndot-1)
      do i = 1, ndot
         h = -1.0d0 + 2.0d0*dble(i-1)/tot1
         theta = acos(h)
         if (i.eq.1 .or. i.eq.ndot) then
            phi = 0.0d0
         else
            phi = mod(phiold+3.6d0/sqrt(tot*(1.0d0-h*h)),2.0d0*pi)
         end if
         dot(1,i) = sin(theta) * cos(phi)
         dot(2,i) = sin(theta) * sin(phi)
         dot(3,i) = cos(theta)
         phiold = phi
      end do
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
      SUBROUTINE CALBAL5
C
C CHANGE RECORD
C **  SUBROUTINES CALBAL CALCULATE GLOBAL VOLUME, MASS, MOMENTUM,
C **  AND ENERGY BALANCES
C
      USE GLOBAL
	IMPLICIT NONE
	REAL::ENEEND,ENEOUT,VOLBMO,SALBMO,DYEBMO,UMOBMO,VMOBMO,ENEBMO
	REAL::VOLERR,SALERR,DYEERR,UMOERR,VMOERR,ENEERR,RVERDE,RSERDE
	REAL::RDERDE,RUERDE,REERDE,RVERDO,RSERDO,RUERDO,ENEBEG,RDERDO
	REAL::REERDO,RUMERDE,RVMERDE,RUMERDO,UUEBMO,VVEBMO,PPEBMO
	REAL::BBEBMO,RVMERDO,PPEEND,BBEEND,AMOEND,VOLEND,SALEND,DYEEND
	REAL::VMOEND,UUEEND,VVEEND,UMOEND
	INTEGER::K,L,LN
C
C **  CHECK FOR END OF BALANCE PERIOD
C
      IF(NBAL.EQ.NTSMMT)THEN
C
C **  CALCULATE ENDING VOLUME, SALT MASS, DYE MASS, MOMENTUM, KINETIC
C **  ENERGY AND POTENTIAL ENERGY, AND ASSOCIATED FLUXES
C
        VOLEND=0.
        SALEND=0.
        DYEEND=0.
        UMOEND=0.
        VMOEND=0.
        UUEEND=0.
        VVEEND=0.
        PPEEND=0.
        BBEEND=0.
        DO L=2,LA
          LN=LNC(L)
          VOLEND=VOLEND+SPB(L)*DXYP(L)*HP(L)
         UMOEND=UMOEND+SPB(L)*0.5*DXYP(L)*HP(L)*(DYIU(L)*HUI(L)*UHDYE(L)
     &        +DYIU(LEAST(L))*HUI(LEAST(L))*UHDYE(LEAST(L)))
         VMOEND=VMOEND+SPB(L)*0.5*DXYP(L)*HP(L)*(DXIV(L)*HVI(L)*VHDXE(L)
     &        +DXIV(LN)*HVI(LN)*VHDXE(LN))
          PPEEND=PPEEND+SPB(L)*0.5*DXYP(L)
     &        *(GI*P(L)*P(L)-G*BELV(L)*BELV(L))
        ENDDO
        AMOEND=SQRT(UMOEND*UMOEND+VMOEND*VMOEND)
        DO K=1,KC
          DO L=2,LA
            LN=LNC(L)
            SALEND=SALEND+SCB(L)*DXYP(L)*HP(L)*SAL(L,K)*DZC(K)
            DYEEND=DYEEND+SCB(L)*DXYP(L)*HP(L)*DYE(L,K)*DZC(K)
            UUEEND=UUEEND+SPB(L)*0.125*DXYP(L)*HP(L)*DZC(K)
     &          *( (U(L,K)+U(LEAST(L),K))*(U(L,K)+U(LEAST(L),K)) )
            VVEEND=VVEEND+SPB(L)*0.125*DXYP(L)*HP(L)*DZC(K)
     &          *( (V(L,K)+V(LN,K))*(V(L,K)+V(LN,K)) )
            BBEEND=BBEEND+SPB(L)*GP*DXYP(L)*HP(L)*DZC(K)*( BELV(L)
     &          +0.5*HP(L)*(Z(K)+Z(K-1)) )*B(L,K)
          ENDDO
        ENDDO
        UUEOUT=DT*UUEOUT
        VVEOUT=DT*VVEOUT
        PPEOUT=DT*PPEOUT
        BBEOUT=DT*BBEOUT
        VOLOUT=DT*VOLOUT
        SALOUT=DT*SALOUT
        DYEOUT=DT*DYEOUT
        UMOOUT=DT*UMOOUT
        VMOOUT=DT*VMOOUT
        ENEBEG=UUEBEG+VVEBEG+PPEBEG+BBEBEG
        ENEEND=UUEEND+VVEEND+PPEEND+BBEEND
        ENEOUT=UUEOUT+VVEOUT+PPEOUT+BBEOUT
        VOLBMO=VOLBEG-VOLOUT
        SALBMO=SALBEG-SALOUT
        DYEBMO=DYEBEG-DYEOUT
        UMOBMO=UMOBEG-DYEOUT
        VMOBMO=VMOBEG-DYEOUT
        ENEBMO=ENEBEG-ENEOUT
        VOLERR=VOLEND-VOLBMO
        SALERR=SALEND-SALBMO
        DYEERR=DYEEND-DYEBMO
        UMOERR=UMOEND-UMOBMO
        VMOERR=VMOEND-VMOBMO
        ENEERR=ENEEND-ENEBMO
        RVERDE=-9999.
        RSERDE=-9999.
        RDERDE=-9999.
        RUERDE=-9999.
        RVERDE=-9999.
        REERDE=-9999.
        RVERDO=-9999.
        RSERDO=-9999.
        RDERDO=-9999.
        RUERDO=-9999.
        RVERDO=-9999.
        REERDO=-9999.
        IF(VOLEND.NE.0.) RVERDE=VOLERR/VOLEND
        IF(SALEND.NE.0.) RSERDE=SALERR/SALEND
        IF(DYEEND.NE.0.) RDERDE=DYEERR/DYEEND
        IF(UMOEND.NE.0.) RUMERDE=UMOERR/UMOEND
        IF(VMOEND.NE.0.) RVMERDE=VMOERR/VMOEND
        IF(ENEEND.NE.0.) REERDE=ENEERR/ENEEND
        IF(VOLOUT.NE.0.) RVERDO=VOLERR/VOLOUT
        IF(SALOUT.NE.0.) RSERDO=SALERR/SALOUT
        IF(DYEOUT.NE.0.) RDERDO=DYEERR/DYEOUT
        IF(UMOOUT.NE.0.) RUMERDO=UMOERR/UMOOUT
        IF(VMOOUT.NE.0.) RVMERDO=VMOERR/VMOOUT
        IF(ENEOUT.NE.0.) REERDO=ENEERR/ENEOUT
C
C **  OUTPUT BALANCE RESULTS TO FILE BAL.OUT
C
        IF(JSBAL.EQ.1)THEN
          OPEN(89,FILE='BAL.OUT',STATUS='UNKNOWN')
          CLOSE(89,STATUS='DELETE')
          OPEN(89,FILE='BAL.OUT',STATUS='UNKNOWN')
          JSBAL=0
        ELSE
          OPEN(89,FILE='BAL.OUT',POSITION='APPEND',STATUS='UNKNOWN')
        ENDIF
        WRITE(89,890)NTSMMT,N
        WRITE(89,891)
        WRITE(89,892)VOLBEG,SALBEG,DYEBEG,ENEBEG,UMOBEG,VMOBEG,AMOBEG
        WRITE(89,900)
        WRITE(89,893)
        WRITE(89,892)VOLOUT,SALOUT,DYEOUT,ENEOUT,UMOOUT,VMOOUT
        WRITE(89,900)
        WRITE(89,894)
        WRITE(89,892)VOLBMO,SALBMO,DYEBMO,ENEBMO,UMOBMO,VMOBMO
        WRITE(89,900)
        WRITE(89,895)
        WRITE(89,892)VOLEND,SALEND,DYEEND,ENEEND,UMOEND,VMOEND,AMOEND
        WRITE(89,900)
        WRITE(89,896)
        WRITE(89,892)VOLERR,SALERR,DYEERR,ENEERR,UMOERR,VMOERR
        WRITE(89,900)
        WRITE(89,897)
        WRITE(89,892)RVERDE,RSERDE,RDERDE,REERDE,RUMERDE,RVMERDE
        WRITE(89,900)
        WRITE(89,898)
        WRITE(89,892)RVERDO,RSERDO,RDERDO,REERDO,RUMERDO,RVMERDO
        WRITE(89,899)
        UUEBMO=UUEBEG-UUEOUT
        VVEBMO=VVEBEG-VVEOUT
        PPEBMO=PPEBEG-PPEOUT
        BBEBMO=BBEBEG-BBEOUT
        WRITE(89,901)UUEBEG
        WRITE(89,902)UUEOUT
        WRITE(89,903)UUEBMO
        WRITE(89,904)UUEEND
        WRITE(89,900)
        WRITE(89,905)VVEBEG
        WRITE(89,906)VVEOUT
        WRITE(89,907)VVEBMO
        WRITE(89,908)VVEEND
        WRITE(89,900)
        WRITE(89,909)PPEBEG
        WRITE(89,910)PPEOUT
        WRITE(89,911)PPEBMO
        WRITE(89,912)PPEEND
        WRITE(89,900)
        WRITE(89,913)BBEBEG
        WRITE(89,914)BBEOUT
        WRITE(89,915)BBEBMO
        WRITE(89,916)BBEEND
        WRITE(89,900)
        WRITE(89,899)
        CLOSE(89)
  890 FORMAT (' VOLUME, MASS, AND ENERGY BALANCE OVER',I5,' TIME STEPS'
     &    ,' ENDING AT TIME STEP',I5,//)
  891 FORMAT (' INITIAL VOLUME    INITIAL SALT    INITIAL DYE     '
     &    ,'INITIAL ENER    INITIAL UMO     INITIAL VMO     '
     &    ,'INITIAL AMO',/)
  892 FORMAT (1X,7(E14.6,2X))
  893 FORMAT (' VOLUME OUT        SALT OUT        DYE OUT         '
     &    ,'ENERGY OUT      UMO OUT         VMO OUT',/)
  894 FORMAT (' INITIAL-OUT VOL   INIT-OUT SALT   INIT-OUT DYE    '
     &    ,'INIT-OUT ENER   INIT-OUT UMO    INIT-OUT VMO',/)
  895 FORMAT (' FINAL VOLUME      FINAL SALT      FINAL DYE       '
     &    ,'FINAL ENERGY    FINAL UMO       FINAL VMO       '
     &    ,'FINAL AMO',/)
  896 FORMAT (' VOLUME ERR        SALT ERR        DYE ERR         '
     &    ,'ENERGY ERR      UMO ERR         VMO ERR',/)
  897 FORMAT (' R VOL/END ER      R SAL/END ER    R DYE/END ER    '
     &    ,'R ENE/END ER    R UMO/END ER    R VMO/END ER',/)
  898 FORMAT (' R VOL/OUT ER      R SAL/OUT ER    R DYE/OUT ER    '
     &    ,'R ENE/OUT ER    R UMO/OUT ER    R VMO/OUT ER',/)
  899 FORMAT (////)
  900 FORMAT (//)
  901 FORMAT(' UUEBEG =  ',E14.6)
  902 FORMAT(' UUEOUT =  ',E14.6)
  903 FORMAT(' UUEBMO =  ',E14.6)
  904 FORMAT(' UUEEND =  ',E14.6)
  905 FORMAT(' VVEBEG =  ',E14.6)
  906 FORMAT(' VVEOUT =  ',E14.6)
  907 FORMAT(' VVEBMO =  ',E14.6)
  908 FORMAT(' VVEEND =  ',E14.6)
  909 FORMAT(' PPEBEG =  ',E14.6)
  910 FORMAT(' PPEOUT =  ',E14.6)
  911 FORMAT(' PPEBMO =  ',E14.6)
  912 FORMAT(' PPEEND =  ',E14.6)
  913 FORMAT(' BBEBEG =  ',E14.6)
  914 FORMAT(' BBEOUT =  ',E14.6)
  915 FORMAT(' BBEBMO =  ',E14.6)
  916 FORMAT(' BBEEND =  ',E14.6)
        NBAL=0
      ENDIF
      NBAL=NBAL+1
      RETURN
      END