Index: trunk/LMDZ.GENERIC/libf/phygeneric/dtridgl.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/dtridgl.F	(revision 4083)
+++ 	(revision )
@@ -1,36 +1,0 @@
-      SUBROUTINE DTRIDGL(L,AF,BF,CF,DF,XK)
-
-C  GCM2.0  Feb 2003
-
-C     DOUBLE PRECISION VERSION OF TRIDGL
-
-!!      PARAMETER (NMAX=201)
-      IMPLICIT REAL*8  (A-H,O-Z)
-      DIMENSION AF(L),BF(L),CF(L),DF(L),XK(L)
-      DIMENSION AS(2*L),DS(2*L)
-
-C*    THIS SUBROUTINE SOLVES A SYSTEM OF TRIDIAGIONAL MATRIX
-C*    EQUATIONS. THE FORM OF THE EQUATIONS ARE:
-C*    A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = D(I)
-C*    WHERE I=1,L  LESS THAN 103.
-C* ..............REVIEWED -CP........
-
-C======================================================================C
-
-      AS(L) = AF(L)/BF(L)
-      DS(L) = DF(L)/BF(L)
-
-      DO I=2,L
-        X         = 1./(BF(L+1-I) - CF(L+1-I)*AS(L+2-I))
-        AS(L+1-I) = AF(L+1-I)*X
-        DS(L+1-I) = (DF(L+1-I)-CF(L+1-I)*DS(L+2-I))*X
-      END DO
- 
-      XK(1)=DS(1)
-      DO I=2,L
-        XKB   = XK(I-1)
-        XK(I) = DS(I)-AS(I)*XKB
-      END DO
-
-      RETURN
-      END
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_fluxes_solver_stellar.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_fluxes_solver_stellar.F	(revision 4083)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_fluxes_solver_stellar.F	(revision 4097)
@@ -45,4 +45,6 @@
 
       use radinc_h, only: L_TAUMAX, L_NLAYRAD, L_NLEVRAD, L_LEVELS
+      use rad_tridiagonal_matrix_solver_mod, only:
+     &                    rad_tridiagonal_matrix_solver
 
       implicit none
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_fluxes_solver_thermal.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_fluxes_solver_thermal.F	(revision 4083)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_correlatedk_fluxes_solver_thermal.F	(revision 4097)
@@ -13,4 +13,6 @@
       use radcommon_h, only: planckir
       use comcstfi_mod, only: pi
+      use rad_tridiagonal_matrix_solver_mod, only:
+     &                        rad_tridiagonal_matrix_solver
       
       IMPLICIT NONE
Index: trunk/LMDZ.GENERIC/libf/phygeneric/rad_tridiagonal_matrix_solver.F
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/rad_tridiagonal_matrix_solver.F	(revision 4083)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/rad_tridiagonal_matrix_solver.F	(revision 4097)
@@ -1,15 +1,27 @@
+      MODULE rad_tridiagonal_matrix_solver_mod
+      
+      IMPLICIT NONE
+      
+      CONTAINS
+      
       SUBROUTINE rad_tridiagonal_matrix_solver(NL,
      *           GAMA,CP,CM,CPM1,CMM1,E1,E2,E3,E4,
      *           BTOP,BSURF,RSF,XK1,XK2)
 
+      USE util_tridiagonal_matrix_solver_mod, ONLY: DTRIDGL
 C  GCM2.0  Feb 2003
 C
 C DOUBLE PRECISION VERSION OF SOLVER
 
