Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zagg.F90	(revision 7180)
@@ -94,5 +94,5 @@
                zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
                ! tranfer of DOC to POC due to brownian motion
-               zaggdoc3 =  ( 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
+               zaggdoc3 =  114. * 0.3 * trb(ji,jj,jk,jpdoc) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
 
                !  Update the trends
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 7180)
@@ -445,9 +445,9 @@
                ! Liu and Millero (1999) only valid 5 - 50 degC
                ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16
-               fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254.0/ztkel1))
-               fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320.0/ztkel1) )
-               fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(zis**0.5)) - (1996.0/ztkel1) )
-               fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086.0/ztkel1) )
-               fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980.0/ztkel1) )
+               fesol(ji,jj,jk,1) = 10**(-13.486 - 0.1856*zis**0.5 + 0.3073*zis + 5254.0/ztkel1)
+               fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139*zis - 1320.0/ztkel1 )
+               fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 )
+               fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 )
+               fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**0.5 - 7980.0/ztkel1 )
             END DO
          END DO
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 7180)
@@ -69,5 +69,5 @@
       INTEGER  ::   ji, jj, jk, jic, jn
       REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
-      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol
+      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll, fe3sol, fe3sol1
       REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
       REAL(wp) ::   ztrc, zdust
@@ -279,16 +279,20 @@
                   zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
                ELSE
+                  zfeequi = zFe3(ji,jj,jk) * 1E-9
                   IF (ln_fecolloid) THEN
-                     zfeequi = zFe3(ji,jj,jk) * 1E-9
                      zhplus   = max( rtrn, hi(ji,jj,jk) )
-                     fe3sol  = fesol(ji,jj,jk,1) * ( fesol(ji,jj,jk,2) * zhplus**2  &
+                     fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
                      &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
                      &         + fesol(ji,jj,jk,5) / zhplus )
                      zfecoll = max( ( 0.1 * zFeL1(ji,jj,jk) * 1E-9 ), ( zFeL1(ji,jj,jk) * 1E-9 -fe3sol ) )
+#if defined key_ligand
+                     zligco  = max( ( 0.1 * trn(ji,jj,jk,jplgw) ), ( trn(ji,jj,jk,jplgw) - fe3sol ) )
+#endif 
                   ELSE
-                     zfeequi = zFe3(ji,jj,jk) * 1E-9 
                      zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
+#if defined key_ligand
+                     zligco  = 0.5 * trn(ji,jj,jk,jplgw)
+#endif 
                      fe3sol  = 0.
-                     kfep    = 0.
                   ENDIF
                ENDIF
@@ -324,29 +328,13 @@
                !  ----------------------------------------------------------------
                zlam1a  = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
-                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) + 0. * trb(ji,jj,jk,jppoc) )
+                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) )
                zaggdfea = zlam1a * zstep * zfecoll
-#if defined key_ligand
-               zligco = max( ( 0.1 * trn(ji,jj,jk,jplgw) ), ( trn(ji,jj,jk,jplgw) - fe3sol ) )
-               zaggliga = zlam1a * zstep * zligco
-#endif
                !
 #if defined key_kriest
                zaggdfeb = 0.
-               !
-#  if defined key_ligand
-               zaggligb = 0.
-#  endif
-               !
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfea - zaggdfeb - zcoag
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea + zaggdfeb
-               !
 #else
-               !
                zlam1b = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
                zaggdfeb = zlam1b * zstep * zfecoll
-               !
-#  if defined key_ligand
-               zaggligb = zlam1b * zstep * zligco
-#  endif
+#endif
                ! precipitation of Fe3+, creation of nanoparticles
                precip(ji,jj,jk) = max( 0., (zfeequi - fe3sol) ) * kfep * zstep
@@ -355,7 +343,14 @@
                &                     - zcoag - precip(ji,jj,jk)
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea
+#if ! defined key_kriest 
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
 #endif
 #if defined key_ligand
+               zaggliga = zlam1a * zstep * zligco
+#   if defined key_kriest
+               zaggligb = 0.
+#   else
+               zaggligb = zlam1b * zstep * zligco
+#   endif
                tra(ji,jj,jk,jpfep) = tra(ji,jj,jk,jpfep) + precip(ji,jj,jk)
                tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) - zaggliga - zaggligb
@@ -404,6 +399,9 @@
       ENDIF
       !
-                       CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
-      IF( ln_fechem )  CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
+      CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zFeL1, zTL1, ztotlig, precip )
+      IF( ln_fechem )  THEN
+         CALL wrk_dealloc( jpi, jpj,      zstrn, zstrn2 )
+         CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2, zTL2, zFeP )
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_fechem')
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 7180)
@@ -75,9 +75,8 @@
                ! ---------------------------------------------------------
                zstep   = xstep
-               zstep2  = zstep / 365 ! per year
 # if defined key_degrad
                zstep = zstep * facvol(ji,jj,jk)
-               zstep2 = zstep2 * facvol(ji,jj,jk)
 # endif
+               zstep2  = zstep / 365. ! per year
                ! production from remineralisation of organic matter
                zlgwp  = orem(ji,jj,jk) * rlig
@@ -88,5 +87,4 @@
                zlgwr = 1. / zlgwr * tgfunc(ji,jj,jk) * zstep2 * trn(ji,jj,jk,jplgw)
                ! photochem loss of weak ligand
-               zlablgw = MAX( 0., trn(ji,jj,jk, jpfer) * plig(ji,jj,jk) )
                zlgwpr = prlgw * zstep * etot(ji,jj,jk) * trn(ji,jj,jk,jplgw) * (1. - fr_i(ji,jj))
                tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zlgwp - zlgwr - zlgwpr
@@ -105,5 +103,5 @@
                ! lower in the aphotici zone
                ! ! 25 Wm-2 constant
-               zrfepa = rfep * (1- EXP(-1*etot(ji,jj,jk) / 25. ) ) * (1.- fr_i(ji,jj))
+               zrfepa = rfep * (1- EXP(-etot(ji,jj,jk) / 25. ) ) * (1.- fr_i(ji,jj))
                zrfepa = MAX( (zrfepa / 10.0), zrfepa ) ! min of 10 days lifetime
                zfepr = rfep * zstep * trn(ji,jj,jk,jpfep)
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 7180)
@@ -187,5 +187,5 @@
                   &   / ( concnno3 * concnnh4 + concnnh4 * trb(ji,jj,jk,jpno3) + concnno3 * trb(ji,jj,jk,jpnh4) ) 
                zlim2  = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnnh4 )
