Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zbio.F90	(revision 3446)
@@ -121,5 +121,5 @@
       ENDIF
 
-      grazing(:,:,1) = 0.e0        ! zooplakton closure ( fbod)
+      xksi(:,:) = 0.e0        ! zooplakton closure ( fbod)
       IF( ln_diatrc ) THEN
          zw2d  (:,:,:) = 0.e0
@@ -186,5 +186,5 @@
                !    closure : flux grazing is redistributed below level jpkbio
                zzoobod = tmminz * zzoo * zzoo
-               grazing(ji,jj,1) = grazing(ji,jj,1) + (1-fdbod) * zzoobod * fse3t(ji,jj,jk)
+               xksi(ji,jj) = xksi(ji,jj) + (1-fdbod) * zzoobod * fse3t(ji,jj,jk)
                zboddet = fdbod * zzoobod
 
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zexp.F90	(revision 3446)
@@ -95,5 +95,5 @@
             DO ji = fs_2, fs_jpim1
                ze3t = 1. / fse3t(ji,jj,jk)
-               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + ze3t * dmin3(ji,jj,jk) * grazing(ji,jj,1)
+               tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) + ze3t * dmin3(ji,jj,jk) * xksi(ji,jj)
             END DO
          END DO
@@ -112,5 +112,5 @@
             ! Deposition of organic matter in the sediment
             zwork = vsed * trn(ji,jj,ikt,jpdet)
-            zsedpoca(ji,jj) = ( zwork + dminl(ji,jj) * grazing(ji,jj,1)   &
+            zsedpoca(ji,jj) = ( zwork + dminl(ji,jj) * xksi(ji,jj)   &
                &           - sedlam * sedpocn(ji,jj) - sedlostpoc * sedpocn(ji,jj) ) * rdt
             zgeolpoc = zgeolpoc + sedlostpoc * sedpocn(ji,jj) * e1e2t(ji,jj)
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 3446)
@@ -83,5 +83,5 @@
       CALL p4z_opt  ( kt, jnt )     ! Optic: PAR in the water column
       CALL p4z_fechem(kt, jnt )     ! Iron chemistry/scavenging
-      CALL p4z_lim  ( kt      )     ! co-limitations by the various nutrients
+      CALL p4z_lim  ( kt, jnt )     ! co-limitations by the various nutrients
       CALL p4z_prod ( kt, jnt )     ! phytoplankton growth rate over the global ocean. 
       !                             ! (for each element : C, Si, Fe, Chl )
@@ -89,5 +89,5 @@
       CALL p4z_mort ( kt      )     ! phytoplankton mortality
       !                             ! zooplankton sources/sinks routines 
-      CALL p4z_micro( kt      )           ! microzooplankton
+      CALL p4z_micro( kt, jnt )           ! microzooplankton
       CALL p4z_meso ( kt, jnt )           ! mesozooplankton
 
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 3446)
@@ -37,5 +37,5 @@
    REAL(wp), PUBLIC ::  ligand     = 0.6E-9_wp  !: ligand concentration in the ocean 
 
-   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con
+   REAL(wp) :: kl1, kl2, kb1, kb2, ks, kpr, spd, con, kth
 
    !!* Substitution
@@ -65,7 +65,7 @@
       !
       INTEGER  ::   ji, jj, jk, jic
-      REAL(wp) ::   zdep, zlam1b, zlamfac
-      REAL(wp) ::   zkeq, zfeequi, zfesatur
-      REAL(wp) ::   zdenom1, zscave, zaggdfe, zcoag
+      REAL(wp) ::   zdep, zlam1a, zlam1b, zlamfac
+      REAL(wp) ::   zkeq, zfeequi, zfesatur, zfecoll
+      REAL(wp) ::   zdenom1, zscave, zaggdfea, zaggdfeb, zcoag
       REAL(wp) ::   ztrc, zdust
 #if ! defined key_kriest
@@ -74,5 +74,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zTL1, zFe3, ztotlig
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zFeL1, zFeL2, zTL2, zFe2, zFeP
-      REAL(wp) :: zkox, zkph1, zkph2, zph, zionic
+      REAL(wp) :: zkox, zkph1, zkph2, zph, zionic, ztligand
       REAL(wp) :: za, zb, zc, zkappa1, zkappa2, za0, za1, za2
       REAL(wp) :: zxs, zfunc, zp, zq, zd, zr, zphi, zfff, zp3, zq2
@@ -85,18 +85,10 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj, jpk, zFe3, zTL1, ztotlig )
-      zFe3   (:,:,:) = 0._wp
-      zTL1   (:,:,:) = 0._wp
-      ztotlig(:,:,:) = 0._wp
+                                          CALL wrk_alloc( jpi, jpj, jpk, zFe3, zTL1, ztotlig )
+      IF( ln_diatrc .AND. lk_iomput )     CALL wrk_alloc( jpi, jpj, jpk, zFeL1 )
+      !
       IF( ln_fechem ) THEN
-         CALL wrk_alloc( jpi, jpj, jpk, zTL2, zFeP )
-         zTL2(:,:,:) = 0._wp
-         zFeP(:,:,:) = 0._wp
-         IF( ln_diatrc .AND. lk_iomput ) THEN
-            CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL1, zFeL2 )
-            zFe2 (:,:,:) = 0._wp
-            zFeL1(:,:,:) = 0._wp
-            zFeL2(:,:,:) = 0._wp
-         ENDIF
+                                          CALL wrk_alloc( jpi, jpj, jpk, zTL2, zFeP )
+         IF( ln_diatrc .AND. lk_iomput )  CALL wrk_alloc( jpi, jpj, jpk, zFe2, zFeL2 )
       ENDIF
       ! Total ligand concentration : Ligands can be chosen to be constant or variable
@@ -104,6 +96,6 @@
       ! -------------------------------------------------
       IF( ln_ligvar ) THEN
-         ztotlig(:,:,:) =  0.09 * ( trn(:,:,:,jpdoc) * 1E6 + 40. ) - 3.2
-         ztotlig(:,:,:) = MAX( MIN( ztotlig(:,:,:), 10. ), 0.4 )
+         ztotlig(:,:,:) =  0.09 * trn(:,:,:,jpdoc) * 1E6 + ligand * 1E9
+         ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
       ELSE
          ztotlig(:,:,:) = ligand * 1E9
@@ -120,6 +112,7 @@
                DO ji = 1, jpi
                   ! Calculate ligand concentrations : assume 2/3rd of excess goes to L2 and 1/3rd to L1
-                  zTL2(ji,jj,jk) = 0.000001       + 0.67 * ( ztotlig(ji,jj,jk) - 0.4 )
-                  zTL1(ji,jj,jk) = 0.4 - 0.000001 + 0.33 * ( ztotlig(ji,jj,jk) - 0.4 )
+                  ztligand       = ztotlig(ji,jj,jk) - ligand * 1E9
+                  zTL2(ji,jj,jk) =                0.000001 + 0.67 * ztligand
+                  zTL1(ji,jj,jk) = ligand * 1E9 - 0.000001 + 0.33 * ztligand
                   ! ionic strength from Millero et al. 1987
                   zionic = 19.9201 * tsn(ji,jj,jk,jp_sal) / ( 1000. - 1.00488 * tsn(ji,jj,jk,jp_sal) + rtrn )
