anchor.diff:          real*8, allocatable :: Estiff(:,:,:)
   c
   c.... ------------------->   layer base elements   <-------------------
  +c
  +      if (numgc .gt. 0) then
   c
   c.... calculate the normal of each growth curve based on new boundary positions
   c
  
positive.code:           call dbl_to_iBCelas(flag(i), iBC(i))
            call setBLbc(disp(i,:), iBC(i), BC(i, ndof+2:ndof+5))
          enddo
  c
  c.... -------------------->   interior elements   <--------------------
  c
  
negative.code: c
  c.... end calculation of growth curve normal
  c
        endif ! end if numgc greater than 0
  c
  c.... ---------------->   Re-position layered mesh   <-----------------
  c
  c.... if no growth curve, skip whole re-positioning method
  c
        if (numgc .gt. 0) then
  c
  c.... loop over growth curves
  c
          do i = 1, numnp
            call iBCelas_to_dbl(iBC(i), flag(i))
          enddo
          listcounter = 0
          ioffset = 1 ! the ID starts from 1 in phasta
          do ngc = 1, numgc
  