-               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  1.E-10   )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  1.E-10  )
                ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
                ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 7180)
@@ -85,5 +85,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-9 ), 0.e0 )
+               zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
                zstep    = xstep
 # if defined key_degrad
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90	(revision 7180)
@@ -64,6 +64,5 @@
       INTEGER  ::   ji, jj, jk, jn
       REAL(wp) ::   zremip, zremig, zdep, zorem, zorem2, zofer
-      REAL(wp) ::   zsizek, zsizek1, alphat, remint
-      REAL(wp) ::   solgoc, zpoc
+      REAL(wp) ::   zsizek, zsizek1, alphat, remint, solgoc, zpoc
 #if ! defined key_kriest
       REAL(wp) ::   zofer2, zofer3
@@ -71,6 +70,6 @@
       REAL(wp) ::   zstep, zrfact2
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: totprod, totthick, totcons 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zremipoc, zremigoc, zorem3
+      REAL(wp), POINTER, DIMENSION(:,:  )   :: totprod, totthick, totcons 
+      REAL(wp), POINTER, DIMENSION(:,:,:)   :: zremipoc, zremigoc, zorem3
       REAL(wp), POINTER, DIMENSION(:,:,:,:) :: alphag
       !!---------------------------------------------------------------------
@@ -79,5 +78,5 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      totprod, totthick, totcons )
+      CALL wrk_alloc( jpi, jpj,      totprod,  totthick, totcons )
       CALL wrk_alloc( jpi, jpj, jpk, zremipoc, zremigoc, zorem3 )
       ALLOCATE( alphag(jpi,jpj,jpk,jcpoc) )
@@ -124,4 +123,7 @@
                   alphat = 0.
                   remint = 0.
+                  !
+                  zsizek1  = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
+                  zsizek = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
                   !
                   IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
@@ -147,14 +149,11 @@
                        ! -----------------------------------------------------
                        !
-                       zsizek  = zdep / (wsbio2 + rtrn) * tgfunc(ji,jj,jk-1)
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
                        ! the concentration of each lability class is calculated
                        ! as the sum of the different sources and sinks
                        ! Please note that production of new GOC experiences
                        ! degradation 
-                       alphag(ji,jj,jk,jn) = alphan(jn) / (reminp(jn) * tgfunc(ji,jj,jk-1) )       &
-                       &   * (1. - exp( -reminp(jn) * zsizek ) ) * exp( -reminp(jn) * zsizek1 )     &
-                       &   * zpoc + prodgoc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) )        &
-                       &   * rday / rfact2 * alphan(jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
+                       &   + prodgoc(ji,jj,jk) * alphan(jn) / tgfunc(ji,jj,jk) / reminp(jn)             &
+                       &   * ( 1. - exp( -reminp(jn) * zsizek ) ) * rday / rfact2 
                        alphat = alphat + alphag(ji,jj,jk,jn)
                        remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
@@ -172,12 +171,8 @@
                     !
                     DO jn = 1, jcpoc
-                       zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn)                &
-                       &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1)  &
-                       &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek ) )                 &
-                       &   * exp( -reminp(jn) * zsizek1 ) + prodgoc(ji,jj,jk) / tgfunc(ji,jj,jk)       &
-                       &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) ) * rday / rfact2 * alphan(jn)     &
-                       &   / reminp(jn)
+                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * ( zsizek              &
+                       &   + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1.           &
+                       &   - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
+                       &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn) 
                        alphat = alphat + alphag(ji,jj,jk,jn)
                        remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
@@ -208,4 +203,5 @@
                zremig = zremigoc(ji,jj,jk) * zstep * tgfunc(ji,jj,jk)
                zorem2  = zremig * trb(ji,jj,jk,jpgoc)
+               orem(ji,jj,jk)      = zorem2
                zorem3(ji,jj,jk) = zremig * solgoc * trb(ji,jj,jk,jpgoc)
                zofer2 = zremig * trb(ji,jj,jk,jpbfe)
@@ -220,5 +216,4 @@
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zofer2 - zofer3
                tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem2
-               orem(ji,jj,jk)      = zorem2
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer2
             END DO
@@ -249,13 +244,11 @@
         DO jj = 1, jpj
            DO ji = 1, jpi
-              IF (tmask(ji,jj,jk) == 1.) THEN
-                zdep = hmld(ji,jj)
-                IF( fsdept(ji,jj,jk) <= zdep ) THEN
-                  totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2
-                  ! The temperature effect is included here
-                  totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk)
-                  totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2    &
-                  &                / ( trb(ji,jj,jk,jppoc) + rtrn )
-                ENDIF
+              zdep = hmld(ji,jj)
+              IF (tmask(ji,jj,jk) == 1. .AND. fsdept(ji,jj,jk) <= zdep ) THEN
+                totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2
+                ! The temperature effect is included here
+                totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk)
+                totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2    &
+                &                / ( trb(ji,jj,jk,jppoc) + rtrn )
               ENDIF
            END DO
@@ -310,4 +303,8 @@
                   remint = 0.
                   !
+                  ! the scale factors are corrected with temperature
+                  zsizek1  = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
+                  zsizek = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
+                  !
                   ! Special treatment of the level just below the MXL
                   ! See the comments in the GOC section
@@ -323,11 +320,10 @@
                     ! 
                     DO jn = 1, jcpoc
-                       ! the scale factor is corrected with temperature
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
                        ! computation of the lability spectrum applying the 
                        ! different sources and sinks
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) + zorem3(ji,jj,jk)                 &
-                       &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday / rfact2                   &
-                       &   * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc  &
+                       &   + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) * alphag(ji,jj,jk,jn) ) &
+                       &   / tgfunc(ji,jj,jk) / reminp(jn) * rday / rfact2 * ( 1. - exp( -reminp(jn)     &
+                       &   * zsizek ) )
                        alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
@@ -346,17 +342,11 @@
                     !
                     DO jn = 1, jcpoc
-                       zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)             &
-                       &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1)  &
-                       &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek ) )              &
-                       &   * exp( -reminp(jn) * zsizek1 ) + prodpoc(ji,jj,jk) / tgfunc(ji,jj,jk)    &
-                       &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) ) * rday / rfact2 * alphan(jn)  &
-                       &   / reminp(jn)
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) + zorem3(ji,jj,jk-1)              &
-                       &  * alphag(ji,jj,jk-1,jn) / tgfunc(ji,jj,jk-1) * rday / rfact2 * ( 1.       &
-                       &  - exp( -reminp(jn) * zsizek ) ) * exp( -reminp(jn) * zsizek1 )                  &
-                       &  + zorem3(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday         &
-                       &  / rfact2 * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)                       &
+                       &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) * alphan(jn)             & 
+                       &   + zorem3(ji,jj,jk-1) * alphag(ji,jj,jk-1,jn) ) * rday / rfact2 / reminp(jn)      &
+                       &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn)  &
+                       &   * zsizek ) + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk)                 &
+                       &   * alphag(ji,jj,jk,jn) ) * rday / rfact2 / reminp(jn) / tgfunc(ji,jj,jk) * ( 1.   &
+                       &   - exp( -reminp(jn) * zsizek ) )
                        alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 7180)