@@ -150,8 +143,8 @@
                   ! calculate some parameters
                   za = 1 + ks / kpr
-                  zb = 1 + zkph1 / ( zkox + rtrn )
+                  zb = 1 + ( zkph1 + kth ) / ( zkox + rtrn )
                   zc = 1 + zkph2 / ( zkox + rtrn )
-                  zkappa1 = ( kb1 + zkph1) / kl1
-                  zkappa2 = ( kb2 + zkph2) / kl2
+                  zkappa1 = ( kb1 + zkph1 + kth ) / kl1
+                  zkappa2 = ( kb2 + zkph2 ) / kl2
                   za2 = zTL1(ji,jj,jk) * zb / za + zTL2(ji,jj,jk) * zc / za + zkappa1 + zkappa2 - ztfe / za
                   za1 = zkappa2 * zTL1(ji,jj,jk) * zb / za + zkappa1 * zTL2(ji,jj,jk) * zc / za &
@@ -216,4 +209,6 @@
             END DO
          END DO
+         IF( ln_diatrc .AND. lk_iomput )   &
+            &  zFeL1(:,:,1:jpkm1) = MAX( 0., trn(:,:,1:jpkm1,jpfer) * 1E9 - zFe3(:,:,1:jpkm1) )
          !
       ENDIF
@@ -234,6 +229,11 @@
                ! Scavenging onto dust is also included as evidenced from the DUNE experiments.
                ! --------------------------------------------------------------------------------------
-               zfeequi = zFe3(ji,jj,jk) * 1E-9 
-               IF( ln_fechem ) zfeequi = zfeequi + zFep(ji,jj,jk) * 1E-9
+               IF( ln_fechem ) THEN
+                  zfeequi = ( zFe3(ji,jj,jk) + zFe2(ji,jj,jk) ) * 1E-9
+                  zfecoll = ( 0.3 * zFeL1(ji,jj,jk) + 0.5 * zFeL2(ji,jj,jk) ) * 1E-9
+               ELSE
+                  zfeequi = zFe3(ji,jj,jk) * 1E-9 
+                  zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
+               ENDIF
 #if defined key_kriest
                ztrc   = ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpgsi) ) * 1.e6 
@@ -258,7 +258,5 @@
                !  due to increased lithogenic particles and let say it is unknown processes (precipitation, ...)
                !  -----------------------------------------------------------
