Changeset 10222

Timestamp:

2018-10-25T11:42:23+02:00 (6 years ago)

Author:

aumont

Message:

update of the sediment module + enhancement, bug correction in PISCES generating O2 negative values

Location:

NEMO/trunk/src

Files:

• OCE/IOM/iom.F90 (modified) (3 diffs)
• OFF/dtadyn.F90 (modified) (3 diffs)
• OFF/nemogcm.F90 (modified) (2 diffs)
• TOP/PISCES/P4Z/p4zche.F90 (modified) (3 diffs)
• TOP/PISCES/P4Z/p4zmeso.F90 (modified) (1 diff)
• TOP/PISCES/P4Z/p4zrem.F90 (modified) (2 diffs)
• TOP/PISCES/P4Z/p4zsms.F90 (modified) (3 diffs)
• TOP/PISCES/SED/par_sed.F90 (modified) (4 diffs)
• TOP/PISCES/SED/sed.F90 (modified) (7 diffs)
• TOP/PISCES/SED/sedadv.F90 (modified) (14 diffs)
• TOP/PISCES/SED/sedarr.F90 (modified) (11 diffs)
• TOP/PISCES/SED/sedbtb.F90 (modified) (5 diffs)
• TOP/PISCES/SED/sedchem.F90 (modified) (6 diffs)
• TOP/PISCES/SED/sedco3.F90 (modified) (4 diffs)
• TOP/PISCES/SED/seddsr.F90 (modified) (7 diffs)
• TOP/PISCES/SED/seddta.F90 (modified) (5 diffs)
• TOP/PISCES/SED/sedini.F90 (modified) (22 diffs)
• TOP/PISCES/SED/sedmat.F90 (modified) (11 diffs)
• TOP/PISCES/SED/sedmbc.F90 (modified) (10 diffs)
• TOP/PISCES/SED/sedmodel.F90 (modified) (3 diffs)
• TOP/PISCES/SED/sedrst.F90 (modified) (7 diffs)
• TOP/PISCES/SED/sedsfc.F90 (modified) (4 diffs)
• TOP/PISCES/SED/sedstp.F90 (modified) (4 diffs)
• TOP/PISCES/SED/sedwri.F90 (modified) (9 diffs)
• TOP/PISCES/sms_pisces.F90 (modified) (2 diffs)
• TOP/PISCES/trcini_pisces.F90 (modified) (2 diffs)
• TOP/PISCES/trcnam_pisces.F90 (modified) (2 diffs)
• TOP/trc.F90 (modified) (1 diff)
• TOP/trcdta.F90 (modified) (1 diff)
• TOP/trcini.F90 (modified) (1 diff)
• TOP/trcrst.F90 (modified) (1 diff)

NEMO/trunk/src/OCE/IOM/iom.F90

@@ -43 +43 @@
    USE ioipsl, ONLY :  ju2ymds    ! for calendar
    USE crs             ! Grid coarsening
+#if defined key_top
+   USE trc, ONLY    :  profsed
+#endif
    USE lib_fortran 
    USE diurnal_bulk, ONLY : ln_diurnal_only, ln_diurnal
@@ -193 +196 @@
       ! vertical grid definition
       IF(.NOT.llrst_context) THEN
-          CALL iom_set_axis_attr( "deptht", paxis = gdept_1d )
-          CALL iom_set_axis_attr( "depthu", paxis = gdept_1d )
-          CALL iom_set_axis_attr( "depthv", paxis = gdept_1d )
-          CALL iom_set_axis_attr( "depthw", paxis = gdepw_1d )
+          CALL iom_set_axis_attr( "deptht",  paxis = gdept_1d )
+          CALL iom_set_axis_attr( "depthu",  paxis = gdept_1d )
+          CALL iom_set_axis_attr( "depthv",  paxis = gdept_1d )
+          CALL iom_set_axis_attr( "depthw",  paxis = gdepw_1d )
 
           ! Add vertical grid bounds
@@ -219 +222 @@
           CALL iom_set_axis_attr( "nstrait", (/ (REAL(ji,wp), ji=1,4) /) )
 # endif
+#if defined key_top
+          CALL iom_set_axis_attr( "profsed", paxis = profsed )
+#endif
           CALL iom_set_axis_attr( "icbcla", class_num )
           CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )

NEMO/trunk/src/OFF/dtadyn.F90

@@ -46 +46 @@
    PRIVATE
 
-   PUBLIC   dta_dyn_init   ! called by opa.F90
-   PUBLIC   dta_dyn        ! called by step.F90
-   PUBLIC   dta_dyn_swp   ! called by step.F90
+   PUBLIC   dta_dyn_init       ! called by opa.F90
+   PUBLIC   dta_dyn            ! called by step.F90
+   PUBLIC   dta_dyn_sed_init   ! called by opa.F90
+   PUBLIC   dta_dyn_sed        ! called by step.F90
+   PUBLIC   dta_dyn_swp        ! called by step.F90
 
    CHARACTER(len=100) ::   cn_dir          !: Root directory for location of ssr files
@@ -183 +185 @@
          CALL prt_ctl(tab3d_1=uslp             , clinfo1=' slp  - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
          CALL prt_ctl(tab3d_1=wslpi            , clinfo1=' slp  - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
-!         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i    - : ', mask1=tmask )
-!         CALL prt_ctl(tab2d_1=hmld             , clinfo1=' hmld    - : ', mask1=tmask )
-!         CALL prt_ctl(tab2d_1=fmmflx           , clinfo1=' fmmflx  - : ', mask1=tmask )
-!         CALL prt_ctl(tab2d_1=emp              , clinfo1=' emp     - : ', mask1=tmask )
-!         CALL prt_ctl(tab2d_1=wndm             , clinfo1=' wspd    - : ', mask1=tmask )
-!         CALL prt_ctl(tab2d_1=qsr              , clinfo1=' qsr     - : ', mask1=tmask )
       ENDIF
       !
@@ -419 +415 @@
    END SUBROUTINE dta_dyn_init
 
+   SUBROUTINE dta_dyn_sed( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dta_dyn  ***
+      !!
+      !! ** Purpose :  Prepares dynamics and physics fields from a NEMO run
+      !!               for an off-line simulation of passive tracers
+      !!
+      !! ** Method : calculates the position of data
+      !!             - computes slopes if needed
+      !!             - interpolates data if needed
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start( 'dta_dyn_sed')
+      !
+      nsecdyn = nsec_year + nsec1jan000   ! number of seconds between Jan. 1st 00h of nit000 year and the middle of time step
+      !
+      IF( kt == nit000 ) THEN    ;    nprevrec = 0
+      ELSE                       ;    nprevrec = sf_dyn(jf_tem)%nrec_a(2)
+      ENDIF
+      CALL fld_read( kt, 1, sf_dyn )      !=  read data at kt time step   ==!
+      !
+      tsn(:,:,:,jp_tem) = sf_dyn(jf_tem)%fnow(:,:,:)  * tmask(:,:,:)    ! temperature
+      tsn(:,:,:,jp_sal) = sf_dyn(jf_sal)%fnow(:,:,:)  * tmask(:,:,:)    ! salinity
+      !
+      CALL eos    ( tsn, rhd, rhop, gdept_0(:,:,:) ) ! In any case, we need rhop
+
+      IF(ln_ctl) THEN                  ! print control
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_tem), clinfo1=' tn      - : ', mask1=tmask,  kdim=jpk   )
+         CALL prt_ctl(tab3d_1=tsn(:,:,:,jp_sal), clinfo1=' sn      - : ', mask1=tmask,  kdim=jpk   )
+      ENDIF
+      !
+      IF( ln_timing )   CALL timing_stop( 'dta_dyn_sed')
+      !
+   END SUBROUTINE dta_dyn_sed
+
+
+   SUBROUTINE dta_dyn_sed_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dta_dyn_init  ***
+      !!
+      !! ** Purpose :   Initialisation of the dynamical data
+      !! ** Method  : - read the data namdta_dyn namelist
+      !!----------------------------------------------------------------------
+      INTEGER  :: ierr, ierr0, ierr1, ierr2, ierr3   ! return error code
+      INTEGER  :: ifpr                               ! dummy loop indice
+      INTEGER  :: jfld                               ! dummy loop arguments
+      INTEGER  :: inum, idv, idimv                   ! local integer
+      INTEGER  :: ios                                ! Local integer output status for namelist read
+      !!
+      CHARACTER(len=100)            ::  cn_dir        !   Root directory for location of core files
+      TYPE(FLD_N), DIMENSION(2) ::  slf_d         ! array of namelist informations on the fields to read
+      TYPE(FLD_N) :: sn_tem , sn_sal   !   "                 "
+      !!
+      NAMELIST/namdta_dyn/cn_dir, ln_dynrnf, ln_dynrnf_depth,  fwbcorr, &
+         &                sn_tem, sn_sal
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )              ! Namelist namdta_dyn in reference namelist : Offline: init. of dynamical data
+      READ  ( numnam_ref, namdta_dyn, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namdta_dyn in reference namelist', lwp )
+      REWIND( numnam_cfg )              ! Namelist namdta_dyn in configuration namelist : Offline: init. of dynamical data
+      READ  ( numnam_cfg, namdta_dyn, IOSTAT = ios, ERR = 902 )
+902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namdta_dyn in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namdta_dyn )
+      !                                         ! store namelist information in an array
+      !                                         ! Control print
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'dta_dyn : offline dynamics '
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namdta_dyn'
+         WRITE(numout,*) '      runoffs option enabled (T) or not (F)            ln_dynrnf        = ', ln_dynrnf
+         WRITE(numout,*) '      runoffs is spread in vertical                    ln_dynrnf_depth  = ', ln_dynrnf_depth
+         WRITE(numout,*) '      annual global mean of empmr for ssh correction   fwbcorr          = ', fwbcorr
+         WRITE(numout,*)
+      ENDIF
+      !
+      jf_tem  = 1     ;   jf_sal  = 2    ;   jfld   = jf_sal
+      !
+      slf_d(jf_tem)  = sn_tem    ;   slf_d(jf_sal)  = sn_sal
+      !
+      ALLOCATE( sf_dyn(jfld), STAT=ierr )         ! set sf structure
+      IF( ierr > 0 )  THEN
+         CALL ctl_stop( 'dta_dyn: unable to allocate sf structure' )   ;   RETURN
+      ENDIF
+      !                                         ! fill sf with slf_i and control print
+      CALL fld_fill( sf_dyn, slf_d, cn_dir, 'dta_dyn_init', 'Data in file', 'namdta_dyn' )
+      !
+      ! Open file for each variable to get his number of dimension
+      DO ifpr = 1, jfld
+         CALL fld_clopn( sf_dyn(ifpr), nyear, nmonth, nday )
+         idv   = iom_varid( sf_dyn(ifpr)%num , slf_d(ifpr)%clvar )        ! id of the variable sdjf%clvar
+         idimv = iom_file ( sf_dyn(ifpr)%num )%ndims(idv)                 ! number of dimension for variable sdjf%clvar
+         IF( sf_dyn(ifpr)%num /= 0 )   CALL iom_close( sf_dyn(ifpr)%num ) ! close file if already open
+         ierr1=0
+         IF( idimv == 3 ) THEN    ! 2D variable
+                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,1)    , STAT=ierr0 )
+            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,1,2)  , STAT=ierr1 )
+         ELSE                     ! 3D variable
+                                      ALLOCATE( sf_dyn(ifpr)%fnow(jpi,jpj,jpk)  , STAT=ierr0 )
+            IF( slf_d(ifpr)%ln_tint ) ALLOCATE( sf_dyn(ifpr)%fdta(jpi,jpj,jpk,2), STAT=ierr1 )
+         ENDIF
+         IF( ierr0 + ierr1 > 0 ) THEN
+            CALL ctl_stop( 'dta_dyn_init : unable to allocate sf_dyn array structure' )   ;   RETURN
+         ENDIF
+      END DO
+      !
+      CALL dta_dyn_sed( nit000 )
+      !
+   END SUBROUTINE dta_dyn_sed_init
 
    SUBROUTINE dta_dyn_swp( kt )

NEMO/trunk/src/OFF/nemogcm.F90

@@ -110 +110 @@
          IF( istp /= nit000 )   CALL day        ( istp )         ! Calendar (day was already called at nit000 in day_init)
                                 CALL iom_setkt  ( istp - nit000 + 1, cxios_context )   ! say to iom that we are at time step kstp
+#if defined key_sed_off
+                                CALL dta_dyn_sed( istp )         ! Interpolation of the dynamical fields
+#else
                                 CALL dta_dyn    ( istp )         ! Interpolation of the dynamical fields
          IF( .NOT.ln_linssh )   CALL dta_dyn_swp( istp )         ! swap of sea  surface height and vertical scale factors
-
+#endif
                                 CALL trc_stp    ( istp )         ! time-stepping
                                 CALL stp_ctl    ( istp, indic )  ! Time loop: control and print
@@ -287 +290 @@
                            CALL trc_nam_run    ! Needed to get restart parameters for passive tracers
                            CALL trc_rst_cal( nit000, 'READ' )   ! calendar
+#if defined key_sed_off
+                           CALL dta_dyn_sed_init ! Initialization for the dynamics
+#else
                            CALL dta_dyn_init   ! Initialization for the dynamics
+#endif
 
                            CALL     trc_init   ! Passive tracers initialization

NEMO/trunk/src/TOP/PISCES/P4Z/p4zche.F90

@@ -37 +37 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akf3       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aks3       !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak1p3      !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak2p3      !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak3p3      !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksi3      !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fluorid    !: ???
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sulfat     !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akw3       !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akf3       !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aks3       !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak1p3      !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak2p3      !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ak3p3      !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   aksi3      !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   borat      !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   fluorid    !: ???
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sulfat     !: ???
 
    !!* Variable for chemistry of the CO2 cycle
@@ -233 +233 @@
         END DO
 
-
-
       ! CHEMICAL CONSTANTS - DEEP OCEAN
       ! -------------------------------
@@ -449 +447 @@
       !
    END SUBROUTINE p4z_che
-
 
    SUBROUTINE ahini_for_at(p_hini)

NEMO/trunk/src/TOP/PISCES/P4Z/p4zmeso.F90

@@ -108 +108 @@
                zdenom    = zfoodlim / ( xkgraz2 + zfoodlim )
                zdenom2   = zdenom / ( zfood + rtrn )
-               zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) 
+               zgraze2   = grazrat2 * xstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) * (1. - nitrfac(ji,jj,jk)) 
 
                zgrazd    = zgraze2  * xprefc   * zcompadi  * zdenom2 

NEMO/trunk/src/TOP/PISCES/P4Z/p4zrem.F90

@@ -116 +116 @@
                   zammonic = zremik * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)
                   denitr(ji,jj,jk)  = zammonic * ( 1. - nitrfac2(ji,jj,jk) )
-                  zoxyrem           = zammonic *        nitrfac2(ji,jj,jk)
+                  denitr(ji,jj,jk)  = MIN( ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit, denitr(ji,jj,jk) )
+                  zoxyrem           = zammonic - denitr(ji,jj,jk)
                   !
                   zolimi (ji,jj,jk) = MAX( 0.e0, zolimi (ji,jj,jk) )
@@ -189 +190 @@
                &         / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) ) 
                zdenitnh4 = nitrif * xstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk)
+               zdenitnh4 = MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenita, zdenitnh4 ) 
                ! Update of the tracers trends
                ! ----------------------------

NEMO/trunk/src/TOP/PISCES/P4Z/p4zsms.F90

@@ -101 +101 @@
       ENDIF
       !
+      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients 
+      !
+#if ! defined key_sed_off
       CALL p4z_che              ! computation of chemical constants
       CALL p4z_int( kt )        ! computation of various rates for biogeochemistry
       !
-      IF( ll_sbc ) CALL p4z_sbc( kt )   ! external sources of nutrients 
-
       DO jnt = 1, nrdttrc          ! Potential time splitting if requested
          !
@@ -149 +150 @@
          END DO
       END IF
-      !
-      IF( lk_sed ) THEN 
+#endif
+      !
+      IF( ln_sediment ) THEN 
          !
          CALL sed_model( kt )     !  Main program of Sediment model
+         !
+         IF( ln_top_euler ) THEN
+            DO jn = jp_pcs0, jp_pcs1
+               trn(:,:,:,jn) = trb(:,:,:,jn)
+            END DO
+         ENDIF
          !
       ENDIF
@@ -352 +360 @@
       IF(lwp)  WRITE(numout,*)
 
-      IF( cn_cfg == "ORCA" .OR. cn_cfg == "orca" ) THEN
+      IF( cn_cfg == "ORCA" .OR. cn_cfg == "orca") THEN
          IF( .NOT. lk_c1d ) THEN      ! ORCA configuration (not 1D) !
             !                                                ! --------------------------- !

NEMO/trunk/src/TOP/PISCES/SED/par_sed.F90

@@ -7 +7 @@
    !!        !  06-12  (C. Ethe)  Orignal
    !!----------------------------------------------------------------------
-#if defined key_sed
-   !!----------------------------------------------------------------------
-   !!   'key_sed'                                                  Sediment 
-   !!----------------------------------------------------------------------
+   !! $Id$
 
    !! Domain characteristics
@@ -19 +16 @@
       jpim1    =>   jpim1 ,  & !: jpi - 1
       jpjm1    =>   jpjm1 ,  & !: jpj - 1
-      jpij     =>   jpij       !: jpi x jpj
-      jp_tem   =>   jp_tem     !: indice of temperature
+      jpij     =>   jpij  ,  & !: jpi x jpj
+      jp_tem   =>   jp_tem,  & !: indice of temperature
       jp_sal   =>   jp_sal     !: indice of salintity
 
-#if ! defined key_sed_off
-   USE par_pisces
-#endif
-
-   INTEGER, PARAMETER :: jpdta = 12
-
+   INTEGER, PARAMETER :: jpdta = 17
 
    ! Vertical sediment geometry
-   INTEGER, PARAMETER :: &
-      jpksed   = 11         ,   & !: number of sediment layers
-      jpksedm1 = jpksed - 1 
-      
+   INTEGER, PUBLIC   ::      &
+      jpksed   = 11 ,        &
+      jpksedm1 = 10
+
    ! sediment tracer species
-   INTEGER, PARAMETER :: &
-      jpsol = 4,         &  !: number of solid component
-      jpwat = 7,         &   !: number of pore water component
-      jpwatp1 = jpwat + 1
+   INTEGER, PARAMETER ::    &
+      jpsol =  8,           &  !: number of solid component
+      jpwat = 10,           &   !: number of pore water component
+      jpwatp1 = jpwat +1,   &
+      jpsol1 = jpsol - 1
+
    
    ! pore water components       
@@ -49 +43 @@
       jwpo4  = 5,        & !: phosphate
       jwalk  = 6,        & !: alkalinity
-      jwc13  = 7           !: dissolved inorganic carbon 13
+      jwnh4  = 7,        & !: Ammonium
+      jwh2s  = 8,        & !: Sulfate
+      jwso4  = 9,        & !: H2S
+      jwfe2  = 10          !: Fe2+
 
    ! solid components       
@@ -56 +53 @@
       jsclay  = 2,        & !: clay
       jspoc   = 3,        & !: organic carbon
-      jscal   = 4           !: calcite
+      jscal   = 4,        & !: calcite
+      jspos   = 5,        & !: semi-ref POC
+      jspor   = 6,        & !: refractory POC
+      jsfeo   = 7,        & !: iron hydroxides
+      jsfes   = 8           !: FeS
 
    INTEGER, PARAMETER ::  &
       jptrased   = jpsol + jpwat , &
-      jpdia3dsed = 3          , &
-      jpdia2dsed = 7
-#endif
+      jpdia3dsed = 2             , &
+      jpdia2dsed = 12
 
-   !!----------------------------------------------------------------------
-   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
-   !! $Id$ 
-   !! Software governed by the CeCILL license (see ./LICENSE)
-   !!======================================================================
 END MODULE par_sed

NEMO/trunk/src/TOP/PISCES/SED/sed.F90

@@ -7 +7 @@
    !!        !  06-12  (C. Ethe)  Orignal
    !!----------------------------------------------------------------------
-#if defined key_sed
-   !!----------------------------------------------------------------------
-   !!   'key_sed'                                                  Sediment 
-   !!----------------------------------------------------------------------
    USE par_sed
+   USE oce_sed
    USE in_out_manager
+
 
    IMPLICIT NONE
    PUBLIC
 
-   PUBLIC sed_alloc   
-
-   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_1d    =>   e3t_1d         !: 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
-   !                                !: it is the accumulated duration of previous runs
-   !                                !: that may have been run with different time steps.
-
-#if ! defined key_sed_off
-
-   USE oce     , ONLY :  tsn        =>   tsn             !: pot. temperature (celsius) and salinity (psu)
-
-   USE trc     , ONLY :  trn        =>   trc             !: tracer 
-   USE trc     , ONLY :  nwritetrc  =>   nwritetrc       !: outputs frequency of tracer model
-
-   USE p4zsink , ONLY :  sinking    =>   sinking         !: sinking flux for POC
-   USE p4zsink , ONLY :  sinking2   =>   sinking2        !: sinking flux for GOC
-   USE p4zsink , ONLY :  sinkcal    =>   sinkcal         !: sinking flux for calcite
-   USE p4zsink , ONLY :  sinksil    =>   sinksil         !: sinking flux for opal ( dsi )
-
-   USE sms_pisces, ONLY : akb3      =>   akb3            !: Chemical constants  
-   USE sms_pisces, ONLY : ak13      =>   ak13            !: Chemical constants  
-   USE sms_pisces, ONLY : ak23      =>   ak23            !: Chemical constants  
-   USE sms_pisces, ONLY : akw3      =>   akw3            !: Chemical constants  
-   USE sms_pisces, ONLY : aksp      =>   aksp            !: Chemical constants  
-   USE sms_pisces, ONLY : borat     =>   borat           !: Chemical constants ( borat ) 
-
-#endif   
-
+   PUBLIC sed_alloc
 
    !! Namelist
-   REAL(wp), PUBLIC, DIMENSION(5) ::  reac                !: reactivity rc in  [l.mol-1.s-1]
    REAL(wp), PUBLIC               ::  reac_sil            !: reactivity of silicate in  [l.mol-1.s-1]
    REAL(wp), PUBLIC               ::  reac_clay           !: reactivity of clay in  [l.mol-1.s-1]
