Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedini.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedini.F90	(revision 10165)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/SED/sedini.F90	(revision 10345)
@@ -4,5 +4,4 @@
    !! Sediment : define sediment variables
    !!=====================================================================
-#if defined key_sed
 
    !!----------------------------------------------------------------------
@@ -11,10 +10,12 @@
    !! * Modules used
    USE sed     ! sediment global variable
-   USE seddta
-   USE sedrst
-   USE sedco3
-   USE sedchem
+   USE sed_oce
    USE sedarr
+   USE sedadv
+   USE trc_oce, ONLY : nn_dttrc
+   USE trcdmp_sed
+   USE trcdta
    USE iom
+   USE lib_mpp         ! distribued memory computing library
 
 
@@ -24,31 +25,55 @@
    !! Module variables
    REAL(wp)    ::  &
-      sisat   =  800.  ,  &  !: saturation for si    [ mol.l-1]
-      claysat =    0.  ,  &  !: saturation for clay  [ mol.l-1]
+      sedzmin = 0.3    ,  &  !: Minimum vertical spacing
+      sedhmax = 10.0   ,  &  !: Maximum depth of the sediment
+      sedkth  = 5.0    ,  &  !: Default parameters
+      sedacr  = 3.0          !: Default parameters
+      
+   REAL(wp)    ::  &
+      porsurf =  0.95  ,  &  !: Porosity at the surface
+      porinf  =  0.75  ,  &  !: Porosity at infinite depth
+      rhox    =  2.0         !: Vertical length scale of porosity variation 
+
+   REAL(wp)    ::  &
       rcopal  =   40.  ,  &  !: reactivity for si    [l.mol-1.an-1]
-      rcclay  =    0.  ,  &  !: reactivity for clay  [l.mol-1.an-1]
       dcoef   =  8.e-6       !: diffusion coefficient (*por)   [cm**2/s]
 
+   REAL(wp), PUBLIC    ::  &
+      redO2    =  172.  ,  &  !: Redfield coef for Oxygen
+      redNo3   =   16.  ,  &  !: Redfield coef for Nitrate
+      redPo4   =    1.  ,  &  !: Redfield coef for Phosphate
+      redC     =  122.  ,  &  !: Redfield coef for Carbon
+      redfep   =  0.175 ,  &  !: Ratio for iron bound phosphorus
+      rcorgl   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]    
+      rcorgs   =   0.5  ,  &  !: reactivity of the semi-labile component
+      rcorgr   =   1E-4 ,  &  !: reactivity of the refractory component
+      rcnh4    =   10E6 ,  &  !: reactivity for O2/NH4 [l.mol-1.an-1]  
+      rch2s    =   1.E5 ,  &  !: reactivity for O2/ODU [l.mol-1.an-1] 
+      rcfe2    =   5.E8 ,  &  !: reactivity for O2/Fe2+ [l.mol-1.an-1]
+      rcfeh2s  =   1.E4 ,  &  !: Reactivity for FEOH/H2S [l.mol-1.an-1]
+      rcfes    =   1.E5 ,  &  !: Reactivity for FE2+/H2S [l.mol-1.an-1]
+      rcfeso   =   3.E5 ,  &  !: Reactivity for FES/O2 [l.mol-1.an-1]
+      xksedo2  =   5E-6 ,  &  !: half-sturation constant for oxic remin.
+      xksedno3 =   5E-6 ,  &  !: half-saturation constant for denitrification
+      xksedfeo =   0.6  ,  & !: half-saturation constant for iron remin
+      xksedso4 =   2E-3       !: half-saturation constant for SO4 remin
+
    REAL(wp)    ::  &
-      redO2   =  172.  ,  &  !: Redfield coef for Oxygene
-      redNo3  =   16.  ,  &  !: Redfield coef for Nitrates
-      redPo4  =    1.  ,  &  !: Redfield coef for Phosphates
-      redC    =  122.  ,  &  !: Redfield coef for Carbone
-      redDnit =  97.6  ,  &  !: Redfield coef for denitrification    
-      rcorg   =   50.  ,  &  !: reactivity for POC/O2 [l.mol-1.an-1]    
-      o2seuil =    1.  ,  &  !: threshold O2 concentration for [mol.l-1]     
-      rcorgN  =   50.       !: reactivity for POC/No3 [l.mol-1.an-1]
-
-   REAL(wp)    ::  &
-      rccal   = 1000.        !: reactivity for calcite         [l.mol-1.an-1]
-
-   REAL(wp)    ::  &
-      dbiot   = 15.          !: coefficient for bioturbation    [cm**2.(1000an-1)]
-
-   LOGICAL     :: &
-      ln_rst_sed = .TRUE.    !: initialisation from a restart file or not
-
+      rccal   = 1000.,      & !: reactivity for calcite         [l.mol-1.an-1]
+      rcligc  = 1.E-4         !: L/C ratio in POC
+
+   REAL(wp), PUBLIC    ::  dbiot   = 15. , &  !: coefficient for bioturbation    [cm**2.(n-1)]
+      dbtbzsc =  10.0  ,    &  !: Vertical scale of variation. If no variation, mixed layer in the sed [cm]
+      xirrzsc = 2.0            !: Vertical scale of irrigation variation.
    REAL(wp)    ::  &
       ryear = 365. * 24. * 3600. !:  1 year converted in second
+
+   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff1
+   DATA diff1/4.59E-6, 1.104E-5, 4.81E-6 , 9.78E-6, 3.58E-6, 4.01E-6, 9.8E-6, 9.73E-6, 5.0E-6, 3.31E-6 /
+
+   REAL(wp), DIMENSION(jpwat), PUBLIC  :: diff2
+   DATA diff2/1.74E-7, 4.47E-7, 2.51E-7, 3.89E-7, 1.77E-7, 2.5E-7, 3.89E-7, 3.06E-7, 2.5E-7, 1.5E-7 /
+
+
 
    !! *  Routine accessibility
@@ -76,10 +101,5 @@
       !!        !  06-07  (C. Ethe)  Re-organization
       !!----------------------------------------------------------------------
-      INTEGER :: ji, jj, ikt
-#if defined key_sed_off
-      INTEGER  :: numblt         
-      INTEGER  :: nummsh   
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zdta
-#endif
+      INTEGER :: ji, jj, ikt, ierr
       !!----------------------------------------------------------------------
 
