Index: branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 9451)
+++ branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 9461)
@@ -270,5 +270,5 @@
                !
                ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
-               IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) ! dust in kg/m2/s
+               IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk)
                zlam1b = 3.e-5 + xlamdust * zdust + xlam1 * ztrc
                zscave = zfeequi * zlam1b * xstep
Index: branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90
===================================================================
--- branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 9451)
+++ branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zligand.F90	(revision 9461)
@@ -68,5 +68,5 @@
                ! there is a larger fraction of refractory OM
                zlgwr = max( rlgs , rlgw * exp( -2 * (trb(ji,jj,jk,jplgw)*1e9) ) ) ! years
-               zlgwr = 1. / zlgwr * tgfunc(ji,jj,jk) * ( xstep / nyear_len(1) ) * trb(ji,jj,jk,jplgw)
+               zlgwr = 1. / zlgwr * tgfunc(ji,jj,jk) * ( xstep / nyear_len(1) ) * blim(ji,jj,jk) * trb(ji,jj,jk,jplgw)
                ! photochem loss of weak ligand
                zlgwpr = prlgw * xstep * etot(ji,jj,jk) * trb(ji,jj,jk,jplgw) * (1. - fr_i(ji,jj))
@@ -104,5 +104,5 @@
             IF( iom_use( "LIGREM" ) )   CALL iom_put("LIGREM"   , zligrem(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
             IF( iom_use( "LIGPR" ) )    CALL iom_put("LIGPR"    , zligpr(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
-            IF( iom_use( "LIGPRODR" ) ) CALL iom_put("LPRODR"   , zrligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
+            IF( iom_use( "LPRODR" ) )   CALL iom_put("LPRODR"   , zrligprod(:,:,:) * 1e9 * 1.e+3 * rfact2r * tmask(:,:,:) )
          ENDIF
       ENDIF
Index: branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 9451)
+++ branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 9461)
@@ -228,4 +228,5 @@
                tra(ji,jj,jk,jpbfe) = tra(ji,jj,jk,jpbfe) + zbactfer*0.08
                zfebact(ji,jj,jk)   = zbactfer * 0.33
+               blim(ji,jj,jk)      =  xlimbacl(ji,jj,jk) * zdepbac(ji,jj,jk) / 1.e-6 * zdepprod(ji,jj,jk)
             END DO
          END DO
Index: branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p5zprod.F90
===================================================================
--- branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p5zprod.F90	(revision 9451)
+++ branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p5zprod.F90	(revision 9461)
@@ -87,5 +87,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zprchln, zprchlp, zprchld
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorcan, zprorcap, zprorcad 
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprofed, zprofep, zprofen
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprofed, zprofep, zprofen, zprosid
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpronewn, zpronewp, zpronewd
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zproregn, zproregp, zproregd
@@ -107,5 +107,5 @@
       CALL wrk_alloc( jpi, jpj, jpk, zpislopeadn, zpislopeadp, zpislopeadd, zysopt ) 
       CALL wrk_alloc( jpi, jpj, jpk, zprdia, zprpic, zprbio, zprorcan, zprorcap, zprorcad )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofed, zprofep, zprofen )
+      CALL wrk_alloc( jpi, jpj, jpk, zprofed, zprofep, zprofen, zprosid )
       CALL wrk_alloc( jpi, jpj, jpk, zpronewn, zpronewp, zpronewd, zproregn, zproregp, zproregd )
       CALL wrk_alloc( jpi, jpj, jpk, zpropo4n, zpropo4p, zpropo4d, zrespn, zrespp, zrespd, zprnut )
@@ -122,5 +122,5 @@
       zpropo4n(:,:,:) = 0._wp ; zpropo4p(:,:,:) = 0._wp ; zpropo4d(:,:,:) = 0._wp
       zprdia  (:,:,:) = 0._wp ; zprpic  (:,:,:) = 0._wp ; zprbio  (:,:,:) = 0._wp
-      zysopt  (:,:,:) = 0._wp
+      zysopt  (:,:,:) = 0._wp ; zprosid (:,:,:) = 0._wp
       zrespn  (:,:,:) = 0._wp ; zrespp  (:,:,:) = 0._wp ; zrespd  (:,:,:) = 0._wp 
 
@@ -227,5 +227,4 @@
                   ENDIF
                   zsilim = MIN( zprdia(ji,jj,jk) / ( prmaxd(ji,jj,jk) + rtrn ), xlimsi2(ji,jj,jk) )
-                  zysopt(ji,jj,jk) = grosip * zlim * zsilfac * zsilfac2
                   zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
                   zratiosi = trb(ji,jj,jk,jpndi)
@@ -352,4 +351,6 @@
                   zpropo4d(ji,jj,jk) = zpropmax * xdiatpo4(ji,jj,jk)
                   zprodopd(ji,jj,jk) = zpropmax * xdiatdop(ji,jj,jk)
+                  ! Uptake of Si
+                  zprosid(ji,jj,jk)  = zprnut(ji,jj,jk) * fvduptk(ji,jj,jk) / rno3 * rfact2 * zysopt(ji,jj,jk)
                   ! Uptake of iron
                   zrat = MIN( 1., zratiof / qfdmax )
@@ -440,5 +441,5 @@
               &                     + zprodopd(ji,jj,jk) * texcretd
               tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcretd
-              tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcretd
+              tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprosid(ji,jj,jk)
               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
               &                     + excretp * zprorcap(ji,jj,jk)
@@ -454,5 +455,5 @@
               zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk) + texcretp * zprofep(ji,jj,jk)
               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
-              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcretd * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
+              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - zprosid(ji,jj,jk)
               tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk) - zprorcap(ji,jj,jk)  &
               &                     + zpsino3 * zpronewn(ji,jj,jk) + zpsinh4 * zproregn(ji,jj,jk)   &