-   REAL(wp), PUBLIC               ::  reac_poc            !: reactivity of poc in  [l.mol-1.s-1]
-   REAL(wp), PUBLIC               ::  reac_no3            !: reactivity of no3 in  [l.mol-1.s-1]
+   REAL(wp), PUBLIC               ::  reac_ligc           !: reactivity of Ligands [l.mol-1.s-1]
+   REAL(wp), PUBLIC               ::  reac_pocl           !: reactivity of pocl in  [s-1]
+   REAL(wp), PUBLIC               ::  reac_pocs           !: reactivity of pocs in  [s-1]
+   REAL(wp), PUBLIC               ::  reac_pocr           !: reactivity of pocr in  [s-1]
+   REAL(wp), PUBLIC               ::  reac_nh4            !: reactivity of NH4 in  [l.mol-1.s-1]
+   REAL(wp), PUBLIC               ::  reac_h2s            !: reactivity of ODU in  [l.mol-1.s-1]
+   REAL(wp), PUBLIC               ::  reac_fe2            !: reactivity of Fe2+ in  [l.mol-1.s-1]
+   REAL(wp), PUBLIC               ::  reac_feh2s          !: reactivity of Fe2+ in  [l.mol-1.s-1]
+   REAL(wp), PUBLIC               ::  reac_fes            !: reactivity of Fe with H2S in  [l.mol-1.s-1]
+   REAL(wp), PUBLIC               ::  reac_feso           !: reactivity of FeS with O2 in  [l.mol-1.s-1]
    REAL(wp), PUBLIC               ::  reac_cal            !: reactivity of cal in  [l.mol-1.s-1]
-   REAL(wp), PUBLIC               ::  sat_sil             !: saturation concentration for silicate in [mol.l-1]
-   REAL(wp), PUBLIC               ::  sat_clay            !: saturation concentration for clay in [mol.l-1]
+   REAL(wp), PUBLIC               ::  adsnh4              !: adsorption coefficient of NH4
+   REAL(wp), PUBLIC               ::  ratligc             !: C/L ratio in POC
    REAL(wp), PUBLIC               ::  so2ut 
    REAL(wp), PUBLIC               ::  srno3 
    REAL(wp), PUBLIC               ::  spo4r 
    REAL(wp), PUBLIC               ::  srDnit 
-   REAL(wp), PUBLIC               ::  sthro2              !: threshold O2 concen. in [mol.l-1]
-   REAL(wp), PUBLIC               ::  pdb = 0.0112372     !: 13C/12C in PD Belemnite
-   REAL(wp), PUBLIC               ::  rc13P  = 0.980      !: 13C/12C in POC = rc13P*PDB
-   REAL(wp), PUBLIC               ::  rc13Ca = 1.001      !: 13C/12C in CaCO3 = rc13Ca*PDB
    REAL(wp), PUBLIC               ::  dtsed               !: sedimentation time step
-   REAL(wp), PUBLIC               ::  db                  !: bioturb coefficient in [cm2.s-1]
-
+   REAL(wp), PUBLIC               ::  dtsed2              !: sedimentation time step
    INTEGER , PUBLIC               ::  nitsed000
    INTEGER , PUBLIC               ::  nitsedend
-   INTEGER , PUBLIC               ::  nwrised
-   INTEGER , PUBLIC               ::  nfreq
-   REAL(wp), PUBLIC               ::  dens                !: density of solid material
+   INTEGER, PUBLIC                ::  nrseddt
+   REAL    , PUBLIC               ::  sedmask
+   REAL(wp), PUBLIC               ::  denssol                !: density of solid material
+   INTEGER , PUBLIC               ::  numrsr, numrsw   !: logical unit for sed restart (read and write)
+   LOGICAL , PUBLIC               ::  lrst_sed       !: logical to control the trc restart write
+   LOGICAL , PUBLIC               ::  ln_rst_sed  = .TRUE.     !: initialisation from a restart file or not
+   LOGICAL , PUBLIC               ::  ln_btbz     = .FALSE.    !: Depth variation of the bioturbation coefficient
+   LOGICAL , PUBLIC               ::  ln_irrig    = .FALSE.    !: iActivation of the bioirrigation
+   LOGICAL , PUBLIC               ::  ln_sed_2way = .FALSE.    !: 2 way coupling with PISCES
+   LOGICAL , PUBLIC               ::  ln_sediment_offline = .FALSE. !: Offline mode for sediment module
+   INTEGER             , PUBLIC   ::  nn_rstsed      !: control of the time step ( 0 or 1 ) for pass. tr.
+   INTEGER , PUBLIC               ::  nn_dtsed = 1   !: frequency of step on passive tracers
+   CHARACTER(len = 80) , PUBLIC   ::  cn_sedrst_in   !: suffix of pass. tracer restart name (input)
+   CHARACTER(len = 256), PUBLIC   ::  cn_sedrst_indir  !: restart input directory
+   CHARACTER(len = 80) , PUBLIC   ::  cn_sedrst_out  !: suffix of pass. tracer restart name (output)
+   CHARACTER(len = 256), PUBLIC   ::  cn_sedrst_outdir  !: restart output directory
+
    !
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  pwcp       !: pore water sediment data at given time-step
@@ -87 +63 @@
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  solcp      !: solid sediment data at given time-step
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  solcp0     !: solid sediment data at initial time
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  trc_dta
+   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  diff
 
    !! * Shared module variables
@@ -102 +80 @@
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  press      !: pressure
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  raintg     !: total massic flux rained in each cell (sum of sol. comp.)
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  fecratio   !: Fe/C ratio in falling particles to the sediments
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  dzdep      !: total thickness of solid material rained [cm] in each cell
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  zkbot      !: total thickness of solid material rained [cm] in each cell
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  wacc       !: total thickness of solid material rained [cm] in each cell
    !
    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  hipor      !: [h+] in mol/kg*densSW 
@@ -127 +108 @@
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksis 
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aksps 
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  sieqs
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  aks3s
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  akf3s
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  sulfats
+   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  fluorids
 
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  mol_wgt    !: molecular weight of solid sediment data
@@ -133 +119 @@
    !! Geometry
    INTEGER , PUBLIC, SAVE                          ::  jpoce, indoce !: Ocean points ( number/indices )
-   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  iarroce       !: Computation of 1D array of sediments points
+   INTEGER , PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  iarroce       !: Computation of 1D array of sediments points
    REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  epkbot        !: ocean bottom layer thickness
+   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  gdepbot       !: Depth of the sediment
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  dzkbot        !: ocean bottom layer thickness in meters
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  tmasksed      !: sediment mask
-   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  sbathy        !: bathymetry
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  dz            !: sediment layers thickness
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  por           !: porosity profile     
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  por1          !: 1-por 
-   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  profsed       !: depth of middle of each layer
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  volw          !: volume of pore water cell fraction
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  vols          !: volume of solid cell fraction
-   REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  diff          !: diffusion ceofficient
+   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  db            !: bioturbation ceofficient
+   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE ::  irrig        !: bioturbation ceofficient
    REAL(wp), PUBLIC, DIMENSION(:    ), ALLOCATABLE ::  rdtsed        !:  sediment model time-step
+   REAL(wp), PUBLIC, DIMENSION(:,:  ), ALLOCATABLE :: sedligand
    REAL(wp)  ::   dens               !: density of solid material
    !! Inputs / Outputs
@@ -171 +157 @@
       !!-------------------------------------------------------------------
       !
-      ALLOCATE( trc_dta(jpi,jpj,jdta)                                     ,   &
-         &      epkbot(jpi,jpj), sbathy(jpi,jpj)                          ,   &
-         &      tmasksed(jpi,jpj,jpksed)                                  ,   &
-         &      dz(jpksed)  , por(jpksed) , por1(jpksed), profsed(jpksed) ,   &
-         &      volw(jpksed), vols(jpksed), diff(jpksed), rdtsed(jpksed)  ,   &
+      ALLOCATE( trc_data(jpi,jpj,jpdta)                                   ,   &
+         &      epkbot(jpi,jpj), gdepbot(jpi,jpj)        ,   &
+         &      dz(jpksed)  , por(jpksed) , por1(jpksed)                  ,   &
+         &      volw(jpksed), vols(jpksed), rdtsed(jpksed)  ,   &
          &      trcsedi  (jpi,jpj,jpksed,jptrased)                        ,   &
          &      flxsedi3d(jpi,jpj,jpksed,jpdia3dsed)                      ,   &
-         &      flxsedi2d(jpi,jpj,jpksed,jpdia2dsed)                      ,   &
+         &      flxsedi2d(jpi,jpj,jpdia2dsed)                             ,   &
          &      mol_wgt(jpsol),                                           STAT=sed_alloc )
 
@@ -185 +170 @@
    END FUNCTION sed_alloc
 
-#else
-   !!======================================================================
-   !! No Sediment model
-   !!======================================================================
-#endif
-
 END MODULE sed

NEMO/trunk/src/TOP/PISCES/SED/sedadv.F90

@@ -4 +4 @@
    !!    Sediment : vertical advection and burial
    !!=====================================================================
-#if defined key_sed
-   !!----------------------------------------------------------------------
-   !!   'key_sed'                                                  Sediment 
+   !! * Modules used
    !!----------------------------------------------------------------------
    !!   sed_adv :
    !!----------------------------------------------------------------------
    USE sed     ! sediment global variable
+   USE lib_mpp         ! distribued memory computing library
+
+   IMPLICIT NONE
+   PRIVATE
 
    PUBLIC sed_adv
-
-   !! * Module variable
-   INTEGER, PARAMETER  ::  nztime = jpksed         ! number of time step between sunrise and sunset
-
-   REAL(wp), DIMENSION(jpksed), SAVE :: dvolsp, dvolsm
-   REAL(wp), DIMENSION(jpksed), SAVE :: c2por, ckpor
+   PUBLIC sed_adv_alloc
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: dvolsp, dvolsm
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: c2por, ckpor
 
    REAL(wp) :: cpor
@@ -49 +49 @@
       ! * local variables
       INTEGER :: ji, jk, js 
-      INTEGER :: jn, ntimes, ikwneg
+      INTEGER :: jn, ntimes, nztime, ikwneg
       
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsolcpno
-      REAL(wp), DIMENSION(:  ), ALLOCATABLE :: zfilled, zfull, zfromup, zempty
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zgap, zwb
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zrainrf
-      REAL(wp), DIMENSION(nztime) ::  zraipush
-
-      REAL(wp) :: zkwnup, zkwnlo, zfrac,  zfromce, zrest
+      REAL(wp), DIMENSION(jpksed,jpsol) :: zsolcpno
+      REAL(wp), DIMENSION(jpksed)       :: zfilled, zfull, zfromup, zempty
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zgap, zwb
+      REAL(wp), DIMENSION(jpoce,jpsol) :: zrainrf
+      REAL(wp), DIMENSION(:  ), ALLOCATABLE :: zraipush
+
+      REAL(wp) :: zkwnup, zkwnlo, zfrac,  zfromce, zrest, sumtot, zsumtot1
 
       !------------------------------------------------------------------------
 
 
+      IF( ln_timing )  CALL timing_start('sed_adv')
+!
       IF( kt == nitsed000 ) THEN
-         WRITE(numsed,*) ' '
-         WRITE(numsed,*) ' sed_adv : vertical sediment advection  '
-         WRITE(numsed,*) ' '
-         por1clay = dens * por1(jpksed) * dz(jpksed) / mol_wgt(jsclay)
+         IF (lwp) THEN
+            WRITE(numsed,*) ' '
+            WRITE(numsed,*) ' sed_adv : vertical sediment advection  '
+            WRITE(numsed,*) ' '
+         ENDIF
+         por1clay = denssol * por1(jpksed) * dz(jpksed)
          cpor     = por1(jpksed) / por1(2)
          DO jk = 2, jpksed
@@ -80 +84 @@
       ENDIF
 
-      ALLOCATE( zsolcpno(jpksed,jpsol), zrainrf(jpoce,jpsol) )
-      ALLOCATE( zfilled(jpksed), zfull(jpksed), zfromup(jpksed), zempty(jpksed) )
-      ALLOCATE( zgap (jpoce,jpksed)   , zwb(jpoce,jpksed)  )  
-
       ! Initialization of data for mass balance calculation
       !---------------------------------------------------
@@ -89 +89 @@
       tosed  (:,:) = 0. 
       rloss  (:,:) = 0.
-
+      ikwneg = 1
+      nztime = jpksed
+
+      ALLOCATE( zraipush(nztime) )
 
       ! Initiate gap 
@@ -104 +107 @@
       zgap(1:jpoce,1:jpksed) = 1. - zgap(1:jpoce,1:jpksed)   
 
-
       ! Initiate burial rates
       !-----------------------
       zwb(:,:) = 0.
       DO jk = 2, jpksed
-         zfrac =  dtsed / ( dens * por1(jk) )     
+         zfrac =  dtsed / ( denssol * por1(jk) )     
          DO ji = 1, jpoce
             zwb(ji,jk) = zfrac * raintg(ji)
@@ -127 +129 @@
       ENDDO
 
-
       zrainrf(:,:) = 0.
       DO ji = 1, jpoce
@@ -206 +207 @@
                   ! quantities to push in deeper sediment
                   tosed  (ji,js) = zsolcpno(jpksed,js) &
-                     &           * zwb(ji,jpksed) * dens * por1(jpksed) / mol_wgt(js)
-               ENDDO
-               
+                     &           * zwb(ji,jpksed) * denssol * por1(jpksed)
+               ENDDO
+
             ELSE ! what is remaining is great than dz(2)
 
                ntimes = INT( zrest / dz(2) ) + 1
-               IF( ntimes > nztime ) THEN
-                  WRITE( numsed,* ) ' sedadv : rest too large at sediment point ji = ', ji
-                  STOP
-               ENDIF
+               IF( ntimes > nztime ) CALL ctl_stop( 'STOP', 'sed_adv : rest too large ' )
                zraipush(1) = dz(2)
                zrest = zrest - zraipush(1)
@@ -249 +247 @@
                      fromsed(ji,js) = 0.
                      tosed  (ji,js) = tosed(ji,js) + zsolcpno(jpksed,js) * zraipush(jn) &
-                        &             * dens * por1(2) / mol_wgt(js)
+                        &             * denssol * por1(2) 
                   ENDDO
                ENDDO
@@ -279 +277 @@
             ! for the last layer, one make go up clay 
             solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
-            !! C. Heinze      fromsed(ji,jsclay) = zempty(jpksed) * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
             fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
-
          ELSE  ! rain > 0 and rain < total reaction loss
 
@@ -323 +319 @@
                ENDDO
                solcp(ji,jpksed,jsclay) =  solcp(ji,jpksed,jsclay) + zempty(jpksed) * 1.
-               !! C. Heinze  fromsed(ji,jsclay) = zempty(jpksed) * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
+               !! C. Heinze  fromsed(ji,jsclay) = zempty(jpksed) * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)
                fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
                
@@ -339 +335 @@
                   solcp(ji,2,js) = zfull(2) * zsolcpno(2,js) + zempty(2) * zrainrf(ji,js)
                ENDDO
-               DO  js = 1, jpsol              
+               DO  js = 1, jpsol
                   DO jk = jpksedm1, 3, -1
                      solcp(ji,jk,js) = zfull(jk) * zsolcpno(jk,js) + zempty(jk) * zsolcpno(jk-1,js)
@@ -349 +345 @@
                ENDDO
                solcp(ji,jpksed,jsclay) = solcp(ji,jpksed,jsclay) + zkwnlo * 1.
-               ! Heinze  fromsed(ji,jsclay) = zkwnlo * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
+               ! Heinze  fromsed(ji,jsclay) = zkwnlo * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)
                fromsed(ji,jsclay) = zkwnlo * 1.* por1clay
             ELSE   ! 2 < ikwneg(ji) <= jpksedm1
@@ -415 +411 @@
                   fromsed(ji,js)   = 0.
                ENDDO
-               ! Heinze  fromsed(ji,jsclay) = zempty * 1. * dens * por1(jpksed) / mol_wgt(jsclay)
+               ! Heinze  fromsed(ji,jsclay) = zempty * 1. * denssol * por1(jpksed) / mol_wgt(jsclay)
                fromsed(ji,jsclay) = zempty(jpksed) * 1. * por1clay
+
             ENDIF ! ikwneg(ji) = 2
          ENDIF  ! zwb > 0
@@ -425 +422 @@
       raintg(:)   = 0.
 
-      DEALLOCATE( zsolcpno )   
-      DEALLOCATE( zrainrf  )
-      DEALLOCATE( zfilled  )
-      DEALLOCATE( zfull    )
-      DEALLOCATE( zfromup  )
-      DEALLOCATE( zempty   )
-      DEALLOCATE( zgap     )  
-      DEALLOCATE( zwb      )
-
+      DEALLOCATE( zraipush )
+
+      IF( ln_timing )  CALL timing_stop('sed_adv')
 
    END SUBROUTINE sed_adv
-#else
-   !!======================================================================
-   !! MODULE sedbtb  :   Dummy module 
-   !!======================================================================
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_adv( kt )         ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'sed_adv: You should not have seen this print! error?', kt
-   END SUBROUTINE sed_adv
-
-   !!======================================================================
-
-#endif
+
+
+   INTEGER FUNCTION sed_adv_alloc()
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_prod_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( dvolsp(jpksed), dvolsm(jpksed), c2por(jpksed),         &
+      &         ckpor(jpksed) ,           STAT = sed_adv_alloc )
+      !
+      IF( sed_adv_alloc /= 0 ) CALL ctl_warn('sed_adv_alloc : failed to allocate arrays.')
+      !
+   END FUNCTION sed_adv_alloc
+
 END MODULE sedadv

NEMO/trunk/src/TOP/PISCES/SED/sedarr.F90

@@ -4 +4 @@
    !!              transform 1D (2D) array to a 2D (1D) table
    !!======================================================================
-#if defined key_sed
+
    !!----------------------------------------------------------------------
    !!   arr_2d_1d  : 2-D to 1-D
@@ -11 +11 @@
    !! * Modules used
    USE par_sed
+   USE dom_oce
+   USE sed
 
    IMPLICIT NONE
@@ -28 +30 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL license (see ./LICENSE)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -42 +44 @@
 
       INTEGER ::  jn, jid, jjd
+
+      IF( ln_timing )   CALL timing_start('pack_arr_2d_1d')
         
       DO jn = 1, ndim1d
@@ -48 +52 @@
          tab1d(jn)  = tab2d(jid, jjd)
       END DO 
+
+      IF( ln_timing )   CALL timing_stop('pack_arr_2d_1d')
 
    END SUBROUTINE pack_arr_2d_1d
@@ -59 +65 @@
       INTEGER ::  jn, jid, jjd
 
+      IF( ln_timing )   CALL timing_start('unpack_arr_1d_2d')
+
       DO jn = 1, ndim1d
          jid             = MOD( tab_ind(jn) - 1, jpi) + 1
@@ -64 +72 @@
          tab2d(jid, jjd) = tab1d(jn)
       END DO
+
+      IF( ln_timing )   CALL timing_stop('unpack_arr_1d_2d')
 
    END SUBROUTINE unpack_arr_1d_2d
@@ -75 +85 @@
       INTEGER, DIMENSION(ndim1d) ::  jid, jjd        
       INTEGER ::    jk, jn , ji, jj
+
+      IF( ln_timing )   CALL timing_start('pack_arr_2d_3d')
       
       DO jn = 1, ndim1d
@@ -88 +100 @@
          ENDDO
       ENDDO
+
+      IF( ln_timing )   CALL timing_stop('pack_arr_2d_3d')
       
    END SUBROUTINE pack_arr_3d_2d
@@ -100 +114 @@
       INTEGER, DIMENSION(ndim1d) ::  jid, jjd        
       INTEGER ::   jk, jn , ji, jj
-
-     DO jn = 1, ndim1d
+      !
+      IF( ln_timing )   CALL timing_start('unpack_arr_2d_3d')
+      !
+      DO jn = 1, ndim1d
          jid(jn) = MOD( tab_ind(jn) - 1, jpi ) + 1
          jjd(jn) = ( tab_ind(jn) - 1 ) / jpi + 1
@@ -114 +130 @@
       ENDDO
 
+      IF( ln_timing )   CALL timing_stop('unpack_arr_2d_3d')
+
    END SUBROUTINE unpack_arr_2d_3d
 
-#else
-   !!======================================================================
-   !! MODULE sedarr  :   Dummy module
-   !!======================================================================
-CONTAINS
-   SUBROUTINE pack_arr         ! Empty routine
-   END SUBROUTINE  pack_arr
-   SUBROUTINE unpack_arr         ! Empty routine
-   END SUBROUTINE  unpack_arr
-   !!======================================================================
-#endif
 END MODULE sedarr

NEMO/trunk/src/TOP/PISCES/SED/sedbtb.F90

@@ -1 +1 @@
 MODULE sedbtb
-#if defined key_sed
    !!======================================================================
    !!              ***  MODULE  sedbtb  ***
@@ -8 +7 @@
    USE sed     ! sediment global variable
    USE sedmat  ! linear system of equations
+   USE lib_mpp         ! distribued memory computing library
+
+   IMPLICIT NONE
+   PRIVATE
 
    PUBLIC sed_btb
@@ -33 +36 @@
       ! * local variables
       INTEGER :: ji, jk, js
-      REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE ::  zsol  !   solution
+      REAL(wp), DIMENSION(jpoce,jpksedm1,jpsol) ::  zsol  !   solution
       !------------------------------------------------------------------------
 
+      IF( ln_timing )  CALL timing_start('sed_btb')
+
       IF( kt == nitsed000 ) THEN
-         WRITE(numsed,*) ' sed_btb : Bioturbation  '
-         WRITE(numsed,*) ' '
+         IF (lwp) WRITE(numsed,*) ' sed_btb : Bioturbation  '
+         IF (lwp) WRITE(numsed,*) ' '
       ENDIF
 
       ! Initializations
       !----------------
-      ALLOCATE( zsol(jpoce,jpksedm1,jpsol) )
-
       zsol(:,:,:) = 0.
-
 
       ! right hand side of coefficient matrix
@@ -58 +60 @@
       ENDDO
 
-      CALL sed_mat( jpsol, jpoce, jpksedm1, zsol )
+      CALL sed_mat( jpsol, jpoce, jpksedm1, zsol, dtsed / 2.0 )
 
 
@@ -70 +72 @@
          ENDDO
       ENDDO
-     
-      DEALLOCATE( zsol )
+
+      IF( ln_timing )  CALL timing_stop('sed_btb')
 
    END SUBROUTINE sed_btb
-#else
-   !!======================================================================
-   !! MODULE sedbtb  :   Dummy module 
-   !!======================================================================
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_btb( kt )         ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'sed_btb: You should not have seen this print! error?', kt
-   END SUBROUTINE sed_btb
 
-   !!======================================================================
-
-#endif
 END MODULE sedbtb

NEMO/trunk/src/TOP/PISCES/SED/sedchem.F90

@@ -1 +1 @@
 MODULE sedchem
 
-#if defined key_sed 
    !!======================================================================
    !!                        ***  Module sedchem  ***
@@ -9 +8 @@
    USE sed     ! sediment global variable
    USE sedarr