@@ -81,8 +81,9 @@
       REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld
       REAL(wp) ::   zmaxday, zdocprod, zpislopen, zpisloped
-      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
+      REAL(wp) ::   zmxltst, zmxlday
+      REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup, chlcnm_n, chlcdm_n
       REAL(wp) ::   zfact
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zw2d
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zw2d, zmixnano, zmixdiat
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopeadn, zpislopeadd, zprdia, zprbio, zprdch, zprnch, zysopt, zw3d   
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorcan, zprorcad, zprofed, zprofen, zpronewn, zpronewd
@@ -93,5 +94,5 @@
       !
       !  Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      zstrn )
+      CALL wrk_alloc( jpi, jpj,      zmixnano, zmixdiat, zstrn )
       CALL wrk_alloc( jpi, jpj, jpk, zpislopeadn, zpislopeadd, zprdia, zprbio, zprdch, zprnch, zysopt ) 
       CALL wrk_alloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
@@ -104,5 +105,5 @@
 
       ! Computation of the optimal production
-      prmax(:,:,:) = 0.6_wp * r1_rday * tgfunc(:,:,:) 
+      prmax(:,:,:) = 0.8_wp * r1_rday * tgfunc(:,:,:) 
       IF( lk_degrad )  prmax(:,:,:) = prmax(:,:,:) * facvol(:,:,:) 
 
@@ -127,13 +128,9 @@
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   zval = MAX( 1., zstrn(ji,jj) )
-                  zmxl_fac(ji,jj,jk) = zval
+                  IF( fsdept(ji,jj,jk) <= hmld(ji,jj) ) THEN
+                     zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
+                  ENDIF
                   zmxl_chl(ji,jj,jk) = zval / 24.
-                  IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                     zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
-                     zmxl_fac(ji,jj,jk) = zmxl_fac(ji,jj,jk) * zval
-                     zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk) * zval
-                  ENDIF
-                  zmxl_fac(ji,jj,jk) = ( 1. - exp( -0.2 * zmxl_fac(ji,jj,jk) ) )
-                  zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk)
+                  zmxl_fac(ji,jj,jk) = 1.5 * zval / ( 12. + zval )
                ENDIF
             END DO
@@ -143,4 +140,7 @@
       zprbio(:,:,:) = prmax(:,:,:) * zmxl_fac(:,:,:)
       zprdia(:,:,:) = zprbio(:,:,:)
+
+      ! Maximum light intensity
+      WHERE( zstrn(:,:) < 1.e0 ) zstrn(:,:) = 24.
 
       ! Computation of the P-I slope for nanos and diatoms
@@ -253,11 +253,15 @@
       END DO
 
+      !  Mixed-layer effect on production 
       !  Sea-ice effect on production
 
       DO jk = 1, jpkm1
-         zprbio(:,:,jk) = zprbio(:,:,jk) * ( 1. - fr_i(:,:) )
-         zprdia(:,:,jk) = zprdia(:,:,jk) * ( 1. - fr_i(:,:) )
-      END DO
-
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+               zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+            END DO
+         END DO
+      END DO
 
       ! Computation of the various production terms 
@@ -303,5 +307,6 @@
                   zprod    = rday * zprorcan(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk)
                   zprochln = chlcmin * 12. * zprorcan (ji,jj,jk)
-                  zprochln = zprochln + (chlcnm-chlcmin) * 12. * zprod / &
+                  chlcnm_n   = MIN ( chlcnm, ( chlcnm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
+                  zprochln = zprochln + (chlcnm_n-chlcmin) * 12. * zprod / &
                                         & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
                   !  production terms for diatoms ( chlorophyll )
@@ -309,5 +314,6 @@
                   zprod    = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk)
                   zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
-                  zprochld = zprochld + (chlcdm-chlcmin) * 12. * zprod / &
+                  chlcdm_n   = MIN ( chlcdm, ( chlcdm / (1. - 1.14 / 43.4 * tsn(ji,jj,jk,jp_tem))) * (1. - 1.14 / 43.4 * 20.))
+                  zprochld = zprochld + (chlcdm_n-chlcmin) * 12. * zprod / &
                                         & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
                   !   Update the arrays TRA which contain the Chla sources and sinks
@@ -485,5 +491,5 @@
      ENDIF
      !
-     CALL wrk_dealloc( jpi, jpj,      zstrn )
+     CALL wrk_dealloc( jpi, jpj,  zmixnano, zmixdiat,    zstrn )
      CALL wrk_dealloc( jpi, jpj, jpk, zpislopeadn, zpislopeadd, zprdia, zprbio, zprdch, zprnch, zysopt ) 
      CALL wrk_dealloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 7180)
@@ -25,4 +25,5 @@
    USE p4zlim          ! Phytoplankton limitation factors
    USE p4zprod         !  Growth rate of the 2 phyto groups
+   USE p4zsink         !  Sinking of particles
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -43,7 +44,9 @@
    REAL(wp), PUBLIC ::  xsilab     !: fraction of labile biogenic silica 
    REAL(wp), PUBLIC ::  oxymin     !: halk saturation constant for anoxia 
+   REAL(wp), PUBLIC ::  oxymin2    !: Minimum O2 concentration for oxic remin.
+   REAL(wp), PUBLIC ::  feratb     !: Fe/C quota in bacteria
+   REAL(wp), PUBLIC ::  xkferb     !: Half-saturation constant for bacteria Fe/C
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitr     !: denitrification array
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   denitnh4   !: -    -    -    -   -
 
 
@@ -71,10 +74,10 @@
       REAL(wp) ::   zremik, zsiremin 
       REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
-      REAL(wp) ::   zbactfer, zolimit
-      REAL(wp) ::   zosil, ztem
+      REAL(wp) ::   zbactfer, zolimit, zdenitnh4
+      REAL(wp) ::   zosil, ztem,ztoto,zpuis
       REAL(wp) ::   zonitr, zstep, zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod, zfacsi, zfacsib
       !!---------------------------------------------------------------------
       !