@@ -90,53 +110,27 @@
       CALL ctl_opn( numsed, 'sediment.output', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
 
-      WRITE(numsed,*)
-      WRITE(numsed,*) '                 L S C E - C E A'
-      WRITE(numsed,*) '                 SEDIMENT model'
-      WRITE(numsed,*) '                version 2.0  (2007) '
-      WRITE(numsed,*)
-      WRITE(numsed,*)
-
-    
-      WRITE(numsed,*) ' sed_init : Initialization of sediment module  '
-      WRITE(numsed,*) ' '
-
-      !                                ! Allocate LOBSTER arrays
-      IF( sed_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' )
-
+      IF (lwp) THEN
+         WRITE(numsed,*)
+         WRITE(numsed,*) '                 PISCES framework'
+         WRITE(numsed,*) '                 SEDIMENT model'
+         WRITE(numsed,*) '                version 3.0  (2018) '
+         WRITE(numsed,*)
+         WRITE(numsed,*)
+      ENDIF
+
+      IF(lwp) WRITE(numsed,*) ' sed_init : Initialization of sediment module  '
+      IF(lwp) WRITE(numsed,*) ' '
+
+      ! Read sediment Namelist
+      !-------------------------
+      CALL sed_init_nam
+
+      ! Allocate SEDIMENT arrays
+      ierr =        sed_alloc()
+      ierr = ierr + sed_oce_alloc()
+      ierr = ierr + sed_adv_alloc() 
+      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'sed_ini: unable to allocate sediment model arrays' )
 
       ! Determination of sediments number of points and allocate global variables
-#if defined key_sed_off
-      ! Reading Netcdf Pisces file for deepest water layer thickness [m]
-      ! This data will be used as mask to merdge space in 1D vector
-      !----------------------------------------------------------------
-
-      CALL iom_open ( 'mesh_mask'  , nummsh )   
-      CALL iom_open ( 'e3tbot'     , numblt )
-
-      ! mask sediment points for outputs
-      CALL iom_get( nummsh, jpdom_data, 'tmask' , tmask )
-      CALL iom_get( nummsh, jpdom_data, 'mbathy', sbathy )
-      
-      ! longitude/latitude values
-      CALL iom_get ( nummsh, jpdom_data,'nav_lon', glamt(:,:) )
-      CALL iom_get ( nummsh, jpdom_data,'nav_lat', gphit(:,:) )
-    
-      ! bottom layer thickness
-      ALLOCATE( zdta(jpi,jpj) )
-      CALL iom_get  ( numblt, jpdom_data, 'E3TBOT', zdta(:,:) )
-      epkbot(:,:) = 0.
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ikt = MAX( INT( sbathy(ji,jj) ), 1 )
-            IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = zdta(ji,jj)
-         ENDDO
-      ENDDO
-
-      CALL iom_close( nummsh )  
-      CALL iom_close( numblt ) 
-
-      DEALLOCATE( zdta )
-#else
-
       epkbot(:,:) = 0.
       DO jj = 1, jpj
@@ -144,13 +138,17 @@
             ikt = mbkt(ji,jj) 
             IF( tmask(ji,jj,ikt) == 1 ) epkbot(ji,jj) = e3t_1d(ikt)
+            gdepbot(ji,jj) = gdepw_0(ji,jj,ikt)
          ENDDO
       ENDDO
-#endif     
-
 
       ! computation of total number of ocean points
       !--------------------------------------------
-      jpoce  = COUNT( epkbot(:,:) > 0. )
-
+      sedmask = 0.
+      IF ( COUNT( epkbot(:,:) > 0. ) == 0 ) THEN 
+          sedmask = 0.
+      ELSE
+          sedmask = 1.
+      ENDIF 
+      jpoce  = MAX( COUNT( epkbot(:,:) > 0. ) , 1 )
 
       ! Allocate memory size of global variables
@@ -159,18 +157,24 @@
       ALLOCATE( rainrm(jpoce,jpsol) )        ;  ALLOCATE( rainrg(jpoce,jpsol) )        ;  ALLOCATE( raintg(jpoce) ) 
       ALLOCATE( dzdep(jpoce) )               ;  ALLOCATE( iarroce(jpoce) )             ;  ALLOCATE( dzkbot(jpoce) )
-      ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) )                
+      ALLOCATE( zkbot(jpoce) )               ;  ALLOCATE( db(jpoce,jpksed) )
+      ALLOCATE( temp(jpoce) )                ;  ALLOCATE( salt(jpoce) )  
+      ALLOCATE( diff(jpoce,jpksed,jpwat ) )  ;  ALLOCATE( irrig(jpoce, jpksed) )
+      ALLOCATE( wacc(jpoce) )                ;  ALLOCATE( fecratio(jpoce) )
       ALLOCATE( press(jpoce) )               ;  ALLOCATE( densSW(jpoce) ) 
-      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) ) 
+      ALLOCATE( hipor(jpoce,jpksed) )        ;  ALLOCATE( co3por(jpoce,jpksed) )
       ALLOCATE( dz3d(jpoce,jpksed) )         ;  ALLOCATE( volw3d(jpoce,jpksed) )       ;  ALLOCATE( vols3d(jpoce,jpksed) )
+      ALLOCATE( sedligand(jpoce, jpksed) )
 
       ! Initialization of global variables
-      pwcp  (:,:,:) = 0.  ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0.  
-      solcp (:,:,:) = 0.  ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0.
-      rainrm(:,:  ) = 0.  ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0. 
-      dzdep (:    ) = 0.  ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0.
-      temp  (:    ) = 0.  ;  salt   (:   ) = 0.  
-      press (:    ) = 0.  ;  densSW (:   ) = 0. 
-      hipor (:,:  ) = 0.  ;  co3por (:,: ) = 0.  
-      dz3d  (:,:  ) = 0.  ;  volw3d (:,: ) = 0.  ;  vols3d   (:,:) = 0. 
+      pwcp  (:,:,:) = 0.   ;  pwcp0 (:,:,:) = 0.  ; pwcp_dta  (:,:) = 0.  
+      solcp (:,:,:) = 0.   ;  solcp0(:,:,:) = 0.  ; rainrm_dta(:,:) = 0.
+      rainrm(:,:  ) = 0.   ;  rainrg(:,:  ) = 0.  ; raintg    (:  ) = 0. 
+      dzdep (:    ) = 0.   ;  iarroce(:   ) = 0   ; dzkbot    (:  ) = 0.
+      temp  (:    ) = 0.   ;  salt   (:   ) = 0.  ; zkbot     (:  ) = 0.
+      press (:    ) = 0.   ;  densSW (:   ) = 0.  ; db        (:,:) = 0. 
+      hipor (:,:  ) = 0.   ;  co3por (:,: ) = 0.  ; irrig     (:,:) = 0. 
+      dz3d  (:,:  ) = 0.   ;  volw3d (:,: ) = 0.  ; vols3d    (:,:) = 0. 
+      fecratio(:)   = 1E-5 
+      sedligand(:,:) = 0.6E-9
 
       ! Chemical variables      
@@ -178,10 +182,13 @@
       ALLOCATE( ak1ps (jpoce) )  ;  ALLOCATE( ak2ps  (jpoce) )  ;  ALLOCATE( ak3ps (jpoce) ) ;  ALLOCATE( aksis (jpoce) )    
       ALLOCATE( aksps (jpoce) )  ;  ALLOCATE( ak12s  (jpoce) )  ;  ALLOCATE( ak12ps(jpoce) ) ;  ALLOCATE( ak123ps(jpoce) )    
-      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )    
+      ALLOCATE( borats(jpoce) )  ;  ALLOCATE( calcon2(jpoce) )  ;  ALLOCATE( sieqs (jpoce) ) 
+      ALLOCATE( aks3s(jpoce) )   ;  ALLOCATE( akf3s(jpoce) )    ;  ALLOCATE( sulfats(jpoce) )
+      ALLOCATE( fluorids(jpoce) )
 
       akbs  (:) = 0. ;   ak1s   (:) = 0. ;  ak2s  (:) = 0. ;   akws   (:) = 0.
       ak1ps (:) = 0. ;   ak2ps  (:) = 0. ;  ak3ps (:) = 0. ;   aksis  (:) = 0.
       aksps (:) = 0. ;   ak12s  (:) = 0. ;  ak12ps(:) = 0. ;   ak123ps(:) = 0.