+   USE eosbn2, ONLY : neos
+   USE lib_mpp         ! distribued memory computing library
+
+   IMPLICIT NONE
+   PRIVATE
 
    !! * Accessibility
-   PUBLIC sed_chem   
+   PUBLIC sed_chem
+   PUBLIC ahini_for_at_sed     !
+   PUBLIC solve_at_general_sed !
+
+   ! Maximum number of iterations for each method
+   INTEGER, PARAMETER :: jp_maxniter_atgen    = 20
+   REAL(wp), PARAMETER :: pp_rdel_ah_target = 1.E-4_wp
 
    !! * Module variables
    REAL(wp) :: &
-      salchl = 1./1.80655 , & ! conversion factor for salinity --> chlorinity (Wooster et al. 1969)
       calcon = 1.03E-2        ! mean calcite concentration [Ca2+] in sea water [mole/kg solution] 
 
-   REAL(wp) :: &              ! coeff. for 1. dissoc. of carbonic acid (Millero, 1995)   
-      c10 = 3670.7   , &
-      c11 = 62.008   , &
-      c12 = 9.7944   , &
-      c13 = 0.0118   , &
-      c14 = 0.000116
-
-   REAL(wp) :: &              ! coeff. for 2. dissoc. of carbonic acid (Millero, 1995)   
-      c20 = 1394.7   , &
-      c21 = 4.777    , &
-      c22 = 0.       , &
-      c23 = 0.0184   , &
-      c24 = 0.000118
-
-   REAL(wp) :: &              ! coeff. for 1. dissoc. of boric acid (Dickson and Goyet, 1994)
-      cb0  = -8966.90, &
-      cb1  = -2890.53, &
-      cb2  = -77.942 , &
-      cb3  = 1.728   , &
-      cb4  = -0.0996 , &
-      cb5  = 148.0248, &
-      cb6  = 137.1942, &
-      cb7  = 1.62142 , &
-      cb8  = -24.4344, &
-      cb9  = -25.085 , &
-      cb10 = -0.2474 , &
-      cb11 = 0.053105
-
-   REAL(wp) :: &             ! borat constants
-      bor1 = 0.000232, &
-      bor2 = 1./10.811
-
-   REAL(wp) :: &             ! constants for calculate concentrations 
-      st1  = 0.14     , &    ! for sulfate (Morris & Riley 1966)
-      st2  = 1./96.062, &
-      ks0  = 141.328  , &
-      ks1  = -4276.1  , &
-      ks2  = -23.093  , &
-      ks3  = -13856.  , &
-      ks4  = 324.57   , &
-      ks5  = -47.986  , &
-      ks6  = 35474.   , &
-      ks7  = -771.54  , &
-      ks8  = 114.723  , &
-      ks9  = -2698.   , &
-      ks10 = 1776.    , &
-      ks11 = 1.       , &
-      ks12 = -0.001005 
-
-   REAL(wp) :: &             ! constants for calculate concentrations 
-      ft1  = 0.000067   , &  ! fluorides (Dickson & Riley 1979 )
-      ft2  = 1./18.9984 , &
-      kf0  = -12.641    , &
-      kf1  = 1590.2     , &
-      kf2  = 1.525      , &
-      kf3  = 1.0        , &
-      kf4  =-0.001005
-
-   REAL(wp) :: &             ! coeff. for dissoc. of water (Dickson and Riley, 1979 )
-      cw0 = -13847.26  , &
-      cw1 = 148.9802   , &
-      cw2 = -23.6521   , &
-      cw3 = 118.67     , &
-      cw4 = -5.977     , &
-      cw5 = 1.0495     , &
-      cw6 = -0.01615
- 
-   REAL(wp) :: &             ! coeff. for dissoc. of phosphate (Millero (1974)
-      cp10 = 115.54    , &
-      cp11 = -4576.752 , &
-      cp12 = -18.453   , &
-      cp13 = -106.736  , &
-      cp14 = 0.69171   , &
-      cp15 = -0.65643  , &
-      cp16 = -0.01844  , &
-      cp20 = 172.1033  , &
-      cp21 = -8814.715 , &
-      cp22 = -27.927   , &
-      cp23 = -160.340  , &
-      cp24 = 1.3566    , &
-      cp25 = 0.37335   , &
-      cp26 = -0.05778  , &
-      cp30 = -18.126   , &
-      cp31 = -3070.75  , &
-      cp32 = 17.27039  , &
-      cp33 = 2.81197   , &
-      cp34 = -44.99486 , &
-      cp35 = -0.09984
-
-   REAL(wp) :: &             ! coeff. for dissoc. of silicates (Millero (1974)  
-      cs10 = 117.40    , &  
-      cs11 = -8904.2   , & 
-      cs12 = -19.334   , & 
-      cs13 = -458.79   , & 
-      cs14 = 3.5913    , & 
-      cs15 = 188.74    , & 
-      cs16 = -1.5998   , & 
-      cs17 = -12.1652  , & 
-      cs18 = 0.07871   , & 
-      cs19 = 0.        , & 
-      cs20 = 1.        , & 
-      cs21 = -0.001005
-
-   REAL(wp) :: &            ! coeff. for apparent solubility equilibrium 
-      ! akcc1 = -34.452 , & ! of calcite (Ingle,1975,  1800, eq. 6)
-      ! akcc2 = -39.866 , &  
-      ! akcc3 = 110.21  , &  
-      ! akcc4 = -7.5752E-6  
-      akcc1 = -171.9065 , &    ! Millero et al. 1995 from Mucci 1983
-      akcc2 = -0.077993 , &  
-      akcc3 = 2839.319  , &  
-      akcc4 = 71.595    , &  
-      akcc5 = -0.77712  , &  
-      akcc6 = 0.0028426 , &  
-      akcc7 = 178.34    , &  
-      akcc8 = -0.07711  , &  
-      akcc9 = 0.0041249
-
-   REAL(wp) :: &                 ! universal gas constant
-      rgas = 83.145              ! bar.cm3/(mol.K) or R=8.31451 J/(g.mol.K)
-
-
-   ! coeff. for seawater pressure correction (millero 95)      
-   REAL(wp), DIMENSION(8)  :: & 
-      devk1, devk2, devk3, devk4, devk5
-
-   DATA devk1/ -25.5   , -15.82    , -29.48  , -25.60    , -48.76    , -14.51   , -23.12   , -26.57   /   
-   DATA devk2/ 0.1271  , -0.0219   , 0.1622  , 0.2324    , -0.5304   , 0.1211   ,  0.1758  , 0.2020   /   
-   DATA devk3/ 0.      , 0.        , 2.608E-3, -3.6246E-3, 0.        , -0.321E-3, -2.647E-3, -3.042E-3/   
-   DATA devk4/-3.08E-3 , 1.13E-3   , -2.84E-3, -5.13E-3  , -11.76E-3 , -2.67E-3 , -5.15E-3 , -4.08E-3 /   
-   DATA devk5/0.0877E-3, -0.1475E-3, 0.      , 0.0794E-3 , -0.3692E-3, 0.0427E-3,  0.09E-3 , 0.0714E-3/
-
+   REAL(wp) ::   rgas   = 83.14472      ! universal gas constants
 
    ! coeff. for density of sea water (Millero & Poisson 1981) 
@@ -162 +40 @@
    REAL(wp), DIMENSION(6)  :: Ddsw                    
    DATA Ddsw / 999.842594 , 6.793952E-2 , -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9/
+
+  REAL(wp) :: devk10  = -25.5
+   REAL(wp) :: devk11  = -15.82
+   REAL(wp) :: devk12  = -29.48
+   REAL(wp) :: devk13  = -20.02
+   REAL(wp) :: devk14  = -18.03
+   REAL(wp) :: devk15  = -9.78
+   REAL(wp) :: devk16  = -48.76
+   REAL(wp) :: devk17  = -14.51
+   REAL(wp) :: devk18  = -23.12
+   REAL(wp) :: devk19  = -26.57
+   REAL(wp) :: devk110  = -29.48
+   !
+   REAL(wp) :: devk20  = 0.1271
+   REAL(wp) :: devk21  = -0.0219
+   REAL(wp) :: devk22  = 0.1622
+   REAL(wp) :: devk23  = 0.1119
+   REAL(wp) :: devk24  = 0.0466
+   REAL(wp) :: devk25  = -0.0090
+   REAL(wp) :: devk26  = 0.5304
+   REAL(wp) :: devk27  = 0.1211
+   REAL(wp) :: devk28  = 0.1758
+   REAL(wp) :: devk29  = 0.2020
+   REAL(wp) :: devk210  = 0.1622
+   !
+   REAL(wp) :: devk30  = 0.
+   REAL(wp) :: devk31  = 0.
+   REAL(wp) :: devk32  = 2.608E-3
+   REAL(wp) :: devk33  = -1.409e-3
+   REAL(wp) :: devk34  = 0.316e-3
+   REAL(wp) :: devk35  = -0.942e-3
+   REAL(wp) :: devk36  = 0.
+   REAL(wp) :: devk37  = -0.321e-3
+   REAL(wp) :: devk38  = -2.647e-3
+   REAL(wp) :: devk39  = -3.042e-3
+   REAL(wp) :: devk310  = -2.6080e-3
+   !
+   REAL(wp) :: devk40  = -3.08E-3
+   REAL(wp) :: devk41  = 1.13E-3
+   REAL(wp) :: devk42  = -2.84E-3
+   REAL(wp) :: devk43  = -5.13E-3
+   REAL(wp) :: devk44  = -4.53e-3
+   REAL(wp) :: devk45  = -3.91e-3
+   REAL(wp) :: devk46  = -11.76e-3
+   REAL(wp) :: devk47  = -2.67e-3
+   REAL(wp) :: devk48  = -5.15e-3
+   REAL(wp) :: devk49  = -4.08e-3
+   REAL(wp) :: devk410  = -2.84e-3
+   !
+   REAL(wp) :: devk50  = 0.0877E-3
+   REAL(wp) :: devk51  = -0.1475E-3
+   REAL(wp) :: devk52  = 0.
+   REAL(wp) :: devk53  = 0.0794E-3
+   REAL(wp) :: devk54  = 0.09e-3
+   REAL(wp) :: devk55  = 0.054e-3
+   REAL(wp) :: devk56  = 0.3692E-3
+   REAL(wp) :: devk57  = 0.0427e-3
+   REAL(wp) :: devk58  = 0.09e-3
+   REAL(wp) :: devk59  = 0.0714e-3
+   REAL(wp) :: devk510  = 0.0
 
    !! $Id$
@@ -179 +117 @@
       INTEGER, INTENT(in) :: kt                     ! time step
 
-#if ! defined key_sed_off
       INTEGER  :: ji, jj, ikt
-      REAL(wp) :: ztkel, ztc, ztc2, zpres, ztr 
-      REAL(wp) :: zsal, zsal2, zsqrt, zsal15  
-      REAL(wp) :: zis, zis2, zisqrt          
+      REAL(wp) :: ztc, ztc2
+      REAL(wp) :: zsal, zsal15  
       REAL(wp) :: zdens0, zaw, zbw, zcw    
-      REAL(wp) :: zbuf1, zbuf2 
-      REAL(wp) :: zcpexp, zcpexp2
-      REAL(wp) :: zck1p, zck2p, zck3p, zcksi   
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zchem_data
-
-#endif
+      REAL(wp), DIMENSION(jpi,jpj,15) ::   zchem_data
       !!----------------------------------------------------------------------
 
-      IF( MOD( kt - 1, nfreq ) /= 0 ) RETURN
-
-      WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt
-      WRITE(numsed,*) ' '
-
-
-#if defined key_sed_off
-      CALL sed_chem_off
-#else
+
+      IF( ln_timing )  CALL timing_start('sed_chem')
+
+      IF (lwp) WRITE(numsed,*) ' Getting Chemical constants from tracer model at time kt = ', kt
+      IF (lwp) WRITE(numsed,*) ' '
+
       ! reading variables
-      ALLOCATE( zchem_data(jpi,jpj,6) )   ;   zchem_data(:,:,:) = 0.
-
-      DO jj = 1,jpj
-         DO ji = 1, jpi
-            ikt = mbkt(ji,jj) 
-            IF ( tmask(ji,jj,ikt) == 1 ) THEN
-               zchem_data(ji,jj,1) = ak13 (ji,jj,ikt)
-               zchem_data(ji,jj,2) = ak23 (ji,jj,ikt)
-               zchem_data(ji,jj,3) = akb3 (ji,jj,ikt)
-               zchem_data(ji,jj,4) = akw3 (ji,jj,ikt)
-               zchem_data(ji,jj,5) = aksp (ji,jj,ikt)
-               zchem_data(ji,jj,6) = borat(ji,jj,ikt)
-            ENDIF
+      zchem_data(:,:,:) = rtrn
+
+      IF (ln_sediment_offline) THEN
+         CALL sed_chem_cst
+      ELSE
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               ikt = mbkt(ji,jj) 
+               IF ( tmask(ji,jj,ikt) == 1 ) THEN
+                  zchem_data(ji,jj,1) = ak13  (ji,jj,ikt)
+                  zchem_data(ji,jj,2) = ak23  (ji,jj,ikt)
+                  zchem_data(ji,jj,3) = akb3  (ji,jj,ikt)
+                  zchem_data(ji,jj,4) = akw3  (ji,jj,ikt)
+                  zchem_data(ji,jj,5) = aksp  (ji,jj,ikt)
+                  zchem_data(ji,jj,6) = borat (ji,jj,ikt)
+                  zchem_data(ji,jj,7) = ak1p3 (ji,jj,ikt)
+                  zchem_data(ji,jj,8) = ak2p3 (ji,jj,ikt)
+                  zchem_data(ji,jj,9) = ak3p3 (ji,jj,ikt)
+                  zchem_data(ji,jj,10)= aksi3 (ji,jj,ikt)
+                  zchem_data(ji,jj,11)= sio3eq(ji,jj,ikt)
+                  zchem_data(ji,jj,12)= aks3  (ji,jj,ikt)
+                  zchem_data(ji,jj,13)= akf3  (ji,jj,ikt)
+                  zchem_data(ji,jj,14)= sulfat(ji,jj,ikt)
+                  zchem_data(ji,jj,15)= fluorid(ji,jj,ikt)
+               ENDIF
+            ENDDO
          ENDDO
-      ENDDO
-
-      CALL pack_arr ( jpoce, ak1s  (1:jpoce), zchem_data(1:jpi,1:jpj,1), iarroce(1:jpoce) )
-      CALL pack_arr ( jpoce, ak2s  (1:jpoce), zchem_data(1:jpi,1:jpj,2), iarroce(1:jpoce) )
-      CALL pack_arr ( jpoce, akbs  (1:jpoce), zchem_data(1:jpi,1:jpj,3), iarroce(1:jpoce) )
-      CALL pack_arr ( jpoce, akws  (1:jpoce), zchem_data(1:jpi,1:jpj,4), iarroce(1:jpoce) )
-      CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5), iarroce(1:jpoce) )
-      CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6), iarroce(1:jpoce) )
-
-      DEALLOCATE( zchem_data )
+
+         CALL pack_arr ( jpoce, ak1s  (1:jpoce), zchem_data(1:jpi,1:jpj,1) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, ak2s  (1:jpoce), zchem_data(1:jpi,1:jpj,2) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, akbs  (1:jpoce), zchem_data(1:jpi,1:jpj,3) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, akws  (1:jpoce), zchem_data(1:jpi,1:jpj,4) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, aksps (1:jpoce), zchem_data(1:jpi,1:jpj,5) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, borats(1:jpoce), zchem_data(1:jpi,1:jpj,6) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, ak1ps (1:jpoce), zchem_data(1:jpi,1:jpj,7) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, ak2ps (1:jpoce), zchem_data(1:jpi,1:jpj,8) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, ak3ps (1:jpoce), zchem_data(1:jpi,1:jpj,9) , iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, aksis (1:jpoce), zchem_data(1:jpi,1:jpj,10), iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, sieqs (1:jpoce), zchem_data(1:jpi,1:jpj,11), iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, aks3s (1:jpoce), zchem_data(1:jpi,1:jpj,12), iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, akf3s (1:jpoce), zchem_data(1:jpi,1:jpj,13), iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, sulfats(1:jpoce), zchem_data(1:jpi,1:jpj,14), iarroce(1:jpoce) )
+         CALL pack_arr ( jpoce, fluorids(1:jpoce), zchem_data(1:jpi,1:jpj,15), iarroce(1:jpoce) )
+      ENDIF
 
       DO ji = 1, jpoce
-         ztkel   = temp(ji) + 273.16
          ztc     = temp(ji)
          ztc2    = ztc * ztc
-         zpres   = press(ji)
          ! zqtt    = ztkel * 0.01
          zsal    = salt(ji)
-         zsal2   = zsal * zsal 
-         zsqrt   = SQRT( zsal )
-         zsal15  = zsqrt * zsal
-         zlogt   = LOG( ztkel )
-         ztr     = 1./ ztkel
-         ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet)
-         zis     = 19.924 * zsal / ( 1000. - 1.005 * zsal )
-         zis2    = zis * zis
-         zisqrt  = SQRT( zis )
+         zsal15  = SQRT( zsal ) * zsal
 
          ! Density of Sea Water - F(temp,sal) [kg/m3]
@@ -256 +194 @@
          densSW(ji) = densSW(ji) * 1E-3   ! to get dens in [kg/l]
 
-         ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970) 
-         ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982) 
-         ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
-         ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.) 
-         ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP 
-         ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286 
-         ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)
-         !-----------------------------------------------------------------------------------------
-         zcpexp  = zpres / ( rgas*ztkel )
-         zcpexp2 = zpres * zcpexp
-
-         ! For Phodphates (phosphoric acid) (DOE 1994)
-         !----------------------------------------------
-         zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt &
-            &      + ( cp15*ztr + cp16 ) * zsal
-         zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt &
-            &      + ( cp25*ztr + cp26 ) * zsal
-         zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt &
-            &      + ( cp34*ztr + cp35 ) * zsal
-
-         ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP)
-         !--------------------------------------------------------
-         zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt &
-            &      + ( cs15*ztr + cs16 ) * zis &
-            &      + ( cs17*ztr + cs18 ) * zis2 &
-            &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal )
-
-
-         !K1, K2 of carbonic acid, KB of boric acid, KW (H2O)
-         !---------------------------------------------------
-         zak1p  = EXP ( zck1p  )
-         zak2p  = EXP ( zck2p  )
-         zak3p  = EXP ( zck3p  )
-         zaksil = EXP ( zcksi  )
-
-         zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(3) + devk5(3)*ztc )
-         aksis(ji)     = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         zbuf1       = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(6) + devk5(6)*ztc )
-         ak1ps(ji)   = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
- 
-         zbuf1       = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(7) + devk5(7)*ztc )
-         ak2ps(ji)   = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         zbuf1       = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(8) + devk5(8)*ztc )
-         ak3ps(ji)   = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
          ak12s  (ji) = ak1s (ji) * ak2s (ji)
          ak12ps (ji) = ak1ps(ji) * ak2ps(ji)
@@ -313 +200 @@
          calcon2(ji) = 0.01028 * ( salt(ji) / 35. ) * densSW(ji)
       ENDDO
-
        
-#endif
+      IF( ln_timing )  CALL timing_stop('sed_chem')
 
    END SUBROUTINE sed_chem
 
-#if defined key_sed_off
-
-   SUBROUTINE sed_chem_off
-      !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE sed_chem_off  ***
-      !!                    
-      !! ** Purpose :   compute chemical constants
+   SUBROUTINE ahini_for_at_sed(p_hini)
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE ahini_for_at  ***
       !!