@@ -83,5 +86,5 @@
       ! Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj,      ztempbac )
-      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
+      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zfacsi, zfacsib )
 
       ! Initialization of local variables
@@ -91,4 +94,6 @@
       zdepprod(:,:,:) = 1._wp
       ztempbac(:,:)   = 0._wp
+      zfacsib(:,:,:)  = xsilab / ( 1.0 - xsilab )
+      zfacsi(:,:,:)   = xsilab 
 
       ! Computation of the mean phytoplankton concentration as
@@ -117,5 +122,6 @@
             DO ji = 1, jpi
                ! denitrification factor computed from O2 levels
-               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
+               ! ----------------------------------------------
+               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( oxymin2 - trb(ji,jj,jk,jpoxy) )    &
                   &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
@@ -147,5 +153,14 @@
                zolimi (ji,jj,jk) = MAX( 0.e0, zolimi (ji,jj,jk) )
                denitr (ji,jj,jk) = MAX( 0.e0, denitr (ji,jj,jk) )
-               !
+               ! Update of the tracers trends
+               ! ----------------------------
+               tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk)
+               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) + zolimi (ji,jj,jk) + denitr(ji,jj,jk)
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - denitr (ji,jj,jk) * rdenit
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) - zolimi (ji,jj,jk) - denitr(ji,jj,jk)
+               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - zolimi (ji,jj,jk) * o2ut
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zolimi (ji,jj,jk) + denitr(ji,jj,jk)
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zolimi(ji,jj,jk)    &
+               &                     + ( rdenit + 1.) * denitr(ji,jj,jk) )
             END DO
          END DO
@@ -165,20 +180,20 @@
                zonitr  = nitrif * zstep * trb(ji,jj,jk,jpnh4) * ( 1.- nitrfac(ji,jj,jk) )  &
                &         / ( 1.+ emoy(ji,jj,jk) ) * ( 1. + fr_i(ji,jj) * emoy(ji,jj,jk) )
-               denitnh4(ji,jj,jk) = nitrif * zstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 
+               zdenitnh4 = nitrif * zstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 
                ! Update of the tracers trends
                ! ----------------------------
-               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zonitr - denitnh4(ji,jj,jk)
-               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zonitr - rdenita * denitnh4(ji,jj,jk)
+               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zonitr - zdenitnh4
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + zonitr - rdenita * zdenitnh4
                tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) - o2nit * zonitr
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2 * rno3 * zonitr + rno3 * ( rdenita - 1. ) * denitnh4(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-
-       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2 * rno3 * zonitr + rno3 * ( rdenita - 1. ) * zdenitnh4
+            END DO
+         END DO
+      END DO
+
+      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('rem1')")
          CALL prt_ctl_trc_info(charout)
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-       ENDIF
+      ENDIF
 
       DO jk = 1, jpkm1
@@ -190,6 +205,6 @@
                ! studies (especially at Papa) have shown this uptake to be significant
                ! ----------------------------------------------------------
-               zbactfer = 10.e-6 *  rfact2 * prmax(ji,jj,jk) * xlimbacl(ji,jj,jk)             &
-                  &              * trb(ji,jj,jk,jpfer) / ( 2.5E-10 + trb(ji,jj,jk,jpfer) )    &
+               zbactfer = feratb *  rfact2 * prmax(ji,jj,jk) * xlimbacl(ji,jj,jk)             &
+                  &              * trb(ji,jj,jk,jpfer) / ( xkferb + trb(ji,jj,jk,jpfer) )    &
                   &              * zdepprod(ji,jj,jk) * zdepbac(ji,jj,jk)
 #if defined key_kriest
@@ -205,9 +220,13 @@
       END DO
 
-       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
+      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('rem2')")
          CALL prt_ctl_trc_info(charout)
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-       ENDIF
+      ENDIF
+
+      ! Initialization of the array which contains the labile fraction
+      ! of bSi. Set to a constant in the upper ocean
+      ! ---------------------------------------------------------------
 
       DO jk = 1, jpkm1
@@ -218,29 +237,27 @@
                zstep = zstep * facvol(ji,jj,jk)
 # endif
+               zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) )
+               zsatur   = MAX( rtrn, ( sio3eq(ji,jj,jk) - trb(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn ) )
+               zsatur2  = ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**37
+               znusil   = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25
+
                ! Remineralization rate of BSi depedant on T and saturation
                ! ---------------------------------------------------------
-               zsatur   = ( sio3eq(ji,jj,jk) - trb(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
-               zsatur   = MAX( rtrn, zsatur )
-               zsatur2  = ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**37
-               znusil   = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25
-               znusil2  = 0.225  * ( 1. + tsn(ji,jj,1,jp_tem) / 15.) + 0.775 * zsatur2
-
-               ! Two classes of BSi are considered : a labile fraction and 
-               ! a more refractory one. The ratio between both fractions is
-               ! constant and specified in the namelist.
-               ! ----------------------------------------------------------
-               zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) ) 
-               zdep     = MAX( 0., fsdept(ji,jj,jk) - zdep )
-               ztem     = MAX( tsn(ji,jj,1,jp_tem), 0. )
-               zfactdep = xsilab * EXP(-( xsiremlab - xsirem ) * znusil2 * zdep / wsbio2 ) * ztem / ( ztem + 10. )
-               zsiremin = ( xsiremlab * zfactdep + xsirem * ( 1. - zfactdep ) ) * zstep * znusil
+               IF ( fsdept(ji,jj,jk) > zdep ) THEN
+                  zfacsib(ji,jj,jk) = zfacsib(ji,jj,jk-1) * EXP( -0.5 * ( xsiremlab - xsirem )  &
+                  &                   * znusil * fse3t(ji,jj,jk) / wsbio4(ji,jj,jk) )
+                  zfacsi(ji,jj,jk)  = zfacsib(ji,jj,jk) / ( 1.0 + zfacsib(ji,jj,jk) )
+                  zfacsib(ji,jj,jk) = zfacsib(ji,jj,jk) * EXP( -0.5 * ( xsiremlab - xsirem )    &
+                  &                   * znusil * fse3t(ji,jj,jk) / wsbio4(ji,jj,jk) )
+               ENDIF
+               zsiremin = ( xsiremlab * zfacsi(ji,jj,jk) + xsirem * ( 1. - zfacsi(ji,jj,jk) ) ) * zstep * znusil
                zosil    = zsiremin * trb(ji,jj,jk,jpgsi)
                !
                tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) - zosil
                tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + zosil