-!!      PARAMETER (NMAX=201)
-      IMPLICIT REAL*8  (A-H,O-Z)
-      DIMENSION GAMA(NL),CP(NL),CM(NL),CPM1(NL),CMM1(NL),XK1(NL),
-     *          XK2(NL),E1(NL),E2(NL),E3(NL),E4(NL)
-      DIMENSION AF(2*NL),BF(2*NL),CF(2*NL),DF(2*NL),XK(2*NL)
+      IMPLICIT NONE
+
+      INTEGER,INTENT(IN) :: NL
+      REAL*8,INTENT(IN) :: GAMA(NL),CP(NL),CM(NL)
+      REAL*8,INTENT(IN) :: CPM1(NL),CMM1(NL)
+      REAL*8,INTENT(OUT) :: XK1(NL), XK2(NL)
+      REAL*8,INTENT(IN) :: E1(NL),E2(NL),E3(NL),E4(NL)
+      REAL*8,INTENT(IN) :: BTOP, BSURF, RSF
+      REAL*8 :: AF(2*NL),BF(2*NL),CF(2*NL),DF(2*NL),XK(2*NL)
+
 C*********************************************************
 C* THIS SUBROUTINE SOLVES FOR THE COEFFICIENTS OF THE    *
@@ -38,4 +50,7 @@
 C*********************************************************
 
+      INTEGER :: L,LM1,LM2
+      INTEGER :: N
+      INTEGER :: I
 C======================================================================C
 
@@ -83,5 +98,5 @@
 C     ***UNMIX THE COEFFICIENTS****
 
-      DO 28 N=1,NL
+      DO N=1,NL
         XK1(N) = XK(2*N-1)+XK(2*N)
         XK2(N) = XK(2*N-1)-XK(2*N)
@@ -92,12 +107,12 @@
 C       BE TREATED CAREFULLY
 
-        IF(XK2(N) .EQ. 0.0) GO TO 28
+        IF(XK2(N) .EQ. 0.0) CONTINUE
 c        IF (ABS (XK2(N)/XK(2*N-1)) .LT. 1.E-30) XK2(N)=0.0
 
         IF (ABS (XK2(N)/(XK(2*N-1)+1.e-20)) .LT. 1.E-30) XK2(N)=0.0   ! For debug only (with -Ktrap=fp option)
 
+      ENDDO
+ 
+      END SUBROUTINE rad_tridiagonal_matrix_solver
 
-   28 CONTINUE
- 
-      RETURN
-      END
+      END MODULE rad_tridiagonal_matrix_solver_mod
Index: trunk/LMDZ.GENERIC/libf/phygeneric/util_tridiagonal_matrix_solver.F90
===================================================================
--- trunk/LMDZ.GENERIC/libf/phygeneric/util_tridiagonal_matrix_solver.F90	(revision 4097)
+++ trunk/LMDZ.GENERIC/libf/phygeneric/util_tridiagonal_matrix_solver.F90	(revision 4097)
@@ -0,0 +1,45 @@
+module util_tridiagonal_matrix_solver_mod
+
+implicit none
+
+contains
+
+      SUBROUTINE DTRIDGL(L,AF,BF,CF,DF,XK)
+
+!  GCM2.0  Feb 2003
+
+!     DOUBLE PRECISION VERSION OF TRIDGL
+
+      REAL*8,INTENT(IN) :: AF(L),BF(L),CF(L),DF(L)
+      REAL*8,INTENT(OUT) :: XK(L)
+      INTEGER,INTENT(IN) :: L
+
+!*    THIS SUBROUTINE SOLVES A SYSTEM OF TRIDIAGIONAL MATRIX
+!*    EQUATIONS. THE FORM OF THE EQUATIONS ARE:
+!*    A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = D(I)
+!* ..............REVIEWED -CP........
+
+      REAL*8 :: AS(2*L),DS(2*L)
+      REAL*8 :: XKB,X
+      INTEGER :: I
+
+!======================================================================C
+
+      AS(L) = AF(L)/BF(L)
+      DS(L) = DF(L)/BF(L)
+
+      DO I=2,L
+        X         = 1./(BF(L+1-I) - CF(L+1-I)*AS(L+2-I))
+        AS(L+1-I) = AF(L+1-I)*X
+        DS(L+1-I) = (DF(L+1-I)-CF(L+1-I)*DS(L+2-I))*X
+      END DO
+ 
+      XK(1)=DS(1)
+      DO I=2,L
+        XKB   = XK(I-1)
+        XK(I) = DS(I)-AS(I)*XKB
+      END DO
+
+      END SUBROUTINE DTRIDGL
+
+end module util_tridiagonal_matrix_solver_mod