-               ztfe = trn(ji,jj,jk,jpfer) * 1.e9 - 1. 
-               IF( ln_fechem )  ztfe = ztfe +  1. - ztotlig(ji,jj,jk) 
-               zlam1b  = xlam1 * MAX( 0.e0, ztfe )
+               zlam1b  = xlam1 * MAX( 0.e0, ( trn(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
                zcoag   = zfeequi * zlam1b * zstep
 
@@ -269,21 +267,22 @@
                zlamfac = MAX( 0.e0, ( gphit(ji,jj) + 55.) / 30. )
                zlamfac = MIN( 1.  , zlamfac )
-               zdep    =  MIN(1., 1000. / fsdept(ji,jj,jk) )
-#if ! defined key_kriest
-               zlam1b = (  80.* ( trn(ji,jj,jk,jpdoc) + 35.e-6 )                           &
-                  &     + 698.*   trn(ji,jj,jk,jppoc) + 1.05e4 * trn(ji,jj,jk,jpgoc)  )    &
-                  &     * xdiss(ji,jj,jk) + 1E-4 * ( 1. - zlamfac ) * zdep
+               zdep    = MIN( 1., 1000. / fsdept(ji,jj,jk) )
+               zlam1a  = ( 0.369  * 0.3 * trn(ji,jj,jk,jpdoc) + 102.4  * trn(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
+                   &   + ( 114.   * 0.3 * trn(ji,jj,jk,jpdoc) + 5.09E3 * trn(ji,jj,jk,jppoc) )
+#if defined key_kriest
+               zlam1a   = zlam1a + 1E-4 * ( 1. - zlamfac ) * zdep
+               zaggdfea = zlam1a * zstep * zfecoll
+               zaggdfeb = 0.
+               !
+               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 = (  80.* (trn(ji,jj,jk,jpdoc) + 35E-6)              &
-                  &     + 698.*  trn(ji,jj,jk,jppoc)  )                    &
-                  &     * xdiss(ji,jj,jk) + 1E-4 * ( 1. - zlamfac ) * zdep
-#endif
-               zaggdfe = zlam1b * zstep * 0.5 * ( trn(ji,jj,jk,jpfer) - zfeequi )
-               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zscave - zaggdfe - zcoag
-#if defined key_kriest
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1
-#else
-               tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1
-               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2
+               zlam1b = 3.53E3 *   trn(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk) + 1E-4 * ( 1. - zlamfac ) * zdep 
+               zaggdfea = zlam1a * zstep * zfecoll
+               zaggdfeb = zlam1b * zstep * zfecoll
+               !
+               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
+               tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zscave * zdenom2 + zaggdfeb
 #endif
             END DO
@@ -291,4 +290,11 @@
       END DO
       !
+      !  Define the bioavailable fraction of iron
+      !  ----------------------------------------
+      IF( ln_fechem ) THEN
+          biron(:,:,:) = MAX( 0., trn(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
+      ELSE
+          biron(:,:,:) = trn(:,:,:,jpfer) 
+      ENDIF
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
@@ -302,13 +308,14 @@
       IF( ln_diatrc .AND. lk_iomput ) THEN
          IF( jnt == nrdttrc ) THEN
-            CALL iom_put("Fe3"    , zFe3   (:,:,:) * tmask(:,:,:) )   ! Fe3+
-            CALL iom_put("TL1"    , zTL1   (:,:,:) * tmask(:,:,:) )   ! TL1
-            CALL iom_put("Totlig" , ztotlig(:,:,:) * tmask(:,:,:) )   ! TL
+            CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
+            CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
+            CALL iom_put("TL1"    , zTL1   (:,:,:)       * tmask(:,:,:) )   ! TL1
+            CALL iom_put("Totlig" , ztotlig(:,:,:)       * tmask(:,:,:) )   ! TL
+            CALL iom_put("Biron"  , biron  (:,:,:) * 1e9 * tmask(:,:,:) )   ! TL
             IF( ln_fechem ) THEN
-               CALL iom_put("Fe2" , zFe2   (:,:,:) * tmask(:,:,:) )   ! Fe2+
-               CALL iom_put("FeL1", zFeL1  (:,:,:) * tmask(:,:,:) )   ! FeL1
-               CALL iom_put("FeL2", zFeL2  (:,:,:) * tmask(:,:,:) )   ! FeL2
-               CALL iom_put("FeP" , zFeP   (:,:,:) * tmask(:,:,:) )   ! FeP
-               CALL iom_put("TL2" , zTL2   (:,:,:) * tmask(:,:,:) )   ! TL2
+               CALL iom_put("Fe2" , zFe2   (:,:,:)       * tmask(:,:,:) )   ! Fe2+
+               CALL iom_put("FeL2", zFeL2  (:,:,:)       * tmask(:,:,:) )   ! FeL2
+               CALL iom_put("FeP" , zFeP   (:,:,:)       * tmask(:,:,:) )   ! FeP
+               CALL iom_put("TL2" , zTL2   (:,:,:)       * tmask(:,:,:) )   ! TL2
             ENDIF
          ENDIF
@@ -321,10 +328,9 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zTL1, ztotlig )
+                                         CALL wrk_dealloc( jpi, jpj, jpk, zFe3, zTL1, ztotlig )
+      IF( ln_diatrc .AND. lk_iomput )    CALL wrk_dealloc( jpi, jpj, jpk, zFeL1 )
       IF( ln_fechem ) THEN
-         CALL wrk_dealloc( jpi, jpj, jpk, zTL2, zFeP )
-         IF( ln_diatrc .AND. lk_iomput ) THEN
-            CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL1, zFeL2 )
-         ENDIF
+                                         CALL wrk_dealloc( jpi, jpj, jpk, zTL2, zFeP )
+         IF( ln_diatrc .AND. lk_iomput ) CALL wrk_dealloc( jpi, jpj, jpk, zFe2, zFeL2 )
       ENDIF
       !
@@ -365,5 +371,5 @@
          ! initialization of some constants used by the complexe chemistry scheme
          ! ----------------------------------------------------------------------
-         spd = 3600.*24.
+         spd = 3600. * 24.
          con = 1.E9
          ! LIGAND KINETICS (values from Witter et al. 2000)
@@ -377,4 +383,6 @@
          ks  = 0.075
          kpr = 0.05
+         ! thermal reduction of Fe3
+         kth = 0.0048 * 24.
       ENDIF
       !
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90	(revision 3446)
@@ -52,6 +52,5 @@
       tgfunc2(:,:,:) = EXP( 0.07608  * tsn(:,:,:,jp_tem) )
 
-      ! Computation of the silicon dependant half saturation
-      ! constant for silica uptake
+      ! Computation of the silicon dependant half saturation  constant for silica uptake
       ! ---------------------------------------------------
       DO ji = 1, jpi
@@ -63,6 +62,6 @@
       !
       IF( nday_year == nyear_len(1) ) THEN
-         xksi    = xksimax
-         xksimax = 0._wp
+         xksi   (:,:) = xksimax(:,:)
+         xksimax(:,:) = 0._wp
       ENDIF
       !
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 3446)
@@ -19,4 +19,5 @@
    USE sms_pisces      ! PISCES variables
    USE p4zopt          ! Optical
+   USE iom             !  I/O manager
 
    IMPLICIT NONE
@@ -33,4 +34,6 @@
    REAL(wp), PUBLIC ::  concnfer  = 1.e-9_wp      !:  Iron half saturation for nanophyto 
    REAL(wp), PUBLIC ::  concdfer  = 2.e-9_wp      !:  Iron half saturation for diatoms  
+   REAL(wp), PUBLIC ::  concbno3  = 2.5e-7_wp      !:  NO3 half saturation  for bacteria 
+   REAL(wp), PUBLIC ::  concbnh4  = 2.5e-8_wp      !:  NH4 half saturation for bacteria
    REAL(wp), PUBLIC ::  xsizedia  = 5.e-7_wp      !:  Minimum size criteria for diatoms
    REAL(wp), PUBLIC ::  xsizephy  = 1.e-6_wp      !:  Minimum size criteria for nanophyto
@@ -59,5 +62,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_lim( kt )
+   SUBROUTINE p4z_lim( kt, jnt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_lim  ***
@@ -69,5 +72,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in)  :: kt
+      INTEGER, INTENT(in)  :: kt, jnt
       !
       INTEGER  ::   ji, jj, jk
@@ -89,5 +92,5 @@
                zno3    = trn(ji,jj,jk,jpno3) / 40.e-6
                zferlim = MAX( 3e-11 * zno3 * zno3, 5e-12 )
-               zferlim = MIN( zferlim, 3e-11 )
+               zferlim = MIN( zferlim, 5e-11 )
                trn(ji,jj,jk,jpfer) = MAX( trn(ji,jj,jk,jpfer), zferlim )
 
@@ -110,4 +113,17 @@
                zconc0nnh4        = MAX( concnnh4, ( zconcn2 * concnnh4 + concnnh4 * xsizern * zconcn ) * z1_trnphy )
 
+               ! Michaelis-Menten Limitation term for nutrients Small bacteria
+               ! -------------------------------------------------------------
+               zdenom = 1. /  ( concbno3 * concbnh4 + concbnh4 * trn(ji,jj,jk,jpno3) + concbno3 * trn(ji,jj,jk,jpnh4) )
+               xnanono3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * concbnh4 * zdenom
+               xnanonh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * concbno3 * zdenom
+               !
+               zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
+               zlim2    = trn(ji,jj,jk,jppo4) / ( trn(ji,jj,jk,jppo4) + concbnh4 )
+               zlim3    = trn(ji,jj,jk,jpfer) / ( concbfe + trn(ji,jj,jk,jpfer) )
+               zlim4    = trn(ji,jj,jk,jpdoc) / ( xkdoc   + trn(ji,jj,jk,jpdoc) )
+               xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
+               xlimbac (ji,jj,jk) = MIN( zlim1, zlim2, zlim3 ) * zlim4
+
                ! Michaelis-Menten Limitation term for nutrients Small flagellates
                ! -----------------------------------------------
@@ -124,8 +140,4 @@
                xlimphy(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
                !
-               zlim3    = trn(ji,jj,jk,jpfer) / ( concbfe + trn(ji,jj,jk,jpfer) )
-               zlim4    = trn(ji,jj,jk,jpdoc) / ( xkdoc   + trn(ji,jj,jk,jpdoc) )
-               xlimbac(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 ) * zlim4
-
                !   Michaelis-Menten Limitation term for nutrients Diatoms
                !   ----------------------------------------------
@@ -137,5 +149,5 @@
                zlim2    = trn(ji,jj,jk,jppo4) / ( trn(ji,jj,jk,jppo4) + zconc1dnh4  )
                zlim3    = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi(ji,jj) )
-               zratio   = trn(ji,jj,jk,jpdfe)/(trn(ji,jj,jk,jpdia)+rtrn)
+               zratio   = trn(ji,jj,jk,jpdfe) * z1_trndia
                zironmin = xcoef1 * trn(ji,jj,jk,jpdch) * z1_trndia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
                zlim4    = MAX( 0., ( zratio - zironmin ) / qdfelim )
@@ -155,14 +167,14 @@
                   &   / ( concnno3 * concnnh4 + concnnh4 * trn(ji,jj,jk,jpno3) + concnno3 * trn(ji,jj,jk,jpnh4) ) 
                zlim2  = trn(ji,jj,jk,jppo4) / ( trn(ji,jj,jk,jppo4) + concnnh4 )
-               zlim3  = trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concbfe  )
+               zlim3  = trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) +  5.E-11   )
                ztem1  = MAX( 0., tsn(ji,jj,jk,jp_tem) )
                ztem2  = tsn(ji,jj,jk,jp_tem) - 10.
                zetot1 = MAX( 0., etot(ji,jj,jk) - 1.) / ( 4. + etot(ji,jj,jk) ) 
-               zetot2 = 1. / ( 30. + etot(ji,jj,jk) ) 
+               zetot2 = 30. / ( 30. + etot(ji,jj,jk) ) 
 
                xfracal(ji,jj,jk) = caco3r * MIN( zlim1, zlim2, zlim3 )                  &
                   &                       * ztem1 / ( 0.1 + ztem1 )                     &
                   &                       * MAX( 1., trn(ji,jj,jk,jpphy) * 1.e6 / 2. )  &
-                  &                       * 2.325 * zetot1 * 30. * zetot2               &
+                  &                       * zetot1 * zetot2               &
                   &                       * ( 1. + EXP(-ztem2 * ztem2 / 25. ) )         &
                   &                       * MIN( 1., 50. / ( hmld(ji,jj) + rtrn ) )
@@ -173,4 +185,15 @@
       END DO
       !
+      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN        ! save output diagnostics
+        !
+        CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
+        CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        CALL iom_put( "LDnut"  , xlimdia(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
+        CALL iom_put( "LNFe"   , xlimnfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        CALL iom_put( "LDFe"   , xlimdfe(:,:,:) * tmask(:,:,:) )  ! Iron limitation term
+        !
+      ENDIF
+
+      !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lim')
       !
@@ -192,6 +215,6 @@
 
       NAMELIST/nampislim/ concnno3, concdno3, concnnh4, concdnh4, concnfer, concdfer, concbfe,   &
-         &                xsizedia, xsizephy, xsizern, xsizerd, xksi1, xksi2, xkdoc,         & 
-         &                qnfelim, qdfelim, caco3r
+         &                concbno3, concbnh4, xsizedia, xsizephy, xsizern, xsizerd,          & 
+         &                xksi1, xksi2, xkdoc, qnfelim, qdfelim, caco3r
 
       REWIND( numnatp )                     ! read numnat
@@ -203,6 +226,6 @@
          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '    mean rainratio                           caco3r    = ', caco3r
-         WRITE(numout,*) '    NO3 half saturation of nanophyto         concnno3  =  ', concnno3
-         WRITE(numout,*) '    NO3 half saturation of diatoms           concdno3  =  ', concdno3
+         WRITE(numout,*) '    NO3 half saturation of nanophyto         concnno3  = ', concnno3
+         WRITE(numout,*) '    NO3 half saturation of diatoms           concdno3  = ', concdno3
          WRITE(numout,*) '    NH4 half saturation for phyto            concnnh4  = ', concnnh4
          WRITE(numout,*) '    NH4 half saturation for diatoms          concdnh4  = ', concdnh4
@@ -214,4 +237,6 @@
          WRITE(numout,*) '    size ratio for nanophytoplankton         xsizern   = ', xsizern
          WRITE(numout,*) '    size ratio for diatoms                   xsizerd   = ', xsizerd
+         WRITE(numout,*) '    NO3 half saturation of bacteria          concbno3  = ', concbno3
+         WRITE(numout,*) '    NH4 half saturation for bacteria         concbnh4  = ', concbnh4
          WRITE(numout,*) '    Minimum size criteria for diatoms        xsizedia  = ', xsizedia
          WRITE(numout,*) '    Minimum size criteria for nanophyto      xsizephy  = ', xsizephy
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 3446)
@@ -81,11 +81,13 @@
       REAL(wp) :: zgrazn, zgrazpoc, zgraznf, zgrazf
       REAL(wp) :: zgrazfff, zgrazffe
+      REAL(wp) :: zrfact2
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing
       CHARACTER (len=25) :: charout
-      REAL(wp) :: zrfact2
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_meso')
       !
-
+      IF( ln_diatrc .AND. lk_iomput )  CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
+      !
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -145,5 +147,6 @@
 
               ! Total grazing ( grazing by microzoo is already computed in p4zmicro )
-              grazing(ji,jj,jk) = grazing(ji,jj,jk) + zgraztot
+              IF( ln_diatrc .AND. lk_iomput )  zgrazing(ji,jj,jk) = zgraztot
+
               !    Mesozooplankton efficiency
               !    --------------------------
@@ -213,12 +216,8 @@
       END DO
       !
-      IF( ln_diatrc .AND. lk_iomput ) THEN
+      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
          zrfact2 = 1.e3 * rfact2r
-         grazing(:,:,:) = grazing(:,:,:) * zrfact2 * tmask(:,:,:)   ! Total grazing of phyto by zoo
-         prodcal(:,:,:) = prodcal(:,:,:) * zrfact2 * tmask(:,:,:)   ! Calcite production
-         IF( jnt == nrdttrc ) THEN
-            CALL iom_put( "GRAZ" , grazing  )  ! Total grazing of phyto by zooplankton
-            CALL iom_put( "PCAL" , prodcal  )  ! Calcite production
-         ENDIF
+         CALL iom_put( "GRAZ2", zgrazing(:,:,:) * zrfact2 * tmask(:,:,:) )  ! Total grazing of phyto by zooplankton
+         CALL iom_put( "PCAL" , prodcal(:,:,:)  * zrfact2 * tmask(:,:,:) )  ! Calcite production
       ENDIF
       !
@@ -228,4 +227,6 @@
         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
       ENDIF
+      !
+      IF( ln_diatrc .AND. lk_iomput )  CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_meso')
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 3446)
@@ -22,4 +22,5 @@
    USE p4zint          !  interpolation and computation of various fields
    USE p4zprod         !  production
+   USE iom             !  I/O manager
    USE prtctl_trc      !  print control for debugging
 
@@ -29,5 +30,4 @@
    PUBLIC   p4z_micro         ! called in p4zbio.F90
    PUBLIC   p4z_micro_init    ! called in trcsms_pisces.F90
-   PUBLIC   p4z_micro_alloc    ! called in trcsms_pisces.F90
 
    !! * Shared module variables
@@ -59,5 +59,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_micro( kt )
+   SUBROUTINE p4z_micro( kt, jnt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_micro  ***
@@ -67,5 +67,7 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt ! ocean time step
+      INTEGER, INTENT(in) ::  kt  ! ocean time step
+      INTEGER, INTENT(in) ::  jnt 
+      !
       INTEGER  :: ji, jj, jk
       REAL(wp) :: zcompadi, zcompaz , zcompaph, zcompapoc
@@ -77,4 +79,6 @@
       REAL(wp) :: zgrazp, zgrazm, zgrazsd
       REAL(wp) :: zgrazmf, zgrazsf, zgrazpf
+      REAL(wp) :: zrfact2
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zgrazing
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -82,5 +86,6 @@
       IF( nn_timing == 1 )  CALL timing_start('p4z_micro')
       !
-      grazing(:,:,:) = 0.  !: grazing set to zero
+      IF( ln_diatrc .AND. lk_iomput )  CALL wrk_alloc( jpi, jpj, jpk, zgrazing )
+      !
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -127,5 +132,5 @@
 
                ! Grazing by microzooplankton
-               grazing(ji,jj,jk) = grazing(ji,jj,jk) + zgraztot
+               IF( ln_diatrc .AND. lk_iomput )  zgrazing(ji,jj,jk) = zgraztot
 
                !    Various remineralization and excretion terms
@@ -185,4 +190,9 @@
       END DO
       !
+      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
+         zrfact2 = 1.e3 * rfact2r
+         CALL iom_put( "GRAZ1" , zgrazing(:,:,:) * zrfact2 * tmask(:,:,:) )  ! Total grazing of phyto by zooplankton
+      ENDIF
+      !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('micro')")
@@ -190,4 +200,6 @@
          CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
       ENDIF
+      !
+      IF( ln_diatrc .AND. lk_iomput )  CALL wrk_dealloc( jpi, jpj, jpk, zgrazing )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_micro')
@@ -240,13 +252,4 @@
    END SUBROUTINE p4z_micro_init
 
-   INTEGER FUNCTION p4z_micro_alloc()
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE p4z_micro_alloc  ***
-      !!----------------------------------------------------------------------
-      ALLOCATE( grazing(jpi,jpj,jpk), STAT=p4z_micro_alloc )
-      IF( p4z_micro_alloc /= 0 ) CALL ctl_warn('p4z_micro_alloc : failed to allocate arrays.')
-
-   END FUNCTION p4z_micro_alloc
-
 #else
    !!======================================================================
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 3446)
@@ -29,4 +29,5 @@
    REAL(wp), PUBLIC :: wchl   = 0.001_wp  !:
    REAL(wp), PUBLIC :: wchld  = 0.02_wp   !:
+   REAL(wp), PUBLIC :: wchldm = 0.05_wp   !:
    REAL(wp), PUBLIC :: mprat  = 0.01_wp   !:
    REAL(wp), PUBLIC :: mprat2 = 0.01_wp   !:
@@ -150,4 +151,5 @@
       REAL(wp) ::  zfactfe,zfactsi,zfactch, zcompadi
       REAL(wp) ::  zrespp2, ztortp2, zmortp2, zstep
+      REAL(wp) ::  zlim2, zlim1
       CHARACTER (len=25) :: charout
       !!---------------------------------------------------------------------
@@ -177,11 +179,13 @@
                !  Phytoplankton respiration 
                !     ------------------------
-!               zrespp2  = 1.e6 * zstep * (  wchl + wchld * max(0.,2. * ( 0.5 - xlimdia(ji,jj,jk) ) )  )    &
-               zrespp2  = 1.e6 * zstep * (  wchl + wchld * ( 1.0 - xlimdia(ji,jj,jk) ) )      &
-                  &                   * xdiss(ji,jj,jk) * zcompadi * trn(ji,jj,jk,jpdia)
+               zlim2   = xlimdia(ji,jj,jk) * xlimdia(ji,jj,jk)
+               zlim1   = 0.25 * ( 1. - zlim2 ) / ( 0.25 + zlim2 ) 
+               zrespp2 = 1.e6 * zstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trn(ji,jj,jk,jpdia)
+              ! zlim1   = 1.0 - xlimdia(ji,jj,jk) 
+              ! zrespp2 = 1.e6 * zstep * (  wchl + wchld * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trn(ji,jj,jk,jpdia)
 
                !     Phytoplankton mortality. 
                !     ------------------------
-               ztortp2  = mprat2 * zstep * trn(ji,jj,jk,jpdia)  / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi 
+               ztortp2 = mprat2 * zstep * trn(ji,jj,jk,jpdia)  / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi 
 
                zmortp2 = zrespp2 + ztortp2
@@ -235,5 +239,5 @@
       !!----------------------------------------------------------------------
 
-      NAMELIST/nampismort/ wchl, wchld, mprat, mprat2, mpratm
+      NAMELIST/nampismort/ wchl, wchld, wchldm, mprat, mprat2, mpratm
 
       REWIND( numnatp )                     ! read numnatp
@@ -246,4 +250,5 @@
          WRITE(numout,*) '    quadratic mortality of phytoplankton      wchl      =', wchl
          WRITE(numout,*) '    maximum quadratic mortality of diatoms    wchld     =', wchld
+         WRITE(numout,*) '    maximum quadratic mortality of diatoms    wchld     =', wchldm
          WRITE(numout,*) '    phytoplankton mortality rate              mprat     =', mprat
          WRITE(numout,*) '    Diatoms mortality rate                    mprat2    =', mprat2
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 3446)
@@ -77,5 +77,5 @@
       !
       INTEGER  ::   ji, jj, jk
-      REAL(wp) ::   zsilfac, zfact, znanotot, zdiattot, zconctemp, zconctemp2
+      REAL(wp) ::   zsilfac, znanotot, zdiattot, zconctemp, zconctemp2
       REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap
       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zproreg, zproreg2
@@ -135,8 +135,8 @@
                DO ji = 1, jpi
                   IF( etot(ji,jj,jk) > 1.E-3 ) THEN
-                  zval = MAX( 1., zstrn(ji,jj) )
-                  zval = 1.5 * zval / ( 12. + zval )
-                  zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
-                  zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
+                     zval = MAX( 1., zstrn(ji,jj) )
+                     zval = 1.5 * zval / ( 12. + zval )
+                     zprbio(ji,jj,jk) = prmax(ji,jj,jk) * zval
+                     zprdia(ji,jj,jk) = zprbio(ji,jj,jk)
                   ENDIF
                END DO
@@ -156,20 +156,16 @@
 !CDIR NOVERRCHK
                DO ji = 1, jpi
-
                   ! Computation of the P-I slope for nanos and diatoms
                   IF( etot(ji,jj,jk) > 1.E-3 ) THEN
-                      ztn    = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
-                      zadap  = ztn / ( 2.+ ztn )
-
+                      ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
+                      zadap       = ztn / ( 2.+ ztn )
                       zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - xsizedia )
                       zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
-
-                      znanotot = enano(ji,jj,jk) * zstrn(ji,jj)
-                      zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
-
-                      zfact  = EXP( -0.21 * znanotot )
-                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )  &
+                      znanotot    = enano(ji,jj,jk) * zstrn(ji,jj)
+                      zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
+                      !
+                      zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap  * EXP( -0.21 * znanotot ) )  &
                          &                   * trn(ji,jj,jk,jpnch) /( trn(ji,jj,jk,jpphy) * 12. + rtrn)