-               !
-            END DO
-         END DO
-      END DO
+            END DO
+         END DO
+      END DO
+               
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
@@ -250,17 +267,4 @@
        ENDIF
 
-      ! Update the arrays TRA which contain the biological sources and sinks
-      ! --------------------------------------------------------------------
-
-      DO jk = 1, jpkm1
-         tra(:,:,jk,jppo4) = tra(:,:,jk,jppo4) + zolimi (:,:,jk) + denitr(:,:,jk)
-         tra(:,:,jk,jpnh4) = tra(:,:,jk,jpnh4) + zolimi (:,:,jk) + denitr(:,:,jk)
-         tra(:,:,jk,jpno3) = tra(:,:,jk,jpno3) - denitr (:,:,jk) * rdenit
-         tra(:,:,jk,jpdoc) = tra(:,:,jk,jpdoc) - zolimi (:,:,jk) - denitr(:,:,jk)
-         tra(:,:,jk,jpoxy) = tra(:,:,jk,jpoxy) - zolimi (:,:,jk) * o2ut
-         tra(:,:,jk,jpdic) = tra(:,:,jk,jpdic) + zolimi (:,:,jk) + denitr(:,:,jk)
-         tra(:,:,jk,jptal) = tra(:,:,jk,jptal) + rno3 * ( zolimi(:,:,jk) + ( rdenit + 1.) * denitr(:,:,jk) )
-      END DO
-
       IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
           zrfact2 = 1.e3 * rfact2r
@@ -268,13 +272,7 @@
           CALL iom_put( "DENIT" , denitr(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2  )  ! Denitrification
       ENDIF
-
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
-         WRITE(charout, FMT="('rem4')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-      ENDIF
       !
       CALL wrk_dealloc( jpi, jpj,      ztempbac )
-      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
+      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi, zfacsi, zfacsib )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_rem')
@@ -296,5 +294,5 @@
       !!----------------------------------------------------------------------
       NAMELIST/nampisrem/ xremik, nitrif, xsirem, xsiremlab, xsilab,   &
-      &                   oxymin
+      &                   oxymin, oxymin2, feratb, xkferb
       INTEGER :: ios                 ! Local integer output status for namelist read
 
@@ -318,9 +316,11 @@
          WRITE(numout,*) '    NH4 nitrification rate                    nitrif    =', nitrif
          WRITE(numout,*) '    halk saturation constant for anoxia       oxymin    =', oxymin
+         WRITE(numout,*) '    Minimum O2 concentration for oxic remin.  oxymin2   =', oxymin2
+         WRITE(numout,*) '    Bacterial Fe/C ratio                      feratb    =', feratb
+         WRITE(numout,*) '    Half-saturation constant for bact. Fe/C   xkferb    =', xkferb
       ENDIF
       !
       nitrfac (:,:,:) = 0._wp
       denitr  (:,:,:) = 0._wp
-      denitnh4(:,:,:) = 0._wp
 
    END SUBROUTINE p4z_rem_init
@@ -331,5 +331,5 @@
       !!                     ***  ROUTINE p4z_rem_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( denitr(jpi,jpj,jpk), denitnh4(jpi,jpj,jpk), STAT=p4z_rem_alloc )
+      ALLOCATE( denitr(jpi,jpj,jpk), STAT=p4z_rem_alloc )
       !
       IF( p4z_rem_alloc /= 0 )   CALL ctl_warn('p4z_rem_alloc: failed to allocate arrays')
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 7180)
@@ -199,5 +199,5 @@
          tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
 #if defined key_ligand
-         trn(:,:,:,jpfep) = trn(:,:,:,jpfep) + (ironsed(:,:,:) * fep_rats ) * rfact2
+         tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + (ironsed(:,:,:) * fep_rats ) * rfact2
 #endif
          !
@@ -211,6 +211,6 @@
          tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
 #if defined key_ligand
-         trn(:,:,:,jpfep) = trn(:,:,:,jpfep) + ( hydrofe(:,:,:) * fep_rath ) * rfact2
-         trn(:,:,:,jplgw) = trn(:,:,:,jplgw) + ( hydrofe(:,:,:) * 0.5 ) * rfact2
+         tra(:,:,:,jpfep) = tra(:,:,:,jpfep) + ( hydrofe(:,:,:) * fep_rath ) * rfact2
+         tra(:,:,:,jplgw) = tra(:,:,:,jplgw) + ( hydrofe(:,:,:) * 0.5 ) * rfact2
 #endif
          !
@@ -343,5 +343,5 @@
             zwstpoc              = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
 #   if defined key_ligand
-            trn(ji,jj,ikt,jpfep) = trn(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep
+            tra(ji,jj,ikt,jpfep) = tra(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep
 #   endif
 # else
@@ -351,5 +351,5 @@
             zwstpoc = trb(ji,jj,ikt,jppoc) * zws3 
 #   if defined key_ligand
-            trn(ji,jj,ikt,jpfep) = trn(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep
+            tra(ji,jj,ikt,jpfep) = tra(ji,jj,ikt,jpfep) - trn(ji,jj,ikt,jpfep) * zwssfep
 #   endif
 # endif
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 7180)
@@ -172,5 +172,5 @@
                  wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * FLOAT( iiter2 ) )
 #if defined key_ligand
-                 wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
+                 wsfep(ji,jj,jk)  = MIN( wsfep(ji,jj,jk), zwsmax * FLOAT( iiter1 ) )
 #endif
                ENDIF
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmeso.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmeso.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmeso.F90	(revision 7180)
@@ -371,6 +371,6 @@
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zgrazpof - zgrazfffp + zfracfe
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) - zgrazfffg + zgrapof - zfracfe
-               zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + rtrn )
-               zgrazcal = ( zgrazffeg + zgrazpoc ) * (1. - part2) * zfracal
+               zfracal = trb(ji,jj,jk,jpcal) / ( trb(ji,jj,jk,jpgoc) + rtrn )
+               zgrazcal = zgrazffeg * (1. - part2) * zfracal
 #endif
                !  calcite production
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmicro.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmicro.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmicro.F90	(revision 7180)
@@ -84,5 +84,5 @@
       REAL(wp) :: zgraztotc, zgraztotn, zgraztotp, zgraztotf, zbasresn, zbasresp, zbasresf
       REAL(wp) :: zgradoc, zgradon, zgradop, zgraref, zgradoct, zgradont, zgradopt, zgrareft
-      REAL(wp) :: zexcess, zgraren, zgrarep, zgrarem, zfracal, zgrazcal
+      REAL(wp) :: zexcess, zgraren, zgrarep, zgrarem
       REAL(wp) :: zgrapoc, zgrapon, zgrapop, zgrapof, zprcaca, zmortz
       REAL(wp) :: zrespz, ztortz, zgrasratf, zgrasratn, zgrasratp