-      borats(:) = 0. ;   calcon2(:) = 0.
+      borats(:) = 0. ;   calcon2(:) = 0. ;  sieqs (:) = 0.
+      aks3s(:)  = 0. ;   akf3s(:)   = 0. ;  sulfats(:) = 0. ;  fluorids(:) = 0.
 
       ! Mass balance calculation  
@@ -191,24 +198,28 @@
       fromsed(:,:) = 0.    ;   tosed(:,:) = 0. ;  rloss(:,:) = 0.  ;   tokbot(:,:) = 0. 
 
-      ! Read sediment Namelist
-      !-------------------------
-      CALL sed_init_nam
-    
       ! Initialization of sediment geometry
       !------------------------------------
       CALL sed_init_geom
 
-
-      ! sets initial sediment composition
-      ! ( only clay or reading restart file )
-      !---------------------------------------
-      CALL sed_init_data
-
-
-      CALL sed_init_wri
-
+      ! Offline specific mode
+      ! ---------------------
+      ln_sediment_offline = .FALSE.
+
+#if defined key_sed_off
+      ln_sediment_offline = .TRUE.
+      IF (lwp) write(numsed,*) 'Sediment module is run in offline mode'
+      IF (lwp) write(numsed,*) 'key_sed_off is activated at compilation time'
+      IF (lwp) write(numsed,*) 'ln_sed_2way is forced to false'
+      IF (lwp) write(numsed,*) '--------------------------------------------'
+      ln_sed_2way = .FALSE.
+#endif
+      ! Initialisation of tracer damping
+      ! --------------------------------
+      IF (ln_sediment_offline) THEN
+         CALL trc_dta_ini(jptra)
+         CALL trc_dmp_sed_ini
+      ENDIF
 
    END SUBROUTINE sed_init
-
 
    SUBROUTINE sed_init_geom
@@ -227,15 +238,11 @@
       !! * Modules used
       !! * local declarations
-      INTEGER ::  &
-        ji, jj, jk
- 
-#if defined key_sed_off
-      REAL(wp) , DIMENSION (jpi,jpj) :: zdta
-      INTEGER  ::  numpres
-#endif
+      INTEGER  :: ji, jj, jk, jn
+      REAL(wp) :: za0, za1, zt, zw, zsum, zsur, zprof, zprofw
+      REAL(wp) :: ztmp1, ztmp2
       !---------------------------------------------------------- 
 
-      WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '
-      WRITE(numsed,*) ' '
+      IF(lwp) WRITE(numsed,*) ' sed_init_geom : Initialization of sediment geometry '
+      IF(lwp) WRITE(numsed,*) ' '
 
       ! Computation of 1D array of sediments points
@@ -250,37 +257,16 @@
       END DO
 
+      IF ( indoce .EQ. 0 ) THEN
+         indoce = 1
+         iarroce(indoce) = 1
+      ENDIF
+
       IF( indoce .NE. jpoce ) THEN
-         WRITE(numsed,*) ' '
-         WRITE(numsed,*) 'Warning : number of ocean points indoce = ', indoce, &
-            &     ' doesn''t match  number of point where epkbot>0  jpoce = ', jpoce
-         WRITE(numsed,*) ' '
-         WRITE(numsed,*) ' '
-         STOP
+         CALL ctl_stop( 'STOP', 'sed_ini: number of ocean points indoce doesn''t match  number of point' )
       ELSE
-         WRITE(numsed,*) ' '
-         WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce
-         WRITE(numsed,*) ' '
-      ENDIF
-
-#if defined key_sed_off
-
-      ! Reading Netcdf Pisces file for deepest water layer thickness [m]
-      ! This data will be used as mask to merdge space in 1D vector
-      !----------------------------------------------------------------
-      CALL iom_open ( 'pressbot'   , numpres )
-      
-      ! pressure  in  bars      
-      CALL iom_get ( numpres, jpdom_data,'BATH', zdta(:,:) )
-      CALL pack_arr ( jpoce, press(1:jpoce), zdta(1:jpi,1:jpj), iarroce(1:jpoce) )
-      press(1:jpoce) = press(1:jpoce) * 1.025e-1 
-
-      CALL iom_close ( numpres )
-#endif
-
-
-      ! mask sediment points for outputs
-      DO jk = 1, jpksed
-         tmasksed(:,:,jk) = tmask(:,:,1)
-      ENDDO
+         IF (lwp) WRITE(numsed,*) ' '
+         IF (lwp) WRITE(numsed,*) ' total number of ocean points jpoce =  ',jpoce
+         IF (lwp) WRITE(numsed,*) ' '
+      ENDIF
 
       ! initialization of dzkbot in [cm]
@@ -288,4 +274,5 @@
       CALL pack_arr ( jpoce, dzkbot(1:jpoce), epkbot(1:jpi,1:jpj), iarroce(1:jpoce) )
       dzkbot(1:jpoce) = dzkbot(1:jpoce) * 1.e+2 
+      CALL pack_arr ( jpoce, zkbot(1:jpoce), gdepbot(1:jpi,1:jpj), iarroce(1:jpoce) )
 
       ! Geometry and  constants 
@@ -293,15 +280,23 @@
       ! (1st layer= diffusive layer = pur water) 
       !------------------------------------------
-      dz(1)  = 0.1
-      dz(2)  = 0.3
-      dz(3)  = 0.3
-      dz(4)  = 0.5
-      dz(5)  = 0.5
-      dz(6)  = 0.5
-      dz(7)  = 1.
-      dz(8)  = 1.
-      dz(9)  = 1.
-      dz(10) = 2.45
-      dz(11) = 2.45
+      za1  = (  sedzmin - sedhmax / FLOAT(jpksed-1)  )                                                      &
+         & / ( TANH((1-sedkth)/sedacr) - sedacr/FLOAT(jpksed-1) * (  LOG( COSH( (jpksed - sedkth) / sedacr) )      &
+         &                                                   - LOG( COSH( ( 1  - sedkth) / sedacr) )  )  )
+      za0  = sedzmin - za1 * TANH( (1-sedkth) / sedacr )
+      zsur = - za0 - za1 * sedacr * LOG( COSH( (1-sedkth) / sedacr )  )
+
+      profsedw(1) = 0.0
+      profsed(1) = -dz(1) / 2.
+      DO jk = 2, jpksed
+         zw = REAL( jk , wp )
+         zt = REAL( jk , wp ) - 0.5_wp
+         profsed(jk)  = ( zsur + za0 * zt + za1 * sedacr * LOG ( COSH( (zt-sedkth) / sedacr ) )  ) 
+         profsedw(jk) = ( zsur + za0 * zw + za1 * sedacr * LOG ( COSH( (zw-sedkth) / sedacr ) )  )
+      END DO
+
+      dz(1) = 0.1
+      DO jk = 2, jpksed
+         dz(jk) = profsedw(jk) - profsedw(jk-1)
+      END DO
 
       DO jk = 1, jpksed
@@ -311,24 +306,10 @@
       ENDDO
 
