subroutine IF_calc(x1, x2, dx1, dx2, dt, imx, t21, tau21l, tau21b)
  ! FORTRAN souce code for computating Liang-Kleeman information flow
  !
  ! Input variables:
  !   imx: dimension of data
  !   dt: time interval
  !   x1: X1 
  !   x2: X2 
  !   dx1: time derivative of X1 using Euler forward scheme
  !   dx2: time derivative of X2 using Euler forward scheme
  !
  ! Output variables:
  !   t21: information flow from X2 to X2 (without normalization) 
  !   tau21l: Normalized information flow from X2 to X2 (Liang 2015)
  !
  ! References:
  !   X. S. Liang, Normalizing the causality between time series, Phys. Rev. E 92, 022126 (2015)
  
    implicit none

    integer :: i, n
    double precision:: c00, c01, c10, c11
    double precision:: dc00, dc01, dc10, dc11
    double precision:: x1_mean, x2_mean, dx1_mean, dx2_mean
    double precision:: detc, a11, a12, f1, Q1, b1, dH1_star, dH1_noise
    integer, intent(in) :: imx
    integer, intent(in) :: dt
    double precision, dimension(imx), intent(in) :: x1, x2, dx1, dx2
    double precision, intent(out) :: t21
    double precision, intent(out) :: tau21l
    double precision, intent(out) :: tau21b
    double precision, parameter :: undef=-1.e+10

    ! -- initialize
    c00 = 0.
    c01 = 0.
    c10 = 0.
    c11 = 0
    dc00 = 0.
    dc01 = 0.
    dc10 = 0.
    dc11 = 0
    x1_mean = 0.
    x2_mean = 0.
    dx1_mean = 0.
    dx2_mean = 0
    b1 = 0.
    q1 = 0.

    ! -- initialize
    n = 0
    do i = 1, imx
       if (.not. isnan(x1(i)) .and. .not. isnan(x2(i)) .and.&
          &.not. isnan(dx1(i)) .and. .not. isnan(dx2(i)))then 
           n = n + 1
           x1_mean = x1_mean + x1(i) 
           x2_mean = x2_mean + x2(i) 
           dx1_mean = dx1_mean + dx1(i) 
           dx2_mean = dx2_mean + dx2(i) 
        end if
    end do
    x1_mean = x1_mean / real(n)
    x2_mean = x2_mean / real(n)
    dx1_mean = dx1_mean / real(n)
    dx2_mean = dx2_mean / real(n)

    do i = 1, imx
       if (.not. isnan(x1(i)) .and. .not. isnan(x2(i)) .and.&
          &.not. isnan(dx1(i)) .and. .not. isnan(dx2(i)))then 
            c00 = c00 + ( (x1(i) - x1_mean) * (x1(i) - x1_mean) ) / real(n)
            c01 = c01 + ( (x1(i) - x1_mean) * (x2(i) - x2_mean) ) / real(n)
            c10 = c10 + ( (x1(i) - x1_mean) * (x2(i) - x2_mean) ) / real(n)
            c11 = c11 + ( (x2(i) - x2_mean) * (x2(i) - x2_mean) ) / real(n)

            dc00 = dc00 + ( (x1(i) - x1_mean) * (dx1(i) - dx1_mean) ) / real(n)
            dc01 = dc01 + ( (x1(i) - x1_mean) * (dx2(i) - dx2_mean) ) / real(n)
            dc10 = dc10 + ( (x2(i) - x2_mean) * (dx1(i) - dx1_mean) ) / real(n)
            dc11 = dc11 + ( (x2(i) - x2_mean) * (dx2(i) - dx2_mean) ) / real(n)
        end if
    end do

    detc = c00 * c11 - c01 * c10
    a11 = (c11 * dc00 - c01 * dc10) / detc
    a12 = (real(-1) * c01 * dc00 + c00 * dc10) / detc

    f1 = dx1_mean - a11 * x1_mean - a12 * x2_mean
    do i = 1, imx
       if (.not. isnan(x1(i)) .and. .not. isnan(x2(i)) .and.&
          &.not. isnan(dx1(i)) .and. .not. isnan(dx2(i)))then 
           Q1 = Q1 + ( dx1(i) - ( f1 + a11 * x1(i) + a12 * x2(i)) )**2.
       end if
    end do
    b1 = sqrt(Q1 * dt / real(n))

    t21 = c01/c00 * (-1. * c10 * dc00 + c00 * dc10) / detc

    dH1_star = a11
    dH1_noise = (b1**2) / (real(2) * c00)

    tau21l = t21 / (abs(t21) + abs(dH1_star) + abs(dh1_noise))

    tau21b = t21 / (abs(dH1_noise) + abs(t21))

end subroutine