Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 4520)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 4521)
@@ -233,5 +233,5 @@
                ! --------------------------------------------------------------------------------------
                IF( ln_fechem ) THEN
-                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) ) * 1E-9
+                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) + zFeP(ji,jj,jk) ) * 1E-9
                   zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
                ELSE
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 4520)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 4521)
@@ -203,5 +203,5 @@
 
                zmortz2 = ztortz2 + zrespz2
-               zmortzgoc = unass2 / ( 1. - zepsherv ) * ztortz2
+               zmortzgoc = unass2 / ( 1. - epsher2 ) * ztortz2 + zrespz2
                tra(ji,jj,jk,jpmes) = tra(ji,jj,jk,jpmes) - zmortz2 + zepsherv * zgraztot 
                tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazd
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 4520)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 4521)
@@ -184,5 +184,6 @@
                tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
 #if defined key_kriest
-               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + ( zmortz - zgrazm ) * xkr_dmicro
+               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zmortz * xkr_dmicro &
+                                                         - zgrazm * trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
 #endif
             END DO
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 4520)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 4521)
@@ -479,5 +479,5 @@
          CALL iom_close( numhydro )
          !
-         hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / fse3t(:,:,:)
+         hydrofe(:,:,:) = ( hydrofe(:,:,:) * hratio ) / ( cvol(:,:,:) * ryyss + rtrn ) / 1000._wp
          !
       ENDIF
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 4520)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 4521)
@@ -70,5 +70,5 @@
       REAL(wp) ::  zlim, zfact, zfactcal
       REAL(wp) ::  zo2, zno3, zflx, zpdenit, z1pdenit, zdenitt, zolimit
-      REAL(wp) ::  zsiloss, zcaloss, zwsbio3, zwsbio4, zwscal, zdep, zwstpoc
+      REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep, zwstpoc
       REAL(wp) ::  ztrfer, ztrpo4, zwdust
       REAL(wp) ::  zrdenittot, zsdenittot, znitrpottot
@@ -77,4 +77,5 @@
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3, zwork4
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zdenit2d, zironice
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsbio3, zwsbio4, zwscal
       REAL(wp), POINTER, DIMENSION(:,:,:) :: znitrpot, zirondep
       !!---------------------------------------------------------------------
@@ -93,5 +94,7 @@
       ! Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zwork4 )
+      CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal )
       CALL wrk_alloc( jpi, jpj, jpk, znitrpot )
+
       zdenit2d(:,:) = 0.e0
 
@@ -196,4 +199,18 @@
       ENDIF
 
+
+      ! OA: Warning, the following part is necessary, especially with Kriest
+      ! to avoid CFL problems above the sediments
+      ! --------------------------------------------------------------------
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ikt  = mbkt(ji,jj)
+            zdep = fse3t(ji,jj,ikt) / xstep
+            zwsbio4(ji,jj) = MIN( 0.99 * zdep, wsbio4(ji,jj,ikt) )
+            zwscal (ji,jj) = MIN( 0.99 * zdep, wscal (ji,jj,ikt) )
+            zwsbio3(ji,jj) = MIN( 0.99 * zdep, wsbio3(ji,jj,ikt) )
+         END DO
+      END DO
+
 #if ! defined key_sed
       ! Computation of the sediment denitrification proportion: The metamodel from midlleburg (2006) is being used
@@ -204,8 +221,8 @@
               ikt = mbkt(ji,jj)
 # if defined key_kriest
-              zflx =    trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)    * 1E3 * 1E6 / 1E4
+              zflx =    trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)    * 1E3 * 1E6 / 1E4
 # else
-              zflx = (  trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt)   &
-                &     + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) )  * 1E3 * 1E6 / 1E4
+              zflx = (  trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
+                &     + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4
 #endif
               zflx  = LOG10( MAX( 1E-3, zflx ) )
@@ -227,14 +244,14 @@
             ikt = mbkt(ji,jj) 
 # if defined key_kriest
-            zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * wscal (ji,jj,ikt)
-            zwork2(ji,jj) = trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt)
+            zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * zwscal (ji,jj)
+            zwork2(ji,jj) = trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)
 # else
