subroutine form_cF_nc(taus,Foo,Fov,Fvv,cFoo,cFov,cFvv)

! Compute F terms in CCSD

  implicit none

! Input variables

  double precision,intent(in)   :: taus(spin_occ_num,spin_occ_num,spin_vir_num,spin_vir_num)

  double precision,intent(in)   :: Foo(spin_occ_num,spin_occ_num)
  double precision,intent(in)   :: Fov(spin_occ_num,spin_vir_num)
  double precision,intent(in)   :: Fvv(spin_vir_num,spin_vir_num)

! Local variables

  integer                       :: i,j,m,n
  integer                       :: a,b,e,f
  double precision,external     :: Kronecker_Delta

! Output variables

  double precision,intent(out)  :: cFoo(spin_occ_num,spin_occ_num)
  double precision,intent(out)  :: cFov(spin_occ_num,spin_vir_num)
  double precision,intent(out)  :: cFvv(spin_vir_num,spin_vir_num)

! Occupied-occupied block

  do m=1,spin_occ_num
    do i=1,spin_occ_num

      cFoo(m,i) = (1d0 - Kronecker_delta(m,i))*Foo(m,i) 

      do e=1,spin_vir_num
        cFoo(m,i) = cFoo(m,i) + 0.5d0*t1_cc(i,e)*Fov(m,e)
      end do

      do e=1,spin_vir_num
        do n=1,spin_occ_num
          cFoo(m,i) = cFoo(m,i) + t1_cc(n,e)*OOOV(m,n,i,e)
        end do
      end do

      do e=1,spin_vir_num
        do n=1,spin_occ_num
          do f=1,spin_vir_num
            cFoo(m,i) = cFoo(m,i) + 0.5d0*taus(i,n,e,f)*OOVV(m,n,e,f)
          end do
        end do
      end do

    end do
  end do

! Occupied-virtual block

  cFov(:,:) = Fov(:,:)

  do m=1,spin_occ_num
    do e=1,spin_vir_num
      
      do n=1,spin_occ_num
        do f=1,spin_vir_num
          cFov(m,e) = cFov(m,e) + t1_cc(n,f)*OOVV(m,n,e,f)
        end do
      end do

    end do
  end do

! Virtual-virtual block

  do a=1,spin_vir_num
    do e=1,spin_vir_num

      cFvv(a,e) = (1d0 - Kronecker_delta(a,e))*Fvv(a,e) 
 
      do m=1,spin_occ_num
        cFvv(a,e) = cFvv(a,e) - 0.5d0*t1_cc(m,a)*Fov(m,e)
      end do

      do m=1,spin_occ_num
        do f=1,spin_vir_num
        cFvv(a,e) = cFvv(a,e) + t1_cc(m,f)*OVVV(m,a,f,e)
        end do
      end do

      do m=1,spin_occ_num
        do n=1,spin_occ_num
          do f=1,spin_vir_num
           cFvv(a,e) = cFvv(a,e) - 0.5d0*taus(m,n,a,f)*OOVV(m,n,e,f)
          end do
        end do
      end do

    end do
  end do

end subroutine form_cF_nc