-
+                      !
                       zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trn(ji,jj,jk,jpdia) + rtrn )   &
                          &                   * trn(ji,jj,jk,jpdch) /( trn(ji,jj,jk,jpdia) * 12. + rtrn)
@@ -201,12 +197,10 @@
                   ! Computation of the P-I slope for nanos and diatoms
                   IF( etot(ji,jj,jk) > 1.E-3 ) THEN
-                      ztn    = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
-                      zadap  = ztn / ( 2.+ ztn )
-
+                      ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
+                      zadap       = ztn / ( 2.+ ztn )
                       zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - xsizedia )
                       zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
-
-                      zfact  = EXP( -0.21 * enano(ji,jj,jk) )
-                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * zfact )
+                      !
+                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -0.21 * enano(ji,jj,jk) ) )
                       zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)
 
@@ -264,8 +258,11 @@
                   zsilim = MIN( zprdia(ji,jj,jk) / ( prmax(ji,jj,jk) + rtrn ), xlimsi(ji,jj,jk) )
                   zsilfac = 4.4 * EXP( -4.23 * zsilim ) * MAX( 0.e0, MIN( 1., 2.2 * ( zlim - 0.5 ) )  ) + 1.e0
-                  zsiborn = MAX( 0.e0, ( trn(ji,jj,jk,jpsil) - 15.e-6 ) )
-                  zsilfac2 = 1.+ 2.* zsiborn / ( zsiborn + xksi2 )
-                  zsilfac = MIN( 5.4, zsilfac * zsilfac2)
-                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac
+                  zsiborn = trn(ji,jj,jk,jpsil) * trn(ji,jj,jk,jpsil) * trn(ji,jj,jk,jpsil)
+                  IF (gphit(ji,jj) < -30 ) THEN
+                    zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
+                  ELSE
+                    zsilfac2 = 1. +      zsiborn / ( zsiborn + xksi2**3 )
+                  ENDIF
+                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
               ENDIF
             END DO