-      !!   History :
-      !!        !  04-10  (N. Emprin, M. Gehlen )  Original code
-      !!        !  06-04  (C. Ethe)  Re-organization 
-      !!----------------------------------------------------------------------  
-      !! * Local declarations
-      INTEGER :: ji
-
-      REAL(wp) ::  ztkel, ztc, ztc2, zpres, ztr 
-      REAL(wp) ::  zsal, zsal2, zsqrt, zsal15   
-      REAL(wp) ::  zis, zis2, zisqrt            
-      REAL(wp) ::  zdens0, zaw, zbw, zcw        
-      REAL(wp) ::  zchl, zst, zft, zbuf1, zbuf2 
-      REAL(wp) ::  zcpexp, zcpexp2              
-      REAL(wp) ::  zckb, zck1, zck2, zckw       
-      REAL(wp) ::  zck1p, zck2p, zck3p, zcksi   
-      REAL(wp) ::  zak1, zak2, zakb, zakw              
-      REAL(wp) ::  zaksp0, zksp, zks, zkf 
-
+      !! Subroutine returns the root for the 2nd order approximation of the
+      !! DIC -- B_T -- A_CB equation for [H+] (reformulated as a cubic 
+      !! polynomial) around the local minimum, if it exists.
+      !! Returns * 1E-03_wp if p_alkcb <= 0
+      !!         * 1E-10_wp if p_alkcb >= 2*p_dictot + p_bortot
+      !!         * 1E-07_wp if 0 < p_alkcb < 2*p_dictot + p_bortot
+      !!                    and the 2nd order approximation does not have 
+      !!                    a solution
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT)  ::  p_hini
+      INTEGER  ::   ji, jk
+      REAL(wp)  ::  zca1, zba1
+      REAL(wp)  ::  zd, zsqrtd, zhmin
+      REAL(wp)  ::  za2, za1, za0
+      REAL(wp)  ::  p_dictot, p_bortot, p_alkcb
+
+      IF( ln_timing )  CALL timing_start('ahini_for_at_sed')
+      !
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            p_alkcb  = pwcp(ji,jk,jwalk) / densSW(ji)
+            p_dictot = pwcp(ji,jk,jwdic) / densSW(ji) 
+            p_bortot = borats(ji) / densSW(ji)
+            IF (p_alkcb <= 0.) THEN
+                p_hini(ji,jk) = 1.e-3
+            ELSEIF (p_alkcb >= (2.*p_dictot + p_bortot)) THEN
+                p_hini(ji,jk) = 1.e-10_wp
+            ELSE
+                zca1 = p_dictot/( p_alkcb + rtrn )
+                zba1 = p_bortot/ (p_alkcb + rtrn )
+           ! Coefficients of the cubic polynomial
+                za2 = aKbs(ji)*(1. - zba1) + ak1s(ji)*(1.-zca1)
+                za1 = ak1s(ji)*akbs(ji)*(1. - zba1 - zca1)    &
+                &     + ak1s(ji)*ak2s(ji)*(1. - (zca1+zca1))
+                za0 = ak1s(ji)*ak2s(ji)*akbs(ji)*(1. - zba1 - (zca1+zca1))
+                                        ! Taylor expansion around the minimum
+                zd = za2*za2 - 3.*za1   ! Discriminant of the quadratic equation
+                                        ! for the minimum close to the root
+
+                IF(zd > 0.) THEN        ! If the discriminant is positive
+                  zsqrtd = SQRT(zd)
+                  IF(za2 < 0) THEN
+                    zhmin = (-za2 + zsqrtd)/3.
+                  ELSE
+                    zhmin = -za1/(za2 + zsqrtd)
+                  ENDIF
+                  p_hini(ji,jk) = zhmin + SQRT(-(za0 + zhmin*(za1 + zhmin*(za2 + zhmin)))/zsqrtd)
+                ELSE
+                  p_hini(ji,jk) = 1.e-7
+                ENDIF
+             !
+            ENDIF
+         END DO
+      END DO
+      !
+      IF( ln_timing )  CALL timing_stop('ahini_for_at_sed')
+      !
+   END SUBROUTINE ahini_for_at_sed
+
+   !===============================================================================
+   SUBROUTINE anw_infsup_sed( p_alknw_inf, p_alknw_sup )
+
+   ! Subroutine returns the lower and upper bounds of "non-water-selfionization"
+   ! contributions to total alkalinity (the infimum and the supremum), i.e
+   ! inf(TA - [OH-] + [H+]) and sup(TA - [OH-] + [H+])
+
+   ! Argument variables
+   INTEGER :: jk
+   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_inf
+   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT) :: p_alknw_sup
+
+   DO jk = 1, jpksed
+      p_alknw_inf(:,jk) =  -pwcp(:,jk,jwpo4) / densSW(:)
+      p_alknw_sup(:,jk) =   (2. * pwcp(:,jk,jwdic) + 2. * pwcp(:,jk,jwpo4) + pwcp(:,jk,jwsil)     &
+   &                        + borats(:) ) / densSW(:)
+   END DO
+
+   END SUBROUTINE anw_infsup_sed
+
+
+   SUBROUTINE solve_at_general_sed( p_hini, zhi )
+
+   ! Universal pH solver that converges from any given initial value,
+   ! determines upper an lower bounds for the solution if required
+
+   ! Argument variables
+   !--------------------
+   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(IN)   :: p_hini
+   REAL(wp), DIMENSION(jpoce,jpksed), INTENT(OUT)  :: zhi
+
+   ! Local variables
+   !-----------------
+   INTEGER   ::  ji, jk, jn
+   REAL(wp)  ::  zh_ini, zh, zh_prev, zh_lnfactor
+   REAL(wp)  ::  zdelta, zh_delta
+   REAL(wp)  ::  zeqn, zdeqndh, zalka
+   REAL(wp)  ::  aphscale
+   REAL(wp)  ::  znumer_dic, zdnumer_dic, zdenom_dic, zalk_dic, zdalk_dic
+   REAL(wp)  ::  znumer_bor, zdnumer_bor, zdenom_bor, zalk_bor, zdalk_bor
+   REAL(wp)  ::  znumer_po4, zdnumer_po4, zdenom_po4, zalk_po4, zdalk_po4
+   REAL(wp)  ::  znumer_sil, zdnumer_sil, zdenom_sil, zalk_sil, zdalk_sil
+   REAL(wp)  ::  znumer_so4, zdnumer_so4, zdenom_so4, zalk_so4, zdalk_so4
+   REAL(wp)  ::  znumer_flu, zdnumer_flu, zdenom_flu, zalk_flu, zdalk_flu
+   REAL(wp)  ::  zalk_wat, zdalk_wat
+   REAL(wp)  ::  zfact, p_alktot, zdic, zbot, zpt, zst, zft, zsit
+   LOGICAL   ::  l_exitnow
+   REAL(wp), PARAMETER :: pz_exp_threshold = 1.0
+   REAL(wp), DIMENSION(jpoce,jpksed) :: zalknw_inf, zalknw_sup, rmask, zh_min, zh_max, zeqn_absmin
+
+   IF( ln_timing )  CALL timing_start('solve_at_general_sed')
+      !  Allocate temporary workspace
+   CALL anw_infsup_sed( zalknw_inf, zalknw_sup )
+
+   rmask(:,:) = 1.0
+   zhi(:,:)   = 0.
+
+   ! TOTAL H+ scale: conversion factor for Htot = aphscale * Hfree
+   DO jk = 1, jpksed
+      DO ji = 1, jpoce
+         IF (rmask(ji,jk) == 1.) THEN
+            p_alktot = pwcp(ji,jk,jwalk) / densSW(ji)
+            aphscale = 1. + sulfats(ji)/aks3s(ji)
+            zh_ini = p_hini(ji,jk)
+
+            zdelta = (p_alktot-zalknw_inf(ji,jk))**2 + 4.*akws(ji) / aphscale
+
+            IF(p_alktot >= zalknw_inf(ji,jk)) THEN
+               zh_min(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_inf(ji,jk) + SQRT(zdelta) )
+            ELSE
+               zh_min(ji,jk) = aphscale * (-(p_alktot-zalknw_inf(ji,jk)) + SQRT(zdelta) ) / 2.
+            ENDIF
+
+            zdelta = (p_alktot-zalknw_sup(ji,jk))**2 + 4.*akws(ji) / aphscale
+
+            IF(p_alktot <= zalknw_sup(ji,jk)) THEN
+               zh_max(ji,jk) = aphscale * (-(p_alktot-zalknw_sup(ji,jk)) + SQRT(zdelta) ) / 2.
+            ELSE
+               zh_max(ji,jk) = 2.*akws(ji) /( p_alktot-zalknw_sup(ji,jk) + SQRT(zdelta) )
+            ENDIF
+
+            zhi(ji,jk) = MAX(MIN(zh_max(ji,jk), zh_ini), zh_min(ji,jk))
+         ENDIF
+      END DO
+   END DO
+
+   zeqn_absmin(:,:) = HUGE(1._wp)
+
+   DO jn = 1, jp_maxniter_atgen
+   DO jk = 1, jpksed
+      DO ji = 1, jpoce
+         IF (rmask(ji,jk) == 1.) THEN
+
+            p_alktot = pwcp(ji,jk,jwalk) / densSW(ji)
+            zdic  = pwcp(ji,jk,jwdic) / densSW(ji)
+            zbot  = borats(ji) / densSW(ji)
+            zpt =  pwcp(ji,jk,jwpo4) / densSW(ji)
+            zsit = pwcp(ji,jk,jwsil) / densSW(ji)
+            zst = sulfats(ji)
+            zft = fluorids(ji)
+            aphscale = 1. + sulfats(ji)/aks3s(ji)
+            zh = zhi(ji,jk)
+            zh_prev = zh
+
+            ! H2CO3 - HCO3 - CO3 : n=2, m=0
+            znumer_dic = 2.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji)
+            zdenom_dic = ak1s(ji)*ak2s(ji) + zh*(ak1s(ji) + zh)
+            zalk_dic   = zdic * (znumer_dic/zdenom_dic)
+            zdnumer_dic = ak1s(ji)*ak1s(ji)*ak2s(ji) + zh     &
+                          *(4.*ak1s(ji)*ak2s(ji) + zh*ak1s(ji))
+            zdalk_dic   = -zdic*(zdnumer_dic/zdenom_dic**2)
+
+
+            ! B(OH)3 - B(OH)4 : n=1, m=0
+            znumer_bor = akbs(ji)
+            zdenom_bor = akbs(ji) + zh
+            zalk_bor   = zbot * (znumer_bor/zdenom_bor)
+            zdnumer_bor = akbs(ji)
+            zdalk_bor   = -zbot*(zdnumer_bor/zdenom_bor**2)
+
+
+            ! H3PO4 - H2PO4 - HPO4 - PO4 : n=3, m=1
+            znumer_po4 = 3.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)  &
+            &            + zh*(2.*ak1ps(ji)*ak2ps(ji) + zh* ak1ps(ji))
+            zdenom_po4 = ak1ps(ji)*ak2ps(ji)*ak3ps(ji)     &
+            &            + zh*( ak1ps(ji)*ak2ps(ji) + zh*(ak1ps(ji) + zh))
+            zalk_po4   = zpt * (znumer_po4/zdenom_po4 - 1.) ! Zero level of H3PO4 = 1
+            zdnumer_po4 = ak1ps(ji)*ak2ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)  &
+            &             + zh*(4.*ak1ps(ji)*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)         &
+            &             + zh*(9.*ak1ps(ji)*ak2ps(ji)*ak3ps(ji)                         &
+            &             + ak1ps(ji)*ak1ps(ji)*ak2ps(ji)                                &
+            &             + zh*(4.*ak1ps(ji)*ak2ps(ji) + zh * ak1ps(ji) ) ) )
+            zdalk_po4   = -zpt * (zdnumer_po4/zdenom_po4**2)
+
+            ! H4SiO4 - H3SiO4 : n=1, m=0
+            znumer_sil = aksis(ji)
+            zdenom_sil = aksis(ji) + zh
+            zalk_sil   = zsit * (znumer_sil/zdenom_sil)
+            zdnumer_sil = aksis(ji)
+            zdalk_sil   = -zsit * (zdnumer_sil/zdenom_sil**2)
+
+            ! HSO4 - SO4 : n=1, m=1
+            aphscale = 1.0 + zst/aks3s(ji)
+            znumer_so4 = aks3s(ji) * aphscale
+            zdenom_so4 = aks3s(ji) * aphscale + zh
+            zalk_so4   = zst * (znumer_so4/zdenom_so4 - 1.)
+            zdnumer_so4 = aks3s(ji) * aphscale
+            zdalk_so4   = -zst * (zdnumer_so4/zdenom_so4**2)
+
+            ! HF - F : n=1, m=1
+            znumer_flu =  akf3s(ji)
+            zdenom_flu =  akf3s(ji) + zh
+            zalk_flu   =  zft * (znumer_flu/zdenom_flu - 1.)
+            zdnumer_flu = akf3s(ji)
+            zdalk_flu   = -zft * (zdnumer_flu/zdenom_flu**2)
+
+            ! H2O - OH
+            zalk_wat   = akws(ji)/zh - zh/aphscale
+            zdalk_wat  = -akws(ji)/zh**2 - 1./aphscale
+
+            ! CALCULATE [ALK]([CO3--], [HCO3-])
+            zeqn = zalk_dic + zalk_bor + zalk_po4 + zalk_sil   &
+            &      + zalk_so4 + zalk_flu                       &
+            &      + zalk_wat - p_alktot
+
+            zalka = p_alktot - (zalk_bor + zalk_po4 + zalk_sil   &
+            &       + zalk_so4 + zalk_flu + zalk_wat)
+
+            zdeqndh = zdalk_dic + zdalk_bor + zdalk_po4 + zdalk_sil &
+            &         + zdalk_so4 + zdalk_flu + zdalk_wat
+
+            ! Adapt bracketing interval
+            IF(zeqn > 0._wp) THEN
+               zh_min(ji,jk) = zh_prev
+            ELSEIF(zeqn < 0._wp) THEN
+               zh_max(ji,jk) = zh_prev
+            ENDIF
+
+            IF(ABS(zeqn) >= 0.5_wp*zeqn_absmin(ji,jk)) THEN
+            ! if the function evaluation at the current point is
+            ! not decreasing faster than with a bisection step (at least linearly)
+            ! in absolute value take one bisection step on [ph_min, ph_max]
+            ! ph_new = (ph_min + ph_max)/2d0
+            !
+            ! In terms of [H]_new:
+            ! [H]_new = 10**(-ph_new)
+            !         = 10**(-(ph_min + ph_max)/2d0)
+            !         = SQRT(10**(-(ph_min + phmax)))
+            !         = SQRT(zh_max * zh_min)
+               zh = SQRT(zh_max(ji,jk) * zh_min(ji,jk))
+               zh_lnfactor = (zh - zh_prev)/zh_prev ! Required to test convergence below
+            ELSE
+            ! dzeqn/dpH = dzeqn/d[H] * d[H]/dpH
+            !           = -zdeqndh * LOG(10) * [H]
+            ! \Delta pH = -zeqn/(zdeqndh*d[H]/dpH) = zeqn/(zdeqndh*[H]*LOG(10))
+            !
+            ! pH_new = pH_old + \deltapH
+            !
+            ! [H]_new = 10**(-pH_new)
+            !         = 10**(-pH_old - \Delta pH)
+            !         = [H]_old * 10**(-zeqn/(zdeqndh*[H]_old*LOG(10)))
+            !         = [H]_old * EXP(-LOG(10)*zeqn/(zdeqndh*[H]_old*LOG(10)))
+            !         = [H]_old * EXP(-zeqn/(zdeqndh*[H]_old))
+
+               zh_lnfactor = -zeqn/(zdeqndh*zh_prev)
+
+               IF(ABS(zh_lnfactor) > pz_exp_threshold) THEN
+                  zh          = zh_prev*EXP(zh_lnfactor)
+               ELSE
+                  zh_delta    = zh_lnfactor*zh_prev
+                  zh          = zh_prev + zh_delta
+               ENDIF
+
+               IF( zh < zh_min(ji,jk) ) THEN
+               ! if [H]_new < [H]_min
+               ! i.e., if ph_new > ph_max then
+               ! take one bisection step on [ph_prev, ph_max]
+               ! ph_new = (ph_prev + ph_max)/2d0
+               ! In terms of [H]_new:
+               ! [H]_new = 10**(-ph_new)
+               !         = 10**(-(ph_prev + ph_max)/2d0)
+               !         = SQRT(10**(-(ph_prev + phmax)))
+               !         = SQRT([H]_old*10**(-ph_max))
+               !         = SQRT([H]_old * zh_min)
+                  zh                = SQRT(zh_prev * zh_min(ji,jk))
+                  zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
+               ENDIF
+
+               IF( zh > zh_max(ji,jk) ) THEN
+               ! if [H]_new > [H]_max
+               ! i.e., if ph_new < ph_min, then
+               ! take one bisection step on [ph_min, ph_prev]
+               ! ph_new = (ph_prev + ph_min)/2d0
+               ! In terms of [H]_new:
+               ! [H]_new = 10**(-ph_new)
+               !         = 10**(-(ph_prev + ph_min)/2d0)
+               !         = SQRT(10**(-(ph_prev + ph_min)))
+               !         = SQRT([H]_old*10**(-ph_min))
+               !         = SQRT([H]_old * zhmax)
+                  zh                = SQRT(zh_prev * zh_max(ji,jk))
+                  zh_lnfactor       = (zh - zh_prev)/zh_prev ! Required to test convergence below
+               ENDIF
+            ENDIF
+
+            zeqn_absmin(ji,jk) = MIN( ABS(zeqn), zeqn_absmin(ji,jk))
+
+            ! Stop iterations once |\delta{[H]}/[H]| < rdel
+            ! <=> |(zh - zh_prev)/zh_prev| = |EXP(-zeqn/(zdeqndh*zh_prev)) -1| < rdel
+            ! |EXP(-zeqn/(zdeqndh*zh_prev)) -1| ~ |zeqn/(zdeqndh*zh_prev)|
+            ! Alternatively:
+            ! |\Delta pH| = |zeqn/(zdeqndh*zh_prev*LOG(10))|
+            !             ~ 1/LOG(10) * |\Delta [H]|/[H]
+            !             < 1/LOG(10) * rdel
+
+            ! Hence |zeqn/(zdeqndh*zh)| < rdel
+
+            ! rdel <-- pp_rdel_ah_target
+            l_exitnow = (ABS(zh_lnfactor) < pp_rdel_ah_target)
+
+            IF(l_exitnow) THEN
+               rmask(ji,jk) = 0.
+            ENDIF
+
+            zhi(ji,jk) =  zh
+
+            IF(jn >= jp_maxniter_atgen) THEN
+               zhi(ji,jk) = -1._wp
+            ENDIF
+
+         ENDIF
+      END DO
+   END DO
+   END DO
+   !
+   IF( ln_timing )  CALL timing_stop('solve_at_general_sed')
+
+   END SUBROUTINE solve_at_general_sed
+
+   SUBROUTINE sed_chem_cst
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sed_chem_cst  ***
+      !!
+      !! ** Purpose :   Sea water chemistry computed following MOCSY protocol
+      !!                Computation is done at the bottom of the ocean only
+      !!
+      !! ** Method  : - ...
+      !!---------------------------------------------------------------------
+      INTEGER  ::   ji
+      REAL(wp), DIMENSION(jpoce) :: saltprac, temps
+      REAL(wp) ::   ztkel, ztkel1, zt , zsal  , zsal2 , zbuf1 , zbuf2
+      REAL(wp) ::   ztgg , ztgg2, ztgg3 , ztgg4 , ztgg5
+      REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
+      REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
+      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1, za2
+      REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
+      REAL(wp) ::   zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi
+      REAL(wp) ::   zst  , zft  , zcks  , zckf  , zaksp1
+      REAL(wp) ::   total2free, free2SWS, total2SWS, SWS2total
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('sed_chem_cst')
+      !
+      ! Computation of chemical constants require practical salinity
+      ! Thus, when TEOS08 is used, absolute salinity is converted to
+      ! practical salinity
+      ! -------------------------------------------------------------
+      IF (neos == -1) THEN
+         saltprac(:) = salt(:) * 35.0 / 35.16504
+      ELSE
+         saltprac(:) = temp(:)
+      ENDIF
+
+      !
+      ! Computations of chemical constants require in situ temperature
+      ! Here a quite simple formulation is used to convert
+      ! potential temperature to in situ temperature. The errors is less than
+      ! 0.04°C relative to an exact computation
+      ! ---------------------------------------------------------------------
+         DO ji = 1, jpoce
+            zpres = zkbot(ji) / 1000.
+            za1 = 0.04 * ( 1.0 + 0.185 * temp(ji) + 0.035 * (saltprac(ji) - 35.0) )
+            za2 = 0.0075 * ( 1.0 - temp(ji) / 30.0 )
+            temps(ji) = temp(ji) - za1 * zpres + za2 * zpres**2
+         END DO
 
       ! CHEMICAL CONSTANTS - DEEP OCEAN
-      !-------------------------------------
-      ! [chem constants]=mol/kg solution (or (mol/kg sol)2 for akws and aksp)
-
+      ! -------------------------------
       DO ji = 1, jpoce
-         ztkel   = temp(ji) + 273.16
-         ztc     = temp(ji)
-         ztc2    = ztc * ztc
-         zpres   = press(ji)
-         ! zqtt    = ztkel * 0.01
-         zsal    = salt(ji)
-         zsal2   = zsal * zsal 
-         zsqrt   = SQRT( zsal )
+         ! SET PRESSION ACCORDING TO SAUNDER (1980)
+         zc1 = 5.92E-3 
+         zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*zkbot(ji)))) / 4.42E-6
+         zpres = zpres / 10.0
+
+         ! SET ABSOLUTE TEMPERATURE
+         ztkel   = temps(ji) + 273.15
+         zsal    = saltprac(ji)
+         zsqrt  = SQRT( zsal )
          zsal15  = zsqrt * zsal
