anchor.diff:       qlim = elim*hlim - flim*glim
        NN = CEILING((LOG10(flim) - JJ) / (2*LOG10(ABS(alpha))))
   
  -     print *, '6) Allocate e, f'
  +     !print *, '6) Allocate e, f'
        IF (ALLOCATED(e)) DEALLOCATE(e)
        ALLOCATE(e(NN))
        e = 0.0d0 
  
positive.code:      w(j) = y(j)-2*y(j+1)+y(j+2)
    END DO
  
    print *, '3) determine lamr, a0, a1, a2, x1, x2'
    lamr = lam*Tcu
    a0 = 6.d0 + lamr*2.d0/3.d0
    a1 = 4.d0 - lamr/6.d0
  
negative.code:      alpha = 0.5d0*(x1 - SQRT(x1**2-4.d0))
       beta  = 0.5d0*(x2 + SQRT(x2**2-4.d0))
    END IF
  ! print *, 'alpha, beta = ', alpha, beta
  
  ! IF (JJ > LOG10(alpha*beta) - (nc-1)*2*LOG10(ABS(alpha))) THEN 
    IF (JJ > LOG10(DBLE(alpha*beta)) - (nc-1)*2*LOG10(ABS(alpha))) THEN 
       !print *, '5) JJ if block: GREATER THAN'
  
       !Use untruncated algorithm since NN > nc-2 
       !Factor coefficient matrix, solve triangular systems, find trace 
  