@@ -312,5 +309,5 @@
                   zprofen(ji,jj,jk) = fecnm * prmax(ji,jj,jk)  &
                   &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    &
-                  &             * trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concnfe(ji,jj,jk) )  &
+                  &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  &
                   &             * zmax * trn(ji,jj,jk,jpphy) * rfact2
                   !  production terms for diatomees
@@ -323,5 +320,5 @@
                   zprofed(ji,jj,jk) = fecdm * prmax(ji,jj,jk)  &
                   &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    &
-                  &             * trn(ji,jj,jk,jpfer) / ( trn(ji,jj,jk,jpfer) + concdfe(ji,jj,jk) )  &
+                  &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concdfe(ji,jj,jk) )  &
                   &             * zmax * trn(ji,jj,jk,jpdia) * rfact2
                ENDIF
@@ -426,20 +423,16 @@
          IF( lk_iomput ) THEN
            IF( jnt == nrdttrc ) THEN
-              CALL iom_put( "PPPHY" , zprorca (:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto
-              CALL iom_put( "PPPHY2", zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by diatom
-              CALL iom_put( "PPNEWN", zpronew (:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by nanophyto
-              CALL iom_put( "PPNEWD", zpronewd(:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by diatom
-              CALL iom_put( "PBSi"  , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:) ) ! biogenic silica production
-              CALL iom_put( "PFeD"  , zprofed (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom
-              CALL iom_put( "PFeN"  , zprofen (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by nanophyto
+              CALL iom_put( "PPPHY"  , zprorca (:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by nanophyto
+              CALL iom_put( "PPPHY2" , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) )  ! primary production by diatom
+              CALL iom_put( "PPNEWN" , zpronew (:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by nanophyto
+              CALL iom_put( "PPNEWD" , zpronewd(:,:,:) * zrfact2 * tmask(:,:,:) )  ! new primary production by diatom
+              CALL iom_put( "PBSi"   , zprorcad(:,:,:) * zrfact2 * tmask(:,:,:) * zysopt(:,:,:) ) ! biogenic silica production
+              CALL iom_put( "PFeD"   , zprofed (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by diatom
+              CALL iom_put( "PFeN"   , zprofen (:,:,:) * zrfact2 * tmask(:,:,:) )  ! biogenic iron production by nanophyto
               CALL iom_put( "Mumax"  , prmax(:,:,:) * tmask(:,:,:) )  ! Maximum growth rate
-              CALL iom_put( "MuN"  , zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for nanophyto
-              CALL iom_put( "MuD"  , zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for diatoms
-              CALL iom_put( "LNnut"  , xlimphy (:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-              CALL iom_put( "LDnut"  , xlimdia (:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
-              CALL iom_put( "LNFe"  , xlimnfe (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-              CALL iom_put( "LDFe"  , xlimdfe (:,:,:) * tmask(:,:,:) )  ! Iron limitation term
-              CALL iom_put( "LNlight"  , zprbio (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
-              CALL iom_put( "LDlight"  , zprdia (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
+              CALL iom_put( "MuN"    , zprbio(:,:,:) * xlimphy(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for nanophyto
+              CALL iom_put( "MuD"    , zprdia(:,:,:) * xlimdia(:,:,:) * tmask(:,:,:) )  ! Realized growth rate for diatoms
+              CALL iom_put( "LNlight", zprbio (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
+              CALL iom_put( "LDlight", zprdia (:,:,:) / (prmax(:,:,:) + rtrn) * tmask(:,:,:) )  ! light limitation term
            ENDIF
          ELSE
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 3446)
@@ -72,14 +72,14 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zremip, zremik, zsiremin 
-      REAL(wp) ::   zsatur, zsatur2, znusil, zdep, zfactdep
+      REAL(wp) ::   zsatur, zsatur2, znusil, znusil2, zdep, zdepmin, zfactdep
       REAL(wp) ::   zbactfer, zorem, zorem2, zofer, zolimit
-      REAL(wp) ::   zosil
+      REAL(wp) ::   zosil, ztem
 #if ! defined key_kriest
       REAL(wp) ::   zofer2
 #endif
-      REAL(wp) ::   zonitr, zstep
+      REAL(wp) ::   zonitr, zstep, zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:  ) :: ztempbac 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zdepbac, zolimi, zdepprod
       !!---------------------------------------------------------------------
       !
@@ -87,10 +87,9 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      ztempbac                 )
-      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zolimi )
+      CALL wrk_alloc( jpi, jpj,      ztempbac                  )
+      CALL wrk_alloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
 
       ! Initialisation of temprary arrys
-      zdepbac (:,:,:) = 0._wp
-      zolimi  (:,:,:) = 0._wp
+      zdepprod(:,:,:) = 1._wp
       ztempbac(:,:)   = 0._wp
 
@@ -108,5 +107,7 @@
                   ztempbac(ji,jj)   = zdepbac(ji,jj,jk)
                ELSE
-                  zdepbac(ji,jj,jk) = MIN( 1., zdep / fsdept(ji,jj,jk) ) * ztempbac(ji,jj)
+                  zdepmin = MIN( 1., zdep / fsdept(ji,jj,jk) )
+                  zdepbac (ji,jj,jk) = zdepmin**0.683 * ztempbac(ji,jj)
+                  zdepprod(ji,jj,jk) = zdepmin**0.273
                ENDIF
             END DO
@@ -190,10 +191,7 @@
                ! studies (especially at Papa) have shown this uptake to be significant
                ! ----------------------------------------------------------
-               zbactfer = 20.e-6 *  rfact2 * prmax(ji,jj,jk)                                 &
-                  &              * trn(ji,jj,jk,jpfer) / ( concbfe + trn(ji,jj,jk,jpfer) )    & 
-                  &              * trn(ji,jj,jk,jpfer) / ( 5E-10 + trn(ji,jj,jk,jpfer) )    &
-                  &              * zdepbac(ji,jj,jk)                &
-                  &              * ( 0.5 + SIGN( 0.5, trn(ji,jj,jk,jpfer) -2.e-11 )  )
-
+               zbactfer = 10.e-6 *  rfact2 * prmax(ji,jj,jk) * xlimbacl(ji,jj,jk)             &
+                  &              * trn(ji,jj,jk,jpfer) / ( 2.5E-10 + trn(ji,jj,jk,jpfer) )    &
+                  &              * zdepprod(ji,jj,jk) * zdepbac(ji,jj,jk)
 #if defined key_kriest
                tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zbactfer*0.05
@@ -276,6 +274,8 @@
                zsatur   = ( sio3eq(ji,jj,jk) - trn(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
                zsatur   = MAX( rtrn, zsatur )
-               zsatur2  = zsatur * ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**4
-               znusil   = 0.225  * ( 1. + tsn(ji,jj,jk,jp_tem) / 15.) * zsatur + 0.775 * zsatur2**9.25
+               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
@@ -284,5 +284,6 @@
                zdep     = MAX( hmld(ji,jj), heup(ji,jj) ) 
                zdep     = MAX( 0., fsdept(ji,jj,jk) - zdep )
-               zfactdep = xsilab * EXP(-( xsiremlab - xsirem ) * zdep / wsbio2 )
+               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
                zosil    = zsiremin * trn(ji,jj,jk,jpgsi)
@@ -314,4 +315,9 @@
       END DO
 
+      IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc ) THEN
+          zrfact2 = 1.e3 * rfact2r
+          CALL iom_put( "REMIN" , zolimi(:,:,:) * tmask(:,:,:) * zrfact2 )  ! Remineralisation rate
+          CALL iom_put( "DENIT" , denitr(:,:,:) * tmask(:,:,:) * zrfact2  )  ! Denitrification
+      ENDIF
 
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
@@ -321,6 +327,6 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      ztempbac        )
-      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zolimi )
+      CALL wrk_dealloc( jpi, jpj,      ztempbac                  )
+      CALL wrk_dealloc( jpi, jpj, jpk, zdepbac, zdepprod, zolimi )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_rem')
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 3446)
@@ -183,4 +183,5 @@
       !
       INTEGER  :: ji, jj, jk, jm, ifpr
+      INTEGER  :: ii0, ii1, ij0, ij1
       INTEGER  :: numdust, numsolub, numriv, numiron, numdepo, numhydro
       INTEGER  :: ierr, ierr1, ierr2, ierr3
@@ -425,4 +426,29 @@
             END DO
          END DO
+         IF( cp_cfg == 'orca' .AND. jp_cfg == 2 ) THEN
+            ii0 = 176   ;   ii1 =  176        ! Southern Island : Kerguelen
+            ij0 =  37   ;   ij1 =   37  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 119   ;   ii1 =  119        ! South Georgia
+            ij0 =  29   ;   ij1 =   29  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 111   ;   ii1 =  111        ! Falklands
+            ij0 =  35   ;   ij1 =   35  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 168   ;   ii1 =  168        ! Crozet
+            ij0 =  40   ;   ij1 =   40  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 119   ;   ii1 =  119        ! South Orkney
+            ij0 =  28   ;   ij1 =   28  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 140   ;   ii1 =  140        ! Bouvet Island
+            ij0 =  33   ;   ij1 =   33  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 = 178   ;   ii1 =  178        ! Prince edwards
+            ij0 =  34   ;   ij1 =   34  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+            !
+            ii0 =  43   ;   ii1 =   43        ! Balleny islands
+            ij0 =  21   ;   ij1 =   21  ;   zcmask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1), 1:jpk ) =  0.3_wp  
+         ENDIF
          CALL lbc_lnk( zcmask , 'T', 1. )      ! lateral boundary conditions on cmask   (sign unchanged)
          DO jk = 1, jpk
@@ -458,13 +484,4 @@
          !
          hydrofe(:,:,:) = hydrofe(:,:,:) * hratio / cvol(:,:,:) * ryyss / fse3t(:,:,:)
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  hydrofe(ji,jj,jk)= ( hydrofe(ji,jj,jk) * hratio )  &
-                    &                / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) * ryyss ) / fse3t(ji,jj,jk)
-!!                hydrofe(ji,jj,jk)= ( hydrofe(ji,jj,jk) * hratio  ) / fse3t(ji,jj,jk)
-               END DO
-            END DO
-         END DO
       ENDIF
       ! 
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 3446)
@@ -177,5 +177,5 @@
          !
          IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   &
-            &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * fse3t(:,:,:) * tmask(:,:,:) ) ! iron inputs from sediments
+            &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments
       ENDIF
 
@@ -186,5 +186,5 @@
          !
          IF( ln_diatrc .AND. lk_iomput .AND. jnt == nrdttrc )   &
-            &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * fse3t(:,:,:) * tmask(:,:,:) ) ! hydrothermal iron input
+            &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
       ENDIF
 
@@ -319,6 +319,6 @@
                zfact = zlim * rfact2
 #endif
-               ztrfer = trn(ji,jj,jk,jpfer) / ( concfediaz + trn(ji,jj,jk,jpfer) )
-               ztrpo4 = trn(ji,jj,jk,jppo4) / ( concnnh4   + trn(ji,jj,jk,jppo4) ) 
+               ztrfer = biron(ji,jj,jk)       / ( concfediaz + biron(ji,jj,jk)       )
+               ztrpo4 = trn  (ji,jj,jk,jppo4) / ( concnnh4   + trn  (ji,jj,jk,jppo4) ) 
                znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   &
                  &         *  zfact * MIN( ztrfer, ztrpo4 ) * ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) )
@@ -343,11 +343,7 @@
       IF( ln_diatrc ) THEN
          zfact = 1.e+3 * rfact2r
-         IF( lk_iomput ) THEN
-            zwork1(:,:)  = znitrpot(:,:,1) * nitrfix * zfact * fse3t(:,:,1) * tmask(:,:,1)
-            zwork4(:,:)  = zwork4(:,:) * zfact * tmask(:,:,1)
-            IF( jnt == nrdttrc ) THEN
-               CALL iom_put( "Nfix"  ,  zwork1 )  ! Nitrate reduction in the sediments
-               CALL iom_put( "Sdenit",  zwork4 )  ! Nitrate reduction in the sediments
-            ENDIF
+         IF( lk_iomput .AND. jnt == nrdttrc ) THEN
+            CALL iom_put( "Nfix"  , znitrpot(:,:,:) * nitrfix * rno3 * zfact * tmask(:,:,:) )  ! nitrogen fixation 
+            CALL iom_put( "Sdenit", zwork4(:,:)               * rno3 * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments
          ELSE
             trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * zfact * fse3t(:,:,1) * tmask(:,:,1)
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 3446)
@@ -200,11 +200,11 @@
                !  Part I : Coagulation dependent on turbulence
                zagg1 = 25.9  * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
-               zagg2 = 4424. * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
+               zagg2 = 4452. * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
 
                ! Part II : Differential settling
 
                !  Aggregation of small into large particles
-               zagg3 =  47.2 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
-               zagg4 =  3.3 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
+               zagg3 =  47.1 * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
+               zagg4 =  3.3  * zstep * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
 
                zagg   = zagg1 + zagg2 + zagg3 + zagg4
@@ -215,14 +215,15 @@
                ! 2nd term is shear aggregation of DOC-POC
                ! 3rd term is differential settling of DOC-POC
-               zaggdoc  = ( ( 34.9 * trn(ji,jj,jk,jpdoc) + 102.4 * trn(ji,jj,jk,jppoc) ) *zfact       &
-               &            + 2.4 * zstep * trn(ji,jj,jk,jppoc) ) * trn(ji,jj,jk,jpdoc)
+               zaggdoc  = ( ( 0.369 * 0.3 * trn(ji,jj,jk,jpdoc) + 102.4 * trn(ji,jj,jk,jppoc) ) * zfact       &
+               &            + 2.4 * zstep * trn(ji,jj,jk,jppoc) ) * 0.3 * trn(ji,jj,jk,jpdoc)
 !               zaggdoc  = ( 0.83 * trn(ji,jj,jk,jpdoc) + 271. * trn(ji,jj,jk,jppoc) ) * zfact * trn(ji,jj,jk,jpdoc)
                ! transfer of DOC to GOC : 
                ! 1st term is shear aggregation
                ! 2nd term is differential settling 
-!               zaggdoc2 = (3.3E3 * zfact + 0.1 * zstep ) * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpdoc)
-               zaggdoc2 = 1.07e4 * zfact * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpdoc)
+               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trn(ji,jj,jk,jpgoc) * 0.3 * trn(ji,jj,jk,jpdoc)
+!               zaggdoc2 = 1.07e4 * zfact * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpdoc)
                ! tranfer of DOC to POC due to brownian motion
-               zaggdoc3 =   0.02 * ( 16706. * trn(ji,jj,jk,jppoc) + 231. * trn(ji,jj,jk,jpdoc) ) * zstep * trn(ji,jj,jk,jpdoc)
+!               zaggdoc3 =   0.02 * ( 16706. * trn(ji,jj,jk,jppoc) + 231. * trn(ji,jj,jk,jpdoc) ) * zstep * trn(ji,jj,jk,jpdoc)
+               zaggdoc3 =  ( 5095. * trn(ji,jj,jk,jppoc) + 114. * 0.3 * trn(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trn(ji,jj,jk,jpdoc)
 
                !  Update the trends
Index: branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 3444)
+++ branches/2012/dev_r3438_LOCEAN15_PISLOB/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 3446)
@@ -19,9 +19,11 @@
    INTEGER ::   numnatp
 
-   !!*  Biological fluxes for light
-   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  neln       !: number of T-levels + 1 in the euphotic layer
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup       !: euphotic layer depth
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  etot       !: par (photosynthetic available radiation)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  grazing    !: Total zooplankton grazing
+   !!*  Biological fluxes for light : variables shared by pisces & lobster
+   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  neln  !: number of T-levels + 1 in the euphotic layer
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  heup  !: euphotic layer depth
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  etot  !: par (photosynthetic available radiation)
+   !
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::  xksi  !:  LOBSTER : zooplakton closure
+   !                                                       !:  PISCES  : silicon dependant half saturation
 
 #if defined key_pisces 
@@ -52,5 +54,4 @@
 
    !!*  Biological fluxes for primary production
-   REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:)  ::   xksi       !: ???
    REAL(wp), ALLOCATABLE, SAVE,   DIMENSION(:,:)  ::   xksimax    !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xnanono3   !: ???
@@ -65,4 +66,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimdfe    !: ???
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   xlimsi     !: ???
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)  ::   biron      !: bioavailable fraction of iron
 
 
@@ -71,4 +73,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   nitrfac    !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xlimbac    !: ??
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xlimbacl   !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xdiss      !: ??
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodcal    !: Calcite production
@@ -112,15 +115,13 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
-      INTEGER ::   ierr(7)        ! Local variables
+      INTEGER ::   ierr(6)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
-      !*  Biological fluxes for light
-      ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj),  STAT=ierr(1) )
-      !
-      ALLOCATE( grazing(jpi,jpj,jpk)                        ,     STAT=ierr(2) )
+      !*  Biological fluxes for light : shared variables for pisces & lobster
+      ALLOCATE( etot(jpi,jpj,jpk), neln(jpi,jpj), heup(jpi,jpj), xksi(jpi,jpj), STAT=ierr(1) )
       !
 #if defined key_pisces
       !*  Biological fluxes for primary production
-      ALLOCATE( xksimax(jpi,jpj)     , xksi(jpi,jpj)        ,       &
+      ALLOCATE( xksimax(jpi,jpj)     , biron   (jpi,jpj,jpk),       &
          &      xnanono3(jpi,jpj,jpk), xdiatno3(jpi,jpj,jpk),       &
          &      xnanonh4(jpi,jpj,jpk), xdiatnh4(jpi,jpj,jpk),       &
@@ -128,22 +129,22 @@
          &      xlimnfe (jpi,jpj,jpk), xlimdfe (jpi,jpj,jpk),       &
          &      xlimsi  (jpi,jpj,jpk), concdfe (jpi,jpj,jpk),       &
-         &      concnfe (jpi,jpj,jpk),                           STAT=ierr(3) ) 
+         &      concnfe (jpi,jpj,jpk),                           STAT=ierr(2) ) 
          !
       !*  SMS for the organic matter
       ALLOCATE( xfracal (jpi,jpj,jpk), nitrfac(jpi,jpj,jpk),       &
          &      xlimbac (jpi,jpj,jpk), xdiss  (jpi,jpj,jpk),       & 
-         &      prodcal(jpi,jpj,jpk)                       ,     STAT=ierr(4) ) 
-         !
+         &      xlimbacl(jpi,jpj,jpk), prodcal(jpi,jpj,jpk),     STAT=ierr(3) )
+
       !* Variable for chemistry of the CO2 cycle
       ALLOCATE( akb3(jpi,jpj,jpk)    , ak13  (jpi,jpj,jpk) ,       &
          &      ak23(jpi,jpj,jpk)    , aksp  (jpi,jpj,jpk) ,       &
          &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       &
-         &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,     STAT=ierr(5) )
+         &      hi  (jpi,jpj,jpk)    , excess(jpi,jpj,jpk) ,     STAT=ierr(4) )
          !
       !* Temperature dependancy of SMS terms
-      ALLOCATE( tgfunc(jpi,jpj,jpk)  , tgfunc2(jpi,jpj,jpk) ,    STAT=ierr(6) )
+      ALLOCATE( tgfunc(jpi,jpj,jpk)  , tgfunc2(jpi,jpj,jpk) ,    STAT=ierr(5) )
          !
       !* Array used to indicate negative tracer values  
-      ALLOCATE( xnegtr(jpi,jpj,jpk)  ,                           STAT=ierr(7) )
+      ALLOCATE( xnegtr(jpi,jpj,jpk)  ,                           STAT=ierr(6) )
 #endif
       !