-            zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * wsbio4(ji,jj,ikt)
-            zwork2(ji,jj) = trn(ji,jj,ikt,jpgoc) * wsbio4(ji,jj,ikt) + trn(ji,jj,ikt,jppoc) * wsbio3(ji,jj,ikt) 
+            zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * zwsbio4(ji,jj)
+            zwork2(ji,jj) = trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj) + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) 
 # endif
             ! For calcite, burial efficiency is made a function of saturation
             zfactcal      = MIN( excess(ji,jj,ikt), 0.2 )
             zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
-            zwork3(ji,jj) = trn(ji,jj,ikt,jpcal) * wscal (ji,jj,ikt) * 2.e0 * zfactcal
+            zwork3(ji,jj) = trn(ji,jj,ikt,jpcal) * zwscal(ji,jj) * 2.e0 * zfactcal
          END DO
       END DO
@@ -256,12 +273,12 @@
             ikt  = mbkt(ji,jj)
             zdep = xstep / fse3t(ji,jj,ikt)
-            zwsbio4 = wsbio4(ji,jj,ikt) * zdep
-            zwscal  = wscal (ji,jj,ikt) * zdep
+            zws4 = zwsbio4(ji,jj) * zdep
+            zwsc = zwscal (ji,jj) * zdep
 # if defined key_kriest
-            zsiloss = trn(ji,jj,ikt,jpgsi) * zwsbio4
+            zsiloss = trn(ji,jj,ikt,jpgsi) * zws4
 # else
-            zsiloss = trn(ji,jj,ikt,jpgsi) * zwscal
+            zsiloss = trn(ji,jj,ikt,jpgsi) * zwsc
 # endif
-            zcaloss = trn(ji,jj,ikt,jpcal) * zwscal
+            zcaloss = trn(ji,jj,ikt,jpcal) * zwsc
             !
             trn(ji,jj,ikt,jpgsi) = trn(ji,jj,ikt,jpgsi) - zsiloss
@@ -282,17 +299,17 @@
             ikt     = mbkt(ji,jj)
             zdep    = xstep / fse3t(ji,jj,ikt)
-            zwsbio4 = wsbio4(ji,jj,ikt) * zdep
-            zwsbio3 = wsbio3(ji,jj,ikt) * zdep
+            zws4 = zwsbio4(ji,jj) * zdep
+            zws3 = zwsbio3(ji,jj) * zdep
 # if ! defined key_kriest
-            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - trn(ji,jj,ikt,jpgoc) * zwsbio4
-            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3
-            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * zwsbio4
-            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3
-            zwstpoc =  trn(ji,jj,ikt,jpgoc) * zwsbio4 + trn(ji,jj,ikt,jppoc) * zwsbio3 
+            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - trn(ji,jj,ikt,jpgoc) * zws4
+            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zws3
+            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * zws4
+            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zws3
+            zwstpoc              =  trn(ji,jj,ikt,jpgoc) * zws4 + trn(ji,jj,ikt,jppoc) * zws3 
 # else
-            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) - trn(ji,jj,ikt,jpnum) * zwsbio4
-            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zwsbio3
-            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zwsbio3
-            zwstpoc = trn(ji,jj,ikt,jppoc) * zwsbio3 
+            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) - trn(ji,jj,ikt,jpnum) * zws4
+            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zws3
+            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zws3
+            zwstpoc = trn(ji,jj,ikt,jppoc) * zws3 
 # endif
 
@@ -383,4 +400,5 @@
       !
       CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zwork4 )
+      CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal )
       CALL wrk_dealloc( jpi, jpj, jpk, znitrpot )
       !
Index: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 4520)
+++ branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 4521)
@@ -121,5 +121,5 @@
       rdenit  = 105._wp /  16._wp
       rdenita =   3._wp /  5._wp
-      o2ut    = 131._wp / 122._wp
+      o2ut    = 133._wp / 122._wp
 
       ! Initialization of tracer concentration in case of  no restart 