@@ -172,9 +172,9 @@
                zgrazpon  = zgrazpoc * trb(ji,jj,jk,jppon) / ( trb(ji,jj,jk,jppoc) + rtrn )
                zgrazpop  = zgrazpoc * trb(ji,jj,jk,jppop) / ( trb(ji,jj,jk,jppoc) + rtrn )
+               zgrazpof  = zgrazpoc* trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
                zgrazdc   = zgraze  * ztmp4  * zdenom
                zgrazdn   = zgrazdc * trb(ji,jj,jk,jpndi) / (trb(ji,jj,jk,jpdia) + rtrn)
                zgrazdp   = zgrazdc * trb(ji,jj,jk,jppdi) / (trb(ji,jj,jk,jpdia) + rtrn)
                zgrazdf   = zgrazdc * trb(ji,jj,jk,jpdfe) / (trb(ji,jj,jk,jpdia) + rtrn)
-               zgrazpof  = zgrazpoc* trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
                !
                zgraztotc = zgraznc + zgrazpoc + zgrazdc + zgrazz + zgrazpc
@@ -296,11 +296,9 @@
                prodcal(ji,jj,jk) = prodcal(ji,jj,jk) + zprcaca  ! prodcal=prodcal(nanophy)+prodcal(microzoo)+prodcal(mesozoo)
                !
-               zfracal = trb(ji,jj,jk,jpcal) / (trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + rtrn )
-               zgrazcal = zgrazpoc * (1. - part ) * zfracal 
                zprcaca = part * zprcaca
-               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarem - zprcaca + zgrazcal
-               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * ( zgrazcal - zprcaca )    &
+               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) + zgrarem - zprcaca
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) - 2. * zprcaca     &
                &                     + rno3 * zgraren
-               tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca - zgrazcal
+               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) + ztortz * xkr_dmicro &
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmort.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmort.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zmort.F90	(revision 7180)
@@ -204,5 +204,4 @@
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zmortp * zfactfe
                prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + zmortp
-               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zmortp
 #endif
             END DO
@@ -285,14 +284,14 @@
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zmortp2 * zfactfe
 #else
-               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zrespp2 + 0.5 * ztortp2
-               tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + (zrespp2 + 0.5 * ztortp2) * zfactn
-               tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + (zrespp2 + 0.5 * ztortp2) * zfactp
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + (zrespp2 + 0.5 * ztortp2) * zfactfe
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + 0.5 * ztortp2
-               tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + 0.5 * ztortp2 * zfactn
-               tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + 0.5 * ztortp2 * zfactp
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + 0.5 * ztortp2 * zfactfe
-               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + 0.5 * ztortp2
-               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zrespp2 + 0.5 * ztortp2
+               tra(ji,jj,jk,jpgoc) = tra(ji,jj,jk,jpgoc) + zrespp2 
+               tra(ji,jj,jk,jpgon) = tra(ji,jj,jk,jpgon) + zrespp2 * zfactn
+               tra(ji,jj,jk,jpgop) = tra(ji,jj,jk,jpgop) + zrespp2 * zfactp
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zrespp2 * zfactfe
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + ztortp2
+               tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) + ztortp2 * zfactn
+               tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) + ztortp2 * zfactp
+               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + ztortp2 * zfactfe
+               prodpoc(ji,jj,jk) = prodpoc(ji,jj,jk) + ztortp2
+               prodgoc(ji,jj,jk) = prodgoc(ji,jj,jk) + zrespp2
 #endif
             END DO
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90	(revision 7180)
@@ -128,4 +128,7 @@
                   remint = 0.
                   !
+                  zsizek1  = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
+                  zsizek = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
+                  !
                   IF ( fsdept(ji,jj,jk-1) <= zdep ) THEN
                     ! 
@@ -150,14 +153,11 @@
                        ! -----------------------------------------------------
                        !
-                       zsizek  = zdep / (wsbio2 + rtrn) * tgfunc(ji,jj,jk-1)
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
                        ! the concentration of each lability class is calculated
                        ! as the sum of the different sources and sinks
                        ! Please note that production of new GOC experiences
                        ! degradation 
-                       alphag(ji,jj,jk,jn) = alphan(jn) / (reminp(jn) * tgfunc(ji,jj,jk-1) )       &
-                       &   * (1. - exp( -reminp(jn) * zsizek ) ) * exp( -reminp(jn) * zsizek1 )     &
-                       &   * zpoc + prodgoc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) )        &
-                       &   * rday / rfact2 * alphan(jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc &
+                       &   + prodgoc(ji,jj,jk) * alphan(jn) / tgfunc(ji,jj,jk) / reminp(jn)             &
+                       &   * ( 1. - exp( -reminp(jn) * zsizek ) ) * rday / rfact2
                        alphat = alphat + alphag(ji,jj,jk,jn)
                        remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
@@ -175,12 +175,8 @@
                     !
                     DO jn = 1, jcpoc
-                       zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio4(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio4(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn)                &
-                       &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1)  &
-                       &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek ) )                 &
-                       &   * exp( -reminp(jn) * zsizek1 ) + prodgoc(ji,jj,jk) / tgfunc(ji,jj,jk)       &
-                       &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) ) * rday / rfact2 * alphan(jn)     &
-                       &   / reminp(jn)
+                       alphag(ji,jj,jk,jn) = alphag(ji,jj,jk-1,jn) * exp( -reminp(jn) * ( zsizek              &
+                       &   + zsizek1 ) ) * zpoc + ( prodgoc(ji,jj,jk-1) / tgfunc(ji,jj,jk-1) * ( 1.           &
+                       &   - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn) * zsizek ) + prodgoc(ji,jj,jk) &
+                       &   / tgfunc(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek ) ) ) * rday / rfact2 / reminp(jn)
                        alphat = alphat + alphag(ji,jj,jk,jn)
                        remint = remint + alphag(ji,jj,jk,jn) * reminp(jn)
@@ -212,8 +208,8 @@
                zopoc2 = zremig  * trb(ji,jj,jk,jpgoc)
                orem(ji,jj,jk) = zopoc2
+               zorem3(ji,jj,jk) = zremig * solgoc * trb(ji,jj,jk,jpgoc)
                zopon2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpgon)
                zopop2 = xremipp / xremipc * zremig * trb(ji,jj,jk,jpgop)
                zofer2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpbfe)