-         zlogt   = LOG( ztkel )
-         ztr     = 1./ ztkel
-         ! zis=ionic strength (ORNL/CDIAC-74, DOE 94,Dickson and Goyet)
-         zis     = 19.924 * zsal / ( 1000. - 1.005 * zsal )
-         zis2    = zis * zis
-         zisqrt  = SQRT( zis )
-
-
-         ! Density of Sea Water - F(temp,sal) [kg/m3]
-         zdens0 =  Ddsw(1) + Ddsw(2) * ztc + Ddsw(3) * ztc2 &
-                  + Ddsw(4) * ztc * ztc2 + Ddsw(5) * ztc2 * ztc2 &
-                  + Ddsw(6) * ztc * ztc2 * ztc2
-         zaw =  Adsw(1) + Adsw(2) * ztc + Adsw(3)* ztc2 + Adsw(4) * ztc * ztc2 &
-              + Adsw(5) * ztc2 * ztc2
-         zbw =  Bdsw(1) + Bdsw(2) * ztc + Bdsw(3) * ztc2
-         zcw =  Cdsw
-         densSW(ji) = zdens0 + zaw * zsal + zbw * zsal15 + zcw * zsal * zsal
-         densSW(ji) = densSW(ji) * 1E-3   ! to get dens in [kg/l]
-
-
-         ! FORMULA FOR CPEXP AFTER EDMOND AND GIESKES (1970) 
-         ! (REFERENCE TO CULBERSON AND PYTKOQICZ (1968) AS MADE IN BROECKER ET AL. (1982) 
-         ! IS INCORRECT; HERE RGAS IS TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
-         ! DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS MULTIPLIED BY LN(10.) 
-         ! TO ALLOW USE OF EXP-FUNCTION WITH BASIS E IN THE FORMULA FOR AKSPP 
-         ! (CF. EDMOND AND GIESKES (1970), P. 1285 AND P. 1286 (THE SMALL FORMULA ON P. 1286 
-         ! IS RIGHT AND CONSISTENT WITH THE SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285)
-         !-----------------------------------------------------------------------------------------
-         zcpexp  = zpres / ( rgas*ztkel )
+         zlogt  = LOG( ztkel )
+         ztr    = 1. / ztkel
+         zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
+         zis2   = zis * zis
+         zisqrt = SQRT( zis )
+         ztc    = temps(ji)
+
+         ! CHLORINITY (WOOSTER ET AL., 1969)
+         zcl     = zsal / 1.80655
+
+         ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
+         zst     = 0.14 * zcl /96.062
+
+         ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
+         zft     = 0.000067 * zcl /18.9984
+
+         ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
+         zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         &
+         &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
+         &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    &
+         &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
+         &         + LOG(1.0 - 0.001005 * zsal))
+
+         ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
+         zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   &
+         &         + LOG(1.0d0 - 0.001005d0*zsal)            &
+         &         + LOG(1.0d0 + zst/zcks))
+
+         ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
+         zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        &
+         &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
+         &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
+         &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      &
+         &      * zlogt + 0.053105*zsqrt*ztkel
+
+         ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO
+         ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
+         zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
+                   - 0.011555*zsal + 0.0001152*zsal*zsal)
+         zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
+                   - 0.01781*zsal + 0.0001122*zsal*zsal)
+
+         ! PKW (H2O) (MILLERO, 1995) from composite data
+         zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
+                   - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
+
+         ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
+         zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
+         &          + (-106.736*ztr + 0.69171) * zsqrt       &
+         &          + (-0.65643*ztr - 0.01844) * zsal
+
+         zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
+         &          + (-160.340*ztr + 1.3566)*zsqrt          &
+         &          + (0.37335*ztr - 0.05778)*zsal
+
+         zck3p    = -3070.75*ztr - 18.126                    &
+         &          + (17.27039*ztr + 2.81197) * zsqrt       &
+         &          + (-44.99486*ztr - 0.09984) * zsal
+
+         ! CONSTANT FOR SILICATE, MILLERO (1995)
+         zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
+         &          + (-458.79*ztr + 3.5913) * zisqrt       &
+         &          + (188.74*ztr - 1.5998) * zis           &
+         &          + (-12.1652*ztr + 0.07871) * zis2       &
+         &          + LOG(1.0 - 0.001005*zsal)
+
+         ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
+         !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
+         zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
+         &         + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
+         &         - 0.07711*zsal + 0.0041249*zsal15
+
+         ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
+         zak1    = 10**(zck1) * total2SWS
+         zak2    = 10**(zck2) * total2SWS
+         zakb    = EXP( zckb ) * total2SWS
+         zakw    = EXP( zckw )
+         zaksp1  = 10**(zaksp0)
+         zak1p   = exp( zck1p )
+         zak2p   = exp( zck2p )
+         zak3p   = exp( zck3p )
+         zaksi   = exp( zcksi )
+         zckf    = zckf * total2SWS
+
+         ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
+         !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
+         !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
+         !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
+         !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
+         !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
+         !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
+         !        & GIESKES (1970), P. 1285-1286 (THE SMALL
+         !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
+         !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
+         zcpexp  = zpres / (rgas*ztkel)
          zcpexp2 = zpres * zcpexp
 
-
-         ! chlorinity (WOOSTER ET AL., 1969)
-         !---------------------------------------
-         zchl = zsal * salchl
-
-         ! total sulfate concentration [mol/kg soln]
-         ! --------------------------------------
-         zst = st1 * zchl * st2
-
-         ! total fluoride concentration [mol/kg soln]
-         ! --------------------------------------
-         zft = ft1 * zchl * ft2
-
-         ! dissociation constant for carbonate (Mehrback 74 - Dickson & Millero 87)
-         !---------------------------------------------------------------------------
-         zck1 = c10*ztr - c11 + c12*zlogt - c13*zsal + c14*zsal2
-         zck2 = c20*ztr + c21 - c22*zlogt - c23*zsal + c24*zsal2
-
-         ! dissociation constant for sulfates (Dickson 1990)
-         !--------------------------------------------------
-         zks = EXP(  ks0 + ks1*ztr + ks2*zlogt &
-            &    + ( ks3*ztr + ks4 + ks5*zlogt ) * zisqrt &
-            &    + ( ks6*ztr + ks7 + ks8*zlogt ) * zis    &
-            &    +   ks9*ztr*zis*zisqrt + ks10*ztr*zis2   &
-            &    +   LOG( ks11 + ks12*zsal ) )
-
-         ! dissociation constant for fluorides (Dickson and Riley 79)
-         !--------------------------------------------------
-         zkf = EXP( kf0 + kf1*ztr + kf2*zisqrt + LOG( kf3 + kf4*zsal ) )
-
-         ! dissociation constant for borates (Doe 94)
-         !--------------------------------------------------
-         zckb = (  cb0 + cb1*zsqrt + cb2*zsal + cb3*zsal15 + cb4*zsal2) * ztr &
-            &  + ( cb5 + cb6*zsqrt + cb7*zsal) &
-            &  + ( cb8 + cb9*zsqrt + cb10*zsal) * zlogt &
-            &  +   cb11*zsqrt*ztkel + LOG( ( 1. + zst/zks + zft/zkf ) / ( 1. + zst/zks ) ) 
-
-         ! PKW (H2O) (DICKSON AND RILEY, 1979)
-         !--------------------------------------
-         zckw =   cw0*ztr + cw1 + cw2*zlogt &
-            & +( cw3*ztr + cw4 + cw5*zlogt )* zsqrt + cw6*zsal
-         
-         ! For Phodphates (phosphoric acid) (DOE 1994)
-         !----------------------------------------------
-         zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt &
-            &      + ( cp15*ztr + cp16 ) * zsal
-         zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt &
-            &      + ( cp25*ztr + cp26 ) * zsal
-         zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt &
-            &      + ( cp34*ztr + cp35 ) * zsal
-
-         ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP)
-         !--------------------------------------------------------
-         zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt &
-            &      + ( cs15*ztr + cs16 ) * zis &
-            &      + ( cs17*ztr + cs18 ) * zis2 &
-            &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal )
-
-         ! apparent solublity product K'SP of calcite in seawater
-         ! (S=27-43, T=2-25 DEG C) AT pres =0 (INGLE, 1975, EQ. 6)
-         ! prob: olivier a log = ln et C. Heize a LOG10(sal)
-         ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log(sal)+akcc4*tkel*tkel)
-         ! aksp0 = 1.E-7*(akcc1+akcc2*sal**(1./3.)+akcc3*log10(sal)+akcc4*tkel*tkel)
-         !--------------------------------------------------------------------
-         zaksp0 = akcc1 + akcc2*ztkel + akcc3*ztr + akcc4 * LOG10(ztkel) &
-            &  + ( akcc5 + akcc6*ztkel+ akcc7*ztr ) * zsqrt &
-            &  +  akcc8*zsal + akcc9*zsal15
-
-         !K1, K2 of carbonic acid, KB of boric acid, KW (H2O)
-         !---------------------------------------------------
-         zak1   = 10**( -zck1  )
-         zak2   = 10**( -zck2  )
-         zakb   = EXP ( zckb   ) 
-         zakw   = EXP ( zckw   )
-         zksp   = 10**( zaksp0 )
-
-
-
-         ! KB of boric acid, K1,K2 of carbonic acid pressure correction 
-         ! after Culberson and  AND Pytkowicz (1968) (CF. BROECKER ET AL., 1982) Millero 95
-         !--------------------------------------------------------------------------------
-         zbuf1       = - ( devk1(1) + devk2(1)*ztc + devk3(1)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(1) + devk5(1)*ztc )
-         ak1s(ji)    = zak1 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         zbuf1       = -( devk1(2) + devk2(2)*ztc + devk3(2)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(2) + devk5(2)*ztc )
-         ak2s(ji)    = zak2 * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(3) + devk5(3) * ztc )
-         akbs(ji)    = zakb * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         zbuf1       = - ( devk1(4) + devk2(4)*ztc + devk3(4)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(4) + devk5(4)*ztc )
-         akws(ji)    = zakw * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-
-         ! APPARENT SOLUBILITY PRODUCT K''SP OF CALCITE (OR ARAGONITE)
-         ! AS FUNCTION OF PRESSURE FOLLWING EDMOND AND GIESKES (1970)
-         ! (P. 1285) AND BERNER (1976)
-         !-----------------------------------------------------------------
-         ! aksp(ji) = aksp0*exp(zcpexp*(devks-devkst*tc))
-         ! or following Mucci
-         zbuf1      = - ( devk1(5) + devk2(5)*ztc + devk3(5)*ztc2 )
-         zbuf2      = 0.5 *( devk4(5) + devk5(5)*ztc )
-         aksps(ji)   = zksp * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         ! For Phodphates (phosphoric acid) (DOE 1994)
-         !----------------------------------------------
-         zck1p = cp10 + cp11*ztr + cp12*zlogt + ( cp13*ztr + cp14 ) * zsqrt &
-            &      + ( cp15*ztr + cp16 ) * zsal
-         zck2p = cp20 + cp21*ztr + cp22*zlogt + ( cp23*ztr + cp24 ) * zsqrt &
-            &      + ( cp25*ztr + cp26 ) * zsal
-         zck3p = cp30 + cp31*ztr + ( cp32*ztr + cp33 ) *  zsqrt &
-            &      + ( cp34*ztr + cp35 ) * zsal
-
-         ! For silicates (DOE 1994) change to mol/kg soln) (OCMIP)
-         !--------------------------------------------------------
-         zcksi = cs10 + cs11*ztr + cs12*zlogt + ( cs13*ztr + cs14) * zisqrt &
-            &      + ( cs15*ztr + cs16 ) * zis &
-            &      + ( cs17*ztr + cs18 ) * zis2 &
-            &      + LOG( 1. + cs19*zsal ) + LOG( cs20 + cs21*zsal )
-
-
-         !K1, K2 of carbonic acid, KB of boric acid, KW (H2O)
-         !---------------------------------------------------
-         zak1p  = EXP ( zck1p  )
-         zak2p  = EXP ( zck2p  )
-         zak3p  = EXP ( zck3p  )
-         zaksil = EXP ( zcksi  )
-
-         zbuf1       = - ( devk1(3) + devk2(3)*ztc + devk3(3)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(3) + devk5(3)*ztc )
-         aksis(ji)     = zaksil * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         zbuf1       = - ( devk1(6) + devk2(6)*ztc + devk3(6)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(6) + devk5(6)*ztc )
-         ak1ps(ji)   = zak1p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
- 
-         zbuf1       = - ( devk1(7) + devk2(7)*ztc + devk3(7)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(7) + devk5(7)*ztc )
-         ak2ps(ji)   = zak2p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         zbuf1       = - ( devk1(8) + devk2(8)*ztc + devk3(8)*ztc2 )
-         zbuf2       = 0.5 * ( devk4(8) + devk5(8)*ztc )
-         ak3ps(ji)   = zak3p * EXP( zbuf1*zcpexp + zbuf2*zcpexp2 )
-
-         ! total borat concentration. [mol/l]
-         ! or from Millero 1995 [mol/l] : borat(l) = 0.000416_8*(sal/35._8)*densSW(l)
-         ! --------------------------------------------------------------------------
-         borats(ji) = bor1 * zchl * bor2 * densSW(ji)
-
-         ak12s  (ji) = ak1s (ji) * ak2s (ji)
-         ak12ps (ji) = ak1ps(ji) * ak2ps(ji)
-         ak123ps(ji) = ak1ps(ji) * ak2ps(ji) * ak3ps(ji)
-
-         calcon2(ji) = 0.01028 * ( zsal / 35. ) * densSW(ji)
-
-      ENDDO
-
-   END SUBROUTINE sed_chem_off
-
-#endif
-
-#else
-   !!======================================================================
-   !! MODULE sedchem  :   Dummy module 
-   !!======================================================================
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_chem( kt )         ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'trc_stp: You should not have seen this print! error?', kt
-   END SUBROUTINE sed_chem
-
-   !!======================================================================
-
-#endif
+         ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
+         !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
+         !        (CF. BROECKER ET AL., 1982)
+
+         zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
+         ak1s(ji) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
+         ak2s(ji) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
+         akbs(ji) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
+         akws(ji) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
+         aks3s(ji) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
+         akf3s(ji) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
+         ak1ps(ji) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
+         ak2ps(ji) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
+         aksis(ji) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         ! CONVERT FROM DIFFERENT PH SCALES
+         total2free  = 1.0/(1.0 + zst/aks3s(ji))
+         free2SWS    = 1. + zst/aks3s(ji) + zft/akf3s(ji)
+         total2SWS   = total2free * free2SWS
+         SWS2total   = 1.0 / total2SWS
+
+         ! Convert to total scale
+         ak1s(ji)  = ak1s(ji)  * SWS2total
+         ak2s(ji)  = ak2s(ji)  * SWS2total
+         akbs(ji)  = akbs(ji)  * SWS2total
+         akws(ji)  = akws(ji)  * SWS2total
+         ak1ps(ji) = ak1ps(ji) * SWS2total
+         ak2ps(ji) = ak2ps(ji) * SWS2total
+         ak3ps(ji) = ak3ps(ji) * SWS2total
+         aksis(ji) = aksis(ji) * SWS2total
+         akf3s(ji) = akf3s(ji) / total2free
+
+         ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE
+         !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
+         !        (P. 1285) AND BERNER (1976)
+         zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
+         zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
+         aksps(ji) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+         ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
+         borats(ji)   = 0.0002414 * zcl / 10.811
+         sulfats(ji)  = zst
+         fluorids(ji) = zft
+
+         ! Iron and SIO3 saturation concentration from ...
+         sieqs(ji) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
+      END DO
+      !
+      IF( ln_timing )  CALL timing_stop('sed_chem_cst')
+      !
+   END SUBROUTINE sed_chem_cst
+
 
 END MODULE sedchem

NEMO/trunk/src/TOP/PISCES/SED/sedco3.F90

@@ -1 +1 @@
 MODULE sedco3
-#if defined key_sed
    !!======================================================================
    !!              ***  MODULE  sedco3  ***
@@ -7 +6 @@
    !! * Modules used
    USE sed     ! sediment global variable
+   USE sedchem
+   USE lib_mpp         ! distribued memory computing library
 
 
@@ -15 +16 @@
    PUBLIC sed_co3     
 
-
-   !! * Module variables
-   REAL(wp) :: epsmax   =  1.e-12      ! convergence limite value
-
    !!----------------------------------------------------------------------
-   !!   OPA 9.0   !   NEMO Consortium   (2003)
+   !!   OPA 9.0   !   LODYC-IPSL   (2003)
    !!----------------------------------------------------------------------
 
@@ -45 +42 @@
       !! * Arguments
       INTEGER, INTENT(in)  :: kt   ! time step
-
       !
       !---Local variables
-      INTEGER  :: jiter, ji, jk, ipt  ! dummy loop indices
+      INTEGER  :: ji, jk           ! dummy loop indices
 
-      INTEGER  :: itermax             ! maximum number of Newton-Raphson iterations
-      INTEGER  :: itime               ! number of time to perform Newton-Raphson algorithm
-      LOGICAL  :: lconv = .FALSE.     ! flag for convergence
-      REAL(wp) :: brems               !  relaxation. parameter
-      REAL(wp) :: zresm, zresm1, zhipor_min 
-      REAL(wp) :: zalk, zbor, zsil, zpo4, zdic
-      REAL(wp) :: zh_old, zh_old2, zh_old3, zh_old4
-      REAL(wp) :: zuu, zvv, zduu, zdvv 
-      REAL(wp) :: zup, zvp, zdup, zdvp
-      REAL(wp) :: zah_old, zah_olds
-      REAL(wp) :: zh_new, zh_new2, zco3
-
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zhinit, zhi
      !!----------------------------------------------------------------------
 
+      IF( ln_timing )  CALL timing_start('sed_co3')
+
       IF( kt == nitsed000 ) THEN
-         WRITE(numsed,*) ' sed_co3 : carbonate ion and proton concentration calculation  '
-         WRITE(numsed,*) ' '
+         IF (lwp) WRITE(numsed,*) ' sed_co3 : carbonate ion and proton concentration calculation  '
+         IF (lwp) WRITE(numsed,*) ' '
       ENDIF
 
-      itermax     = 30
-      brems       = 1.
-      itime       = 0
+      DO jk = 1, jpksed
+         zhinit(:,jk)   = hipor(:,jk) / densSW(:)
+      END DO
 
+      !     -------------------------------------------
+      !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
+      !     -------------------------------------------
+
+      CALL solve_at_general_sed(zhinit, zhi)
 
       DO jk = 1, jpksed
-         DO WHILE( itime <= 2 )
-            lconv  = .FALSE.
-            IF( itime > 0 ) THEN  
-               ! increase max number of iterations and relaxation parameter
-               itermax = 200
-!!               brems   = 0.3
-               IF( itime == 2 ) hipor(1:jpoce,jk) = 3.e-9 ! re-initilazation of [H] values
-            ENDIF
+         DO ji = 1, jpoce
+            co3por(ji,jk) = pwcp(ji,jk,jwdic) * ak1s(ji) * ak2s(ji) / (zhi(ji,jk)**2   &
+            &               + ak1s(ji) * zhi(ji,jk) + ak1s(ji) * ak2s(ji) + rtrn )
+            hipor(ji,jk)  = zhi(ji,jk) * densSW(ji)
+         END DO
+      END DO
 
-iflag:      DO jiter = 1, itermax
-
-                ! Store previous hi field.   
-               zresm = -1.e10
-               ipt = 1
-               DO ji = 1, jpoce
-                  ! dissociation constant are in mol/kg of solution
-                  ! convert pwcp concentration [mol/l] in mol/kg for solution
-                  zalk    = pwcp(ji,jk,jwalk) / densSW(ji)
-                  zh_old  = hipor(ji,jk) / densSW(ji)
-                  zh_old2 = zh_old  * zh_old
-                  zh_old3 = zh_old2 * zh_old
-                  zh_old4 = zh_old3 * zh_old
-                  zbor    = borats(ji) / densSW(ji)
-                  zsil    = pwcp(ji,jk,jwsil) / densSW(ji)
-                  zpo4    = pwcp(ji,jk,jwpo4) / densSW(ji)
-                  zdic    = pwcp(ji,jk,jwdic) / densSW(ji)               
-                  ! intermediate calculation 
-                  zuu     = zdic * ( ak1s(ji) / zh_old + 2.* ak12s(ji) / zh_old2 )
-                  zvv     = 1. + ak1s(ji) / zh_old + ak12s(ji) / zh_old2
-                  zduu    = zdic * ( -ak1s(ji) / zh_old2 - 4. * ak12s(ji) / zh_old3 )
-                  zdvv    = -ak1s(ji) / zh_old2 - 2. * ak12s(ji) / zh_old3
-                  zup     = zpo4 * ( ak12ps(ji) / zh_old2 + 2. * ak123ps(ji) / zh_old3 - 1.)
-                  zvp     = 1. + ak1ps(ji) / zh_old + ak12ps(ji) / zh_old2 + ak123ps(ji) / zh_old3
-                  zdup    = zpo4 * ( -2. * ak12ps(ji) / zh_old3 - 6. * ak123ps(ji) / zh_old4 )
-                  zdvp    = -ak1ps(ji) / zh_old2 - 2.* ak12ps(ji) / zh_old3 - 3. * ak123ps(ji) / zh_old4
-                  
-                  zah_old  = zuu / zvv + zbor / ( 1. + zh_old / akbs(ji) ) + &
-                     &      akws(ji) / zh_old - zh_old + zsil / ( 1. + zh_old / aksis(ji) ) + &
-                     &      zup / zvp
-                  
-                  zah_olds = ( ( zduu * zvv - zdvv * zuu ) / ( zvv * zvv ) )      - &
-                     &        zbor / akbs(ji) * ( 1. + zh_old / akbs(ji) )**(-2) - &
-                     &        akws(ji) / zh_old2 - 1. -                            &
-                     &        zsil / aksis(ji) * ( 1. + zh_old / aksis(ji) )**(-2) + &
-                     &       ( ( zdup * zvp - zdvp * zup ) / ( zvp * zvp ) )
-                  !
-                  zh_new = zh_old - brems * ( zah_old - zalk ) / zah_olds
-                  ! 
-                  zresm1 = ABS( zh_new - zh_old )
-                  IF( zresm1 > zresm ) THEN 
-                     zresm = zresm1  
-                  ENDIF
-                  !
-                  zh_new2  = zh_new * zh_new
-                  zco3   = ( ak12s(ji) * zdic ) / ( ak12s(ji) + ak1s(ji) * zh_new + zh_new2)
-                  ! again in mol/l
-                  hipor (ji,jk) = zh_new * densSW(ji)
-                  co3por(ji,jk) = zco3   * densSW(ji)
-                  
-               ENDDO  ! end loop ji
-               
-               ! convergence test
-               IF( zresm <= epsmax ) THEN
-                  lconv = .TRUE.
-                  !minimum value of hipor
-                  zhipor_min = MINVAL( hipor(1:jpoce,jk ) )
-                  EXIT iflag
-               ENDIF
-
-            ENDDO iflag
-
-            IF( lconv ) THEN
-!               WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ;  res =',zresm  
-               IF( zhipor_min < 0. ) THEN 
-                  IF ( itime == 0 ) THEN
-!                     WRITE(numsed,*) '    but hipor < 0 ; try one more time for jk = ', jk 
-!                     WRITE(numsed,*) '    with re-initialization of initial PH field '       
-                     itime = 2
-                  ELSE
-!                     WRITE(numsed,*) ' convergence after iter =', jiter, ' iterations ;  res =',zresm 
-!                     WRITE(numsed,*) '    but hipor < 0, again for second time for jk = ', jk 
-!                     WRITE(numsed,*) ' We stop : STOP '
-                     STOP
-                  ENDIF
-               ELSE
-!                  WRITE(numsed,*) ' successfull convergence for level jk = ',jk,&
-!                     &               '  after iter =', jiter, ' iterations ;  res =',zresm  
-!                  WRITE(numsed,*) ' '
-                  itime = 3
-               ENDIF
-            ELSE
-               itime = itime + 1
-               WRITE(numsed,*) ' No convergence for jk = ', jk, ' after ', itime, '  try'            
-               IF ( itime == 1 ) THEN
-                  WRITE(numsed,*) ' try one more time with more iterations and higher relax. value'
-               ELSE IF ( itime == 2 ) THEN
-                  WRITE(numsed,*) ' try one more time for with more iterations, higher relax. value'               
-                  WRITE(numsed,*) ' and with re-initialization of initial PH field ' 
-               ELSE       
-                  WRITE(numsed,*) ' No more... we stop '
-                  STOP
-               ENDIF
-            ENDIF
-         ENDDO ! End of WHILE LOOP
-     ENDDO
+     IF( ln_timing )  CALL timing_stop('sed_co3')
 
    END SUBROUTINE sed_co3
-#else
-   !!======================================================================
-   !! MODULE sedco3  :   Dummy module 
-   !!======================================================================
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_co3( kt )         ! Empty routine
-      INTEGER, INTENT(in) :: kt
-      WRITE(*,*) 'sed_co3: You should not have seen this print! error?', kt
-   END SUBROUTINE sed_co3
-
-   !!======================================================================
-
-#endif
 
 END MODULE sedco3

NEMO/trunk/src/TOP/PISCES/SED/seddsr.F90

@@ -1 +1 @@
 MODULE seddsr
-#if defined key_sed
    !!======================================================================
    !!              ***  MODULE  seddsr  ***
-   !!    Sediment : dissolution and reaction in pore water
+   !!    Sediment : dissolution and reaction in pore water related 
+   !!    related to organic matter
    !!=====================================================================
    !! * Modules used
    USE sed     ! sediment global variable
-   USE sedmat  ! linear system of equations
-   USE sedco3  ! carbonate ion and proton concentration 
+   USE sed_oce
+   USE sedini
+   USE lib_mpp         ! distribued memory computing library
+   USE lib_fortran
+
+   IMPLICIT NONE
+   PRIVATE
 
    PUBLIC sed_dsr
@@ -14 +19 @@
    !! * Module variables
 
-   REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: cons_o2
-   REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: cons_no3
-   REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: sour_no3
-   REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC :: sour_c13
-   REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC ::  dens_mol_wgt  ! molecular density 
+   REAL(wp) :: zadsnh4
+   REAL(wp), DIMENSION(jpsol), PUBLIC      :: dens_mol_wgt  ! molecular density 
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc    ! temp. variables
+
 
    !! $Id$
 CONTAINS
    
-   SUBROUTINE sed_dsr( kt ) 
+   SUBROUTINE sed_dsr( kt, knt ) 
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE sed_dsr  ***
@@ -29 +33 @@
       !!  ** 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 
