Index: /NEMO/trunk/src/TOP/PISCES/SED/oce_sed.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/oce_sed.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/oce_sed.F90	(revision 10225)
@@ -0,0 +1,60 @@
+MODULE oce_sed
+   !!======================================================================
+   !!                        ***  sed  ***
+   !! Sediment :   set sediment global variables
+   !!======================================================================
+   !! History :
+   !!        !  06-12  (C. Ethe)  Orignal
+   !!----------------------------------------------------------------------
+   USE par_sed
+   USE timing
+   USE par_trc
+
+   USE dom_oce , ONLY :   nidom     =>   nidom          !:
+   USE dom_oce , ONLY :   glamt     =>   glamt          !: longitude of t-point (degre)
+   USE dom_oce , ONLY :   gphit     =>   gphit          !: latitude  of t-point (degre)
+   USE dom_oce , ONLY :   e3t_n     =>   e3t_n          !: latitude  of t-point (degre)
+   USE dom_oce , ONLY :   e3t_1d    =>   e3t_1d         !: reference depth of t-points (m)
+   USE dom_oce , ONLY :   gdepw_0   =>   gdepw_0        !: reference depth of t-points (m)
+   USE dom_oce , ONLY :   mbkt      =>   mbkt           !: vertical index of the bottom last T- ocean level
+   USE dom_oce , ONLY :   tmask     =>   tmask          !: land/ocean mask at t-points
+   USE dom_oce , ONLY :   rdt       =>   rdt            !: time step for the dynamics
+   USE dom_oce , ONLY :   nyear     =>   nyear          !: Current year
+   USE dom_oce , ONLY :   nmonth    =>   nmonth         !: Current month
+   USE dom_oce , ONLY :   nday      =>   nday           !: Current day
+   USE dom_oce , ONLY :   ndastp    =>   ndastp         !: time step date in year/month/day aammjj
+   USE dom_oce , ONLY :   nday_year =>   nday_year      !: curent day counted from jan 1st of the current year
+   USE dom_oce , ONLY :   adatrj    =>   adatrj         !: number of elapsed days since the begining of the run
+   USE trc     , ONLY :  nittrc000  =>   nittrc000
+   !                                !: it is the accumulated duration of previous runs
+   !                                !: that may have been run with different time steps.
+
+   USE oce     , ONLY :  tsn        =>   tsn             !: pot. temperature (celsius) and salinity (psu)
+   USE trc     , ONLY :  trb        =>   trb             !: pot. temperature (celsius) and salinity (psu)
+
+   USE sms_pisces, ONLY : wsbio4    =>   wsbio4          !: sinking flux for POC
+   USE sms_pisces, ONLY : wsbio3    =>   wsbio3          !: sinking flux for GOC
+   USE sms_pisces, ONLY : wscal     =>   wscal           !: sinking flux for calcite
+   USE sms_pisces, ONLY : wsbio2    =>   wsbio2           !: sinking flux for calcite
+   USE sms_pisces, ONLY : wsbio     =>   wsbio           !: sinking flux for calcite
+   USE p4zche, ONLY     : akb3      =>   akb3            !: Chemical constants  
+   USE sms_pisces, ONLY : ak13      =>   ak13            !: Chemical constants  
+   USE sms_pisces, ONLY : ak23      =>   ak23            !: Chemical constants  
+   USE p4zche, ONLY     : akw3      =>   akw3            !: Chemical constants  
+   USE sms_pisces, ONLY : aksp      =>   aksp            !: Chemical constants  
+   USE p4zche, ONLY     : borat     =>   borat           !: Chemical constants ( borat ) 
+   USE p4zche, ONLY     : ak1p3     =>   ak1p3           !: Chemical constants  
+   USE p4zche, ONLY     : ak2p3     =>   ak2p3           !: Chemical constants  
+   USE p4zche, ONLY     : ak3p3     =>   ak3p3           !: Chemical constants  
+   USE p4zche, ONLY     : aksi3     =>   aksi3           !: Chemical constants  
+   USE p4zche, ONLY     : aks3      =>   aks3            !: Chemical constants  
+   USE p4zche, ONLY     : akf3      =>   akf3            !: Chemical constants  
+   USE p4zche, ONLY     : fluorid   =>   fluorid         !: Chemical constants  
+   USE p4zche, ONLY     : sulfat    =>   sulfat          !: Chemical constants  
+   USE p4zche, ONLY     : sio3eq    =>   sio3eq          !: Chemical constants  
+   USE p4zsbc, ONLY     : dust      =>   dust
+   USE trc       , ONLY : r2dttrc   =>   r2dttrc
+
+END MODULE oce_sed
+
+
Index: /NEMO/trunk/src/TOP/PISCES/SED/p4zsed.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/p4zsed.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/p4zsed.F90	(revision 10225)
@@ -0,0 +1,517 @@
+MODULE p4zsed
+   !!======================================================================
+   !!                         ***  MODULE p4sed  ***
+   !! TOP :   PISCES Compute loss of organic matter in the sediments
+   !!======================================================================
+   !! History :   1.0  !  2004-03 (O. Aumont) Original code
+   !!             2.0  !  2007-12 (C. Ethe, G. Madec)  F90
+   !!             3.4  !  2011-06 (C. Ethe) USE of fldread
+   !!             3.5  !  2012-07 (O. Aumont) improvment of river input of nutrients 
+   !!----------------------------------------------------------------------
+   !!   p4z_sed        :  Compute loss of organic matter in the sediments
+   !!----------------------------------------------------------------------
+   USE oce_trc         !  shared variables between ocean and passive tracers
+   USE trc             !  passive tracers common variables 
+   USE sms_pisces      !  PISCES Source Minus Sink variables
+   USE p4zice          !  Co-limitations of differents nutrients
+   USE p4zsbc          !  External source of nutrients 
+   USE p4zint          !  interpolation and computation of various fields
+   USE sed             !  Sediment module
+   USE iom             !  I/O manager
+   USE prtctl_trc      !  print control for debugging
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   p4z_sed  
+   PUBLIC   p4z_sed_alloc
+ 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrpot    !: Nitrogen fixation 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments
+   REAL(wp) :: r1_rday                  !: inverse of rday
+   LOGICAL, SAVE :: lk_sed
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id: p4zsed.F90 10127 2018-09-13 15:27:54Z cetlod $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE p4z_sed( kt, knt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_sed  ***
+      !!
+      !! ** Purpose :   Compute loss of organic matter in the sediments. This
+      !!              is by no way a sediment model. The loss is simply 
+      !!              computed to balance the inout from rivers and dust
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      !
+      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
+      INTEGER  ::  ji, jj, jk, ikt
+      REAL(wp) ::  zrivalk, zrivsil, zrivno3
+      REAL(wp) ::  zwflux, zfminus, zfplus
+      REAL(wp) ::  zlim, zfact, zfactcal
+      REAL(wp) ::  zo2, zno3, zflx, zpdenit, z1pdenit, zolimit
+      REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep
+      REAL(wp) ::  zwstpoc, zwstpon, zwstpop
+      REAL(wp) ::  ztrfer, ztrpo4s, ztrdp, zwdust, zmudia, ztemp
+      REAL(wp) ::  xdiano3, xdianh4
+      REAL(wp) ::  zwssfep
+      !
+      CHARACTER (len=25) :: charout
+      REAL(wp), DIMENSION(jpi,jpj    ) :: zdenit2d, zbureff, zwork
+      REAL(wp), DIMENSION(jpi,jpj    ) :: zwsbio3, zwsbio4, zwscal
+      REAL(wp), DIMENSION(jpi,jpj    ) :: zsedcal, zsedsi, zsedc
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zsoufer, zlight
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrpo4, ztrdop, zirondep, zpdep
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zsidep, zwsfep, zironice
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('p4z_sed')
+      !
+      IF( kt == nittrc000 .AND. knt == 1 )   THEN 
+          r1_rday  = 1. / rday
+          IF (ln_sediment .AND. ln_sed_2way) THEN
+             lk_sed = .TRUE.
+          ELSE
+             lk_sed = .FALSE.
+          ENDIF
+      ENDIF
+      !
+      ! Allocate temporary workspace
+      IF( ln_p5z )    ALLOCATE( ztrpo4(jpi,jpj,jpk), ztrdop(jpi,jpj,jpk) )
+      IF( ln_ligand ) ALLOCATE( zwsfep(jpi,jpj) )
+
+
+      zdenit2d(:,:) = 0.e0
+      zbureff (:,:) = 0.e0
+      zwork   (:,:) = 0.e0
+      zsedsi  (:,:) = 0.e0
+      zsedcal (:,:) = 0.e0
+      zsedc   (:,:) = 0.e0
+
+
+      ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al.
+      ! ----------------------------------------------------
+      IF( ln_ironice ) THEN  
+         !                                              
+         ALLOCATE( zironice(jpi,jpj) )
+         !                                              
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zdep    = rfact2 / e3t_n(ji,jj,1)
+               zwflux  = fmmflx(ji,jj) / 1000._wp
+               zfminus = MIN( 0._wp, -zwflux ) * trb(ji,jj,1,jpfer) * zdep
+               zfplus  = MAX( 0._wp, -zwflux ) * icefeinput * zdep
+               zironice(ji,jj) =  zfplus + zfminus
+            END DO
+         END DO
+         !
+         tra(:,:,1,jpfer) = tra(:,:,1,jpfer) + zironice(:,:) 
+         ! 
+         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironice" ) )   &
+            &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! iron flux from ice
+         !
+         DEALLOCATE( zironice )
+         !                                              
+      ENDIF
+
+      ! Add the external input of nutrients from dust deposition
+      ! ----------------------------------------------------------
+      IF( ln_dust ) THEN
+         !                                              
+         ALLOCATE( zsidep(jpi,jpj), zpdep(jpi,jpj,jpk), zirondep(jpi,jpj,jpk) )
+         !                                              ! Iron and Si deposition at the surface
+         IF( ln_solub ) THEN
+            zirondep(:,:,1) = solub(:,:) * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
+         ELSE
+            zirondep(:,:,1) = dustsolub  * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 55.85 + 3.e-10 * r1_ryyss 
+         ENDIF
+         zsidep(:,:)   = 8.8 * 0.075 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 28.1 
+         zpdep (:,:,1) = 0.1 * 0.021 * dust(:,:) * mfrac * rfact2 / e3t_n(:,:,1) / 31. / po4r 
+         !                                              ! Iron solubilization of particles in the water column
+         !                                              ! dust in kg/m2/s ---> 1/55.85 to put in mol/Fe ;  wdust in m/j
+         zwdust = 0.03 * rday / ( wdust * 55.85 ) / ( 270. * rday )
+         DO jk = 2, jpkm1
+            zirondep(:,:,jk) = dust(:,:) * mfrac * zwdust * rfact2 * EXP( -gdept_n(:,:,jk) / 540. )
+            zpdep   (:,:,jk) = zirondep(:,:,jk) * 0.023
+         END DO
+         !                                              ! Iron solubilization of particles in the water column
+         tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:)
+         DO jk = 1, jpkm1
+            tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zpdep   (:,:,jk)
+            tra(:,:,jk,jpfer) = tra(:,:,jk,jpfer) + zirondep(:,:,jk) 
+         ENDDO
+         ! 
+         IF( lk_iomput ) THEN
+            IF( knt == nrdttrc ) THEN
+                IF( iom_use( "Irondep" ) )   &
+                &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * e3t_n(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
+                IF( iom_use( "pdust" ) )   &
+                &  CALL iom_put( "pdust"  , dust(:,:) / ( wdust * rday )  * tmask(:,:,1) ) ! dust concentration at surface
+            ENDIF
+         ENDIF
+         DEALLOCATE( zsidep, zpdep, zirondep )
+         !                                              
+      ENDIF
+     
+      ! Add the external input of nutrients from river
+      ! ----------------------------------------------------------
+      IF( ln_river ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 1, nk_rnf(ji,jj)
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2
+                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2
+                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2
+                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2
+               ENDDO
+            ENDDO
+         ENDDO
+         IF( ln_p5z ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  DO jk = 1, nk_rnf(ji,jj)
+                     tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + rivdop(ji,jj) * rfact2
+                     tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + rivdon(ji,jj) * rfact2
+                     tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + rivdoc(ji,jj) * rfact2
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+      ENDIF
+      
+      ! Add the external input of nutrients from nitrogen deposition
+      ! ----------------------------------------------------------
+      IF( ln_ndepo ) THEN
+         tra(:,:,1,jpno3) = tra(:,:,1,jpno3) + nitdep(:,:) * rfact2
+         tra(:,:,1,jptal) = tra(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2
+      ENDIF
+
+      ! Add the external input of iron from hydrothermal vents
+      ! ------------------------------------------------------
+      IF( ln_hydrofe ) THEN
+            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
+         IF( ln_ligand ) THEN
+            tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + ( hydrofe(:,:,:) * fep_rath ) * rfact2
+            tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * lgw_rath ) * rfact2
+         ENDIF
+         !
+         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   &
+            &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
+      ENDIF
+
+      ! OA: Warning, the following part is necessary to avoid CFL problems above the sediments
+      ! --------------------------------------------------------------------
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt  = mbkt(ji,jj)
+            zdep = e3t_n(ji,jj,ikt) / xstep
+            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
+            zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) )
+            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
+         END DO
+      END DO
+      !
+      IF( ln_ligand ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikt  = mbkt(ji,jj)
+               zdep = e3t_n(ji,jj,ikt) / xstep
+               zwsfep(ji,jj)  = MIN( 0.99 * zdep, wsfep(ji,jj,ikt)  )
+            END DO
+         ENDDO
+      ENDIF
+
+      IF( .NOT.lk_sed ) THEN
+
+         ! Add the external input of iron from sediment mobilization
+         ! ------------------------------------------------------
+         IF( ln_ironsed ) THEN
+                            tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
+            IF( ln_ligand ) tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + ( ironsed(:,:,:) * fep_rats ) * rfact2
+            !
+            IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
+               &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments
+         ENDIF
+
+         ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used
+         ! Computation of the fraction of organic matter that is permanently buried from Dunne's model
+         ! -------------------------------------------------------
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+              IF( tmask(ji,jj,1) == 1 ) THEN
+                 ikt = mbkt(ji,jj)
+                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
+                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4
+                 zflx  = LOG10( MAX( 1E-3, zflx ) )
+                 zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) )
+                 zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) )
+                 zdep  = LOG10( gdepw_n(ji,jj,ikt+1) )
+                 zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    &
+                   &                + 0.4721 * zo2 - 0.0996 * zdep + 0.4256 * zflx * zo2
+                 zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) )
+                   !
+                 zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
+                   &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6
+                 zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
+              ENDIF
+            END DO
+         END DO 
+         !
+      ENDIF
+
+      ! This loss is scaled at each bottom grid cell for equilibrating the total budget of silica in the ocean.
+      ! Thus, the amount of silica lost in the sediments equal the supply at the surface (dust+rivers)
+      ! ------------------------------------------------------
+      IF( .NOT.lk_sed )  zrivsil = 1._wp - sedsilfrac
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt  = mbkt(ji,jj)
+            zdep = xstep / e3t_n(ji,jj,ikt) 
+            zwsc = zwscal (ji,jj) * zdep
+            zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc
+            zcaloss = trb(ji,jj,ikt,jpcal) * zwsc
+            !
+            tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss
+            tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss
+         END DO
+      END DO
+      !
+      IF( .NOT.lk_sed ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikt  = mbkt(ji,jj)
+               zdep = xstep / e3t_n(ji,jj,ikt) 
+               zwsc = zwscal (ji,jj) * zdep
+               zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc
+               zcaloss = trb(ji,jj,ikt,jpcal) * zwsc
+               tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
+               !
+               zfactcal = MIN( excess(ji,jj,ikt), 0.2 )
+               zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
+               zrivalk  = sedcalfrac * zfactcal
+               tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
+               tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk
+               zsedcal(ji,jj) = (1.0 - zrivalk) * zcaloss * e3t_n(ji,jj,ikt) 
+               zsedsi (ji,jj) = (1.0 - zrivsil) * zsiloss * e3t_n(ji,jj,ikt) 
+            END DO
+         END DO
+      ENDIF
+      !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt  = mbkt(ji,jj)
+            zdep = xstep / e3t_n(ji,jj,ikt) 
+            zws4 = zwsbio4(ji,jj) * zdep
+            zws3 = zwsbio3(ji,jj) * zdep
+            tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4 
+            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3
+            tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4
+            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3
+         END DO
+      END DO
+      !
+      IF( ln_ligand ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikt     = mbkt(ji,jj)
+               zdep    = xstep / e3t_n(ji,jj,ikt) 
+               zwssfep = zwsfep(ji,jj) * zdep
+               tra(ji,jj,ikt,jpfep) = tra(ji,jj,ikt,jpfep) - trb(ji,jj,ikt,jpfep) * zwssfep
+            END DO
+         END DO
+      ENDIF
+      !
+      IF( ln_p5z ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikt  = mbkt(ji,jj)
+               zdep = xstep / e3t_n(ji,jj,ikt) 
+               zws4 = zwsbio4(ji,jj) * zdep
+               zws3 = zwsbio3(ji,jj) * zdep
+               tra(ji,jj,ikt,jpgon) = tra(ji,jj,ikt,jpgon) - trb(ji,jj,ikt,jpgon) * zws4
+               tra(ji,jj,ikt,jppon) = tra(ji,jj,ikt,jppon) - trb(ji,jj,ikt,jppon) * zws3
+               tra(ji,jj,ikt,jpgop) = tra(ji,jj,ikt,jpgop) - trb(ji,jj,ikt,jpgop) * zws4
+               tra(ji,jj,ikt,jppop) = tra(ji,jj,ikt,jppop) - trb(ji,jj,ikt,jppop) * zws3
+            END DO
+         END DO
+      ENDIF
+
+      IF( .NOT.lk_sed ) THEN
+         ! The 0.5 factor in zpdenit is to avoid negative NO3 concentration after
+         ! denitrification in the sediments. Not very clever, but simpliest option.
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ikt  = mbkt(ji,jj)
+               zdep = xstep / e3t_n(ji,jj,ikt) 
+               zws4 = zwsbio4(ji,jj) * zdep
+               zws3 = zwsbio3(ji,jj) * zdep
+               zrivno3 = 1. - zbureff(ji,jj)
+               zwstpoc = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
+               zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
+               z1pdenit = zwstpoc * zrivno3 - zpdenit
+               zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
+               tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit
+               tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit
+               tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit
+               tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * zpdenit
+               tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut
+               tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * zpdenit )
+               tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit 
+               sdenit(ji,jj) = rdenit * zpdenit * e3t_n(ji,jj,ikt)
+               zsedc(ji,jj)   = (1. - zrivno3) * zwstpoc * e3t_n(ji,jj,ikt)
+               IF( ln_p5z ) THEN
+                  zwstpop              = trb(ji,jj,ikt,jpgop) * zws4 + trb(ji,jj,ikt,jppop) * zws3
+                  zwstpon              = trb(ji,jj,ikt,jpgon) * zws4 + trb(ji,jj,ikt,jppon) * zws3
+                  tra(ji,jj,ikt,jpdon) = tra(ji,jj,ikt,jpdon) + ( z1pdenit - zolimit ) * zwstpon / (zwstpoc + rtrn)
+                  tra(ji,jj,ikt,jpdop) = tra(ji,jj,ikt,jpdop) + ( z1pdenit - zolimit ) * zwstpop / (zwstpoc + rtrn)
+               ENDIF
+            END DO
+         END DO
+       ENDIF
+
+
+      ! Nitrogen fixation process
+      ! Small source iron from particulate inorganic iron
+      !-----------------------------------
+      DO jk = 1, jpkm1
+         zlight (:,:,jk) =  ( 1.- EXP( -etot_ndcy(:,:,jk) / diazolight ) ) * ( 1. - fr_i(:,:) ) 
+         zsoufer(:,:,jk) = zlight(:,:,jk) * 2E-11 / ( 2E-11 + biron(:,:,jk) )
+      ENDDO
+      IF( ln_p4z ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !                      ! Potential nitrogen fixation dependant on temperature and iron
+                  zlim = ( 1.- xnanono3(ji,jj,jk) - xnanonh4(ji,jj,jk) )
+                  IF( zlim <= 0.2 )   zlim = 0.01
+                  zfact = zlim * rfact2
+
+                  ztrfer  = biron(ji,jj,jk)       / ( concfediaz + biron(ji,jj,jk)       )
+                  ztrpo4s = trb  (ji,jj,jk,jppo4) / ( concnnh4   + trb  (ji,jj,jk,jppo4) ) 
+                  nitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday ) &
+                    &                *  zfact * MIN( ztrfer, ztrpo4s ) * zlight(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSE       ! p5z
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  !                      ! Potential nitrogen fixation dependant on temperature and iron
+                  ztemp = tsn(ji,jj,jk,jp_tem)
+                  zmudia = MAX( 0.,-0.001096*ztemp**2 + 0.057*ztemp -0.637 ) * 7.625
+                  !       Potential nitrogen fixation dependant on temperature and iron
+                  xdianh4 = trb(ji,jj,jk,jpnh4) / ( concnnh4 + trb(ji,jj,jk,jpnh4) )
+                  xdiano3 = trb(ji,jj,jk,jpno3) / ( concnno3 + trb(ji,jj,jk,jpno3) ) * (1. - xdianh4)
+                  zlim = ( 1.- xdiano3 - xdianh4 )
+                  IF( zlim <= 0.1 )   zlim = 0.01
+                  zfact = zlim * rfact2
+                  ztrfer = biron(ji,jj,jk) / ( concfediaz + biron(ji,jj,jk) )
+                  ztrpo4(ji,jj,jk) = trb(ji,jj,jk,jppo4) / ( 1E-6 + trb(ji,jj,jk,jppo4) )
+                  ztrdop(ji,jj,jk) = trb(ji,jj,jk,jpdop) / ( 1E-6 + trb(ji,jj,jk,jpdop) ) * (1. - ztrpo4(ji,jj,jk))
+                  ztrdp = ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk)
+                  nitrpot(ji,jj,jk) =  zmudia * r1_rday * zfact * MIN( ztrfer, ztrdp ) * zlight(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      ! Nitrogen change due to nitrogen fixation
+      ! ----------------------------------------
+      IF( ln_p4z ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zfact = nitrpot(ji,jj,jk) * nitrfix
+                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) +             zfact
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3      * zfact
+                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2nit     * zfact 
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) &
+                  &                     * 0.002 * trb(ji,jj,jk,jpdoc) * xstep
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * xstep
+              END DO
+            END DO 
+         END DO
+      ELSE    ! p5z
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zfact = nitrpot(ji,jj,jk) * nitrfix
+                  tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zfact / 3.0
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * zfact / 3.0
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - 16.0 / 46.0 * zfact * ( 1.0 - 1.0 / 3.0 ) &
+                  &                     * ztrpo4(ji,jj,jk) / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
+                  tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zfact * 1.0 / 3.0
+                  tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zfact * 1.0 / 3.0
+                  tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + 16.0 / 46.0 * zfact / 3.0  &
+                  &                     - 16.0 / 46.0 * zfact * ztrdop(ji,jj,jk)   &
+                  &                     / (ztrpo4(ji,jj,jk) + ztrdop(ji,jj,jk) + rtrn)
+                  tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zfact * 1.0 / 3.0 * 2.0 / 3.0
+                  tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + zfact * 1.0 / 3.0 * 2.0 /3.0
+                  tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 2.0 /3.0
+                  tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zfact * 1.0 / 3.0 * 1.0 / 3.0
+                  tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zfact * 1.0 / 3.0 * 1.0 /3.0
+                  tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + 16.0 / 46.0 * zfact * 1.0 / 3.0 * 1.0 /3.0
+                  tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + ( o2ut + o2nit ) * zfact * 2.0 / 3.0 + o2nit * zfact / 3.0
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - 30E-6 * zfact * 1.0 / 3.0 
+                  tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 30E-6 * zfact * 1.0 / 3.0 * 2.0 / 3.0
+                  tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + 30E-6 * zfact * 1.0 / 3.0 * 1.0 / 3.0
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+              END DO
+            END DO 
+         END DO
+         !
+      ENDIF
+
+      IF( lk_iomput ) THEN
+         IF( knt == nrdttrc ) THEN
+            zfact = 1.e+3 * rfact2r !  conversion from molC/l/kt  to molN/m3/s
+            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation 
+            IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated )
+               zwork(:,:) = 0.
+               DO jk = 1, jpkm1
+                 zwork(:,:) = zwork(:,:) + nitrpot(:,:,jk) * nitrfix * rno3 * zfact * e3t_n(:,:,jk) * tmask(:,:,jk)
+               ENDDO
+               CALL iom_put( "INTNFIX" , zwork ) 
+            ENDIF
+            IF( iom_use("SedCal" ) ) CALL iom_put( "SedCal", zsedcal(:,:) * zfact )
+            IF( iom_use("SedSi" ) )  CALL iom_put( "SedSi",  zsedsi (:,:) * zfact )
+            IF( iom_use("SedC" ) )   CALL iom_put( "SedC",   zsedc  (:,:) * zfact )
+            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", sdenit (:,:) * zfact * rno3 )
+         ENDIF
+      ENDIF
+      !
+      IF(ln_ctl) THEN  ! print mean trends (USEd for debugging)
+         WRITE(charout, fmt="('sed ')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+      ENDIF
+      !
+      IF( ln_p5z )    DEALLOCATE( ztrpo4, ztrdop )
+      IF( ln_ligand ) DEALLOCATE( zwsfep )
+      !
+      IF( ln_timing )  CALL timing_stop('p4z_sed')
+      !
+   END SUBROUTINE p4z_sed
+
+
+   INTEGER FUNCTION p4z_sed_alloc()
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_sed_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( nitrpot(jpi,jpj,jpk), sdenit(jpi,jpj), STAT=p4z_sed_alloc )
+      !
+      IF( p4z_sed_alloc /= 0 )   CALL ctl_warn('p4z_sed_alloc: failed to allocate arrays')
+      !
+   END FUNCTION p4z_sed_alloc
+
+   !!======================================================================
+END MODULE p4zsed
Index: /NEMO/trunk/src/TOP/PISCES/SED/sed_oce.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/sed_oce.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/sed_oce.F90	(revision 10225)
@@ -0,0 +1,36 @@
+MODULE sed_oce
+   !!======================================================================
+   !!                        ***  sed  ***
+   !! Sediment :   set sediment global variables
+   !!======================================================================
+
+   !! History :
+   !!        !  06-12  (C. Ethe)  Orignal
+   !!----------------------------------------------------------------------
+   USE par_sed
+   USE trc, ONLY : profsed
+
+   IMPLICIT NONE
+   PUBLIC
+
+   PUBLIC sed_oce_alloc
+
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  profsedw       !: depth of middle of each layer
+
+   !! $Id: sed.F90 7646 2017-02-06 09:25:03Z timgraham $
+CONTAINS
+
+   INTEGER FUNCTION sed_oce_alloc()
+      !!-------------------------------------------------------------------
+      !!                    *** ROUTINE sed_alloc ***
+      !!-------------------------------------------------------------------
+      USE lib_mpp, ONLY: ctl_warn
+      !!-------------------------------------------------------------------
+      !
+      ALLOCATE( profsed(jpksed) , profsedw(jpksed) , STAT=sed_oce_alloc )
+
+      IF( sed_oce_alloc /= 0 )   CALL ctl_warn('sed_oce_alloc: failed to allocate arrays')
+      !
+   END FUNCTION sed_oce_alloc
+
+END MODULE sed_oce
Index: /NEMO/trunk/src/TOP/PISCES/SED/seddiff.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/seddiff.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/seddiff.F90	(revision 10225)
@@ -0,0 +1,120 @@
+MODULE seddiff
+   !!======================================================================
+   !!              ***  MODULE  seddsr  ***
+   !!    Sediment : dissolution and reaction in pore water related 
+   !!    related to organic matter
+   !!=====================================================================
+   !! * Modules used
+   USE sed     ! sediment global variable
+   USE sed_oce
+   USE sedmat  ! linear system of equations
+   USE sedini
+   USE lib_mpp         ! distribued memory computing library
+   USE lib_fortran
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC sed_diff
+
+   !! * Module variables
+
+   !! $Id: seddsr.F90 5215 2015-04-15 16:11:56Z nicolasmartin $
+CONTAINS
+   
+   SUBROUTINE sed_diff( kt, knt ) 
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE sed_diff  ***
+      !! 
+      !!  ** Purpose :  computes pore water diffusion
+      !!
+      !!  ** Methode :  implicit computation of undersaturation
+      !!               resulting from diffusive pore water transport.
+      !!
+      !!  ** Remarks :
+      !!              - undersaturation : deviation from saturation concentration
+      !!   History :
+      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
+      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
+      !!        !  06-04 (C. Ethe)  Re-organization
+      !!        !  19-08 (O. Aumont) Debugging and improvement of the model
+      !!----------------------------------------------------------------------
+      !! Arguments
+      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration
+      ! --- local variables
+      INTEGER :: ji, jk, js   ! dummy looop indices
+
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2   ! reaction rate in pore water
+      !!
+      !!----------------------------------------------------------------------
+
+      IF( ln_timing )  CALL timing_start('sed_diff')
+!
+      IF( kt == nitsed000 .AND. knt == 1 ) THEN
+         IF (lwp) THEN
+            WRITE(numsed,*) ' sed_diff : pore-water diffusion '
+            WRITE(numsed,*) ' '
+         ENDIF
+      ENDIF
+
+     ! Initializations
+     !----------------------
+      zrearat1(:,:)   = 0.
+      zrearat2(:,:) = 0.
+
+      !---------------------------
+      ! Solves PO4 diffusion 
+      !----------------------------
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwpo4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwpo4), dtsed2 / 2.0 )
+
+      !---------------------------
+      ! Solves NH4 diffusion 
+      !----------------------------
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwnh4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwnh4), dtsed2 / 2.0 )
+
+      !---------------------------
+      ! Solves Fe2+ diffusion 
+      !----------------------------
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwfe2, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwfe2), dtsed2 / 2.0 )
+
+      !---------------------------
+      ! Solves H2S diffusion 
+      !----------------------------
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwh2s, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwh2s), dtsed2 / 2.0  )
+
+      !---------------------------
+      ! Solves SO4 diffusion 
+      !----------------------------
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwso4, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwso4), dtsed2 / 2.0 )
+
+      !---------------------------
+      ! Solves O2 diffusion
+      !----------------------------
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwoxy, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwoxy), dtsed2 / 2.0 )
+
+      !---------------------------
+      ! Solves NO3 diffusion
+      !----------------------------
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwno3, jpoce, jpksed, zrearat1, zrearat2, pwcp(:,:,jwno3), dtsed2 / 2.0 )
+
+      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, sedligand(:,:), dtsed2 / 2.0 )
+
+      IF( ln_timing )  CALL timing_stop('sed_diff')
+!      
+   END SUBROUTINE sed_diff
+
+END MODULE seddiff
Index: /NEMO/trunk/src/TOP/PISCES/SED/sedinitrc.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/sedinitrc.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/sedinitrc.F90	(revision 10225)
@@ -0,0 +1,192 @@
+MODULE sedinitrc
+   !!======================================================================
+   !!              ***  MODULE  sedinitrc  ***
+   !! Sediment : define sediment variables
+   !!=====================================================================
+
+   !!----------------------------------------------------------------------
+   !!   sed_init    : initialization, namelist read, and parameters control
+   !!----------------------------------------------------------------------
+   !! * Modules used
+   USE sed     ! sediment global variable
+   USE sed_oce
+   USE sedini
+   USE seddta
+   USE sedrst
+   USE sedco3
+   USE sedchem
+   USE sedarr
+   USE lib_mpp         ! distribued memory computing library
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   REAL(wp)    ::  &
+      ryear = 365. * 24. * 3600. !:  1 year converted in second
+
+   !! *  Routine accessibility
+   PUBLIC sed_initrc          ! routine called by opa.F90
+
+   !! $Id: sedini.F90 5215 2015-04-15 16:11:56Z nicolasmartin $
+CONTAINS
+
+
+   SUBROUTINE sed_initrc
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE sed_init  ***
+      !!
+      !! ** Purpose :  Initialization of sediment module
+      !!               - Reading namelist
+      !!               - Read the deepest water layer thickness
+      !!                 ( using as mask ) in Netcdf file
+      !!               - Convert unity if necessary
+      !!               - sets initial sediment composition
+      !!                 ( only clay or reading restart file )
+      !!               - sets sediment grid, porosity and others constants
+      !!
+      !!   History :
+      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
+      !!        !  06-07  (C. Ethe)  Re-organization
+      !!----------------------------------------------------------------------
+      INTEGER :: ji, jj, ikt
+      !!----------------------------------------------------------------------
+
+
+      ! Initialize the sediment tracers concentrations
+      !------------------------------------------------
+
+      IF(lwp) WRITE(numsed,*) ' sed_initrc : Initialization of sediment concentration '
+      IF(lwp) WRITE(numsed,*) ' '
+
+      ! Determination of sediments number of points and allocate global variables
+
+      ! sets initial sediment composition
+      ! ( only clay or reading restart file )
+      !---------------------------------------
+      CALL sed_init_data
+
+
+      CALL sed_init_wri
+
+
+   END SUBROUTINE sed_initrc
+
+
+   SUBROUTINE sed_init_data
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE sed_init_data  ***
+      !!
+      !! ** Purpose :  Initialization of sediment module
+      !!               - sets initial sediment composition
+      !!                 ( only clay or reading restart file )
+      !!
+      !!   History :
+      !!        !  06-07  (C. Ethe)  original
+      !!----------------------------------------------------------------------
+ 
+      ! local variables
+      INTEGER :: &
+         ji, jk, zhipor
+
+      !--------------------------------------------------------------------
+ 
+
+      IF( .NOT. ln_rst_sed ) THEN
+
+         IF (lwp) WRITE(numsed,*) ' Initilization of default values of sediment components'
+
+         ! default values for initial pore water concentrations [mol/l]
+         pwcp(:,:,:) = 0.
+         ! default value for initial solid component (fraction of dry weight dim=[0])
+         ! clay
+         solcp(:,:,:) = 0.
+         solcp(:,2:jpksed,jsclay) = 1.0 * 0.965
+         solcp(:,2:jpksed,jsfeo)  = 1.0 * 0.035
+
+         ! Initialization of [h+] and [co3--]
+
+         zhipor = 8.0
+         ! Initialization of [h+] in mol/kg
+         DO jk = 1, jpksed
+            DO ji = 1, jpoce
+               hipor (ji,jk) = 10.**( -1. * zhipor )
+            ENDDO
+         ENDDO
+
+         co3por(:,:) = 1E-6
+
+      ELSE   
+  
+         IF (lwp) WRITE(numsed,*) ' Initilization of Sediment components from restart'
+
+         CALL sed_rst_cal( nitsed000, 'READ' )
+         CALL sed_rst_read
+
+      ENDIF
+
+
+      ! Load initial Pisces Data for bot. wat. Chem and fluxes
+      CALL sed_dta ( nitsed000 ) 
+
+      ! Initialization of chemical constants
+      CALL sed_chem ( nitsed000 )
+
+      ! Stores initial sediment data for mass balance calculation
+      pwcp0 (1:jpoce,1:jpksed,1:jpwat ) = pwcp (1:jpoce,1:jpksed,1:jpwat ) 
+      solcp0(1:jpoce,1:jpksed,1:jpsol ) = solcp(1:jpoce,1:jpksed,1:jpsol) 
+
+      ! Conversion of [h+] in mol/Kg to get it in mol/l ( multiplication by density)
+      DO jk = 1, jpksed
+         hipor(1:jpoce,jk) = hipor(1:jpoce,jk) * densSW(1:jpoce)
+      ENDDO
+
+
+      ! In default case - no restart - sedco3 is run to initiate [h+] and [co32-]
+      ! Otherwise initiate values of pH and co3 read in restart
+      IF( .NOT. ln_rst_sed ) THEN
+         ! sedco3 is run to initiate[h+] [co32-] in mol/l of solution
+         CALL sed_co3 ( nitsed000 )
+
+      ENDIF
+            
+   END SUBROUTINE sed_init_data
+
+   SUBROUTINE sed_init_wri
+
+      INTEGER :: jk
+
+      IF (lwp) THEN
+         WRITE(numsed,*)' '
+         WRITE(numsed,*)'======== Write summary of sediment char.  ============'
+         WRITE(numsed,*)' '
+         WRITE(numsed,*)' '
+         WRITE(numsed,*)'-------------------------------------------------------------------'
+         WRITE(numsed,*)' Initial Conditions '
+         WRITE(numsed,*)'-------------------------------------------------------------------'
+         WRITE(numsed,*)'dzm = dzkbot minimum to calculate ', 0.
+         WRITE(numsed,*)'Local zone : jpi, jpj, jpksed : ',jpi, jpj, jpksed
+         WRITE(numsed,*)'jpoce = ',jpoce,' nbtot pts = ',jpij,' nb earth pts = ',jpij - jpoce
+         WRITE(numsed,*)'sublayer thickness dz(1) [cm] : ', dz(1)
+         WRITE(numsed,*)'Vertical domain of the sediment'
+         WRITE(numsed,*)'-------------------------------'
+         WRITE(numsed,*)' Indice, profsed, dz'
+         DO jk = 2, jpksed
+            WRITE(numsed,*) jk,profsed(jk),dz(jk) 
+         END DO
+         WRITE(numsed,*)' nb solid comp : ',jpsol
+         WRITE(numsed,*)'(1=opal,2=clay,3=POC,4=CaCO3), 5=POS, 6=POR, 7=FEO, 8=FeS'
+         WRITE(numsed,*)'weight mol 1,2,3,4,5,6,7'
+         WRITE(numsed,'(8(F0.2,3X))')mol_wgt(jsopal),mol_wgt(jsclay),mol_wgt(jspoc),mol_wgt(jscal),mol_wgt(jspos),mol_wgt(jspor),mol_wgt(jsfeo),mol_wgt(jsfes)
+         WRITE(numsed,*)'nb dissolved comp',jpwat
+         WRITE(numsed,*)'1=silicic acid,,2=O2,3=DIC,4=NO3,5=PO4,6=Alk,7=NH4,8=ODU'
+         WRITE(numsed,*)'redfield coef C,O,N P Dit '
+         WRITE(numsed,'(5(F0.2,3X))')1./spo4r,so2ut/spo4r,srno3/spo4r,spo4r/spo4r,srDnit/spo4r
+         WRITE(numsed,*) ' '
+         WRITE(numsed,*) ' End Of Initialization '
+         WRITE(numsed,*) ' '
+      ENDIF
+!
+   END SUBROUTINE sed_init_wri
+
+END MODULE sedinitrc
Index: /NEMO/trunk/src/TOP/PISCES/SED/sedinorg.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/sedinorg.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/sedinorg.F90	(revision 10225)
@@ -0,0 +1,300 @@
+MODULE sedinorg
+   !!======================================================================
+   !!              ***  MODULE  sedinorg  ***
+   !!    Sediment : dissolution and reaction in pore water of 
+   !!               inorganic species
+   !!=====================================================================
+   !! * Modules used
+   USE sed     ! sediment global variable
+   USE sed_oce
+   USE sedmat  ! linear system of equations
+   USE sedco3  ! carbonate ion and proton concentration 
+   USE sedini
+   USE seddsr
+   USE lib_mpp         ! distribued memory computing library
+   USE lib_fortran
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC sed_inorg
+
+   !! $Id: seddsr.F90 5215 2015-04-15 16:11:56Z nicolasmartin $
+CONTAINS
+   
+   SUBROUTINE sed_inorg( kt ) 
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE sed_inorg  ***
+      !! 
+      !!  ** Purpose :  computes pore water dissolution and reaction
+      !!
+      !!  ** Methode :  implicit simultaneous computation of undersaturation
+      !!               resulting from diffusive pore water transport and chemical
+      !!               pore water reactions. Solid material is consumed according
+      !!               to redissolution and remineralisation
+      !!
+      !!  ** Remarks :
+      !!              - undersaturation : deviation from saturation concentration
+      !!              - reaction rate   : sink of undersaturation from dissolution
+      !!                                 of solid material 
+      !!
+      !!   History :
+      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
+      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
+      !!        !  06-04 (C. Ethe)  Re-organization
+      !!        !  19-08 (O. Aumont) Debugging and improvement of the model
+      !!----------------------------------------------------------------------
+      !! Arguments
+      INTEGER, INTENT(in) ::   kt       ! number of iteration
+      ! --- local variables
+      INTEGER :: ji, jk, js, jw         ! dummy looop indices
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2    ! reaction rate in pore water
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite   
+      REAL(wp), DIMENSION(jpoce) :: zco3eq
+      REAL(wp), DIMENSION(jpoce,jpksed,jpsol) :: zvolc    ! temp. variables
+      REAL(wp), DIMENSION(jpoce) :: zsieq
+      REAL(wp)  ::  zsolid1, zvolw, zreasat
+      REAL(wp)  ::  zsatur, zsatur2, znusil, zsolcpcl, zsolcpsi
+      !!
+      !!----------------------------------------------------------------------
+
+      IF( ln_timing )  CALL timing_start('sed_inorg')
+!
+      IF( kt == nitsed000 ) THEN
+         IF (lwp) THEN
+            WRITE(numsed,*) ' sed_inorg : Dissolution reaction '
+            WRITE(numsed,*) ' '
+         ENDIF
+!         ! 
+      ENDIF
+
+     ! Initializations
+     !----------------------
+      
+      zrearat1(:,:) = 0.    ;   zundsat(:,:)  = 0. 
+      zrearat2(:,:) = 0.    ;   zrearat2(:,:) = 0.
+      zco3eq(:)     = rtrn
+      zvolc(:,:,:)  = 0.
+
+      ! -----------------------------------------------
+      ! Computation of Si solubility
+      ! Param of Ridgwell et al. 2002
+      ! -----------------------------------------------
+
+      DO ji = 1, jpoce
+         zsolcpcl = 0.0
+         zsolcpsi = 0.0
+         DO jk = 1, jpksed
+            zsolcpsi = zsolcpsi + solcp(ji,jk,jsopal) * dz(jk)
+            zsolcpcl = zsolcpcl + solcp(ji,jk,jsclay) * dz(jk)
+         END DO
+         zsieq(ji) = sieqs(ji) * MAX(0.25, 1.0 - (0.045 * zsolcpcl / zsolcpsi )**0.58 )
+         zsieq(ji) = MAX( rtrn, sieqs(ji) )
+      END DO
+
+      DO js = 1, jpsol
+         DO jk = 1, jpksed
+            DO ji = 1, jpoce    
+               zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  &
+                  &              ( volw3d(ji,jk) * 1.e-3  )     
+            ENDDO
+         ENDDO
+      ENDDO
+
+      !----------------------------------------------------------
+      ! 5.  Beginning of  Pore Water diffusion and solid reaction
+      !---------------------------------------------------------
+      
+      !-----------------------------------------------------------------------------
+      ! For jk=2,jpksed, and for couple 
+      !  1 : jwsil/jsopal  ( SI/Opal )
+      !  2 : jsclay/jsclay ( clay/clay ) 
+      !  3 : jwoxy/jspoc   ( O2/POC )
+      !  reaction rate is a function of solid=concentration in solid reactif in [mol/l] 
+      !  and undersaturation in [mol/l].
+      !  Solid weight fractions should be in ie [mol/l])
+      !  second member and solution are in zundsat variable
+      !-------------------------------------------------------------------------
+
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            ! For Silicic Acid and clay
+            zundsat(ji,jk) = zsieq(ji) - pwcp(ji,jk,jwsil)
+         ENDDO
+      ENDDO
+      
+      ! Definition of reaction rates [rearat]=sans dim 
+      ! For jk=1 no reaction (pure water without solid) for each solid compo
+      DO ji = 1, jpoce
+         zrearat1(ji,:) = 0.
+         zrearat2(ji,:) = 0.
+      ENDDO
+
+      ! left hand side of coefficient matrix
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
+            zsatur = MAX(0., zundsat(ji,jk) / zsieq(ji) )
+            zsatur2 = (1.0 + temp(ji) / 400.0 )**37
+            znusil = ( 0.225 * ( 1.0 + temp(ji) / 15.) + 0.775 * zsatur2 * zsatur**2.25 ) / zsieq(ji)
+            zrearat1(ji,jk)  = ( reac_sil * znusil * dtsed * zsolid1 ) / &
+               &                ( 1. + reac_sil * znusil * dtsed * zundsat(ji,jk) )
+         ENDDO
+      ENDDO
+
+      CALL sed_mat( jwsil, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
+
+      ! New solid concentration values (jk=2 to jksed) for each couple 
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / ( zvolc(ji,jk,jsopal) )
+            solcp(ji,jk,jsopal) = solcp(ji,jk,jsopal) - zreasat
+         ENDDO
+      ENDDO
+
+      ! New pore water concentrations    
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            ! Acid Silicic 
+            pwcp(ji,jk,jwsil)  = zsieq(ji) - zundsat(ji,jk)
+         ENDDO
+      ENDDO
+
+      !---------------------------------------------------------------
+      ! Performs CaCO3 particle deposition and redissolution (indice 9)
+      !--------------------------------------------------------------
+
+      ! computes co3por from the updated pwcp concentrations (note [co3por] = mol/l)
+
+      CALL sed_co3( kt )
+
+      ! *densSW(l)**2 converts aksps [mol2/kg sol2] into [mol2/l2] to get [undsat] in [mol/l]
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            zco3eq(ji)     = aksps(ji) * densSW(ji) * densSW(ji) / ( calcon2(ji) + rtrn )
+            zco3eq(ji)     = MAX( rtrn, zco3eq(ji) ) 
+            zundsat(ji,jk) = MAX(0., zco3eq(ji) - co3por(ji,jk) )
+         ENDDO
+      ENDDO
+
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
+            zrearat1(ji,jk) = ( reac_cal * dtsed * zsolid1 / zco3eq(ji) ) / &
+                  &               ( 1. + reac_cal * dtsed * zundsat(ji,jk) / zco3eq(ji) )
+         END DO
+      END DO
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
+
+      ! New solid concentration values (jk=2 to jksed) for cacO3
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zreasat = zrearat1(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jscal)
+            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat
+         ENDDO
+      ENDDO
+
+      ! New dissolved concentrations
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            zreasat = zrearat1(ji,jk) * zundsat(ji,jk)    
+            ! For DIC
+            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
+            ! For alkalinity
+            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.0 * zreasat 
+         ENDDO
+      ENDDO
+
+      !-------------------------------------------------
+      ! Beginning DIC, Alkalinity 
+      !-------------------------------------------------
+      
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce      
+            zundsat(ji,jk)   = pwcp(ji,jk,jwdic)
+            zrearat1(ji,jk)  = 0.
+         ENDDO
+      ENDDO
+
+      ! solves tridiagonal system
+      CALL sed_mat( jwdic, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
+
+     ! New dissolved concentrations      
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce                      
+            pwcp(ji,jk,jwdic) = zundsat(ji,jk)
+         ENDDO
+      ENDDO            
+
+      !-------------------------------------------------
+      ! Beginning DIC, Alkalinity 
+      !-------------------------------------------------
+
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            zundsat(ji,jk) = pwcp(ji,jk,jwalk)
+            zrearat1(ji,jk) = 0.
+         ENDDO
+      ENDDO
+!
+!      ! solves tridiagonal system
+      CALL sed_mat( jwalk, jpoce, jpksed, zrearat1, zrearat2, zundsat, dtsed )
+!
+!      ! New dissolved concentrations      
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            pwcp(ji,jk,jwalk) = zundsat(ji,jk)
+         ENDDO
+      ENDDO
+      
+      !----------------------------------
+      !   Back to initial geometry
+      !-----------------------------
+      
+      !---------------------------------------------------------------------
+      !   1/ Compensation for ajustement of the bottom water concentrations
+      !      (see note n° 1 about *por(2))
+      !--------------------------------------------------------------------
+      DO jw = 1, jpwat
+         DO ji = 1, jpoce
+            pwcp(ji,1,jw) = pwcp(ji,1,jw) + &
+               &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / dzkbot(ji)
+         END DO
+      ENDDO
+      
+      !-----------------------------------------------------------------------
+      !    2/ Det of new rainrg taking account of the new weight fraction obtained 
+      !      in dz3d(2) after diffusion/reaction (react/diffu are also in dzdep!)
+      !      This new rain (rgntg rm) will be used in advection/burial routine
+      !------------------------------------------------------------------------
+      DO js = 1, jpsol
+         DO ji = 1, jpoce
+            rainrg(ji,js) = raintg(ji) * solcp(ji,2,js)
+            rainrm(ji,js) = rainrg(ji,js) / mol_wgt(js)
+         END DO
+      ENDDO
+
+      !  New raintg
+      raintg(:) = 0.
+      DO js = 1, jpsol
+         DO ji = 1, jpoce
+            raintg(ji) = raintg(ji) + rainrg(ji,js)
+         END DO
+      ENDDO
+      
+      !--------------------------------
+      !    3/ back to initial geometry
+      !--------------------------------
+      DO ji = 1, jpoce
+         dz3d  (ji,2) = dz(2)
+         volw3d(ji,2) = dz3d(ji,2) * por(2)
+         vols3d(ji,2) = dz3d(ji,2) * por1(2)
+      ENDDO
+
+      IF( ln_timing )  CALL timing_stop('sed_inorg')
+!      
+   END SUBROUTINE sed_inorg
+
+END MODULE sedinorg
Index: /NEMO/trunk/src/TOP/PISCES/SED/sedorg.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/sedorg.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/sedorg.F90	(revision 10225)
@@ -0,0 +1,159 @@
+MODULE sedorg
+   !!======================================================================
+   !!              ***  MODULE  seddsr  ***
+   !!    Sediment : dissolution and reaction in pore water related 
+   !!    related to organic matter
+   !!=====================================================================
+   !! * Modules used
+   USE sed     ! sediment global variable
+   USE sed_oce
+   USE sedini
+   USE seddiff
+   USE seddsr
+   USE lib_mpp         ! distribued memory computing library
+   USE lib_fortran
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC sed_org
+
+   !! * Module variables
+
+   REAL(wp) :: zadsnh4
+
+   !! $Id: seddsr.F90 5215 2015-04-15 16:11:56Z nicolasmartin $
+CONTAINS
+   
+   SUBROUTINE sed_org( kt ) 
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE sed_org  ***
+      !! 
+      !!  ** Purpose :  computes pore water diffusion and reaction
+      !!
+      !!  ** Methode :  Computation of the redox reactions in sediment.
+      !!                The main redox reactions are solved in sed_dsr whereas
+      !!                the secondary reactions are solved in sed_dsr_redoxb.
+      !!                A strand spliting approach is being used here (see 
+      !!                sed_dsr_redoxb for more information). 
+      !!                Diffusive fluxes are computed in sed_diff
+      !!
+      !!   History :
+      !!        !  98-08 (E. Maier-Reimer, Christoph Heinze )  Original code
+      !!        !  04-10 (N. Emprin, M. Gehlen ) f90
+      !!        !  06-04 (C. Ethe)  Re-organization
+      !!        !  19-08 (O. Aumont) Debugging and improvement of the model.
+      !!                             The original method is replaced by a 
+      !!                             Strand splitting method which deals 
+      !!                             well with stiff reactions.
+      !!----------------------------------------------------------------------
+      !! Arguments
+      INTEGER, INTENT(in) ::   kt
+      ! --- local variables
+      INTEGER  :: ji, jk, js, jw, jnt   ! dummy looop indices
+      REAL(wp) :: zadsnh4
+      !!
+      !!----------------------------------------------------------------------
+
+      IF( ln_timing )  CALL timing_start('sed_org')
+!
+      IF( kt == nitsed000 ) THEN
+         IF (lwp) THEN
+            WRITE(numsed,*) ' sed_org : Organic degradation related reactions and diffusion'
+            WRITE(numsed,*) ' '
+         ENDIF
+!         ! 
+         dens_mol_wgt(1:jpsol) = denssol / mol_wgt(1:jpsol)
+         ! 
+      ENDIF
+
+      dtsed2 = dtsed / REAL( nrseddt, wp )
+
+      ! 1. Change of geometry
+      !    Increase of dz3d(2) thickness : dz3d(2) = dz3d(2)+dzdep
+      !    Warning : no change for dz(2)
+      !---------------------------------------------------------
+      dz3d(1:jpoce,2) = dz3d(1:jpoce,2) + dzdep(1:jpoce)
+
+      ! New values for volw3d(:,2) and vols3d(:,2)
+      ! Warning : no change neither for volw(2) nor  vols(2)
+      !------------------------------------------------------
+      volw3d(1:jpoce,2) = dz3d(1:jpoce,2) * por(2)
+      vols3d(1:jpoce,2) = dz3d(1:jpoce,2) * por1(2)
+
+      ! 2. Change of previous solid fractions (due to volum changes) for k=2
+      !---------------------------------------------------------------------
+
+      DO js = 1, jpsol
+         DO ji = 1, jpoce
+            solcp(ji,2,js) = solcp(ji,2,js) * dz(2) / dz3d(ji,2)
+         ENDDO
+      END DO
+
+      ! 3. New solid fractions (including solid rain fractions) for k=2
+      !------------------------------------------------------------------
+      DO js = 1, jpsol
+         DO ji = 1, jpoce
+            IF (raintg(ji) .ne. 0) THEN
+               solcp(ji,2,js) = solcp(ji,2,js) + &
+               &           ( rainrg(ji,js) / raintg(ji) ) * ( dzdep(ji) / dz3d(ji,2) )
+               ! rainrm are temporary cancel
+               rainrm(ji,js) = 0.
+            ENDIF
+         END DO
+      ENDDO
+
+      ! 4.  Adjustment of bottom water concen.(pwcp(1)):
+      !     We impose that pwcp(2) is constant. Including dzdep in dz3d(:,2) we assume
+      !     that dzdep has got a porosity of por(2). So pore water volum of jk=2 increase.
+      !     To keep pwcp(2) cste we must compensate this "increase" by a slight adjusment
+      !     of bottom water concentration.
+      !     This adjustment is compensate at the end of routine
+      !-------------------------------------------------------------
+      DO jw = 1, jpwat
+         DO ji = 1, jpoce
+            pwcp(ji,1,jw) = pwcp(ji,1,jw) - &
+               &            pwcp(ji,2,jw) * dzdep(ji) * por(2) / ( dzkbot(ji) + rtrn )
+         END DO
+      ENDDO
+
+      zadsnh4 = 1.0 / ( 1.0 + adsnh4 )
+
+      ! --------------------------------------------------
+      ! Computation of the diffusivities
+      ! --------------------------------------------------
+
+      DO js = 1, jpwat
+         DO jk = 1, jpksed
+            DO ji = 1, jpoce
+               diff(ji,jk,js) = ( diff1(js) + diff2(js) * temp(ji) ) / ( 1.0 - 2.0 * log( por(jk) ) )
+            END DO
+         END DO
+      END DO
+
+      ! Impact of bioirrigation and adsorption on diffusion
+      ! ---------------------------------------------------
+
+      diff(:,:,jwnh4) = diff(:,:,jwnh4) * ( 1.0 + irrig(:,:) ) * zadsnh4
+      diff(:,:,jwsil) = diff(:,:,jwsil) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwoxy) = diff(:,:,jwoxy) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwdic) = diff(:,:,jwdic) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwno3) = diff(:,:,jwno3) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwpo4) = diff(:,:,jwpo4) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwalk) = diff(:,:,jwalk) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwh2s) = diff(:,:,jwh2s) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwso4) = diff(:,:,jwso4) * ( 1.0 + irrig(:,:) )
+      diff(:,:,jwfe2) = diff(:,:,jwfe2) * ( 1.0 + 0.2 * irrig(:,:) )
+
+      DO jnt = 1, nrseddt
+         CALL sed_diff( kt, jnt )        ! 1st pass in diffusion to get values at t+1/2
+         CALL sed_dsr ( kt, jnt )        ! Dissolution reaction
+         CALL sed_diff( kt, jnt )        ! 2nd pass in diffusion to get values at t+1
+      END DO
+
+
+      IF( ln_timing )  CALL timing_stop('sed_org')
+!      
+   END SUBROUTINE sed_org
+
+END MODULE sedorg
Index: /NEMO/trunk/src/TOP/PISCES/SED/trcdmp_sed.F90
===================================================================
--- /NEMO/trunk/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 10225)
+++ /NEMO/trunk/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 10225)
@@ -0,0 +1,157 @@
+MODULE trcdmp_sed
+   !!======================================================================
+   !!                       ***  MODULE  trcdmp  ***
+   !! Ocean physics: internal restoring trend on passive tracers
+   !!======================================================================
+   !! History :  OPA  !  1991-03  (O. Marti, G. Madec)  Original code
+   !!                 !  1996-01  (G. Madec) statement function for e3
+   !!                 !  1997-05  (H. Loukos)  adapted for passive tracers
+   !!    NEMO    9.0  !  2004-03  (C. Ethe)    free form + modules
+   !!            3.2  !  2007-02  (C. Deltel)  Diagnose ML trends for passive tracers
+   !!            3.3  !  2010-06  (C. Ethe, G. Madec) merge TRA-TRC 
+   !!----------------------------------------------------------------------
+#if  defined key_top
+   !!----------------------------------------------------------------------
+   !!   trc_dmp      : update the tracer trend with the internal damping
+   !!   trc_dmp_init : initialization, namlist read, parameters control
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! ocean dynamics and tracers variables
+   USE trc             ! ocean passive tracers variables
+   USE sed , ONLY : dtsed => dtsed      ! ocean dynamics and tracers variables
+   USE trc             ! ocean passive tracers variables
+   USE trcdta
+   USE prtctl_trc      ! Print control for debbuging
+   USE iom
+
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC trc_dmp_sed      
+   PUBLIC trc_dmp_sed_alloc  
+   PUBLIC trc_dmp_sed_ini    
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   restosed   ! restoring coeff. on tracers (s-1)
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcdmp.F90 7646 2017-02-06 09:25:03Z timgraham $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION trc_dmp_sed_alloc()
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE trc_dmp_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( restosed(jpi,jpj,jpk) , STAT=trc_dmp_sed_alloc )
+      !
+      IF( trc_dmp_sed_alloc /= 0 )   CALL ctl_warn('trc_dmp_sed_alloc: failed to allocate array')
+      !
+   END FUNCTION trc_dmp_sed_alloc
+
+
+   SUBROUTINE trc_dmp_sed( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE trc_dmp_sed  ***
+      !!                  
+      !! ** Purpose :   Compute the passive tracer trend due to a newtonian damping
+      !!      of the tracer field towards given data field and add it to the
+      !!      general tracer trends.
+      !!
+      !! ** Method  :   Newtonian damping towards trdta computed 
+      !!      and add to the general tracer trends:
+      !!                     trn = tra + restotr * (trdta - trb)
+      !!         The trend is computed either throughout the water column
+      !!      (nlmdmptr=0) or in area of weak vertical mixing (nlmdmptr=1) or
+      !!      below the well mixed layer (nlmdmptr=2)
+      !!
+      !! ** Action  : - update the tracer trends tra with the newtonian 
+      !!                damping trends.
+      !!              - save the trends ('key_trdmxl_trc')
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER ::   ji, jj, jk, jn, jl, ikt   ! dummy loop indices
+      CHARACTER (len=22) ::   charout
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   ztrcdta   ! 3D  workspace
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('trc_dmp_sed')
+      !
+      !
+      IF( nb_trcdta > 0 ) THEN  ! Initialisation of tracer from a file that may also be used for damping
+         !
+         DO jn = 1, jptra                                           ! tracer loop
+            !                                                       ! ===========
+            IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
+               !
+               jl = n_trc_index(jn) 
+               CALL trc_dta( kt, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
+               !
+               DO jj = 1, jpj
+                  DO ji = 1, jpi   ! vector opt.
+                     ikt = mbkt(ji,jj)
+                     trb(ji,jj,ikt,jn) = ztrcdta(ji,jj,ikt) + ( trb(ji,jj,ikt,jn) -  ztrcdta(ji,jj,ikt) )     &
+                     &                  * exp( -restosed(ji,jj,ikt) * dtsed )
+                  END DO
+               END DO
+               ! 
+            ENDIF
+         END DO                                                     ! tracer loop
+         !                                                          ! ===========
+      ENDIF
+      !
+      !                                          ! print mean trends (used for debugging)
+      IF( ln_ctl ) THEN
+         WRITE(charout, FMT="('dmp ')")
+         CALL prt_ctl_trc_info(charout)
+         CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      ENDIF
+      !
+      IF( ln_timing )  CALL timing_stop('trc_dmp_sed')
+      !
+   END SUBROUTINE trc_dmp_sed
+
+
+   SUBROUTINE trc_dmp_sed_ini
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_dmp_ini  ***
+      !! 
+      !! ** Purpose :   Initialization for the newtonian damping 
+      !!
+      !! ** Method  :   read the nammbf namelist and check the parameters
+      !!              called by trc_dmp at the first timestep (nittrc000)
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )  CALL timing_start('trc_dmp_sed_ini')
+
+      IF (lwp) WRITE(numout,*) '   tracer damping throughout the water column'
+      !
+      IF( trc_dmp_sed_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'trc_dmp_sed_ini: unable to allocate arrays' )
+      !
+      IF( .NOT.lk_c1d ) THEN
+         !Read in mask from file
+          restosed(:,:,:) = 0.5 / rday
+         !
+      ENDIF
+      IF( ln_timing )  CALL timing_stop('trc_dmp_sed_ini')
+     !
+   END SUBROUTINE trc_dmp_sed_ini
+
+#else
+   !!----------------------------------------------------------------------
+   !!  Dummy module :                                     No passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_dmp_sed( kt )        ! Empty routine
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trc_dmp_sed: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_dmp_sed
+#endif
+
+   !!======================================================================
+END MODULE trcdmp_sed