-               zorem3(ji,jj,jk) = zremig * solgoc * trb(ji,jj,jk,jpgoc)
 
                ! Update the appropriate tracers trends
@@ -258,13 +254,11 @@
         DO jj = 1, jpj
            DO ji = 1, jpi
-              IF (tmask(ji,jj,jk) == 1.) THEN
-                zdep = hmld(ji,jj)
-                IF( fsdept(ji,jj,jk) <= zdep ) THEN
-                  totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2
-                  ! The temperature effect is included here
-                  totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk)
-                  totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2    &
-                     &                / ( trb(ji,jj,jk,jppoc) + rtrn )
-                ENDIF
+              zdep = hmld(ji,jj)
+              IF (tmask(ji,jj,jk) == 1. .AND. fsdept(ji,jj,jk) <= zdep ) THEN
+                 totprod(ji,jj) = totprod(ji,jj) + prodpoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2
+                 ! The temperature effect is included here
+                 totthick(ji,jj) = totthick(ji,jj) + fse3t(ji,jj,jk)* tgfunc(ji,jj,jk)
+                 totcons(ji,jj) = totcons(ji,jj) - conspoc(ji,jj,jk) * fse3t(ji,jj,jk) * rday/ rfact2    &
+                 &                / ( trb(ji,jj,jk,jppoc) + rtrn )
               ENDIF
            END DO
@@ -319,4 +313,8 @@
                   remint = 0.
                   !
+                  ! the scale factors are corrected with temperature
+                  zsizek1  = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
+                  zsizek = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
+                  !
                   ! Special treatment of the level just below the MXL
                   ! See the comments in the GOC section
@@ -332,11 +330,10 @@
                     !
                     DO jn = 1, jcpoc
-                       ! the scale factor is corrected with temperature
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
                        ! computation of the lability spectrum applying the 
                        ! different sources and sinks
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) + zorem3(ji,jj,jk)                 &
-                       &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday / rfact2                   &
-                       &   * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                      alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn) * zsizek ) * zpoc  &
+                       &   + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk) * alphag(ji,jj,jk,jn) ) &
+                       &   / tgfunc(ji,jj,jk) / reminp(jn) * rday / rfact2 * ( 1. - exp( -reminp(jn)     &
+                       &   * zsizek ) )
                        alphap(ji,jj,jk,jn) = MAX( 0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
@@ -355,17 +352,11 @@
                     !
                     DO jn = 1, jcpoc
-                       zsizek  = fse3t(ji,jj,jk-1) / 2. / (wsbio3(ji,jj,jk-1) + rtrn) * tgfunc(ji,jj,jk-1)
-                       zsizek1 = fse3t(ji,jj,jk) / 2. / (wsbio3(ji,jj,jk) + rtrn) * tgfunc(ji,jj,jk)
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)             &
-                       &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1)  &
-                       &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek ) )              &
-                       &   * exp( -reminp(jn) * zsizek1 ) + prodpoc(ji,jj,jk) / tgfunc(ji,jj,jk)    &
-                       &   * ( 1. - exp( -reminp(jn) * zsizek1 ) ) ) * rday / rfact2 * alphan(jn)  &
-                       &   / reminp(jn)
-                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk,jn) + zorem3(ji,jj,jk-1)              &
-                       &  * alphag(ji,jj,jk-1,jn) / tgfunc(ji,jj,jk-1) * rday / rfact2 * ( 1.       &
-                       &  - exp( -reminp(jn) * zsizek ) ) * exp( -reminp(jn) * zsizek1 )                  &
-                       &  + zorem3(ji,jj,jk) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * rday         &
-                       &  / rfact2 * alphag(ji,jj,jk,jn) / reminp(jn) / tgfunc(ji,jj,jk)
+                       alphap(ji,jj,jk,jn) = alphap(ji,jj,jk-1,jn) * exp( -reminp(jn)                       &
+                       &   * ( zsizek + zsizek1 ) ) * zpoc + ( prodpoc(ji,jj,jk-1) * alphan(jn)             &
+                       &   + zorem3(ji,jj,jk-1) * alphag(ji,jj,jk-1,jn) ) * rday / rfact2 / reminp(jn)      &
+                       &   / tgfunc(ji,jj,jk-1) * ( 1. - exp( -reminp(jn) * zsizek1 ) ) * exp( -reminp(jn)  &
+                       &   * zsizek ) + ( prodpoc(ji,jj,jk) * alphan(jn) + zorem3(ji,jj,jk)                 &
+                       &   * alphag(ji,jj,jk,jn) ) * rday / rfact2 / reminp(jn) / tgfunc(ji,jj,jk) * ( 1.   &
+                       &   - exp( -reminp(jn) * zsizek ) )
                        alphap(ji,jj,jk,jn) = max(0., alphap(ji,jj,jk,jn) )
                        alphat = alphat + alphap(ji,jj,jk,jn)
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90	(revision 7180)
@@ -88,5 +88,5 @@
       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zprontot, zproptot, zprodtot
       REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
-      REAL(wp) ::   zpislopen, zpislopep, zpisloped
+      REAL(wp) ::   zpislopen, zpislopep, zpisloped, thetannm_n, thetandm_n, thetanpm_n
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
       REAL(wp) ::   zrfact2