-      ! Depth(jk)= depth of middle of each layer
-      !----------------------------------------
-      profsed(1) = -dz(1)/ 2.
-      DO jk = 2, jpksed
-         profsed(jk) = profsed(jk-1) + dz(jk-1) / 2. + dz(jk) / 2.
-      ENDDO
-
       !  Porosity profile [0]
       !---------------------
-      por(1)  = 1.
-      por(2)  = 0.95
-      por(3)  = 0.9
-      por(4)  = 0.85
-      por(5)  = 0.81
-      por(6)  = 0.75
-      por(7)  = 0.75      
-      por(8)  = 0.75      
-      por(9)  = 0.75      
-      por(10) = 0.75
-      por(11) = 0.75
+      por(1) = 1.0
+      DO jk = 2, jpksed
+         por(jk) = porinf + ( porsurf-porinf) * exp(-rhox * (profsed(jk) ) )
+      END DO
  
       ! inverse of  Porosity profile
@@ -353,9 +334,7 @@
       dz3d(1:jpoce,1) = dz(1)
 
-
       !---------------------------------------------
       ! Molecular weight [g/mol] for solid species
       !---------------------------------------------
-
 
       ! opal=sio2*0.4(h20)=28+2*16+0.4*(2+16)
@@ -371,4 +350,8 @@
          &              0.59 * 27. + 10. * 16.
 
+      mol_wgt(jsfeo)  = 55.0 + 3.0 * ( 16.0 + 1.0)
+
+      mol_wgt(jsfes)  = 55.0 + 32.0
+
       ! for chemistry Poc : C(122)H(244)O(86)N(16)P(1)
       ! But den sity of Poc is an Hydrated material (= POC + 30H2O)
@@ -377,4 +360,6 @@
       mol_wgt(jspoc) = ( 122. * 12. + 355. + 120. * 16.+ &
          &                16. * 14. + 31. ) / 122.
+      mol_wgt(jspos) = mol_wgt(jspoc)
+      mol_wgt(jspor) = mol_wgt(jspoc)
 
       ! CaCO3
@@ -384,15 +369,28 @@
       ! Density of solid material in sediment [g/cm**3]
       !------------------------------------------------
-      dens = 2.6
+      denssol = 2.6
 
       ! Initialization of diffusion coefficient as function of porosity [cm**2/s]
       !--------------------------------------------------------------------
-      diff(:) = dcoef * por(:)
+!      DO jn = 1, jpsol
+!         DO jk = 1, jpksed
+!            DO ji = 1, jpoce
+!               diff(ji,jk,jn) = dcoef / ( 1.0 - 2.0 * log(por(jk)) )
+!            END DO
+!         END DO
+!      END DO
+
+      ! Accumulation rate from Burwicz et al. (2011). This is used to
+      ! compute the flux of clays and minerals
+      ! --------------------------------------------------------------
+      DO ji = 1, jpoce
+          ztmp1 = 0.117 / ( 1.0 + ( zkbot(ji) / 200.)**3 )
+          ztmp2 = 0.006 / ( 1.0 + ( zkbot(ji) / 4000.)**10 )
+          wacc(ji) = ztmp1 + ztmp2
+      END DO
 
 
       ! Initialization of time step as function of porosity [cm**2/s]
       !------------------------------------------------------------------
-      rdtsed(2:jpksed) = dtsed / ( dens * por1(2:jpksed) )
-     
    END SUBROUTINE sed_init_geom
 
@@ -408,6 +406,8 @@
       !!----------------------------------------------------------------------
 
-      INTEGER :: &
-         numnamsed = 28
+      INTEGER ::   numnamsed_ref = -1           !! Logical units for namelist sediment
+      INTEGER ::   numnamsed_cfg = -1           !! Logical units for namelist sediment
+      INTEGER :: ios                 ! Local integer output status for namelist read
+      CHARACTER(LEN=20)   ::   clname
 
       TYPE PSED
@@ -422,57 +422,84 @@
       TYPE(PSED) , DIMENSION(jpdia2dsed) :: seddiag2d
 
-      NAMELIST/nam_time/nfreq
+      NAMELIST/nam_run/nrseddt,ln_sed_2way
+      NAMELIST/nam_geom/jpksed, sedzmin, sedhmax, sedkth, sedacr, porsurf, porinf, rhox
       NAMELIST/nam_trased/sedsol, sedwat
       NAMELIST/nam_diased/seddiag3d, seddiag2d
-      NAMELIST/nam_reac/sisat, claysat, rcopal, rcclay, dcoef 
-      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redDnit, &
-         &           rcorg, o2seuil, rcorgN
-      NAMELIST/nam_cal/rccal
-      NAMELIST/nam_dc13/pdb, rc13P, rc13Ca
-      NAMELIST/nam_btb/dbiot
-      NAMELIST/nam_rst/ln_rst_sed
-
-      INTEGER :: jn, jn1
+      NAMELIST/nam_inorg/rcopal, dcoef, rccal, ratligc, rcligc
+      NAMELIST/nam_poc/redO2, redNo3, redPo4, redC, redfep, rcorgl, rcorgs,  &
+         &             rcorgr, rcnh4, rch2s, rcfe2, rcfeh2s, rcfes, rcfeso,  &
+         &             xksedo2, xksedno3, xksedfeo, xksedso4
+      NAMELIST/nam_btb/dbiot, ln_btbz, dbtbzsc, adsnh4, ln_irrig, xirrzsc
+      NAMELIST/nam_rst/ln_rst_sed, cn_sedrst_indir, cn_sedrst_outdir, cn_sedrst_in, cn_sedrst_out
+
+      INTEGER :: ji, jn, jn1
       !-------------------------------------------------------
 
-      WRITE(numsed,*) ' sed_init_nam : Read namelists '
-      WRITE(numsed,*) ' '
+      IF(lwp) WRITE(numsed,*) ' sed_init_nam : Read namelists '
+      IF(lwp) WRITE(numsed,*) ' '
 
       ! ryear = 1 year converted in second
       !------------------------------------
-      WRITE(numsed,*) ' '
-      WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
-      WRITE(numsed,*) ' '     
+      IF (lwp) THEN
+         WRITE(numsed,*) ' '
+         WRITE(numsed,*) 'number of seconds in one year : ryear = ', ryear
+         WRITE(numsed,*) ' '     
+      ENDIF
 
       ! Reading namelist.sed variables
       !---------------------------------
