Reintroduce the difcub function to fleur

It is needed for the Wannier code.

math/CMakeLists.txt

@@ -22,6 +22,7 @@ math/sphbes.f
 math/sphpts.f
 math/util.F
 math/ylm4.f
+math/difcub.f
 )
 set(fleur_F90 ${fleur_F90}
 math/differentiate.f90

math/difcub.f

+!--------------------------------------------------------------------------------
+! Copyright (c) 2017 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
+! This file is part of FLEUR and available as free software under the conditions
+! of the MIT license as expressed in the LICENSE file in more detail.
+!--------------------------------------------------------------------------------
+
+      MODULE m_difcub
+
+      CONTAINS
+
+      REAL FUNCTION difcub(x,f,xi)
+c     **********************************************************
+c     differentiate the function f, given at the
+c     points x0,x1,x2,x3 at the point xi by lagrange
+c     interpolation for polynomial of 3rd order
+c     r.p.
+c     ***********************************************************
+      IMPLICIT NONE
+C     .. Scalar Arguments ..
+      REAL xi
+C     ..
+C     .. Array Arguments ..
+      REAL f(0:3),x(0:3)
+C     ..
+      difcub = ((xi-x(1))* (xi-x(2))+ (xi-x(1))* (xi-x(3))+
+     +         (xi-x(2))* (xi-x(3)))*f(0)/ ((x(0)-x(1))* (x(0)-x(2))*
+     +         (x(0)-x(3))) + ((xi-x(0))* (xi-x(2))+
+     +         (xi-x(0))* (xi-x(3))+ (xi-x(2))* (xi-x(3)))*f(1)/
+     +         ((x(1)-x(0))* (x(1)-x(2))* (x(1)-x(3))) +
+     +         ((xi-x(0))* (xi-x(1))+ (xi-x(0))* (xi-x(3))+
+     +         (xi-x(1))* (xi-x(3)))*f(2)/ ((x(2)-x(0))* (x(2)-x(1))*
+     +         (x(2)-x(3))) + ((xi-x(0))* (xi-x(1))+
+     +         (xi-x(0))* (xi-x(2))+ (xi-x(1))* (xi-x(2)))*f(3)/
+     +         ((x(3)-x(0))* (x(3)-x(1))* (x(3)-x(2)))
+      RETURN
+      END FUNCTION difcub
+
+      END MODULE m_difcub

wannier/uhu/d2fdz2.F

@@ -3,10 +3,9 @@
 
       SUBROUTINE d2fdz2(jmtd,jri,rmsh,dx,f,fac,d2f)
 
-      IMPLICIT NONE
+      USE m_difcub
 
-      REAL difcub
-      EXTERNAL difcub
+      IMPLICIT NONE
       
       REAL, INTENT(INOUT) :: d2f(jmtd)
       REAL, INTENT(IN)    :: fac(jmtd)

wannier/uhu/d2fdz2cmplx.F

@@ -3,10 +3,9 @@
 
       SUBROUTINE d2fdz2cmplx(jmtd,jri,rmsh,dx,f,fac,d2f)
 
-      IMPLICIT NONE
+      USE m_difcub
 
-      REAL difcub
-      EXTERNAL difcub
+      IMPLICIT NONE
       
       COMPLEX, INTENT(INOUT) :: d2f(jmtd)
       COMPLEX, INTENT(IN)    :: fac(jmtd)

wannier/uhu/dujdr.F

@@ -4,10 +4,9 @@
       SUBROUTINE dujdr(jmtd,jri,rmsh,dx,ub,j2,b,l,lmaxd,dub)
 
       USE m_constants
-      IMPLICIT NONE
+      USE m_difcub
 
-      REAL difcub
-      EXTERNAL difcub
+      IMPLICIT NONE
       
       REAL, INTENT(IN) :: rmsh(jmtd),dx
       REAL, INTENT(IN) :: ub(jmtd,2)       ! u(b2)

wannier/uhu/wann_uHu_radintsra.F

@@ -5,10 +5,9 @@
 
       USE m_intgr, ONLY : intgr3
       USE m_constants
-      IMPLICIT NONE
+      USE m_difcub
 
-      REAL difcub
-      EXTERNAL difcub
+      IMPLICIT NONE
       
       REAL, INTENT(IN) :: f(jmtd,2),g(jmtd,2)
       REAL, INTENT(IN) :: vr(jmtd)

wannier/uhu/wann_uHu_radintsra2.F

@@ -9,10 +9,9 @@
 
       USE m_intgr, ONLY : intgr3
       USE m_constants
-      IMPLICIT NONE
+      USE m_difcub
 
-      REAL difcub
-      EXTERNAL difcub
+      IMPLICIT NONE
       
       REAL, INTENT(IN) :: p(jmtd,2)
       REAL, INTENT(IN) :: q(jmtd,2)