--- trunk/libf/filtrez/inifgn.f	2010/03/05 16:43:45	25
+++ trunk/Sources/filtrez/inifgn.f	2015/08/24 16:30:33	167
@@ -1,102 +1,80 @@
-!
-! $Header: /home/cvsroot/LMDZ4/libf/filtrez/inifgn.F,v 1.1.1.1 2004/05/19 12:53:09 lmdzadmin Exp $
-!
-      SUBROUTINE inifgn(dv)
-c  
-c    ...  H.Upadyaya , O.Sharma  ... 
-c
-      use dimens_m
-      use paramet_m
-      use comgeom
-      use serre
-            use coefils
-      IMPLICIT NONE
-c
-
-c
-      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
-C
-c
-      EXTERNAL SSUM, acc, jacobi
-CC      EXTERNAL eigen
-      REAL SSUM
-c
-
-      imm1  = iim -1
-      pi = 2.* ASIN(1.)
-C
-      DO 5 i=1,iim
-       dlonu(i)=  xprimu( i )
-       dlonv(i)=  xprimv( i )
-   5  CONTINUE
-
-      DO 12 i=1,iim
-      sddv(i)   = SQRT(dlonv(i))
-      sddu(i)   = SQRT(dlonu(i))
-      unsddu(i) = 1./sddu(i)
-      unsddv(i) = 1./sddv(i)
-  12  CONTINUE
-C
-      DO 17 j=1,iim
-      DO 17 i=1,iim
-      vec(i,j)     = 0.
-      vec1(i,j)    = 0.
-      eignfnv(i,j) = 0.
-      eignfnu(i,j) = 0.
-  17  CONTINUE
-c
-c
-      eignfnv(1,1)    = -1.
-      eignfnv(iim,1)  =  1.
-      DO 20 i=1,imm1
-      eignfnv(i+1,i+1)= -1.
-      eignfnv(i,i+1)  =  1.
-  20  CONTINUE
-      DO 25 j=1,iim
-      DO 25 i=1,iim
-      eignfnv(i,j) = eignfnv(i,j)/(sddu(i)*sddv(j))
-  25  CONTINUE
-      DO 30 j=1,iim
-      DO 30 i=1,iim
-      eignfnu(i,j) = -eignfnv(j,i)
-  30  CONTINUE
-c
-      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)
-       ENDDO
-      ENDDO
-      ENDDO
-
-c
-      CALL jacobi(vec,iim,iim,dv,eignfnv,nrot)
-      CALL acc(eignfnv,d,iim)
-      CALL eigen_sort(dv,eignfnv,iim,iim)
-c
-      CALL jacobi(vec1,iim,iim,du,eignfnu,nrot)
-      CALL acc(eignfnu,d,iim)
-      CALL eigen_sort(du,eignfnu,iim,iim)
-
-cc   ancienne version avec appels IMSL
-c
-c     CALL MXM(eignfnu,iim,eignfnv,iim,vec,iim)
-c     CALL MXM(eignfnv,iim,eignfnu,iim,vec1,iim)
-c     CALL EVCSF(iim,vec,iim,dv,eignfnv,iim)
-c     CALL acc(eignfnv,d,iim)
-c     CALL eigen(eignfnv,dv)
-c
-c     CALL EVCSF(iim,vec1,iim,du,eignfnu,iim)
-c     CALL acc(eignfnu,d,iim)
-c     CALL eigen(eignfnu,du)
+module inifgn_m
 
-      RETURN
-      END
+  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(eignval_v, eignfnu, eignfnv)
+
+    ! From LMDZ4/libf/filtrez/inifgn.F, v 1.1.1.1 2004/05/19 12:53:09
+
+    ! 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 dynetat0_m, ONLY: xprimu, xprimv
+    use jumble, only: new_unit
+    use numer_rec_95, only: jacobi, eigsrt
+
+    real, intent(out):: eignval_v(:) ! (iim)
+    ! eigenvalues sorted in descending order
+
+    real, intent(out):: eignfnu(:, :), eignfnv(:, :) ! (iim, iim) eigenvectors
+
+    ! Local:
+
+    REAL delta(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, unit
+
+    !----------------------------------------------------------------
+
+    print *, "Call sequence information: inifgn"
+
+    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)
+
+    delta = matmul(deriv_v, deriv_u) ! second derivative at v longitudes
+    CALL jacobi(delta, eignval_v, eignfnv)
+    CALL acc(eignfnv)
+    CALL eigsrt(eignval_v, eignfnv)
+
+    delta = matmul(deriv_u, deriv_v) ! second derivative at u longitudes
+    CALL jacobi(delta, eignval_u, eignfnu)
+    CALL acc(eignfnu)
+    CALL eigsrt(eignval_u, eignfnu)
+
+    call new_unit(unit)
+    open(unit, file = "inifgn_out.txt", status = "replace", action = "write") 
+    write(unit, fmt = *) '"eignval_v"', eignval_v
+    close(unit)
+
+  END SUBROUTINE inifgn
+
+end module inifgn_m