--- trunk/filtrez/inifgn.f90	2014/03/05 14:38:41	81
+++ trunk/Sources/filtrez/inifgn.f	2015/06/23 18:26:18	153
@@ -1,106 +1,63 @@
+module inifgn_m
 
-! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifgn.F,v 1.1.1.1 2004/05/19
-! 12:53:09 lmdzadmin Exp $
+  use dimens_m, only: iim
 
-SUBROUTINE inifgn(dv)
+  IMPLICIT NONE
 
-  ! ...  H.Upadyaya , O.Sharma  ...
+  private iim
 
-  USE dimens_m
-  USE paramet_m
-  USE comgeom
-  USE serre
-  USE coefils
-  IMPLICIT NONE
+  real sddu(iim), sddv(iim) ! SQRT(dx / di)
+  real unsddu(iim), unsddv(iim)
+
+  real eignfnu(iim, iim), eignfnv(iim, iim)
+  ! eigenfunctions of the discrete laplacian
+
+contains
+
+  SUBROUTINE inifgn(dv)
+
+    ! From LMDZ4/libf/filtrez/inifgn.F, v 1.1.1.1 2004/05/19 12:53:09
+
+    ! H. Upadyaya, O. Sharma 
+
+    use acc_m, only: acc 
+    USE dimens_m, ONLY: iim
+    USE dynetat0_m, ONLY: xprimu, xprimv
+    use numer_rec_95, only: jacobi, eigsrt
+
+    real, intent(out):: dv(:) ! (iim) eigenvalues sorted in descending order
+
+    ! Local:
+    REAL, dimension(iim, iim):: a, b, c
+    REAL du(iim)
+    INTEGER i
+
+    !----------------------------------------------------------------
+
+    print *, "Call sequence information: inifgn"
+
+    sddv = sqrt(xprimv(:iim))
+    sddu = sqrt(xprimu(:iim))
+    unsddu = 1. / sddu
+    unsddv = 1. / sddv
+
+    b = 0.
+    b(iim, 1) = unsddu(iim) * unsddv(1)
+    forall (i = 1:iim) b(i, i) = - unsddu(i) * unsddv(i)
+    forall (i = 1:iim - 1) b(i, i + 1) = unsddu(i) * unsddv(i + 1)
 
+    c = - transpose(b)
 
+    a = matmul(c, b)
+    CALL jacobi(a, dv, eignfnv)
+    CALL acc(eignfnv)
+    CALL eigsrt(dv, eignfnv)
 
-  REAL vec(iim, iim), vec1(iim, iim)
-  REAL dlonu(iim), dlonv(iim)
-  REAL du(iim), dv(iim), d(iim)
-  REAL pi
-  INTEGER i, j, k, imm1, nrot
-
-
-  EXTERNAL ssum, acc, jacobi
-  ! C      EXTERNAL eigen
-  REAL ssum
-
-
-  imm1 = iim - 1
-  pi = 2.*asin(1.)
-
-  DO i = 1, iim
-    dlonu(i) = xprimu(i)
-    dlonv(i) = xprimv(i)
-  END DO
-
-  DO i = 1, iim
-    sddv(i) = sqrt(dlonv(i))
-    sddu(i) = sqrt(dlonu(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)
-
-  ! c   ancienne version avec appels IMSL
-
-  ! CALL MXM(eignfnu,iim,eignfnv,iim,vec,iim)
-  ! CALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim)
-  ! CALL EVCSF(iim,vec,iim,dv,eignfnv,iim)
-  ! CALL acc(eignfnv,d,iim)
-  ! CALL eigen(eignfnv,dv)
-
-  ! CALL EVCSF(iim,vec1,iim,du,eignfnu,iim)
-  ! CALL acc(eignfnu,d,iim)
-  ! CALL eigen(eignfnu,du)
+    a = matmul(b, c)
+    CALL jacobi(a, du, eignfnu)
+    CALL acc(eignfnu)
+    CALL eigsrt(du, eignfnu)
 
-  RETURN
-END SUBROUTINE inifgn
+  END SUBROUTINE inifgn
 
+end module inifgn_m