+      !!  ** 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). 
       !!
       !!   History :
@@ -43 +43 @@
       !!        !  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       ! number of iteration
+      INTEGER, INTENT(in) ::   kt, knt       ! number of iteration
       ! --- local variables
-      INTEGER :: ji, jk, js, jw   ! dummy looop indices
-      INTEGER :: nv               ! number of variables in linear tridiagonal eq
-
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zrearat    ! reaction rate in pore water
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite   
-      REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zmo2_0, zmo2_1  ! temp. array for mass balance calculation
-      REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zmno3_0, zmno3_1, zmno3_2
-      REAL(wp), DIMENSION(:    ), ALLOCATABLE :: zmc13_0, zmc13_1, zmc13_2, zmc13_3
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zvolc    ! temp. variables
+      INTEGER :: ji, jk, js, jw, jn   ! dummy looop indices
+
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zrearat1, zrearat2, zrearat3    ! reaction rate in pore water
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zundsat    ! undersaturation ; indice jpwatp1 is for calcite   
+      REAL(wp), DIMENSION(jpoce,jpksed) :: zkpoc, zkpos, zkpor, zlimo2, zlimno3, zlimso4, zlimfeo    ! undersaturation ; indice jpwatp1 is for calcite   
+      REAL(wp), DIMENSION(jpoce)        :: zsumtot
       REAL(wp)  ::  zsolid1, zsolid2, zsolid3, zvolw, zreasat
-
+      REAL(wp)  ::  zsatur, zsatur2, znusil, zkpoca, zkpocb, zkpocc
+      REAL(wp)  ::  zratio, zgamma, zbeta, zlimtmp, zundsat2
       !!
       !!----------------------------------------------------------------------
 
-      IF( kt == nitsed000 ) THEN
-         WRITE(numsed,*) ' sed_dsr : Dissolution reaction '
-         WRITE(numsed,*) ' '
-         ! 
-         ALLOCATE( dens_mol_wgt((jpoce) ) 
-         dens_mol_wgt(1:jpsol) = dens / mol_wgt(1:jpsol)
-         ! 
-         ALLOCATE( cons_o2 (jpoce) ) ;   ALLOCATE( cons_no3(jpoce) ) 
-         ALLOCATE( sour_no3(jpoce) ) ;   ALLOCATE( sour_c13(jpoce) )  
+      IF( ln_timing )  CALL timing_start('sed_dsr')
+!
+      IF( kt == nitsed000 .AND. knt == 1 ) THEN
+         IF (lwp) THEN
+            WRITE(numsed,*) ' sed_dsr : Dissolution reaction '
+            WRITE(numsed,*) ' '
+         ENDIF
       ENDIF
 
-      ! Initialization of data for mass balance calculation
-      !---------------------------------------------------
-
-      tokbot(:,:)  = 0.
-      cons_o2 (:)  = 0. 
-      cons_no3(:)  = 0.  
-      sour_no3(:)  = 0.  
-      sour_c13(:)  = 0.       
-   
-      ! Initializations
-      !----------------------
-      ALLOCATE( zmo2_0 (jpoce) ) ;  ALLOCATE( zmo2_1 (jpoce) ) 
-      ALLOCATE( zmno3_0(jpoce) ) ;  ALLOCATE( zmno3_1(jpoce) )  ;  ALLOCATE( zmno3_2(jpoce) ) 
-      ALLOCATE( zmc13_0(jpoce) ) ;  ALLOCATE( zmc13_1(jpoce) )  ;  ALLOCATE( zmc13_2(jpoce) ) ;  ALLOCATE( zmc13_3(jpoce) ) 
-
-      zmo2_0 (:) = 0.  ; zmo2_1 (:) = 0.
-      zmno3_0(:) = 0.  ; zmno3_1(:) = 0.  ;  zmno3_2(:) = 0.
-      zmc13_0(:) = 0.  ; zmc13_1(:) = 0.  ;  zmc13_2(:) = 0.  ; zmc13_3(:) = 0.
+     ! Initializations
+     !----------------------
       
-      ALLOCATE( zrearat(jpoce,jpksed,3) ) ;  ALLOCATE( zundsat(jpoce,jpksed,3) )      
-      zrearat(:,:,:)   = 0.    ;   zundsat(:,:,:) = 0. 
-
-
-      ALLOCATE( zvolc(jpoce,jpksed,jpsol) ) 
-      zvolc(:,:,:)   = 0.
-
-      !--------------------------------------------------------------------
-      ! Temporary accomodation to take account of  particule rain deposition
-      !---------------------------------------------------------------------
-      
-      
-      ! 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)
+      zrearat1(:,:)   = 0.    ;   zundsat(:,:) = 0. ; zkpoc(:,:) = 0.
+      zlimo2 (:,:)    = 0.    ;   zlimno3(:,:) = 0. ; zrearat2(:,:) = 0.
+      zlimso4(:,:)    = 0.    ;   zkpor(:,:)   = 0. ; zrearat3(:,:) = 0.
+      zkpos  (:,:)    = 0.
+      zsumtot(:)      = rtrn
+  
+      ALLOCATE( zvolc(jpoce, jpksed, jpsol) )
+      zvolc(:,:,:)    = 0.
+      zadsnh4 = 1.0 / ( 1.0 + adsnh4 )
+
+      ! Inhibition terms for the different redox equations
+      ! --------------------------------------------------
+      DO jk = 1, jpksed
+         DO ji = 1, jpoce
+            zkpoc(ji,jk) = reac_pocl 
+            zkpos(ji,jk) = reac_pocs
+            zkpor(ji,jk) = reac_pocr
+         END DO
+      END DO
 
       ! Conversion of volume units
@@ -120 +99 @@
       DO js = 1, jpsol
          DO jk = 1, jpksed
-            DO ji = 1, jpoce    
+            DO ji = 1, jpoce
                zvolc(ji,jk,js) = ( vols3d(ji,jk) * dens_mol_wgt(js) ) /  &
-                  &              ( volw3d(ji,jk) * 1.e-3 )     
+                  &              ( volw3d(ji,jk) * 1.e-3 )
             ENDDO
          ENDDO
       ENDDO
 
-      ! 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
-            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.
-         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)
-         END DO
-      ENDDO
-
- 
       !----------------------------------------------------------
-      ! 5.  Beginning of  Pore Water diffusion and solid reaction
+      ! 5.  Beginning of 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
-      !-------------------------------------------------------------------------
-
-      !number of variables
-      nv  = 3
-     
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            ! For Silicic Acid and clay
-            zundsat(ji,jk,1) = sat_sil   - pwcp(ji,jk,jwsil)
-            zundsat(ji,jk,2) = sat_clay
-            ! For O2
-            zundsat(ji,jk,3) = pwcp(ji,jk,jwoxy) / so2ut 
-         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
-         zrearat(ji,1,:) = 0.
-      ENDDO
-
+      zrearat1(:,:) = 0.
+      zrearat2(:,:) = 0.
+      zrearat3(:,:) = 0.
+
+      zundsat(:,:) = pwcp(:,:,jwoxy)
+
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 )
+            zsolid1 = zvolc(ji,jk,jspoc)  * solcp(ji,jk,jspoc)
+            zsolid2 = zvolc(ji,jk,jspos)  * solcp(ji,jk,jspos)
+            zsolid3 = zvolc(ji,jk,jspor)  * solcp(ji,jk,jspor)
+            zkpoca  = zkpoc(ji,jk) * zlimo2(ji,jk)
+            zkpocb  = zkpos(ji,jk) * zlimo2(ji,jk)
+            zkpocc  = zkpor(ji,jk) * zlimo2(ji,jk)
+            zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+            &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 )
+            zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+            &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 )
+            zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+            &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 )
+         ENDDO
+      ENDDO
 
       ! left hand side of coefficient matrix
-      DO jk = 2, jpksed
-         DO ji = 1, jpoce
-            zsolid1 = zvolc(ji,jk,jsopal) * solcp(ji,jk,jsopal)
-            zsolid2 = zvolc(ji,jk,jsclay) * solcp(ji,jk,jsclay)
-            zsolid3 = zvolc(ji,jk,jspoc)  * solcp(ji,jk,jspoc)
-
-            zrearat(ji,jk,1)  = ( reac_sil * dtsed * zsolid1 ) / &
-               &                ( 1. + reac_sil * dtsed * zundsat(ji,jk,1 ) )
-            zrearat(ji,jk,2)  = ( reac_clay * dtsed * zsolid2 ) / &
-               &                ( 1. + reac_clay * dtsed * zundsat(ji,jk,2 ) )
-            zrearat(ji,jk,3)  = ( reac_poc  * dtsed * zsolid3 ) / &
-               &                ( 1. + reac_poc  * dtsed * zundsat(ji,jk,3 ) )
-         ENDDO
-      ENDDO
-
-
-      CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat )
-
+!      DO jn = 1, 5
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+jflag1:     DO jn = 1, 10
+               zsolid1 = zvolc(ji,jk,jspoc)  * solcp(ji,jk,jspoc)
+               zsolid2 = zvolc(ji,jk,jspos)  * solcp(ji,jk,jspos)
+               zsolid3 = zvolc(ji,jk,jspor)  * solcp(ji,jk,jspor)
+               zbeta   = xksedo2 - pwcp(ji,jk,jwoxy) + so2ut * ( zrearat1(ji,jk)    &
+               &         + zrearat2(ji,jk) + zrearat3(ji,jk) )
+               zgamma = - xksedo2 * pwcp(ji,jk,jwoxy)
+               zundsat2 = zundsat(ji,jk)
+               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0
+               zlimo2(ji,jk) = 1.0 / ( zundsat(ji,jk) + xksedo2 )
+               zkpoca  = zkpoc(ji,jk) * zlimo2(ji,jk)
+               zkpocb  = zkpos(ji,jk) * zlimo2(ji,jk)
+               zkpocc  = zkpor(ji,jk) * zlimo2(ji,jk)
+               zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+               &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 )
+               zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+               &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 )
+               zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+               &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 )
+               IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN
+                  EXIT jflag1
+               ENDIF
+            END DO jflag1
+         END DO
+      END DO
 
       ! New solid concentration values (jk=2 to jksed) for each couple 
-      DO js = 1, nv
-         DO jk = 2, jpksed
-            DO ji = 1, jpoce
-               zreasat = zrearat(ji,jk,js) * zundsat(ji,jk,js) / zvolc(ji,jk,js)
-               solcp(ji,jk,js) = solcp(ji,jk,js) - zreasat
-            ENDDO
-         ENDDO
-      ENDDO
-      ! mass of O2/NO3 before POC remin. for mass balance check 
-      ! det. of o2 consomation/NO3 production Mc13
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zvolw = volw3d(ji,jk) * 1.e-3
-            zmo2_0 (ji)  = zmo2_0 (ji) + pwcp(ji,jk,jwoxy) * zvolw
-            zmno3_0(ji)  = zmno3_0(ji) + pwcp(ji,jk,jwno3) * zvolw
-            zmc13_0(ji)  = zmc13_0(ji) + pwcp(ji,jk,jwc13) * zvolw
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zreasat = zrearat1(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc)
+            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat
+            zreasat = zrearat2(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos)
+            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat
+            zreasat = zrearat3(ji,jk) * zlimo2(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor)
+            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat
          ENDDO
       ENDDO
 
       ! New pore water concentrations    
-      DO jk = 1, jpksed
+      DO jk = 2, jpksed
          DO ji = 1, jpoce
             ! Acid Silicic 
-            pwcp(ji,jk,jwsil)  = sat_sil - zundsat(ji,jk,1)            
-            ! For O2 (in mol/l)
-            pwcp(ji,jk,jwoxy)  = zundsat(ji,jk,3) * so2ut 
-            zreasat = zrearat(ji,jk,3) * zundsat(ji,jk,3)    ! oxygen         
+            pwcp(ji,jk,jwoxy)  = zundsat(ji,jk)
+            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimo2(ji,jk) * zundsat(ji,jk)    ! oxygen         
             ! For DIC
             pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
-            ! For nitrates
-            pwcp(ji,jk,jwno3)  = pwcp(ji,jk,jwno3) + zreasat * srno3            
+            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3
             ! For Phosphate (in mol/l)
-            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r            
+            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r
+            ! For iron (in mol/l)
+            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat
             ! For alkalinity
-            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) - zreasat * ( srno3 + 2.* spo4r )           
-            ! For DIC13
-            pwcp(ji,jk,jwc13)  = pwcp(ji,jk,jwc13) + zreasat * rc13P * pdb
-         ENDDO
-      ENDDO
-
-
-      ! Mass of O2 for mass balance check and det. of o2 consomation
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zvolw = volw3d(ji,jk) * 1.e-3
-            zmo2_1 (ji) = zmo2_1 (ji) + pwcp(ji,jk,jwoxy) * zvolw
-            zmno3_1(ji) = zmno3_1(ji) + pwcp(ji,jk,jwno3) * zvolw
-            zmc13_1(ji) = zmc13_1(ji) + pwcp(ji,jk,jwc13) * zvolw
-         ENDDO
-      ENDDO
-
-      DO ji = 1, jpoce
-         cons_o2 (ji) = zmo2_0 (ji) - zmo2_1 (ji)
-         sour_no3(ji) = zmno3_1(ji) - zmno3_0(ji)     
-         sour_c13(ji) = zmc13_1(ji) - zmc13_0(ji) 
-      ENDDO
- 
+            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r )
+            ! Ammonium
+            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4
+            ! Ligands
+            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat - reac_ligc * sedligand(ji,jk)
+         ENDDO
+      ENDDO
 
       !--------------------------------------------------------------------
@@ -282 +205 @@
       !--------------------------------------------------------------------
 
-      nv = 1
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zundsat(ji,jk,1) = pwcp(ji,jk,jwno3) / srDnit
-         ENDDO
-      ENDDO
-      DO jk = 2, jpksed
-         DO ji = 1, jpoce
-            IF( pwcp(ji,jk,jwoxy) < sthrO2 ) THEN
+      zrearat1(:,:) = 0.
+      zrearat2(:,:) = 0.
+      zrearat3(:,:) = 0.
+
+      zundsat(:,:) = pwcp(:,:,jwno3)
+
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 )
+            zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc)
+            zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos)
+            zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor)
+            zkpoca = zkpoc(ji,jk) * zlimno3(ji,jk)
+            zkpocb = zkpos(ji,jk) * zlimno3(ji,jk)
+            zkpocc = zkpor(ji,jk) * zlimno3(ji,jk)
+            zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+            &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 )
+            zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+            &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 )
+            zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+            &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 )
+        END DO
+      END DO
+
+!      DO jn = 1, 5
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+jflag2:    DO jn = 1, 10
+               zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) )
                zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc)
-               zrearat(ji,jk,1) = ( reac_no3 * dtsed * zsolid1 ) / &
-                  &                    ( 1. + reac_no3 * dtsed * zundsat(ji,jk,1 ) )
-            ELSE
-               zrearat(ji,jk,1) = 0.
-            ENDIF
-         END DO
-      END DO
-
-
-      ! solves tridiagonal system
-      CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat )
+               zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos)
+               zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor)
+               zbeta   = xksedno3 - pwcp(ji,jk,jwno3) + srDnit * ( zrearat1(ji,jk)    &
+               &         + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimtmp
+               zgamma = - xksedno3 * pwcp(ji,jk,jwno3)
+               zundsat2 = zundsat(ji,jk)
+               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0
+               zlimno3(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) / ( zundsat(ji,jk) + xksedno3 )
+               zkpoca  = zkpoc(ji,jk) * zlimno3(ji,jk)
+               zkpocb  = zkpos(ji,jk) * zlimno3(ji,jk)
+               zkpocc  = zkpor(ji,jk) * zlimno3(ji,jk)
+               zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+               &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 )
+               zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+               &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 )
+               zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+               &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 )
+               IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN
+                  EXIT jflag2
+               ENDIF
+            END DO jflag2
+         END DO
+      END DO
 
 
@@ -308 +263 @@
       DO jk = 2, jpksed
          DO ji = 1, jpoce
-            zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1) / zvolc(ji,jk,jspoc)
+            zreasat = zrearat1(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc)
             solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat
+            zreasat = zrearat2(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos)
+            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat
+            zreasat = zrearat3(ji,jk) * zlimno3(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor)
+            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat
          ENDDO
       ENDDO
 
       ! New dissolved concentrations
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1)    
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
             ! For nitrates
-            pwcp(ji,jk,jwno3)  =  zundsat(ji,jk,1) * srDnit
+            pwcp(ji,jk,jwno3)  =  zundsat(ji,jk)
+            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimno3(ji,jk) * zundsat(ji,jk)
             ! For DIC
             pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
+            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3
             ! For Phosphate (in mol/l)
             pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r            
+            ! Ligands
+            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat
+            ! For iron (in mol/l)
+            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat
             ! For alkalinity
-            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srDnit - 2.* spo4r )           
-            ! For DIC13
-            pwcp(ji,jk,jwc13)  = pwcp(ji,jk,jwc13) + zreasat * rc13P * pdb
-         ENDDO
-      ENDDO
-
-
-      ! Mass of O2 for mass balance check and det. of o2 consomation
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zvolw = volw3d(ji,jk) * 1.e-3
-            zmno3_2(ji) = zmno3_2(ji) + pwcp(ji,jk ,jwno3) * zvolw
-            zmc13_2(ji) = zmc13_2(ji) + pwcp(ji,jk ,jwc13) * zvolw   
-         ENDDO
-      ENDDO
-
-      DO ji = 1, jpoce
-         cons_no3(ji) = zmno3_1(ji) - zmno3_2(ji)  
-         sour_c13(ji) = sour_c13(ji) + zmc13_2(ji) - zmc13_1(ji)     
-      ENDDO
-
-
-      !---------------------------
-      ! Solves PO4 diffusion 
-      !----------------------------
-
-      nv = 1
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zundsat(ji,jk,1) = pwcp(ji,jk,jwpo4)
-            zrearat(ji,jk,1) = 0.
-         ENDDO
-      ENDDO
-
-
-      ! solves tridiagonal system
-      CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat )
-
-
-      ! New undsaturation values and dissolved concentrations
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            pwcp(ji,jk,jwpo4) = zundsat(ji,jk,1)
-         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 )
-
-
-      nv = 1
-      ! *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
-            zundsat(ji,jk,1) = aksps(ji) * densSW(ji) * densSW(ji) / calcon2(ji) &
-               &                     - co3por(ji,jk)
-            ! positive values of undersaturation
-            zundsat(ji,jk,1) = MAX( 0., zundsat(ji,jk,1) )            
-         ENDDO
-      ENDDO
-
-      DO jk = 2, jpksed
-         DO ji = 1, jpoce
-            zsolid1 = zvolc(ji,jk,jscal) * solcp(ji,jk,jscal)
-            zrearat(ji,jk,1) = ( reac_cal * dtsed * zsolid1 ) / &
-                  &               ( 1. + reac_cal * dtsed * zundsat(ji,jk,1) )
-         END DO
-      END DO
-
-
-      ! solves tridiagonal system
-      CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat )
-
-
-      ! New solid concentration values (jk=2 to jksed) for cacO3
-      DO jk = 2, jpksed
-         DO ji = 1, jpoce
-            zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1) / zvolc(ji,jk,jscal)
-            solcp(ji,jk,jscal) = solcp(ji,jk,jscal) - zreasat
-         ENDDO
-      ENDDO
+            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srDnit + srno3 * zadsnh4 - 2.* spo4r )           
+            ! Ammonium
+            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4
+         ENDDO
+      ENDDO
+
+      !--------------------------------------------------------------------
+      ! Begining POC iron reduction
+      ! (indice n�5 for couple POFe(OH)3 ie solcp(:,:,jspoc)/pwcp(:,:,jsfeo))
+      !--------------------------------------------------------------------
+
+      zrearat1(:,:) = 0.
+      zrearat2(:,:) = 0.
+      zrearat3(:,:) = 0.
+
+      zundsat(:,:) = solcp(:,:,jsfeo)
+
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    &
+            &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo )
+            zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc)
+            zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos)
+            zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor)
+            zkpoca = zkpoc(ji,jk) * zlimfeo(ji,jk)
+            zkpocb = zkpos(ji,jk) * zlimfeo(ji,jk)
+            zkpocc = zkpor(ji,jk) * zlimfeo(ji,jk)
+            zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+            &                    ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 )
+            zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+            &                    ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 )
+            zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+            &                    ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 )
+         END DO
+      END DO
+
+!      DO jn = 1, 5
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+jflag3:     DO jn = 1, 10
+               zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    &
+               &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) )
+               zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc)
+               zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos)
+               zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor)
+               zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) / zvolc(ji,jk,jsfeo)
+               zbeta   = xksedfeo - solcp(ji,jk,jsfeo) + 4.0 * zreasat * zlimtmp
+               zgamma  = -xksedfeo * solcp(ji,jk,jsfeo)
+               zundsat2 = zundsat(ji,jk)
+               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0
+               zlimfeo(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    &
+               &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) / ( zundsat(ji,jk) + xksedfeo )
+               zkpoca  = zkpoc(ji,jk) * zlimfeo(ji,jk)
+               zkpocb  = zkpos(ji,jk) * zlimfeo(ji,jk)
+               zkpocc  = zkpor(ji,jk) * zlimfeo(ji,jk)
+               zrearat1(ji,jk) = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+               &                    ( 1. + zkpoca * zundsat(ji,jk) * dtsed2 )
+               zrearat2(ji,jk) = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+               &                    ( 1. + zkpocb * zundsat(ji,jk) * dtsed2 )
+               zrearat3(ji,jk) = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+               &                    ( 1. + zkpocc * zundsat(ji,jk) * dtsed2 )
+               IF ( ABS( (zundsat(ji,jk)-zundsat2)/( MAX(0.,zundsat2)+rtrn)) < 1E-8 ) THEN
+                  EXIT jflag3
+               ENDIF
+            END DO jflag3
+         END DO
+      END DO
+
+
+         ! New solid concentration values (jk=2 to jksed) for each couple 
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zreasat = zrearat1(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc)
+            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat
+            zreasat = zrearat2(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos)
+            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat
+            zreasat = zrearat3(ji,jk) * zlimfeo(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor)
+            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat
+         END DO
+      END DO
 
       ! New dissolved concentrations
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zreasat = zrearat(ji,jk,1) * zundsat(ji,jk,1)    
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimfeo(ji,jk) * zundsat(ji,jk)
+            ! For FEOH
+            solcp(ji,jk,jsfeo) = zundsat(ji,jk)
             ! For DIC
             pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
+            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3
+            ! For Phosphate (in mol/l)
+            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * ( spo4r + 4.0 * redfep )
+            ! Ligands
+            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat
+            ! For iron (in mol/l)
+            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat
             ! For alkalinity