-      CALL ctl_opn( numnamsed, 'namelist.sediment', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
-
-      dtsed = rdt
+      clname = 'namelist_sediment'
+      IF(lwp) WRITE(numsed,*) ' sed_init_nam : read SEDIMENT namelist'
+      IF(lwp) WRITE(numsed,*) ' ~~~~~~~~~~~~~~'
+      CALL ctl_opn( numnamsed_ref, TRIM( clname )//'_ref', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
+      CALL ctl_opn( numnamsed_cfg, TRIM( clname )//'_cfg', 'OLD'    , 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
+
       nitsed000 = nittrc000
       nitsedend = nitend
-#if ! defined key_sed_off
-      nwrised   = nwritetrc
-#else
-      nwrised   = nwrite
-#endif
-     ! Diffraction/reaction parameters
-      !----------------------------------
-      READ( numnamsed, nam_time )
-      WRITE(numsed,*) ' namelist nam_time'
-
-#if ! defined key_sed_off
-      nfreq     = 1
-#endif
-
-      WRITE(numsed,*) ' sedimentation time step              dtsed      = ', dtsed
-      WRITE(numsed,*) ' 1st time step for sediment.          nitsed000  = ', nitsed000
-      WRITE(numsed,*) ' last time step for sediment.         nitsedend  = ', nitsedend
-      WRITE(numsed,*) ' frequency of sediment outputs        nwrised    = ', nwrised
-      WRITE(numsed,*) ' frequency of restoring inputs data   nfreq      = ', nfreq
-      WRITE(numsed,*) ' '
-
-      REWIND( numnamsed )               ! read nattrc
-      READ  ( numnamsed, nam_trased )
+      ! Namelist nam_run
+      REWIND( numnamsed_ref )              ! Namelist nam_run in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_run, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_run in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_run, IOSTAT = ios, ERR = 902)
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_run in configuration namelist', lwp )
+
+      IF (lwp) THEN
+         WRITE(numsed,*) ' namelist nam_run'
+         WRITE(numsed,*) ' Nb of iterations for fast species    nrseddt = ', nrseddt
+         WRITE(numsed,*) ' 2-way coupling between PISCES and Sed ln_sed_2way = ', ln_sed_2way
+      ENDIF
+
+      REWIND( numnamsed_ref )              ! Namelist nam_geom in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_geom, IOSTAT = ios, ERR = 903)
+903   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_geom in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_geom, IOSTAT = ios, ERR = 904)
+904   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_geom in configuration namelist', lwp )
+
+      IF (lwp) THEN 
+         WRITE(numsed,*) ' namelist nam_geom'
+         WRITE(numsed,*) ' Number of vertical layers            jpksed  = ', jpksed
+         WRITE(numsed,*) ' Minimum vertical spacing             sedzmin = ', sedzmin
+         WRITE(numsed,*) ' Maximum depth of the sediment        sedhmax = ', sedhmax
+         WRITE(numsed,*) ' Default parameter                    sedkth  = ', sedkth
+         WRITE(numsed,*) ' Default parameter                    sedacr  = ', sedacr
+         WRITE(numsed,*) ' Sediment porosity at the surface     porsurf = ', porsurf
+         WRITE(numsed,*) ' Sediment porosity at infinite depth  porinf  = ', porinf
+         WRITE(numsed,*) ' Length scale of porosity variation   rhox    = ', rhox
+      ENDIF
+
+      jpksedm1  = jpksed - 1
+      dtsed = r2dttrc
+
+      REWIND( numnamsed_ref )              ! Namelist nam_trased in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_trased, IOSTAT = ios, ERR = 905)
+905   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_trased in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_trased, IOSTAT = ios, ERR = 906)
+906   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_trased in configuration namelist', lwp )
 
       DO jn = 1, jpsol
@@ -489,18 +516,24 @@
       END DO
 
-      WRITE(numsed,*) ' namelist nam_trased'
-      WRITE(numsed,*) ' '
-      DO jn = 1, jptrased
-         WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
-         WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
-         WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
-         WRITE(numsed,*) ' '
-      END DO
-      WRITE(numsed,*) ' '
-
+      IF (lwp) THEN
+         WRITE(numsed,*) ' namelist nam_trased'
+         WRITE(numsed,*) ' '
+         DO jn = 1, jptrased
+            WRITE(numsed,*) 'name of 3d output sediment field number :',jn,' : ',TRIM(sedtrcd(jn))
+            WRITE(numsed,*) 'long name ', TRIM(sedtrcl(jn))
+            WRITE(numsed,*) ' in unit = ', TRIM(sedtrcu(jn))
+            WRITE(numsed,*) ' '
+         END DO
+         WRITE(numsed,*) ' '
+      ENDIF
+
+      REWIND( numnamsed_ref )              ! Namelist nam_diased in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_diased, IOSTAT = ios, ERR = 907)
+907   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_diased in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_diased, IOSTAT = ios, ERR = 908)
+908   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diased in configuration namelist', lwp )
       
-      REWIND( numnamsed )
-      READ( numnamsed, nam_diased )
-
       DO jn = 1, jpdia3dsed
          seddia3d(jn) = seddiag3d(jn)%snamesed
@@ -515,352 +548,144 @@
       END DO
 
-      WRITE(numsed,*) ' namelist nam_diased'
-      WRITE(numsed,*) ' '
-      DO jn = 1, jpdia3dsed
-         WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
-         WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
-         WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
-         WRITE(numsed,*) ' '
-      END DO
-
-      DO jn = 1, jpdia2dsed
-         WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
-         WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
-         WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
-         WRITE(numsed,*) ' '
-      END DO
-
-      WRITE(numsed,*) ' '
-
-
-      ! Diffraction/reaction parameters
+      IF (lwp) THEN
+         WRITE(numsed,*) ' namelist nam_diased'
+         WRITE(numsed,*) ' '
+         DO jn = 1, jpdia3dsed
+            WRITE(numsed,*) 'name of 3D output diag number :',jn, ' : ', TRIM(seddia3d(jn))
+            WRITE(numsed,*) 'long name ', TRIM(seddia3l(jn))
+            WRITE(numsed,*) ' in unit = ',TRIM(seddia3u(jn))
+            WRITE(numsed,*) ' '
+         END DO
+
+         DO jn = 1, jpdia2dsed
+            WRITE(numsed,*) 'name of 2D output diag number :',jn, ' : ', TRIM(seddia2d(jn))
+            WRITE(numsed,*) 'long name ', TRIM(seddia2l(jn))
+            WRITE(numsed,*) ' in unit = ',TRIM(seddia2u(jn))
+            WRITE(numsed,*) ' '
+         END DO
+
+         WRITE(numsed,*) ' '
+      ENDIF
+
+      ! Inorganic chemistry parameters
       !----------------------------------
-      REWIND( numnamsed )
-      READ( numnamsed, nam_reac )
-      WRITE(numsed,*) ' namelist nam_reac'
-      WRITE(numsed,*) ' saturation for si    sisat   = ', sisat
-      WRITE(numsed,*) ' saturation for clay  claysat = ', claysat
-      WRITE(numsed,*) ' reactivity for Si    rcopal  = ', rcopal
-      WRITE(numsed,*) ' reactivity for clay  rcclay  = ', rcclay
-      WRITE(numsed,*) ' diff. coef for por.  dcoef   = ', dcoef
-      WRITE(numsed,*) ' '
-
+      REWIND( numnamsed_ref )              ! Namelist nam_inorg in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_inorg, IOSTAT = ios, ERR = 909)
+909   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_inorg in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_inorg, IOSTAT = ios, ERR = 910)
+910   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_inorg in configuration namelist', lwp )
+
+      IF (lwp) THEN
+         WRITE(numsed,*) ' namelist nam_inorg'
+         WRITE(numsed,*) ' reactivity for Si      rcopal  = ', rcopal
+         WRITE(numsed,*) ' diff. coef for por.    dcoef   = ', dcoef
+         WRITE(numsed,*) ' reactivity for calcite rccal   = ', rccal
+         WRITE(numsed,*) ' L/C ratio in POC       ratligc = ', ratligc
+         WRITE(numsed,*) ' reactivity for ligands rcligc  = ', rcligc
+         WRITE(numsed,*) ' '
+      ENDIF
 
       ! Unity conversion to get saturation conc. psat in [mol.l-1]
       ! and reactivity rc in  [l.mol-1.s-1]
       !----------------------------------------------------------