@@ -382,15 +382,21 @@
                   znanotot = enano(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod = rday * (zpronewn(ji,jj,jk) + zproregn(ji,jj,jk)) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
-                  zprochln = thetannm * zprod / ( zpislopeadn(ji,jj,jk) * znanotot + rtrn )
+                  thetannm_n   = MIN ( thetannm, ( thetannm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
+                  &               * (1. - 1.14 / 43.4 * 20.))
+                  zprochln = thetannm_n * zprod / ( zpislopeadn(ji,jj,jk) * znanotot + rtrn )
                   zprochln = MAX(zprochln, chlcmin * 12. * zprorcan (ji,jj,jk) )
                      !  production terms for picophyto. ( chlorophyll )
                   zpicotot = epico(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod = rday * (zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)) * zprchlp(ji,jj,jk) * xlimpic(ji,jj,jk)
-                  zprochlp = thetanpm * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
+                  thetanpm_n   = MIN ( thetanpm, ( thetanpm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
+                  &               * (1. - 1.14 / 43.4 * 20.))
+                  zprochlp = thetanpm_n * zprod / ( zpislopeadp(ji,jj,jk) * zpicotot + rtrn )
                   zprochlp = MAX(zprochlp, chlcmin * 12. * zprorcap(ji,jj,jk) )
                   !  production terms for diatomees ( chlorophyll )
                   zdiattot = ediat(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod = rday * (zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
-                  zprochld = thetandm * zprod / ( zpislopeadd(ji,jj,jk) * zdiattot + rtrn )
+                  thetandm_n   = MIN ( thetandm, ( thetandm / (1. - 1.14 / 43.4 *tsn(ji,jj,jk,jp_tem)))   &
+                  &               * (1. - 1.14 / 43.4 * 20.))
+                  zprochld = thetandm_n * zprod / ( zpislopeadd(ji,jj,jk) * zdiattot + rtrn )
                   zprochld = MAX(zprochld, chlcmin * 12. * zprorcad(ji,jj,jk) )
                   !   Update the arrays TRA which contain the Chla sources and sinks
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zrem.F90	(revision 7180)
@@ -75,10 +75,9 @@
       REAL(wp) ::   zremik, zremikc, zremikn, zremikp, zsiremin 
       REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
-      REAL(wp) ::   zbactfer, zolimit
+      REAL(wp) ::   zbactfer, zolimit, zonitr, zstep, zrfact2
       REAL(wp) ::   zosil, ztem, zdenitnh4, zolimic, zolimin, zolimip, zdenitrn, zdenitrp
-      REAL(wp) ::   zonitr, zstep, zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zwork1, zdepprod
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zwork1, zdepprod, zfacsi, zfacsib
       !!---------------------------------------------------------------------
       !
@@ -87,9 +86,11 @@
       ! Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj,      ztempbac                  )
-      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zwork1 )
+      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zwork1, zfacsi, zfacsib )
 
       ! Initialisation of temprary arrys
       zdepprod(:,:,:) = 1._wp
       ztempbac(:,:)   = 0._wp
+      zfacsib(:,:,:)  = xsilab / ( 1.0 - xsilab )
+      zfacsi(:,:,:)   = xsilab
 
       ! Computation of the mean phytoplankton concentration as
@@ -236,4 +237,8 @@
       ENDIF
 
+      ! Initialization of the array which contains the labile fraction
+      ! of bSi. Set to a constant in the upper ocean
+      ! ---------------------------------------------------------------
+
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -243,26 +248,24 @@
                zstep = zstep * facvol(ji,jj,jk)
 # endif
+               zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) )
+               zsatur   = MAX( rtrn, ( sio3eq(ji,jj,jk) - trb(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn ) )
+               zsatur2  = ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**37
+               znusil   = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25
+
                ! Remineralization rate of BSi depedant on T and saturation
                ! ---------------------------------------------------------
-               zsatur   = ( sio3eq(ji,jj,jk) - trb(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
-               zsatur   = MAX( rtrn, zsatur )
-               zsatur2  = ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**37
-               znusil   = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2 * zsatur**9.25
-               znusil2  = 0.225  * ( 1. + tsn(ji,jj,1,jp_tem) / 15.) + 0.775 * zsatur2
-
-               ! Two classes of BSi are considered : a labile fraction and 
-               ! a more refractory one. The ratio between both fractions is
-               ! constant and specified in the namelist.
-               ! ----------------------------------------------------------
-               zdep     = MAX( hmld(ji,jj), heup_01(ji,jj) ) 
-               zdep     = MAX( 0., fsdept(ji,jj,jk) - zdep )
-               ztem     = MAX( tsn(ji,jj,1,jp_tem), 0. )
-               zfactdep = xsilab * EXP(-( xsiremlab - xsirem ) * znusil2 * zdep / wsbio2 ) * ztem / ( ztem + 10. )
-               zsiremin = ( xsiremlab * zfactdep + xsirem * ( 1. - zfactdep ) ) * zstep * znusil
+               IF ( fsdept(ji,jj,jk) > zdep ) THEN
+                  zfacsib(ji,jj,jk) = zfacsib(ji,jj,jk-1) * EXP( -0.5 * ( xsiremlab - xsirem )  &
+                  &                   * znusil * fse3t(ji,jj,jk) / wsbio4(ji,jj,jk) )
+                  zfacsi(ji,jj,jk)  = zfacsib(ji,jj,jk) / ( 1.0 + zfacsib(ji,jj,jk) )
+                  zfacsib(ji,jj,jk) = zfacsib(ji,jj,jk) * EXP( -0.5 * ( xsiremlab - xsirem )    &
+                  &                   * znusil * fse3t(ji,jj,jk) / wsbio4(ji,jj,jk) )
+               ENDIF
+               zsiremin = ( xsiremlab * zfacsi(ji,jj,jk) + xsirem * ( 1. - zfacsi(ji,jj,jk) ) ) * zstep * znusil
                zosil    = zsiremin * trb(ji,jj,jk,jpgsi)
                !
                tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) - zosil
                tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) + zosil
-               !
+               ENDIF
             END DO
          END DO
@@ -280,13 +283,7 @@
           CALL iom_put( "DENIT" , denitrc(:,:,:) * rdenit * rno3 * tmask(:,:,:) * zrfact2  )  ! Denitrification
       ENDIF
-
-      IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
-         WRITE(charout, FMT="('rem4')")
-         CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
-      ENDIF
       !
       CALL wrk_dealloc( jpi, jpj,      ztempbac                  )
-      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zwork1 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zwork1, zfacsi, zfacsib )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p5z_rem')
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 7180)
@@ -176,5 +176,5 @@
 #if defined key_ligand
          trn(:,:,:,jplgw) = 0.6E-9
-         trn(:,:,:,jpfep) = 0. * 5.e-6
+         trn(:,:,:,jpfep) = 0.
 #endif
 
@@ -316,5 +316,5 @@
 #if defined key_ligand
          trn(:,:,:,jplgw) = 0.6E-9
-         trn(:,:,:,jpfep) = 0. * 5.e-6
+         trn(:,:,:,jpfep) = 0.
 #endif
 
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl.F90	(revision 7180)
@@ -54,4 +54,5 @@
       CHARACTER (len=22) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd
+      INTEGER :: jn
       !!----------------------------------------------------------------------
       !
Index: /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7179)
+++ /branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 7180)
@@ -102,5 +102,5 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~ '
 
-         IF( ln_rsttr .AND.    &                     ! Restart: read in restart  file
+         IF( ln_rsttr .AND. .NOT. ln_top_euler .AND.   &     ! Restart: read in restart  file
             iom_varid( numrtr, 'sbc_'//TRIM(ctrcnm(1))//'_b', ldstop = .FALSE. ) > 0 ) THEN
             IF(lwp) WRITE(numout,*) '          nittrc000-nn_dttrc surface tracer content forcing fields red in the restart file'