-            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + 2.* zreasat 
-            ! For DIC13
-            pwcp(ji,jk,jwc13)  = pwcp(ji,jk,jwc13) + zreasat * rc13Ca * pdb
-         ENDDO
-      ENDDO
-
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce
-            zmc13_3(ji) = zmc13_3(ji) + pwcp(ji,jk,jwc13) * volw3d(ji,jk) * 1.e-3                  
-         ENDDO
-      ENDDO
-
-      DO ji = 1, jpoce     
-         sour_c13(ji) = sour_c13(ji) + zmc13_3(ji) - zmc13_2(ji)   
-      ENDDO
-      
-      !-------------------------------------------------
-      ! Beginning DIC, Alkalinity and DIC13 diffusion
-      !-------------------------------------------------
-      
-      nv = 3
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce      
-            zundsat(ji,jk,1) = pwcp(ji,jk,jwdic)
-            zundsat(ji,jk,2) = pwcp(ji,jk,jwalk)
-            zundsat(ji,jk,3) = pwcp(ji,jk,jwc13)
-      
-            zrearat(ji,jk,1) = 0.
-            zrearat(ji,jk,2) = 0.
-            zrearat(ji,jk,3) = 0.
-      
-         ENDDO
-      ENDDO
-
-
-      ! solves tridiagonal system
-      CALL sed_mat( nv, jpoce, jpksed, zrearat, zundsat )
-
-
-      ! New dissolved concentrations      
-      DO jk = 1, jpksed
-         DO ji = 1, jpoce                      
-            pwcp(ji,jk,jwdic) = zundsat(ji,jk,1)
-            pwcp(ji,jk,jwalk) = zundsat(ji,jk,2)
-            pwcp(ji,jk,jwc13) = zundsat(ji,jk,3)
-         ENDDO
-      ENDDO            
-      
-      !----------------------------------
-      !   Back to initial geometry
-      !-----------------------------
-      
-      !---------------------------------------------------------------------
-      !   1/ Compensation for ajustement of the bottom water concentrations
-      !      (see note n° 1 about *por(2))
+            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + 8.0 * zreasat
+            ! Ammonium
+            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4
+            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + zreasat * 4.0
+         ENDDO
+      ENDDO
+
       !--------------------------------------------------------------------
-      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
-      !--------------------------------
+      ! Begining POC denitrification and NO3- diffusion
+      ! (indice n�5 for couple POC/NO3- ie solcp(:,:,jspoc)/pwcp(:,:,jwno3))
+      !--------------------------------------------------------------------
+
+      zrearat1(:,:) = 0.
+      zrearat2(:,:) = 0.
+      zrearat3(:,:) = 0.
+
+      zundsat(:,:) = pwcp(:,:,jwso4)
+
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    &
+            &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo)  &
+            &                / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 )
+            zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc)
+            zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos)
+            zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor)
+            zkpoca = zkpoc(ji,jk) * zlimso4(ji,jk)
+            zkpocb = zkpos(ji,jk) * zlimso4(ji,jk)
+            zkpocc = zkpor(ji,jk) * zlimso4(ji,jk)
+            zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+            &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 )
+            zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+            &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 )
+            zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+            &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 )
+        END DO
+      END DO
+!
+!      DO jn = 1, 5 
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+jflag4:     DO jn = 1, 10
+               zlimtmp = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    &
+               &         / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo)  &
+               &         / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) 
+               zsolid1 = zvolc(ji,jk,jspoc) * solcp(ji,jk,jspoc)
+               zsolid2 = zvolc(ji,jk,jspos) * solcp(ji,jk,jspos)
+               zsolid3 = zvolc(ji,jk,jspor) * solcp(ji,jk,jspor)
+               zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) 
+               zbeta   = xksedso4 - pwcp(ji,jk,jwso4) + 0.5 * zreasat * zlimtmp
+               zgamma = - xksedso4 * pwcp(ji,jk,jwso4)
+               zundsat2 = zundsat(ji,jk)
+               zundsat(ji,jk) = ( - zbeta + SQRT( zbeta**2 - 4.0 * zgamma ) ) / 2.0
+               zlimso4(ji,jk) = ( 1.0 - pwcp(ji,jk,jwoxy) * zlimo2(ji,jk) ) * ( 1.0 - pwcp(ji,jk,jwno3)    &
+               &                / ( pwcp(ji,jk,jwno3) + xksedno3 ) ) * ( 1. - solcp(ji,jk,jsfeo)  &
+               &                / ( solcp(ji,jk,jsfeo) + xksedfeo ) ) / ( zundsat(ji,jk) + xksedso4 )
+               zkpoca  = zkpoc(ji,jk) * zlimso4(ji,jk)
+               zkpocb  = zkpos(ji,jk) * zlimso4(ji,jk)
+               zkpocc  = zkpor(ji,jk) * zlimso4(ji,jk)
+               zrearat1(ji,jk)  = ( zkpoc(ji,jk) * dtsed2 * zsolid1 ) / &
+               &                 ( 1. + zkpoca * zundsat(ji,jk ) * dtsed2 )
+               zrearat2(ji,jk)  = ( zkpos(ji,jk) * dtsed2 * zsolid2 ) / &
+               &                 ( 1. + zkpocb * zundsat(ji,jk ) * dtsed2 )
+               zrearat3(ji,jk)  = ( zkpor(ji,jk) * dtsed2 * zsolid3 ) / &
+               &                 ( 1. + zkpocc * zundsat(ji,jk ) * dtsed2 )
+               IF ( ABS( (zundsat(ji,jk)-zundsat2)/(zundsat2+rtrn)) < 1E-8 ) THEN
+                  EXIT jflag4
+               ENDIF
+            END DO jflag4
+         END DO
+      END DO
+
+     ! New solid concentration values (jk=2 to jksed) for each couple 
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            zreasat = zrearat1(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspoc)
+            solcp(ji,jk,jspoc) = solcp(ji,jk,jspoc) - zreasat
+            zreasat = zrearat2(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspos)
+            solcp(ji,jk,jspos) = solcp(ji,jk,jspos) - zreasat
+            zreasat = zrearat3(ji,jk) * zlimso4(ji,jk) * zundsat(ji,jk) / zvolc(ji,jk,jspor)
+            solcp(ji,jk,jspor) = solcp(ji,jk,jspor) - zreasat
+         ENDDO
+      ENDDO
+!
+      ! New dissolved concentrations
+      DO jk = 2, jpksed
+         DO ji = 1, jpoce
+            ! For sulfur
+            pwcp(ji,jk,jwh2s)  = pwcp(ji,jk,jwh2s) - ( zundsat(ji,jk) - pwcp(ji,jk,jwso4) ) 
+            pwcp(ji,jk,jwso4)  =  zundsat(ji,jk)
+            zreasat = ( zrearat1(ji,jk) + zrearat2(ji,jk) + zrearat3(ji,jk) ) * zlimso4(ji,jk) * zundsat(ji,jk)
+            ! For DIC
+            pwcp(ji,jk,jwdic)  = pwcp(ji,jk,jwdic) + zreasat
+            zsumtot(ji) = zsumtot(ji) + zreasat / dtsed2 * volw3d(ji,jk) * 1.e-3 * 86400. * 365. * 1E3
+            ! For Phosphate (in mol/l)
+            pwcp(ji,jk,jwpo4)  = pwcp(ji,jk,jwpo4) + zreasat * spo4r
+            ! Ligands
+            sedligand(ji,jk)   = sedligand(ji,jk) + ratligc * zreasat
+            ! For iron (in mol/l)
+            pwcp(ji,jk,jwfe2)  = pwcp(ji,jk,jwfe2) + fecratio(ji) * zreasat
+            ! For alkalinity
+            pwcp(ji,jk,jwalk)  = pwcp(ji,jk,jwalk) + zreasat * ( srno3 * zadsnh4 - 2.* spo4r ) + zreasat
+            ! Ammonium
+            pwcp(ji,jk,jwnh4)  = pwcp(ji,jk,jwnh4) + zreasat * srno3 * zadsnh4
+         ENDDO
+      ENDDO
+
+      ! Oxydation of the reduced products. Here only ammonium and ODU is accounted for
+      ! There are two options here: A simple time splitting scheme and a modified 
+     ! Patankar scheme
+      ! ------------------------------------------------------------------------------
+
+      call sed_dsr_redoxb
+
+      ! -------------------------------------------------------------- 
+      !    4/ Computation of the bioturbation coefficient
+      !       This parameterization is taken from Archer et al. (2002)
+      ! --------------------------------------------------------------
+
       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
-      
-      !----------------------------------------------------------------------
-      !    4/ Saving new amount of material in dzkbot for mass balance check
-      !       tokbot in [mol] (implicit *1cm*1cm for spacial dim)
-      !----------------------------------------------------------------------
-      DO jw = 1, jpwat 
-         DO ji = 1, jpoce
-            tokbot(ji,jw) = pwcp(ji,1,jw) * 1.e-3 * dzkbot(ji)
-         END DO
-      ENDDO
-
-      DEALLOCATE( zmo2_0  ) ;  DEALLOCATE( zmno3_1 )  ;  DEALLOCATE( zmno3_2 ) 
-      DEALLOCATE( zmc13_0 ) ;  DEALLOCATE( zmc13_1 )  ;  DEALLOCATE( zmc13_2 ) ;  DEALLOCATE( zmc13_3 ) 
-      
-      DEALLOCATE( zrearat ) ;  DEALLOCATE( zundsat )  ;  DEALLOCATE( zvolc )   
-      
+         db(ji,:) = dbiot * zsumtot(ji) * pwcp(ji,1,jwoxy) / (pwcp(ji,1,jwoxy) + 20.E-6)
+      END DO
+
+      ! ------------------------------------------------------
+      !    Vertical variations of the bioturbation coefficient
+      ! ------------------------------------------------------
+      IF (ln_btbz) THEN
+         DO ji = 1, jpoce
+            db(ji,:) = db(ji,:) * exp( -(profsedw(:) / dbtbzsc)**2 ) / (365.0 * 86400.0)
+         END DO
+      ELSE
+         DO jk = 1, jpksed
+            IF (profsedw(jk) > dbtbzsc) THEN
+                db(:,jk) = 0.0
+            ENDIF
+         END DO
+      ENDIF
+
+      IF (ln_irrig) THEN
+         DO jk = 1, jpksed
+            DO ji = 1, jpoce
+               irrig(ji,jk) = ( 7.63752 - 7.4465 * exp( -0.89603 * zsumtot(ji) ) ) * pwcp(ji,1,jwoxy)  &
+               &             / (pwcp(ji,1,jwoxy) + 20.E-6) 
+               irrig(ji,jk) = irrig(ji,jk) * exp( -(profsedw(jk) / xirrzsc) )
+            END DO
+         END DO
+      ELSE
+         irrig(:,:) = 0.0
+      ENDIF
+
+      DEALLOCATE( zvolc )
+
+      IF( ln_timing )  CALL timing_stop('sed_dsr')
+!      
    END SUBROUTINE sed_dsr
-#else
-   !!======================================================================
-   !! MODULE seddsr  :   Dummy module 
-   !!======================================================================
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_dsr ( kt )
-     INTEGER, INTENT(in) :: kt
-     WRITE(*,*) 'sed_dsr: You should not have seen this print! error?', kt
-  END SUBROUTINE sed_dsr
-#endif
-   
+
+   SUBROUTINE sed_dsr_redoxb
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE sed_dsr_redox  ***
+      !! 
+      !!  ** Purpose :  computes secondary redox reactions
+      !!
+      !!  ** Methode :  It uses Strand splitter algorithm proposed by 
+      !!                Nguyen et al. (2009) and modified by Wang et al. (2018)
+      !!                Basically, each equation is solved analytically when 
+      !!                feasible, otherwise numerically at t+1/2. Then 
+      !!                the last equation is solved at t+1. The other equations
+      !!                are then solved at t+1 starting in the reverse order.
+      !!                Ideally, it's better to start from the fastest reaction
+      !!                to the slowest and then reverse the order to finish up
+      !!                with the fastest one. But random order works well also.
+      !!                The scheme is second order, positive and mass 
+      !!                conserving. It works well for stiff systems.
+      !! 
+      !!   History :
+      !!        !  18-08 (O. Aumont)  Original code
+      !!----------------------------------------------------------------------
+      !! Arguments
+      ! --- local variables
+      INTEGER   ::  ji, jk, jn   ! dummy looop indices
+
+      REAL, DIMENSION(6)  :: zsedtrn, zsedtra
+      REAL(wp)  ::  zalpha, zbeta, zgamma, zdelta, zepsi, zsedfer
+      !!
+      !!----------------------------------------------------------------------
+
+      IF( ln_timing )  CALL timing_start('sed_dsr_redoxb')
+
+      DO ji = 1, jpoce
+         DO jk = 2, jpksed
+            zsedtrn(1)  = pwcp(ji,jk,jwoxy)
+            zsedtrn(2)  = MAX(0., pwcp(ji,jk,jwh2s) )
+            zsedtrn(3)  = pwcp(ji,jk,jwnh4)
+            zsedtrn(4)  = MAX(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) )
+            zsedfer     = MIN(0., pwcp(ji,jk,jwfe2) - sedligand(ji,jk) )
+            zsedtrn(5)  = solcp(ji,jk,jsfeo) * zvolc(ji,jk,jsfeo)
+            zsedtrn(6)  = solcp(ji,jk,jsfes) * zvolc(ji,jk,jsfes)
+            zsedtra(:)  = zsedtrn(:) 
+
+            ! First pass of the scheme. At the end, it is 1st order 
+            ! -----------------------------------------------------
+            ! Fe + O2
+            zalpha = zsedtra(1) - 0.25 * zsedtra(4)
+            zbeta  = zsedtra(4) + zsedtra(5) 
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
+            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4)
+            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn )
+            zsedtra(1) = zalpha + 0.25 * zsedtra(4) 
+            zsedtra(5) = zbeta  - zsedtra(4)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
+            pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4)
+            ! H2S + O2
+            zalpha = zsedtra(1) - 2.0 * zsedtra(2)
+            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2)
+            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn )
+            zsedtra(1) = zalpha + 2.0 * zsedtra(2)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2)
+            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2)
+            ! NH4 + O2
+            zalpha = zsedtra(1) - 2.0 * zsedtra(3)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3)
+            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn )
+            zsedtra(1) = zalpha + 2.0 * zsedtra(3)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3)
+            ! FeS - O2
+            zalpha = zsedtra(1) - 2.0 * zsedtra(6)
+            zbeta  = zsedtra(4) + zsedtra(6)
+            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6)
+            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn )
+            zsedtra(1) = zalpha + 2.0 * zsedtra(6)
+            zsedtra(4) = zbeta  - zsedtra(6)
+            pwcp(ji,jk,jwso4) = zgamma - zsedtra(6)
+!            ! Fe - H2S
+            zalpha = zsedtra(2) - zsedtra(4)
+            zbeta  = zsedtra(4) + zsedtra(6)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
+            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn )
+            zsedtra(2) = zalpha + zsedtra(4)
+            zsedtra(6) = zbeta  - zsedtra(4)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
+            ! FEOH + H2S
+            zalpha = zsedtra(5) - 2.0 * zsedtra(2)
+            zbeta  = zsedtra(5) + zsedtra(4)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
+            zdelta = pwcp(ji,jk,jwso4) + zsedtra(2)
+            zepsi  = pwcp(ji,jk,jwpo4) + redfep * zsedtra(5)
+            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_feh2s * zalpha * dtsed2 ) - 1.0 )  &
+            &            + zalpha * exp( reac_feh2s * zalpha * dtsed2 ) + rtrn )
+            zsedtra(5) = zalpha + 2.0 * zsedtra(2)
+            zsedtra(4) = zbeta  - zsedtra(5)
+            pwcp(ji,jk,jwso4) = zdelta - zsedtra(2)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
+            pwcp(ji,jk,jwpo4) = zepsi - redfep * zsedtra(5)
+            ! Fe - H2S
+            zalpha = zsedtra(2) - zsedtra(4)
+            zbeta  = zsedtra(4) + zsedtra(6)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
+            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( zsedtra(4) * ( exp( reac_fes * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_fes * zalpha * dtsed2 /2. ) + rtrn )
+            zsedtra(2) = zalpha + zsedtra(4)
+            zsedtra(6) = zbeta  - zsedtra(4)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
+            ! FeS - O2
+            zalpha = zsedtra(1) - 2.0 * zsedtra(6)
+            zbeta  = zsedtra(4) + zsedtra(6)
+            zgamma = pwcp(ji,jk,jwso4) + zsedtra(6)
+            zsedtra(6) = ( zsedtra(6) * zalpha ) / ( 2.0 * zsedtra(6) * ( exp( reac_feso * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_feso * zalpha * dtsed2 /2. ) + rtrn )
+            zsedtra(1) = zalpha + 2.0 * zsedtra(6)
+            zsedtra(4) = zbeta  - zsedtra(6)
+            pwcp(ji,jk,jwso4) = zgamma - zsedtra(6)
+            ! NH4 + O2
+            zalpha = zsedtra(1) - 2.0 * zsedtra(3)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(3)
+            zsedtra(3) = ( zsedtra(3) * zalpha ) / ( 2.0 * zsedtra(3) * ( exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_nh4 * zadsnh4 * zalpha * dtsed2 /2. ) + rtrn )
+            zsedtra(1) = zalpha + 2.0 * zsedtra(3)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(3)
+            ! H2S + O2
+            zalpha = zsedtra(1) - 2.0 * zsedtra(2)
+            zbeta  = pwcp(ji,jk,jwso4) + zsedtra(2)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(2)
+            zsedtra(2) = ( zsedtra(2) * zalpha ) / ( 2.0 * zsedtra(2) * ( exp( reac_h2s * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_h2s * zalpha * dtsed2 / 2. ) + rtrn )
+            zsedtra(1) = zalpha + 2.0 * zsedtra(2)
+            pwcp(ji,jk,jwso4) = zbeta - zsedtra(2)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(2)
+            ! Fe + O2
+            zalpha = zsedtra(1) - 0.25 * zsedtra(4)
+            zbeta  = zsedtra(4) + zsedtra(5)
+            zgamma = pwcp(ji,jk,jwalk) - 2.0 * zsedtra(4)
+            zdelta = pwcp(ji,jk,jwpo4) - redfep * zsedtra(4)
+            zsedtra(4) = ( zsedtra(4) * zalpha ) / ( 0.25 * zsedtra(4) * ( exp( reac_fe2 * zalpha * dtsed2 / 2. ) - 1.0 )  &
+            &            + zalpha * exp( reac_fe2 * zalpha * dtsed2 / 2. ) + rtrn )
+            zsedtra(1) = zalpha + 0.25 * zsedtra(4)
+            zsedtra(5) = zbeta  - zsedtra(4)
+            pwcp(ji,jk,jwpo4) = zdelta + redfep * zsedtra(4)
+            pwcp(ji,jk,jwalk) = zgamma + 2.0 * zsedtra(4)
+            pwcp(ji,jk,jwoxy)  = zsedtra(1)
+            pwcp(ji,jk,jwh2s)  = zsedtra(2)
+            pwcp(ji,jk,jwnh4)  = zsedtra(3)
+            pwcp(ji,jk,jwfe2)  = zsedtra(4) + sedligand(ji,jk) + zsedfer
+            pwcp(ji,jk,jwno3)  = pwcp(ji,jk,jwno3)  + ( zsedtrn(3) - pwcp(ji,jk,jwnh4) )
+            solcp(ji,jk,jsfeo) = zsedtra(5) / zvolc(ji,jk,jsfeo)
+            solcp(ji,jk,jsfes) = zsedtra(6) / zvolc(ji,jk,jsfes)
+         END DO
+      END DO
+
+      IF( ln_timing )  CALL timing_stop('sed_dsr_redoxb')
+
+  END SUBROUTINE sed_dsr_redoxb
+
 END MODULE seddsr

NEMO/trunk/src/TOP/PISCES/SED/seddta.F90

@@ -4 +4 @@
    !! Sediment data  :  read sediment input data from a file
    !!=====================================================================
-#if defined key_sed
+
    !! * Modules used
    USE sed
    USE sedarr
    USE iom
+   USE lib_mpp         ! distribued memory computing library
 
    IMPLICIT NONE
@@ -17 +18 @@
 
    !! *  Module variables
-   REAL(wp), DIMENSION(:), ALLOCATABLE ::  smask       ! mask for sediments points
    REAL(wp) ::  rsecday  ! number of second per a day
-   REAL(wp) ::  conv1    ! [m/day]--->[cm/s]  
    REAL(wp) ::  conv2    ! [kg/m2/month]-->[g/cm2/s] ( 1 month has 30 days )
-
-   INTEGER ::  numbio
-
-#if defined key_sed_off
-   INTEGER ::  numoce
-#endif
 
    !! $Id$
@@ -54 +47 @@
       INTEGER  ::  ji, jj, js, jw, ikt
 
-      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zdta
-      REAL(wp), DIMENSION(:)  , ALLOCATABLE :: zdtap, zdtag
-
+      REAL(wp), DIMENSION(jpoce) :: zdtap, zdtag
+      REAL(wp), DIMENSION(jpi,jpj) :: zwsbio4, zwsbio3, zwscal
+      REAL(wp) :: zf0, zf1, zf2, zkapp, zratio, zdep
 
       !----------------------------------------------------------------------
@@ -64 +57 @@
       !-------------------------------------------------------------
 
-      WRITE(numsed,*)
-      WRITE(numsed,*) ' sed_dta : Bottom layer fields'
-      WRITE(numsed,*) ' ~~~~~~'
-      WRITE(numsed,*) ' Data from SMS model'
-      WRITE(numsed,*)
+      IF( ln_timing )  CALL timing_start('sed_dta')
+
+      IF (lwp) THEN
+         WRITE(numsed,*)
+         WRITE(numsed,*) ' sed_dta : Bottom layer fields'
+         WRITE(numsed,*) ' ~~~~~~'
+         WRITE(numsed,*) ' Data from SMS model'
+         WRITE(numsed,*)
+      ENDIF
 
 
       ! open file
       IF( kt == nitsed000 ) THEN
-         WRITE(numsed,*) ' sed_dta : Sediment fields'
-         CALL iom_open ( 'data_bio_bot'     , numbio )
-#if defined key_sed_off
-         CALL iom_open( 'data_oce_bot', numoce)
-#endif
+         IF (lwp) WRITE(numsed,*) ' sed_dta : Sediment fields'
+         dtsed = r2dttrc
          rsecday = 60.* 60. * 24.
-         conv1   = 1.0e+2 / rsecday 
-         conv2   = 1.0e+3 / ( 1.0e+4 * rsecday * 30. ) 
-
-         ! Compute sediment mask
-         ALLOCATE( zdta(jpi,jpj) ) 
+!         conv2   = 1.0e+3 / ( 1.0e+4 * rsecday * 30. )
+         conv2 = 1.0e+3 /  1.0e+4 
+         rdtsed(2:jpksed) = dtsed / ( denssol * por1(2:jpksed) )
+      ENDIF
+
+      ! Initialization of temporaries arrays  
+      zdtap(:)    = 0. 
+      zdtag(:)    = 0.  
+
+      ! reading variables
+      IF (lwp) WRITE(numsed,*)
+      IF (lwp) WRITE(numsed,*) ' sed_dta : Bottom layer fields at time  kt = ', kt
+      ! reading variables
+      !
+      !    Sinking speeds of detritus is increased with depth as shown
+      !    by data and from the coagulation theory
+      !    -----------------------------------------------------------
+      IF (ln_sediment_offline) THEN
          DO jj = 1, jpj
             DO ji = 1, jpi
-               ikt = MAX( INT( sbathy(ji,jj) ) - 1, 1 )
-               zdta(ji,jj) = tmask(ji,jj,ikt) 
-            ENDDO
-         ENDDO
-         ALLOCATE( smask(jpoce) )
-         smask(:) = 0.
-         CALL pack_arr( jpoce, smask(1:jpoce), zdta(1:jpi,1:jpj), iarroce(1:jpoce) )
-      ENDIF
-
-      ! Initialization of temporaries arrays  
-      ALLOCATE( zdtap(jpoce) )    ;   zdtap(:)    = 0. 
-      ALLOCATE( zdtag(jpoce) )    ;   zdtag(:)    = 0.  
-
-      IF( MOD( kt - 1, nfreq ) == 0 ) THEN
-         ! reading variables
-         WRITE(numsed,*)
-         WRITE(numsed,*) ' sed_dta : Bottom layer fields at time  kt = ', kt
-         ! reading variables
-         trc_data(:,:,:) = 0.
-#if ! defined key_sed_off
-         DO jj = 1,jpj
+               ikt = mbkt(ji,jj)
+               zwsbio4(ji,jj) = wsbio2 / rday
+               zwsbio3(ji,jj) = wsbio  / rday
+               zwscal (ji,jj) = wsbio2 / rday
+            END DO
+         END DO
+      ELSE
+         DO jj = 1, jpj
             DO ji = 1, jpi
                ikt = mbkt(ji,jj)
-               IF ( tmask(ji,jj,ikt) == 1 ) THEN
-                  trc_data(ji,jj,1)  = trn  (ji,jj,ikt,jptal)
-                  trc_data(ji,jj,2)  = trn  (ji,jj,ikt,jpdic)
-                  trc_data(ji,jj,3)  = trn  (ji,jj,ikt,jpno3) / 7.6
-                  trc_data(ji,jj,4)  = trn  (ji,jj,ikt,jppo4) / 122.
-                  trc_data(ji,jj,5)  = trn  (ji,jj,ikt,jpoxy)
-                  trc_data(ji,jj,6)  = trn  (ji,jj,ikt,jpsil)
-                  trc_data(ji,jj,7 ) = sinksil (ji,jj,ikt)
-                  trc_data(ji,jj,8 ) = sinking (ji,jj,ikt)
-                  trc_data(ji,jj,9 ) = sinking2(ji,jj,ikt)
-                  trc_data(ji,jj,10) = sinkcal (ji,jj,ikt)
-                  trc_data(ji,jj,11) = tsn     (ji,jj,ikt,jp_tem)
-                  trc_data(ji,jj,12) = tsn     (ji,jj,ikt,jp_sal)
-               ENDIF
-            ENDDO
+               zdep = e3t_n(ji,jj,ikt) / r2dttrc
+               zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) / rday )
+               zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) / rday )
+               zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) / rday )
+            END DO
+         END DO
+      ENDIF
+
+      trc_data(:,:,:) = 0.
+      DO jj = 1,jpj
+         DO ji = 1, jpi
+            ikt = mbkt(ji,jj)
+            IF ( tmask(ji,jj,ikt) == 1 ) THEN
+               trc_data(ji,jj,1)   = trb(ji,jj,ikt,jpsil)
+               trc_data(ji,jj,2)   = trb(ji,jj,ikt,jpoxy)
+               trc_data(ji,jj,3)   = trb(ji,jj,ikt,jpdic)
+               trc_data(ji,jj,4)   = trb(ji,jj,ikt,jpno3) / 7.625
+               trc_data(ji,jj,5)   = trb(ji,jj,ikt,jppo4) / 122.
+               trc_data(ji,jj,6)   = trb(ji,jj,ikt,jptal)
+               trc_data(ji,jj,7)   = trb(ji,jj,ikt,jpnh4) / 7.625
+               trc_data(ji,jj,8)   = 0.0
+               trc_data(ji,jj,9)   = 28.0E-3
+               trc_data(ji,jj,10)  = trb(ji,jj,ikt,jpfer)
+               trc_data(ji,jj,11 ) = MIN(trb(ji,jj,ikt,jpgsi), 1E-4) * zwsbio4(ji,jj) * 1E3
+               trc_data(ji,jj,12 ) = MIN(trb(ji,jj,ikt,jppoc), 1E-4) * zwsbio3(ji,jj) * 1E3
+               trc_data(ji,jj,13 ) = MIN(trb(ji,jj,ikt,jpgoc), 1E-4) * zwsbio4(ji,jj) * 1E3
+               trc_data(ji,jj,14)  = MIN(trb(ji,jj,ikt,jpcal), 1E-4) * zwscal (ji,jj) * 1E3
+               trc_data(ji,jj,15)  = tsn(ji,jj,ikt,jp_tem)
+               trc_data(ji,jj,16)  = tsn(ji,jj,ikt,jp_sal)
+               trc_data(ji,jj,17 ) = ( trb(ji,jj,ikt,jpsfe) * zwsbio3(ji,jj) + trb(ji,jj,ikt,jpbfe)  &
+               &                     * zwsbio4(ji,jj)  ) * 1E3 / ( trc_data(ji,jj,12 ) + trc_data(ji,jj,13 ) + rtrn )
+               trc_data(ji,jj,17 ) = MIN(1E-3, trc_data(ji,jj,17 ) )
+            ENDIF
          ENDDO
