anchor.diff:    a2 = SQRT(atmp)
     x1 =  (a1 + a2) / 2.d0
     x2 = -(a2 - a1) / 2.d0
  -  print *,' a1**2 - 4.d0*(a0-2.d0) = ', atmp
  -  print *, 'a0, a1, a2 =', a0, a1, a2
  -  print *, 'x1, x2 =', x1, x2
  -
  -  print *, '4) determine alpha, beta'
  +! print *,' a1**2 - 4.d0*(a0-2.d0) = ', atmp
  +! print *, 'a0, a1, a2 =', a0, a1, a2
  +! print *, 'x1, x2 =', x1, x2
  +
  +! print *, '4) determine alpha, beta'
     IF (lamr > 24.d0) THEN 
        alpha = 0.5d0*(x1 + SQRT(x1**2-4.d0))
        beta  = 0.5d0*(x2 + SQRT(x2**2-4.d0))
  
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:      qlim = elim*hlim - flim*glim
       NN = CEILING((LOG10(flim) - JJ) / (2*LOG10(ABS(alpha))))
  
       !print *, '6) Allocate e, f'
       IF (ALLOCATED(e)) DEALLOCATE(e)
       ALLOCATE(e(NN))
       e = 0.0d0 
  