Index: NEMO/branches/UKMO/dev_r9950_old_tidal_mixing/src/TOP/PISCES/SED/sedmat.F90
===================================================================
--- NEMO/branches/UKMO/dev_r9950_old_tidal_mixing/src/TOP/PISCES/SED/sedmat.F90	(revision 9989)
+++ NEMO/branches/UKMO/dev_r9950_old_tidal_mixing/src/TOP/PISCES/SED/sedmat.F90	(revision 10324)
@@ -1,4 +1,3 @@
 MODULE sedmat
-#if defined key_sed
    !!======================================================================
    !!              ***  MODULE  sedmat  ***
@@ -9,4 +8,6 @@
 
    USE sed     ! sediment global variable
+   USE lib_mpp         ! distribued memory computing library
+
 
    IMPLICIT NONE
@@ -25,5 +26,5 @@
  CONTAINS
 
-    SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psol )
+    SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in )
        !!---------------------------------------------------------------------
        !!                  ***  ROUTINE sed_mat_dsr  ***
@@ -48,10 +49,12 @@
        !!----------------------------------------------------------------------
        !! * Arguments
-       INTEGER , INTENT(in) ::  nvar  ! number of variables
+       INTEGER , INTENT(in) ::  nvar  ! number of variable
        INTEGER , INTENT(in) ::  ndim  ! number of points
        INTEGER , INTENT(in) ::  nlev  ! number of sediment levels
 
-       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(in   ) :: preac  ! reaction rates
-       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: psol   ! solution ( undersaturation values )
+       REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: preac  ! reaction rates
+       REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: psms  ! reaction rates
+       REAL(wp), DIMENSION(ndim,nlev), INTENT(inout) :: psol   ! solution ( undersaturation values )
+       REAL(wp), INTENT(in) ::  dtsed_in
  
        !---Local declarations
@@ -67,4 +70,5 @@
        !----------------------------------------------------------------------
 
+       IF( ln_timing )  CALL timing_start('sed_mat_dsr')
 
        ! Computation left hand side of linear system of 
@@ -73,10 +77,11 @@
 
 
+       jn = nvar
        ! first sediment level          
        DO ji = 1, ndim
-          aplus  = ( ( volw3d(ji,1) / dz3d(ji,1) ) + &
-                     ( volw3d(ji,2) / dz3d(ji,2) ) ) / 2.
+          aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + &
+                        ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2.
           dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
-          rplus  = ( dtsed / volw3d(ji,1) ) * diff(1) * aplus / dxplus 
+          rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 
 
           za(ji,1) = 0.
@@ -87,15 +92,15 @@
        DO jk = 2, nlev - 1
           DO ji = 1, ndim
-             aminus  = ( ( volw3d(ji,jk-1) / dz3d(ji,jk-1) ) + &
-                &        ( volw3d(ji,jk  ) / dz3d(ji,jk  ) ) ) / 2.
+             aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + &
+             &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2.
              dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.
 
-             aplus   = ( ( volw3d(ji,jk  ) / dz3d(ji,jk  ) ) + &
-                &        ( volw3d(ji,jk+1) / dz3d(ji,jk+1) ) ) / 2.
+             aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + &
+             &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2.
              dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
-             !
-             rminus  = ( dtsed / volw3d(ji,jk) ) * diff(jk-1) * aminus / dxminus
-             rplus   = ( dtsed / volw3d(ji,jk) ) * diff(jk)   * aplus / dxplus
-             !     
+                !
+             rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus
+             rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus / dxplus
+                !     
              za(ji,jk) = -rminus
              zb(ji,jk) = 1. + rminus + rplus 
@@ -105,8 +110,8 @@
 
        DO ji = 1, ndim
-          aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1) ) + &
-             &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2.
+          aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1)  ) + &
+          &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2.
           dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2.
-          rminus  = ( dtsed / volw3d(ji,nlev) ) * diff(nlev-1) * aminus / dxminus
+          rminus  = ( dtsed_in / volw3d(ji,nlev) ) * diff(ji,nlev-1,jn) * aminus / dxminus
           !
           za(ji,nlev) = -rminus
@@ -119,31 +124,30 @@
        ! -----------------------------------------------
 
-       DO jn = 1, nvar
-
-          zr  (:,:)    = psol(:,:,jn)
-          zbet(:  )    = zb(:,1) + preac(:,1,jn)
-          psol(:,1,jn) = zr(:,1) / zbet(:)
+       zr  (:,:) = psol(:,:) + psms(:,:)
+       zb  (:,:) = zb(:,:) + preac(:,:)
+       zbet(:  ) = zb(:,1)
+       psol(:,1) = zr(:,1) / zbet(:)
+
           ! 
-          DO jk = 2, nlev
-             DO ji = 1, ndim
-                zgamm(ji,jk)   =  zc(ji,jk-1) / zbet(ji)
-                zbet(ji)       = ( zb(ji,jk) + preac(ji,jk,jn) ) - za(ji,jk) * zgamm(ji,jk)
-                psol(ji,jk,jn) = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1,jn) ) / zbet(ji)
-             END DO
-          ENDDO
+       DO jk = 2, nlev
+          DO ji = 1, ndim
+             zgamm(ji,jk) =  zc(ji,jk-1) / zbet(ji)
+             zbet(ji)     =  zb(ji,jk) - za(ji,jk) * zgamm(ji,jk)
+             psol(ji,jk)  = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1) ) / zbet(ji)
+          END DO
+       ENDDO
           ! 