-      sat_sil   = sisat   * 1.e-6   
-      sat_clay  = claysat * 1.e-6
-
-      reac_sil  = rcopal / ryear     
-      reac_clay  = rcclay / ryear
-
+      reac_sil   = rcopal / ryear     
+      reac_ligc  = rcligc / ryear
 
       ! Additional parameter linked to POC/O2/No3/Po4
       !----------------------------------------------
-      REWIND( numnamsed )
-      READ( numnamsed, nam_poc )
-      WRITE(numsed,*) ' namelist nam_poc'
-      WRITE(numsed,*) ' Redfield coef for oxy     redO2   = ', redO2
-      WRITE(numsed,*) ' Redfield coef for no3     redNo3  = ', redNo3
-      WRITE(numsed,*) ' Redfield coef for po4     redPo4  = ', redPo4
-      WRITE(numsed,*) ' Redfield coef for carbone redC    = ', redC
-      WRITE(numsed,*) ' Redfield coef for denitri redDnit = ', redDnit
-      WRITE(numsed,*) ' reactivity for POC/O2     rcorg   = ', rcorg
-      WRITE(numsed,*) ' threshold O2 concen.      o2seuil = ', o2seuil
-      WRITE(numsed,*) ' reactivity for POC/NO3    rcorgN  = ', rcorgN
-      WRITE(numsed,*) ' '
-
-
-      so2ut  = redO2   / redC
-      srno3  = redNo3  / redC
-      spo4r  = redPo4  / redC
-      srDnit = redDnit / redC
-      sthro2 = o2seuil * 1.e-6 ! threshold O2 concen. in [mol.l-1]
+      REWIND( numnamsed_ref )              ! Namelist nam_poc in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_poc, IOSTAT = ios, ERR = 911)
+911   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_poc in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_poc, IOSTAT = ios, ERR = 912)
+912   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_poc in configuration namelist', lwp )
+
+      IF (lwp) THEN
+         WRITE(numsed,*) ' namelist nam_poc'
+         WRITE(numsed,*) ' Redfield coef for oxy            redO2    = ', redO2
+         WRITE(numsed,*) ' Redfield coef for no3            redNo3   = ', redNo3
+         WRITE(numsed,*) ' Redfield coef for po4            redPo4   = ', redPo4
+         WRITE(numsed,*) ' Redfield coef for carbon         redC     = ', redC
+         WRITE(numsed,*) ' Ration for iron bound P          redfep   = ', redfep
+         WRITE(numsed,*) ' reactivity for labile POC        rcorgl   = ', rcorgl
+         WRITE(numsed,*) ' reactivity for semi-refract. POC rcorgs   = ', rcorgs
+         WRITE(numsed,*) ' reactivity for refractory POC    rcorgr   = ', rcorgr
+         WRITE(numsed,*) ' reactivity for NH4               rcnh4    = ', rcnh4
+         WRITE(numsed,*) ' reactivity for H2S               rch2s    = ', rch2s
+         WRITE(numsed,*) ' reactivity for Fe2+              rcfe2    = ', rcfe2
+         WRITE(numsed,*) ' reactivity for FeOH/H2S          rcfeh2s  = ', rcfeh2s
+         WRITE(numsed,*) ' reactivity for Fe2+/H2S          rcfes    = ', rcfes
+         WRITE(numsed,*) ' reactivity for FeS/O2            rcfeso   = ', rcfeso
+         WRITE(numsed,*) ' Half-sat. cste for oxic remin    xksedo2  = ', xksedo2
+         WRITE(numsed,*) ' Half-sat. cste for denit.        xksedno3 = ', xksedno3
+         WRITE(numsed,*) ' Half-sat. cste for iron remin    xksedfeo = ', xksedfeo
+         WRITE(numsed,*) ' Half-sat. cste for SO4 remin     xksedso4 = ', xksedso4
+         WRITE(numsed,*) ' '
+      ENDIF
+
+
+      so2ut  = redO2    / redC
+      srno3  = redNo3   / redC
+      spo4r  = redPo4   / redC
+      srDnit = ( (redO2 + 32. ) * 0.8 - redNo3 - redNo3 * 0.6 ) / redC
       ! reactivity rc in  [l.mol-1.s-1]
-      reac_poc  = rcorg / ryear
-      reac_no3 = rcorgN  / ryear
-
-
-      ! Carbonate parameters
-      !---------------------
-      READ( numnamsed, nam_cal )
-      WRITE(numsed,*) ' namelist  nam_cal'
-      WRITE(numsed,*) ' reactivity for calcite rccal = ', rccal
-      WRITE(numsed,*) ' '
+      reac_pocl  = rcorgl / ryear
+      reac_pocs  = rcorgs / ryear
+      reac_pocr  = rcorgr / ryear
+      reac_nh4   = rcnh4  / ryear
+      reac_h2s   = rch2s  / ryear
+      reac_fe2   = rcfe2  / ryear
+      reac_feh2s = rcfeh2s/ ryear
+      reac_fes   = rcfes  / ryear
+      reac_feso  = rcfeso / ryear
 
       ! reactivity rc in  [l.mol-1.s-1]      
       reac_cal = rccal / ryear
 
-
-      ! C13  parameters
-      !----------------
-      READ( numnamsed, nam_dc13 )
-      WRITE(numsed,*) ' namelist nam_dc13 ' 
-      WRITE(numsed,*) ' 13C/12C in PD Belemnite        PDB    = ', pdb
-      WRITE(numsed,*) ' 13C/12C in POC = rc13P*PDB     rc13P  = ', rc13P
-      WRITE(numsed,*) ' 13C/12C in CaCO3 = rc13Ca*PDB  rc13Ca = ', rc13Ca
-      WRITE(numsed,*) ' '
-
-      
       ! Bioturbation parameter
       !------------------------
-      READ( numnamsed, nam_btb )
-      WRITE(numsed,*) ' namelist nam_btb ' 
-      WRITE(numsed,*) ' coefficient for bioturbation   dbiot    = ', dbiot
-      WRITE(numsed,*) ' '
-
-      ! Unity convertion to get bioturb coefficient in [cm2.s-1]
-      db = dbiot / ( ryear * 1000. )
+      REWIND( numnamsed_ref )              ! Namelist nam_btb in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_btb, IOSTAT = ios, ERR = 913)
+913   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_btb in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_btb, IOSTAT = ios, ERR = 914)
+914   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_btb in configuration namelist', lwp )
+
+      IF (lwp) THEN
+         WRITE(numsed,*) ' namelist nam_btb ' 
+         WRITE(numsed,*) ' coefficient for bioturbation      dbiot    = ', dbiot
+         WRITE(numsed,*) ' Depth varying bioturbation        ln_btbz  = ', ln_btbz
+         WRITE(numsed,*) ' coefficient for btb attenuation   dbtbzsc  = ', dbtbzsc
+         WRITE(numsed,*) ' Adsorption coefficient of NH4     adsnh4   = ', adsnh4
+         WRITE(numsed,*) ' Bioirrigation in sediment         ln_irrig = ', ln_irrig
+         WRITE(numsed,*) ' coefficient for irrig attenuation xirrzsc  = ', xirrzsc
+         WRITE(numsed,*) ' '
+      ENDIF
 
       ! Initial value (t=0) for sediment pore water and solid components
       !----------------------------------------------------------------
