anchor.diff:    abserr = 1e-6_dp
     relerr = 1e-3_dp
   
  -! Set experimental data
  -  t   = 1_dp
  -  mu  =-0.9_dp*t
  -  x0  = 1e-2_dp
  -  epsf=-0.9_dp*x0
  -  V   = 1e-2_dp
  -  Uc  = 1_dp
  -  Uf  = 1e-3_dp
  -
   ! frequency
     omega = dcmplx(om,del)
  -
  -! Set up matrix as an array of size 16
  -  Gkinv = (/ complex(dp) :: &
  -       0,  V, -Uc,  0, &
  -       V,  0,   0, Uf, &
  -     -Uc,  0,   0,  V, &
  -       0, Uf,   V,  0 /)
  -  Gkqinv = (/ complex(dp) :: &
  -       0,  V, -Uc,  0, &
  -       V,  0,   0, Uf, &
  -     -Uc,  0,   0,  V, &
  -       0, Uf,   V,  0 /)
  -
  -  open(log, file="susy_qpi.log", position="append", status="old")
  -  write(log, *) "t:     ", t
  -  write(log, *) "mu:    ", mu
  -  write(log, *) "x0:    ", x0
  -  write(log, *) "epsf:  ", epsf
  -  write(log, *) "V:     ", V
  -  write(log, *) "Uc:    ", Uc
  -  write(log, *) "Uf:    ", Uf
  -  write(log, *) "omega: ", omega
  -  close(log)
   
   ! Step through 1/8 triangle of BZ
     steps=100
  
positive.code:   close(log)
    open(dat, file=filename, status="unknown")
  
  ! Set constants
    Pi=4.0_dp*datan(1.0_dp)
  
  ! Settings for dcuhre
    a(1) = -1.0_dp*Pi
    b(1) = Pi
  
negative.code:         close(log)
        end if
        write(dat,*)  qx,  qy, result(1)
        write(dat,*)  qx, -qy, result(1)
        write(dat,*) -qx,  qy, result(1)
        write(dat,*) -qx, -qy, result(1)
        write(dat,*)  qy,  qx, result(1)
        write(dat,*)  qy, -qx, result(1)
        write(dat,*) -qy,  qx, result(1)
        write(dat,*) -qy, -qx, result(1)      
        call flush(dat)
      end do
    end do
  
    close(dat)
    open(log, file="susy_qpi.log", position="append", status="old")
    write(log,*) "END: "//filename
    call cpu_time(end)
    write(log,*) "TOTAL TIME: ", end-start
    close(log)
  