--- trunk/filtrez/inifgn.f	2015/01/28 16:10:02	121
+++ trunk/Sources/filtrez/inifgn.f	2015/07/07 17:49:23	154
@@ -1,86 +1,73 @@
 module inifgn_m
 
+  use dimens_m, only: iim
+
   IMPLICIT NONE
 
+  private iim
+
+  real sddu(iim), sddv(iim) 
+  ! sdd[uv] = sqrt(2 pi / iim * (derivative of the longitudinal zoom
+  ! function)(rlon[uv]))
+
+  real unsddu(iim), unsddv(iim)
+
 contains
 
-  SUBROUTINE inifgn(dv)
+  SUBROUTINE inifgn(eignval_v, eignfnu, eignfnv)
 
     ! From LMDZ4/libf/filtrez/inifgn.F, v 1.1.1.1 2004/05/19 12:53:09
 
-    ! H.Upadyaya, O.Sharma 
+    ! Authors: H. Upadyaya, O. Sharma
 
+    ! Computes the eigenvalues and eigenvectors of the discrete analog
+    ! of the second derivative with respect to longitude.
+
+    use acc_m, only: acc 
     USE dimens_m, ONLY: iim
-    USE comgeom, ONLY: xprimu, xprimv
-    USE coefils, ONLY: eignfnu, eignfnv, sddu, sddv, unsddu, unsddv
+    USE dynetat0_m, ONLY: xprimu, xprimv
+    use numer_rec_95, only: jacobi, eigsrt
+
+    real, intent(out):: eignval_v(:) ! (iim)
+    ! eigenvalues sorted in descending order
 
-    real dv(iim)
+    real, intent(out):: eignfnu(:, :), eignfnv(:, :) ! (iim, iim) eigenvectors
 
     ! Local:
-    REAL vec(iim, iim), vec1(iim, iim)
-    REAL du(iim)
-    real d(iim)
-    REAL pi
-    INTEGER i, j, k, imm1, nrot
 
-    EXTERNAL acc, jacobi
+    REAL a(iim, iim) ! second derivative, symmetric, elements are angle^{-2}
+
+    REAL deriv_u(iim, iim), deriv_v(iim, iim) 
+    ! first derivative at u and v longitudes, elements are angle^{-1}
+
+    REAL eignval_u(iim)
+    INTEGER i
 
     !----------------------------------------------------------------
 
-    imm1 = iim - 1
-    pi = 2.*asin(1.)
+    print *, "Call sequence information: inifgn"
 
-    DO i = 1, iim
-       sddv(i) = sqrt(xprimv(i))
-       sddu(i) = sqrt(xprimu(i))
-       unsddu(i) = 1./sddu(i)
-       unsddv(i) = 1./sddv(i)
-    END DO
-
-    DO j = 1, iim
-       DO i = 1, iim
-          vec(i, j) = 0.
-          vec1(i, j) = 0.
-          eignfnv(i, j) = 0.
-          eignfnu(i, j) = 0.
-       END DO
-    END DO
-
-    eignfnv(1, 1) = -1.
-    eignfnv(iim, 1) = 1.
-    DO i = 1, imm1
-       eignfnv(i+1, i+1) = -1.
-       eignfnv(i, i+1) = 1.
-    END DO
-    DO j = 1, iim
-       DO i = 1, iim
-          eignfnv(i, j) = eignfnv(i, j)/(sddu(i)*sddv(j))
-       END DO
-    END DO
-    DO j = 1, iim
-       DO i = 1, iim
-          eignfnu(i, j) = -eignfnv(j, i)
-       END DO
-    END DO
-
-    DO j = 1, iim
-       DO i = 1, iim
-          vec(i, j) = 0.0
-          vec1(i, j) = 0.0
-          DO k = 1, iim
-             vec(i, j) = vec(i, j) + eignfnu(i, k)*eignfnv(k, j)
-             vec1(i, j) = vec1(i, j) + eignfnv(i, k)*eignfnu(k, j)
-          END DO
-       END DO
-    END DO
-
-    CALL jacobi(vec, iim, iim, dv, eignfnv, nrot)
-    CALL acc(eignfnv, d, iim)
-    CALL eigen_sort(dv, eignfnv, iim, iim)
-
-    CALL jacobi(vec1, iim, iim, du, eignfnu, nrot)
-    CALL acc(eignfnu, d, iim)
-    CALL eigen_sort(du, eignfnu, iim, iim)
+    sddv = sqrt(xprimv(:iim))
+    sddu = sqrt(xprimu(:iim))
+    unsddu = 1. / sddu
+    unsddv = 1. / sddv
+
+    deriv_u = 0.
+    deriv_u(iim, 1) = unsddu(iim) * unsddv(1)
+    forall (i = 1:iim) deriv_u(i, i) = - unsddu(i) * unsddv(i)
+    forall (i = 1:iim - 1) deriv_u(i, i + 1) = unsddu(i) * unsddv(i + 1)
+
+    deriv_v = - transpose(deriv_u)
+
+    a = matmul(deriv_v, deriv_u) ! second derivative at v longitudes
+    CALL jacobi(a, eignval_v, eignfnv)
+    CALL acc(eignfnv)
+    CALL eigsrt(eignval_v, eignfnv)
+
+    a = matmul(deriv_u, deriv_v) ! second derivative at u longitudes
+    CALL jacobi(a, eignval_u, eignfnu)
+    CALL acc(eignfnu)
+    CALL eigsrt(eignval_u, eignfnu)
 
   END SUBROUTINE inifgn