-      READ( numnamsed, nam_rst )
-      WRITE(numsed,*) ' namelist  nam_rst ' 
-      WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 
-      WRITE(numsed,*) ' '
-
-      CLOSE( numnamsed )
+      REWIND( numnamsed_ref )              ! Namelist nam_rst in reference namelist : Pisces variables
+      READ  ( numnamsed_ref, nam_rst, IOSTAT = ios, ERR = 915)
+915   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in reference namelist', lwp )
+
+      REWIND( numnamsed_cfg )              ! Namelist nam_rst in reference namelist : Pisces variables
+      READ  ( numnamsed_cfg, nam_rst, IOSTAT = ios, ERR = 916)
+916   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_rst in configuration namelist', lwp )
+
+      IF (lwp) THEN
+         WRITE(numsed,*) ' namelist  nam_rst ' 
+         WRITE(numsed,*) '  boolean term for restart (T or F) ln_rst_sed = ', ln_rst_sed 
+         WRITE(numsed,*) ' '
+      ENDIF
+      nn_dtsed = nn_dttrc
+
+      CLOSE( numnamsed_cfg )
+      CLOSE( numnamsed_ref )
 
    END SUBROUTINE sed_init_nam
 
-   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
-
-          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.
-
-         ! Initialization of [h+] and [co3--]
-
-         zhipor = 8.
-         ! Initialization of [h+] in mol/kg
-         DO jk = 1, jpksed
-            DO ji = 1, jpoce
-               hipor (ji,jk) = 10.**( -1. * zhipor )
-            ENDDO
-         ENDDO
-
-         co3por(:,:) = 0.
-
-      ELSE   
-  
-         WRITE(numsed,*) ' Initilization of Sediment components from restart'
-
-         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
-
-      WRITE(numsed,*)' '
-      WRITE(numsed,*)'======== Write summary of initial state   ============'
-      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 : ',jpi, jpj
-      WRITE(numsed,*)'jpoce = ',jpoce,' nbtot pts = ',jpij,' nb earth pts = ',jpij - jpoce
-      WRITE(numsed,*)'sublayer thickness dz(1) [cm] : ', dz(1)
-      WRITE(numsed,*)'Coeff diff for k=1 (cm2/s) : ',diff(1)
-      WRITE(numsed,*)' nb solid comp : ',jpsol
-      WRITE(numsed,*)'(1=opal,2=clay,3=POC,4=CaCO3)'
-      WRITE(numsed,*)'weight mol 1,2,3,4'
-      WRITE(numsed,'(4(F0.2,3X))')mol_wgt(jsopal),mol_wgt(jsclay),mol_wgt(jspoc),mol_wgt(jscal)
-      WRITE(numsed,*)'nb dissolved comp',jpwat
-      WRITE(numsed,*)'(1=silicic acid,2="dissolved" clay,3=O2,4=DIC,5=Nitrate,&
-         &6=Phosphates,7=Alk))'
-      WRITE(numsed,*)'Psat (umol/l) for silicic Acid and "dissolved" clay'
-      WRITE(numsed,'(2(F0.2,3X))') sat_sil * 1e+6,  sat_clay * 1e+6
-      WRITE(numsed,*)'reaction rate rc for Op/si,Clay,POC/O2,caco3, POC/No3 (an-1)'
-      WRITE(numsed,'(5(F0.2,3X))') reac_sil * ryear, reac_clay * ryear, reac_poc * ryear, &
-                                   reac_cal * ryear, reac_no3 * ryear
-      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,*)'threshold for stating denitrification [mol/l]'
-      WRITE(numsed,'(1PE8.2)') sthrO2
-      WRITE(numsed,*)'-------------------------------------------------------------------'
-      WRITE(numsed,*)'Min-Max-Mean'
-      WRITE(numsed,*)'For each variable : min, max, moy value observed on selected local zone'
-      WRITE(numsed,*)'-------------------------------------------------------------------'
-      WRITE(numsed,*)'thickness of the last wet layer dzkbot [m]'
-      WRITE(numsed,'(3(F0.2,3X))') MINVAL(dzkbot(1:jpoce))/100.,MAXVAL(dzkbot(1:jpoce))/100.,&
-         &SUM(dzkbot(1:jpoce))/jpoce/100.
-      WRITE(numsed,*)'temp [°C]'
-      WRITE(numsed,'(3(F0.2,3X))') MINVAL(temp(1:jpoce)),MAXVAL(temp(1:jpoce)),&
-         &                         SUM(temp(1:jpoce))/jpoce
-      WRITE(numsed,*)'salt o/oo'
-      WRITE(numsed,'(3(F0.2,3X))')MINVAL(salt(1:jpoce)),MAXVAL(salt(1:jpoce)),&
-         &                        SUM(salt(1:jpoce))/jpoce
-#if defined key_sed_off
-      WRITE(numsed,*)'pressure [bar] (depth in m is about 10*pressure)'
-      WRITE(numsed,'(3(F0.2,3X))') MINVAL(press(1:jpoce)),MAXVAL(press(1:jpoce)),&
-         &                         SUM(press(1:jpoce))/jpoce
-#endif
-      WRITE(numsed,*)'density of Sea Water'
-      WRITE(numsed,'(3(F0.2,3X))') MINVAL(densSW(1:jpoce)), MAXVAL(densSW(1:jpoce)),&
-         &                         SUM(densSW(1:jpoce))/jpoce
-      WRITE(numsed,*)''
-      WRITE(numsed,*)'     Dissolved Components '
-      WRITE(numsed,*)'     ====================='
-      WRITE(numsed,*)'[Si(OH)4] dissolved (1)(k=1)(µmol/l)(and min value in mol/kg  of solution)'
-      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwsil))*1.e+6,&
-         &                         SUM(pwcp(1:jpoce,1,jwsil))*1.e+6/jpoce,&
-         &                         MINVAL(pwcp(1:jpoce,1,jwsil)*1.e+6/densSW(1:jpoce))
-      WRITE(numsed,*)'[O2]    dissolved (3) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
-      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwoxy))*1.e+6,&
-         &SUM(pwcp(1:jpoce,1,jwoxy))*1.e+6/jpoce,&
-         &MINVAL(pwcp(1:jpoce,1,jwoxy)*1.e+6/densSW(1:jpoce))
-      WRITE(numsed,*)'[DIC]    dissolved (4) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
-      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwdic))*1.e+6,&
-         &SUM(pwcp(1:jpoce,1,jwdic))*1.e+6/jpoce,&
-         &MINVAL(pwcp(1:jpoce,1,jwdic)*1.e+6/densSW(1:jpoce))
-      WRITE(numsed,*)'[NO3]    dissolved (5) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
-      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwno3))*1.e+6,&
-         &SUM(pwcp(1:jpoce,1,jwno3))*1.e+6/jpoce,&
-         &MINVAL(pwcp(1:jpoce,1,jwno3)*1.e+6/densSW(1:jpoce))
-      WRITE(numsed,*)'[PO4]    dissolved (6) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
-      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwpo4))*1.e+6,&
-         &SUM(pwcp(1:jpoce,1,jwpo4))*1.e+6/jpoce,&
-         &MINVAL(pwcp(1:jpoce,1,jwpo4)*1.e+6/densSW(1:jpoce))
-      WRITE(numsed,*)'[Alk]    dissolved (7) (k=1)(µequi)(and min value in mol/kg  of solution)'
-      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwalk))*1.e+6,&
-         &SUM(pwcp(1:jpoce,1,jwalk))*1.e+6/jpoce,&
-         &MINVAL(pwcp(1:jpoce,1,jwalk)*1.e+6/densSW(1:jpoce))
-      WRITE(numsed,*)'[DIC13]  dissolved (8) (k=1)(µmol/l)(and min value in mol/kg  of solution)'
-      WRITE(numsed,'(4(F0.3,2X))') MINVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,MAXVAL(pwcp(1:jpoce,1,jwc13))*1.e+6,&
-         &SUM(pwcp(1:jpoce,1,jwc13))*1.e+6/jpoce,&
-         &MINVAL(pwcp(1:jpoce,1,jwc13)*1.e+6/densSW(1:jpoce))
-      WRITE(numsed,*)''
-      WRITE(numsed,*)'     Solid Components '
-      WRITE(numsed,*)'     ====================='
-      WRITE(numsed,*)'nmole of Opale rained per dt'
-      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsopal))*dtsed,MAXVAL(rainrm(1:jpoce,jsopal))*dtsed,&
-         &SUM(rainrm(1:jpoce,1))*dtsed/jpoce
-      WRITE(numsed,*)'nmole of Clay rained per dt'
-      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jsclay))*dtsed,MAXVAL(rainrm(1:jpoce,jsclay))*dtsed,&
-         &SUM(rainrm(1:jpoce,jsclay))*dtsed/jpoce
-      WRITE(numsed,*)'nmole of POC rained per dt'
-      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jspoc))*dtsed,MAXVAL(rainrm(1:jpoce,jspoc))*dtsed,&
-         &SUM(rainrm(1:jpoce,jspoc))*dtsed/jpoce
-      WRITE(numsed,*)'nmole of CaCO3 rained per dt'
-      WRITE(numsed,'(3(1PE9.3,2X))') MINVAL(rainrm(1:jpoce,jscal))*dtsed,MAXVAL(rainrm(1:jpoce,jscal))*dtsed,&
-         &SUM(rainrm(1:jpoce,jscal))*dtsed/jpoce
-      WRITE(numsed,*)' '
-      WRITE(numsed,*)'Weight frac of opal rained (%) '
-      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&
-         &MAXVAL(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100.,&
-         & SUM(rainrg(1:jpoce,jsopal)/raintg(1:jpoce))*100./jpoce
-      WRITE(numsed,*)'Weight frac of clay  rained (%) '
-      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&
-         &MAXVAL(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100.,&
-         &SUM(rainrg(1:jpoce,jsclay)/raintg(1:jpoce))*100./jpoce
-      WRITE(numsed,*)'Weight frac of POC rained (%)'
-      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&
-         &MAXVAL(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100.,&
-         &SUM(rainrg(1:jpoce,jspoc)/raintg(1:jpoce))*100./jpoce
-      WRITE(numsed,*)'Weight frac of CaCO3  rained (%)'
-      WRITE(numsed,'(3(F0.3,7X))') MINVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&
-         &MAXVAL(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100.,&
-         &SUM(rainrg(1:jpoce,jscal)/raintg(1:jpoce))*100./jpoce
-      WRITE(numsed,*)''
-      WRITE(numsed,*)' Thickness of sediment layer added by particule rain, dzdep cm '
-      WRITE(numsed,*)' =============================================================='
-      WRITE(numsed,'(3(F0.5,2X))') MINVAL(dzdep(1:jpoce)),MAXVAL(dzdep(1:jpoce)),SUM(dzdep(1:jpoce))/jpoce
-      WRITE(numsed,*)''
-      WRITE(numsed,*)' chemical constants K1,pK1,K2,pK2,Kw,pKw and Kb pKb (min max) [mol/kgsol] or [mol/kgsol]2 '
-      WRITE(numsed,*)' ========================================================================================='
-      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak1s(1:jpoce)),MAXVAL(ak1s(1:jpoce)),-LOG10(MINVAL(ak1s(1:jpoce))),&
-         &-LOG10(MAXVAL(ak1s(1:jpoce)))
-      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(ak2s(1:jpoce)),MAXVAL(ak2s(1:jpoce)),-LOG10(MINVAL(ak2s(1:jpoce))),&
-         &-LOG10(MAXVAL(ak2s(1:jpoce)))
-      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akws(1:jpoce)),MAXVAL(akws(1:jpoce)),-LOG10(MINVAL(akws(1:jpoce))),&
-         &-LOG10(MAXVAL(akws(1:jpoce)))
-      WRITE(numsed,'(4(1PE10.3,2X))')MINVAL(akbs(1:jpoce)),MAXVAL(akbs(1:jpoce)),-LOG10(MINVAL(akbs(1:jpoce))),&
-         &-LOG10(MAXVAL(akbs(1:jpoce)))
-      WRITE(numsed,*)'and boron'
-      WRITE(numsed,'(2(1PE10.3,2X))')MINVAL(borats(1:jpoce)),MAXVAL(borats(1:jpoce))
-      WRITE(numsed,*)''
-      WRITE(numsed,*)' Compo of initial sediment for point jpoce'
-      WRITE(numsed,*)' ========================================='
-      WRITE(numsed,*)'solcp(1),   solcp(2),   solcp(3),   solcp(4),   hipor,      pH,         co3por'
-      DO jk = 1,jpksed
-         WRITE(numsed,'(7(1PE10.3,2X))')solcp(jpoce,jk,jsopal),solcp(jpoce,jk,jsclay),solcp(jpoce,jk,jspoc),solcp(jpoce,jk,jscal),&
-            &       hipor(jpoce,jk),-LOG10(hipor(jpoce,jk)/densSW(jpoce)),co3por(jpoce,jk)
-      ENDDO
-      WRITE(numsed,'(A82)')'pwcp(1),    pwcp(2),    pwcp(3),    pwcp(4),    pwcp(5),    pwcp(6),    pwcp(7)'
-      DO jk = 1, jpksed
-         WRITE(numsed,'(7(1PE10.3,2X))')pwcp(jpoce,jk,jwsil),pwcp(jpoce,jk,jwoxy),pwcp(jpoce,jk,jwdic),&
-            & pwcp(jpoce,jk,jwno3),pwcp(jpoce,jk,jwpo4),pwcp(jpoce,jk,jwalk),pwcp(jpoce,jk,jwc13)
-      ENDDO
-      WRITE(numsed,*) ' '
-      WRITE(numsed,*) ' End Of Initialization '
-      WRITE(numsed,*) ' '
-!
-   END SUBROUTINE sed_init_wri
-#else
-   !!----------------------------------------------------------------------
-   !!   Dummy module :                      NO Sediment model
-   !!----------------------------------------------------------------------
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_ini              ! Empty routine
-   END SUBROUTINE sed_ini
-#endif
-
-
 END MODULE sedini