-          DO jk = nlev - 1, 1, -1
-             DO ji = 1,ndim
-                psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(ji,jk+1) * psol(ji,jk+1,jn)
-             END DO
-          ENDDO
-
+       DO jk = nlev - 1, 1, -1
+          DO ji = 1,ndim
+             psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1)
+          END DO
        ENDDO
 
+       IF( ln_timing )  CALL timing_stop('sed_mat_dsr')
+
 
     END SUBROUTINE sed_mat_dsr
 
-    
-    SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol )
+    SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol, dtsed_in )
        !!---------------------------------------------------------------------
        !!                  ***  ROUTINE sed_mat_btb  ***
@@ -170,4 +174,6 @@
           psol      ! solution
 
+      REAL(wp), INTENT(in) :: dtsed_in
+
        !---Local declarations
        INTEGER  ::  &
@@ -192,76 +198,68 @@
 
 
+      IF( ln_timing )  CALL timing_start('sed_mat_btb')
+
        ! first sediment level          
-       aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2.
-       dxplus = ( dz(2) + dz(3) ) / 2.
-       rplus  = ( dtsed / vols(2) ) * db * aplus / dxplus 
-
-       za(1) = 0.
-       zb(1) = 1. + rplus
-       zc(1) = -rplus
+      DO ji = 1, ndim
+         aplus  = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2.
+         dxplus = ( dz(2) + dz(3) ) / 2.
+         rplus  = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 
+
+         za(1) = 0.
+         zb(1) = 1. + rplus
+         zc(1) = -rplus
 
              
-       DO jk = 2, nlev - 1
-          aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2.
-          dxminus = ( dz(jk) + dz(jk+1) ) / 2.
-          rminus  = ( dtsed / vols(jk+1) ) * db * aminus / dxminus
-          !
-          aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2.
-          dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2.
-          rplus   = ( dtsed / vols(jk+1) ) * db * aplus / dxplus
-          !     
-          za(jk) = -rminus
-          zb(jk) = 1. + rminus + rplus 
-          zc(jk) = -rplus
-       ENDDO
- 
-       aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2.
-       dxminus = ( dz(nlev) + dz(nlev+1) ) / 2.
-       rminus  = ( dtsed / vols(nlev+1) ) * db * aminus / dxminus
+         DO jk = 2, nlev - 1
+            aminus  = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2.
+            dxminus = ( dz(jk) + dz(jk+1) ) / 2.
+            rminus  = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus
+            !
+            aplus   = ( ( vols(jk+1) / dz(jk+1  ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2.
+            dxplus  = ( dz(jk+1) + dz(jk+2) ) / 2.
+            rplus   = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus
+            !     
+            za(jk) = -rminus
+            zb(jk) = 1. + rminus + rplus 
+            zc(jk) = -rplus
+         ENDDO
+ 
+         aminus  = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2.
+         dxminus = ( dz(nlev) + dz(nlev+1) ) / 2.
+         rminus  = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus
+         !
+         za(nlev) = -rminus
+         zb(nlev) = 1. + rminus
+         zc(nlev) = 0.
+
+
+         ! solves tridiagonal system of linear equations 
+         ! -----------------------------------------------    
+         DO jn = 1, nvar
+          
+            DO jk = 1, nlev
+               zr  (ji,jk)    = psol(ji,jk,jn)
+            END DO
+            zbet          = zb(1)
+            psol(ji,1,jn) = zr(ji,1) / zbet
+            ! 
+            DO jk = 2, nlev
+               zgamm(jk) =  zc(jk-1) / zbet
+               zbet      =  zb(jk) - za(jk) * zgamm(jk)
+               psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet
+            ENDDO
+            ! 
+            DO jk = nlev - 1, 1, -1
+                psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn)
+            ENDDO
+
+         ENDDO
+
+       END DO
        !
- 
-       za(nlev) = -rminus
-       zb(nlev) = 1. + rminus
-       zc(nlev) = 0.
-
-
-       ! solves tridiagonal system of linear equations 
-       ! -----------------------------------------------    
-       DO jn = 1, nvar
-          
-          zr  (:,:)    = psol(:,:,jn)
-          zbet         = zb(1)
-          psol(:,1,jn) = zr(:,1) / zbet
-          ! 
-          DO jk = 2, nlev
-             zgamm(jk) =  zc(jk-1) / zbet
-             zbet      =  zb(jk) - za(jk) * zgamm(jk)
-             DO ji = 1, ndim
-                psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet
-             END DO
-          ENDDO
-          ! 
-          DO jk = nlev - 1, 1, -1
-             DO ji = 1,ndim
-                psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn)
-             END DO
-          ENDDO
-
-       ENDDO
+       IF( ln_timing )  CALL timing_stop('sed_mat_btb')
 
        
     END SUBROUTINE sed_mat_btb
 
-
-#else
-   !!======================================================================
-   !! MODULE sedmat  :   Dummy module 
-   !!======================================================================
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_mat         ! Empty routine
-   END SUBROUTINE sed_mat
-   !!======================================================================
-#endif
-
  END MODULE sedmat