@@ -515,5 +516,5 @@
           ENDIF
           IF( iom_use( "PBSi" ) )  THEN
-              zw3d(:,:,:) = zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:) ! biogenic silica production
+              zw3d(:,:,:) = zprosid(:,:,:) * zfact * tmask(:,:,:)   ! biogenic silica production
               CALL iom_put( "PBSi"  , zw3d )
           ENDIF
@@ -589,5 +590,5 @@
       CALL wrk_dealloc( jpi, jpj, jpk, zpislopeadn, zpislopeadp, zpislopeadd, zysopt )                           
       CALL wrk_dealloc( jpi, jpj, jpk, zprdia, zprpic, zprbio, zprorcan, zprorcap, zprorcad )
-      CALL wrk_dealloc( jpi, jpj, jpk, zprofed, zprofep, zprofen ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, zprofed, zprofep, zprofen, zprosid ) 
       CALL wrk_dealloc( jpi, jpj, jpk, zpronewn, zpronewp, zpronewd, zproregn, zproregp, zproregd )
       CALL wrk_dealloc( jpi, jpj, jpk, zpropo4n, zpropo4p, zpropo4d, zrespn, zrespp, zrespd, zprnut )
Index: branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 9451)
+++ branches/CNRS/dev_r8832_PISCO/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 9461)
@@ -106,5 +106,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   prodgoc    !: Calcite production
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   consgoc    !: Calcite production
-
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   blim       !: bacterial production factor
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizen      !: size of diatoms 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sizep      !: size of diatoms 
@@ -163,8 +163,8 @@
          ALLOCATE( xfracal (jpi,jpj,jpk), orem(jpi,jpj,jpk)    ,    &
             &      nitrfac(jpi,jpj,jpk), nitrfac2(jpi,jpj,jpk) ,    &
-            &      prodcal(jpi,jpj,jpk) , xdiss   (jpi,jpj,jpk),    &
+            &      prodcal(jpi,jpj,jpk) , xdiss  (jpi,jpj,jpk),    &
             &      prodpoc(jpi,jpj,jpk) , conspoc(jpi,jpj,jpk) ,    &
             &      prodgoc(jpi,jpj,jpk) , consgoc(jpi,jpj,jpk) ,    &
-            &      prmax  (jpi,jpj,jpk) , STAT=ierr(4) )
+            &      prmax  (jpi,jpj,jpk) , blim   (jpi,jpj,jpk) ,  STAT=ierr(4) )
 
          !* Variable for chemistry of the CO2 cycle