-
-#else
-         CALL iom_get( numbio, jpdom_data, 'ALKBOT'     , trc_data(:,:,1 ) )
-         CALL iom_get( numbio, jpdom_data, 'DICBOT'     , trc_data(:,:,2 ) )
-         CALL iom_get( numbio, jpdom_data, 'NO3BOT'     , trc_data(:,:,3 ) )
-         CALL iom_get( numbio, jpdom_data, 'PO4BOT'     , trc_data(:,:,4 ) )
-         CALL iom_get( numbio, jpdom_data, 'O2BOT'      , trc_data(:,:,5 ) )
-         CALL iom_get( numbio, jpdom_data, 'SIBOT'      , trc_data(:,:,6 ) )
-         CALL iom_get( numbio, jpdom_data, 'OPALFLXBOT' , trc_data(:,:,7 ) ) 
-         CALL iom_get( numbio, jpdom_data, 'POCFLXBOT'  , trc_data(:,:,8 ) ) 
-         CALL iom_get( numbio, jpdom_data, 'GOCFLXBOT'  , trc_data(:,:,9 ) ) 
-         CALL iom_get( numbio, jpdom_data, 'CACO3FLXBOT', trc_data(:,:,10) ) 
-         CALL iom_get( numoce, jpdom_data, 'TBOT'       , trc_data(:,:,11) ) 
-         CALL iom_get( numoce, jpdom_data, 'SBOT'       , trc_data(:,:,12) ) 
-#endif
-
-         ! Pore water initial concentration [mol/l] in  k=1
-         !-------------------------------------------------
-
-          ! Alkalinity ( 1 umol = 10-6equivalent )
-         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwalk), trc_data(1:jpi,1:jpj,1), iarroce(1:jpoce) )
-         ! DIC 
-         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwdic), trc_data(1:jpi,1:jpj,2), iarroce(1:jpoce) )
-         ! Nitrates (1 umol/l = 10-6 mol/l)
-         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwno3), trc_data(1:jpi,1:jpj,3), iarroce(1:jpoce) )
-         ! Phosphates (1 umol/l = 10-6 mol/l)
-         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwpo4), trc_data(1:jpi,1:jpj,4), iarroce(1:jpoce) )
-         ! Oxygen (1 umol/l = 10-6 mol/l)
-         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwoxy), trc_data(1:jpi,1:jpj,5), iarroce(1:jpoce) )        
-         ! Silicic Acid [mol.l-1]
-         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwsil), trc_data(1:jpi,1:jpj,6), iarroce(1:jpoce) )                  
-         ! DIC13 (mol/l)obtained from dc13 and DIC (12) and PDB 
-         CALL iom_get ( numbio,jpdom_data,'DC13',zdta(:,:) )
-         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jwc13), zdta(1:jpi,1:jpj), iarroce(1:jpoce) )
-         pwcp_dta(1:jpoce,jwc13) = pdb * ( pwcp_dta(1:jpoce,jwc13) * 1.0e-3 + 1.0 )  &
-            &                          *   pwcp_dta(1:jpoce,jwdic)         
-         
-         !  Solid components : 
-         !-----------------------
-         !  Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
-         CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,7), iarroce(1:jpoce) ) 
-         rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4
-         !  Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
-         CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,8) , iarroce(1:jpoce) )      
-         CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,9) , iarroce(1:jpoce) )
-         rainrm_dta(1:jpoce,jspoc) =   ( zdtap(1:jpoce) +  zdtag(1:jpoce) ) * 1e-4
-         !  Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
-         CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,10), iarroce(1:jpoce) )
-         rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4
-         ! vector temperature [°C] and salinity 
-         CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,11), iarroce(1:jpoce) )
-         CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,12), iarroce(1:jpoce) )
-        
-         ! Clay rain rate in [mol/(cm**2.s)] 
-         ! inputs data in [kg.m-2.mois-1] ---> 1e+3/(1e+4*60*24*60*60) [g.cm-2.s-1]   
-         ! divided after by molecular weight g.mol-1      
-         zdta(:,:)   = 0.
-         CALL iom_get( numbio, jpdom_data, 'CLAY', zdta(:,:) ) 
-         CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsclay) , zdta(1:jpi,1:jpj), iarroce(1:jpoce) )      
-         rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay)
-
-      ENDIF
+      ENDDO
+
+      ! Pore water initial concentration [mol/l] in  k=1
+      !-------------------------------------------------
+      DO jw = 1, jpwat
+         CALL pack_arr ( jpoce,  pwcp_dta(1:jpoce,jw), trc_data(1:jpi,1:jpj,jw), iarroce(1:jpoce) )
+      END DO
+      !  Solid components : 
+      !-----------------------
+      !  Sinking fluxes for OPAL in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
+      CALL pack_arr ( jpoce, rainrm_dta(1:jpoce,jsopal), trc_data(1:jpi,1:jpj,11), iarroce(1:jpoce) ) 
+      rainrm_dta(1:jpoce,jsopal) = rainrm_dta(1:jpoce,jsopal) * 1e-4
+      !  Sinking fluxes for POC in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
+      CALL pack_arr ( jpoce, zdtap(1:jpoce), trc_data(1:jpi,1:jpj,12) , iarroce(1:jpoce) )      
+      CALL pack_arr ( jpoce, zdtag(1:jpoce), trc_data(1:jpi,1:jpj,13) , iarroce(1:jpoce) )
+      DO ji = 1, jpoce
+!        zkapp  = MIN( (1.0 - 0.02 ) * reac_poc, 3731.0 * max(100.0, zkbot(ji) )**(-1.011) / ( 365.0 * 24.0 * 3600.0 ) )
+!        zkapp   = MIN( 0.98 * reac_poc, 100.0 * max(100.0, zkbot(ji) )**(-0.6) / ( 365.0 * 24.0 * 3600.0 ) )
+!        zratio = ( ( 1.0 - 0.02 ) * reac_poc + 0.02 * reac_poc * 0. - zkapp) / ( ( 0.02 - 1.0 ) * reac_poc / 100. - 0.02 * reac_poc * 0. + zkapp )
+!        zf1    = ( 0.02 * (reac_poc - reac_poc * 0.) + zkapp - reac_poc ) / ( reac_poc / 100. - reac_poc )
+!        zf1    = MIN(0.98, MAX(0., zf1 ) )
+         zf1    = 0.48
+         zf0    = 1.0 - 0.02 - zf1
+         zf2    = 0.02
+         rainrm_dta(ji,jspoc) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf0
+         rainrm_dta(ji,jspos) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf1
+         rainrm_dta(ji,jspor) =   ( zdtap(ji) +  zdtag(ji) ) * 1e-4 * zf2
+      END DO
+      !  Sinking fluxes for Calcite in mol.m-2.s-1 ; conversion in mol.cm-2.s-1
+      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jscal), trc_data(1:jpi,1:jpj,14), iarroce(1:jpoce) )
+      rainrm_dta(1:jpoce,jscal) = rainrm_dta(1:jpoce,jscal) * 1e-4
+      ! vector temperature [°C] and salinity 
+      CALL pack_arr ( jpoce,  temp(1:jpoce), trc_data(1:jpi,1:jpj,15), iarroce(1:jpoce) )
+      CALL pack_arr ( jpoce,  salt(1:jpoce), trc_data(1:jpi,1:jpj,16), iarroce(1:jpoce) )
+      
+      ! Clay rain rate in [mol/(cm**2.s)] 
+      ! inputs data in [kg.m-2.sec-1] ---> 1e+3/(1e+4) [g.cm-2.s-1]   
+      ! divided after by molecular weight g.mol-1      
+      CALL pack_arr ( jpoce,  rainrm_dta(1:jpoce,jsclay), dust(1:jpi,1:jpj), iarroce(1:jpoce) )
+      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * conv2 / mol_wgt(jsclay)   &
+      &                            + wacc(1:jpoce) * por1(2) * denssol / mol_wgt(jsclay) / ( rsecday * 365.0 )
+      rainrm_dta(1:jpoce,jsclay) = rainrm_dta(1:jpoce,jsclay) * 0.965
+      rainrm_dta(1:jpoce,jsfeo)  = rainrm_dta(1:jpoce,jsclay) * mol_wgt(jsclay) / mol_wgt(jsfeo) * 0.035 / 0.965
+!    rainrm_dta(1:jpoce,jsclay) = 1.0E-4 * conv2 / mol_wgt(jsclay)
+
+      ! Iron monosulphide rain rates. Set to 0
+      rainrm_dta(1:jpoce,jsfes)  = 0. 
+
+      ! Fe/C ratio in sinking particles that fall to the sediments
+      CALL pack_arr ( jpoce,  fecratio(1:jpoce), trc_data(1:jpi,1:jpj,17), iarroce(1:jpoce) )
+
+      sedligand(:,1) = 1.E-9
 
       ! sediment pore water at 1st layer (k=1)
       DO jw = 1, jpwat
-         pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw) * smask(1:jpoce)
+         pwcp(1:jpoce,1,jw) = pwcp_dta(1:jpoce,jw)
       ENDDO
 
       !  rain
       DO js = 1, jpsol
-         rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js) * smask(1:jpoce)
+         rainrm(1:jpoce,js) = rainrm_dta(1:jpoce,js)
       ENDDO
 
@@ -212 +214 @@
       dzdep(1:jpoce) = raintg(1:jpoce) * rdtsed(2) 
 
-
-      DEALLOCATE( zdta ) 
-      DEALLOCATE( zdtap    ) ;  DEALLOCATE( zdtag    ) 
-
-      IF( kt == nitsedend )   THEN
-         CALL iom_close ( numbio )
-#if defined key_sed_off
-         CALL iom_close ( numoce )
-#endif
-      ENDIF
+      IF( lk_iomput ) THEN
+          IF( iom_use("sflxclay" ) ) CALL iom_put( "sflxclay", dust(:,:) * conv2 * 1E4 )
+          IF( iom_use("sflxcal" ) )  CALL iom_put( "sflxcal", trc_data(:,:,13) )
+          IF( iom_use("sflxbsi" ) )  CALL iom_put( "sflxbsi", trc_data(:,:,10) )
+          IF( iom_use("sflxpoc" ) )  CALL iom_put( "sflxpoc", trc_data(:,:,11) + trc_data(:,:,12) )
+      ENDIF
+
+      IF( ln_timing )  CALL timing_stop('sed_dta')
       
    END SUBROUTINE sed_dta
 
-#else
-   !!======================================================================
-   !! MODULE seddta  :   Dummy module 
-   !!======================================================================
-   !! $Id$
-CONTAINS
-   SUBROUTINE sed_dta ( kt )
-     INTEGER, INTENT(in) :: kt
-     WRITE(*,*) 'sed_stp: You should not have seen this print! error?', kt 
-  END SUBROUTINE sed_dta
-#endif
-
 END MODULE seddta

NEMO/trunk/src/TOP/PISCES/SED/sedini.F90

@@ -4 +4 @@
    !! Sediment : define sediment variables
    !!=====================================================================
-#if defined key_sed
 
    !!----------------------------------------------------------------------
@@ -11 +10 @@
    !! * Modules used
    USE sed     ! sediment global variable
-   USE seddta
-   USE sedrst
-   USE sedco3
-   USE sedchem
+   USE sed_oce
    USE sedarr
+   USE sedadv
+   USE trcdmp_sed
+   USE trcdta
    USE iom
+   USE lib_mpp         ! distribued memory computing library
 
 
@@ -24 +24 @@
    !! 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 +100 @@
       !!        !  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 +109 @@
       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 +137 @@
             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 +156 @@
       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 +181 @@
       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 +197 @@
       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 +237 @@
       !! * 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 +256 @@
       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 +273 @@
       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 +279 @@
       ! (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 +305 @@
       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 +333 @@
       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 +349 @@
          &              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 +359 @@
       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 +368 @@
       ! 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 +405 @@
       !!----------------------------------------------------------------------
 
-      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 +421 @@
       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 +515 @@
       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 +547 @@
       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

NEMO/trunk/src/TOP/PISCES/SED/sedmat.F90

@@ -1 +1 @@
 MODULE sedmat
-#if defined key_sed
    !!======================================================================
    !!              ***  MODULE  sedmat  ***
@@ -9 +8 @@
 
    USE sed     ! sediment global variable
+   USE lib_mpp         ! distribued memory computing library
+
 
    IMPLICIT NONE
@@ -25 +26 @@
  CONTAINS
 
-    SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psol )
+    SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in )
        !!---------------------------------------------------------------------
        !!                  ***  ROUTINE sed_mat_dsr  ***
@@ -48 +49 @@
        !!----------------------------------------------------------------------
        !! * Arguments
-       INTEGER , INTENT(in) ::  nvar  ! number of variables
+       INTEGER , INTENT(in) ::  nvar  ! number of variable
        INTEGER , INTENT(in) ::  ndim  ! number of points
        INTEGER , INTENT(in) ::  nlev  ! number of sediment levels
 
-       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(in   ) :: preac  ! reaction rates
-       REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: psol   ! solution ( undersaturation values )
+       REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: preac  ! reaction rates
+       REAL(wp), DIMENSION(ndim,nlev), INTENT(in   ) :: psms  ! reaction rates
+       REAL(wp), DIMENSION(ndim,nlev), INTENT(inout) :: psol   ! solution ( undersaturation values )
+       REAL(wp), INTENT(in) ::  dtsed_in
  
        !---Local declarations
@@ -67 +70 @@
        !----------------------------------------------------------------------
 
+       IF( ln_timing )  CALL timing_start('sed_mat_dsr')
 
        ! Computation left hand side of linear system of 
@@ -73 +77 @@
 
 
+       jn = nvar
        ! first sediment level          
        DO ji = 1, ndim
-          aplus  = ( ( volw3d(ji,1) / dz3d(ji,1) ) + &
-                     ( volw3d(ji,2) / dz3d(ji,2) ) ) / 2.
+          aplus  = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + &
+                        ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2.
           dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2.
-          rplus  = ( dtsed / volw3d(ji,1) ) * diff(1) * aplus / dxplus 
+          rplus  = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 
 
           za(ji,1) = 0.
@@ -87 +92 @@
        DO jk = 2, nlev - 1
           DO ji = 1, ndim
-             aminus  = ( ( volw3d(ji,jk-1) / dz3d(ji,jk-1) ) + &
-                &        ( volw3d(ji,jk  ) / dz3d(ji,jk  ) ) ) / 2.
+             aminus  = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + &
+             &        ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) ) / 2.
              dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2.
 
-             aplus   = ( ( volw3d(ji,jk  ) / dz3d(ji,jk  ) ) + &
-                &        ( volw3d(ji,jk+1) / dz3d(ji,jk+1) ) ) / 2.
+             aplus   = ( ( volw3d(ji,jk  ) / ( dz3d(ji,jk  ) ) ) + &
+             &        ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2.
              dxplus  = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2
-             !
-             rminus  = ( dtsed / volw3d(ji,jk) ) * diff(jk-1) * aminus / dxminus
-             rplus   = ( dtsed / volw3d(ji,jk) ) * diff(jk)   * aplus / dxplus
-             !     
+                !
+             rminus  = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus
+             rplus   = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn)   * aplus / dxplus
+                !     
              za(ji,jk) = -rminus
              zb(ji,jk) = 1. + rminus + rplus 
@@ -105 +110 @@
 
        DO ji = 1, ndim
-          aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1) ) + &
-             &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2.
+          aminus  = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1)  ) + &
+          &        ( volw3d(ji,nlev)  / dz3d(ji,nlev) ) ) / 2.
           dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2.
-          rminus  = ( dtsed / volw3d(ji,nlev) ) * diff(nlev-1) * aminus / dxminus
+          rminus  = ( dtsed_in / volw3d(ji,nlev) ) * diff(ji,nlev-1,jn) * aminus / dxminus
           !
           za(ji,nlev) = -rminus
@@ -119 +124 @@
        ! -----------------------------------------------
 
-       DO jn = 1, nvar
-
-          zr  (:,:)    = psol(:,:,jn)
-          zbet(:  )    = zb(:,1) + preac(:,1,jn)
-          psol(:,1,jn) = zr(:,1) / zbet(:)
+       zr  (:,:) = psol(:,:) + psms(:,:)
+       zb  (:,:) = zb(:,:) + preac(:,:)
+       zbet(:  ) = zb(:,1)
+       psol(:,1) = zr(:,1) / zbet(:)
+
           ! 
-          DO jk = 2, nlev
-             DO ji = 1, ndim
-                zgamm(ji,jk)   =  zc(ji,jk-1) / zbet(ji)
-                zbet(ji)       = ( zb(ji,jk) + preac(ji,jk,jn) ) - za(ji,jk) * zgamm(ji,jk)
-                psol(ji,jk,jn) = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1,jn) ) / zbet(ji)