Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 6966)
@@ -27,5 +27,4 @@
    USE p4zpoc          !  Remineralization of organic particles
    USE p4zagg          !  Aggregation of particles
-   USE p4zlys          !  Dissolution of calcite
    USE p4zfechem       !  Iron chemistry
    USE p4zligand
@@ -82,5 +81,4 @@
       CALL p4z_opt   ( kt, knt )     ! Optic: PAR in the water column
       CALL p4z_sink  ( kt, knt )     ! vertical flux of particulate organic matter
-      CALL p4z_lys   (kt, knt )      ! Dissolution of calcite
       CALL p4z_fechem(kt, knt )      ! Iron chemistry/scavenging
       CALL p4z_lim   ( kt, knt )     ! co-limitations by the various nutrients
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 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zche.F90	(revision 6966)
@@ -23,4 +23,5 @@
    USE sms_pisces    !  PISCES Source Minus Sink variables
    USE lib_mpp       !  MPP library
+   USE eosbn2, ONLY : nn_eos
 
    IMPLICIT NONE
@@ -37,4 +38,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   :: chemo2    ! Solubilities of O2 and CO2
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: fesol    ! solubility of Fe
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tempis   ! In situ temperature
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   salinprac  ! Practical salinity
 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   akb3       !: ???
@@ -154,5 +157,5 @@
       REAL(wp) ::   zpres, ztc  , zcl   , zcpexp, zoxy  , zcpexp2
       REAL(wp) ::   zsqrt, ztr  , zlogt , zcek1, zc1, zplat
-      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt
+      REAL(wp) ::   zis  , zis2 , zsal15, zisqrt, za1, za2
       REAL(wp) ::   zckb , zck1 , zck2  , zckw  , zak1 , zak2  , zakb , zaksp0, zakw
       REAL(wp) ::   zck1p, zck2p, zck3p, zcksi, zak1p, zak2p, zak3p, zaksi
@@ -163,4 +166,31 @@
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_che')
+      !
+      ! Computation of chemical constants require practical salinity
+      ! Thus, when TEOS08 is used, absolute salinity is converted to 
+      ! practical salinity
+      ! -------------------------------------------------------------
+      IF (nn_eos == -1) THEN
+         salinprac(:,:,:) = tsn(:,:,:,jp_sal) * 35.0 / 35.16504
+      ELSE
+         salinprac(:,:,:) = tsn(:,:,:,jp_sal)
+      ENDIF
+
+      !
+      ! Computations of chemical constants require in situ temperature
+      ! Here a quite simple formulation is used to convert 
+      ! potential temperature to in situ temperature. The errors is less than 
+      ! 0.04°C relative to an exact computation
+      ! ---------------------------------------------------------------------
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zpres = fsdept(ji,jj,jk) / 1000.
+               za1 = 0.04 * ( 1.0 + 0.185 * tsn(ji,jj,jk,jp_tem) + 0.035 * (salinprac(ji,jj,jk) - 35.0) )
+               za2 = 0.0075 * ( 1.0 - tsn(ji,jj,jk,jp_tem) / 30.0 )
+               tempis(ji,jj,jk) = tsn(ji,jj,jk,jp_tem) - za1 * zpres + za2 * zpres**2
+            END DO
+         END DO
+      END DO
       !
       ! CHEMICAL CONSTANTS - SURFACE LAYER
@@ -171,7 +201,7 @@
          DO ji = 1, jpi
             !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tsn(ji,jj,1,jp_tem) + 273.15
+            ztkel = tempis(ji,jj,1) + 273.15
             zt    = ztkel * 0.01
-            zsal  = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
+            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
             !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
             !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
@@ -193,8 +223,8 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-              ztkel = tsn(ji,jj,jk,jp_tem) + 273.15
-              zsal  = tsn(ji,jj,jk,jp_sal) + ( 1.- tmask(ji,jj,jk) ) * 35.
+              ztkel = tempis(ji,jj,jk) + 273.15
+              zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
               zsal2 = zsal * zsal
-              ztgg  = LOG( ( 298.15 - tsn(ji,jj,jk,jp_tem) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
+              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
               ztgg2 = ztgg  * ztgg
               ztgg3 = ztgg2 * ztgg
@@ -229,6 +259,6 @@
 
                ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tsn(ji,jj,jk,jp_tem) + 273.15
-               zsal    = tsn(ji,jj,jk,jp_sal) + ( 1.-tmask(ji,jj,jk) ) * 35.
+               ztkel   = tempis(ji,jj,jk) + 273.15
+               zsal    = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
                zsqrt  = SQRT( zsal )
                zsal15  = zsqrt * zsal
@@ -238,5 +268,5 @@
                zis2   = zis * zis
                zisqrt = SQRT( zis )
-               ztc     = tsn(ji,jj,jk,jp_tem) + ( 1.- tmask(ji,jj,jk) ) * 20.
+               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
 
                ! CHLORINITY (WOOSTER ET AL., 1969)
@@ -414,10 +444,10 @@
 
                ! Liu and Millero (1999) only valid 5 - 50 degC
-               ztkel1 = MAX( 5. , tsn(ji,jj,jk,jp_tem) ) + 273.16
-               fesol(ji,jj,jk,1) = 10**((-13.486) - (0.1856* (zis**0.5)) + (0.3073*zis) + (5254/ztkel1))
-               fesol(ji,jj,jk,2) = 10**(2.517 - (0.885*(zis**0.5)) + (0.2139 * zis) - (1320/ztkel1) )
-               fesol(ji,jj,jk,3) = 10**(0.4511 - (0.3305*(ZIS**0.5)) - (1996/ztkel1) )
-               fesol(ji,jj,jk,4) = 10**(-0.2965 - (0.7881*(zis**0.5)) - (4086/ztkel1) )
-               fesol(ji,jj,jk,5) = 10**(4.4466 - (0.8505*(zis**0.5)) - (7980/ztkel1) )
+               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) )
             END DO
          END DO
@@ -794,10 +824,11 @@
       ALLOCATE( sio3eq(jpi,jpj,jpk), fekeq(jpi,jpj,jpk), chemc(jpi,jpj,3), chemo2(jpi,jpj,jpk), STAT=ierr(1) )
 
-      ALLOCATE( akb3(jpi,jpj,jpk)    ,                             &
-         &      akw3(jpi,jpj,jpk)    , borat (jpi,jpj,jpk) ,       &
-         &      aks3(jpi,jpj,jpk)    , akf3(jpi,jpj,jpk)   ,       &
-         &      ak1p3(jpi,jpj,jpk)   , ak2p3(jpi,jpj,jpk)  ,       &
-         &      ak3p3(jpi,jpj,jpk)   , aksi3(jpi,jpj,jpk)  ,       &
-         &      fluorid(jpi,jpj,jpk) , sulfat(jpi,jpj,jpk) ,     STAT=ierr(2) )
+      ALLOCATE( akb3(jpi,jpj,jpk)     , tempis(jpi, jpj, jpk),       &
+         &      akw3(jpi,jpj,jpk)     , borat (jpi,jpj,jpk)  ,       &
+         &      aks3(jpi,jpj,jpk)     , akf3(jpi,jpj,jpk)    ,       &
+         &      ak1p3(jpi,jpj,jpk)    , ak2p3(jpi,jpj,jpk)   ,       &
+         &      ak3p3(jpi,jpj,jpk)    , aksi3(jpi,jpj,jpk)   ,       &
+         &      fluorid(jpi,jpj,jpk)  , sulfat(jpi,jpj,jpk)  ,       &
+         &      salinprac(jpi,jpj,jpk),                 STAT=ierr(2) )
 
       ALLOCATE( fesol(jpi,jpj,jpk,5), STAT=ierr(3) )
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 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 6966)
@@ -163,10 +163,9 @@
                   zTL2(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 )
                   zph    = -LOG10( MAX( hi(ji,jj,jk), rtrn) )
-                  zoxy   = trb(ji,jj,jk,jpoxy) * ( rhop(ji,jj,jk) / 1.e3 )
+                  zoxy   = trb(ji,jj,jk,jpoxy)
                   ! Fe2+ oxydation rate from Santana-Casiano et al. (2005)
-                  zkox   = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tsn(ji,jj,jk,jp_tem) + 273.15 )  &
-                    &    - 0.04406 * SQRT( tsn(ji,jj,jk,jp_sal) ) - 0.002847 * tsn(ji,jj,jk,jp_sal)
+                  zkox   = 35.407 - 6.7109 * zph + 0.5342 * zph * zph - 5362.6 / ( tempis(ji,jj,jk) + 273.15 )  &
+                    &    - 0.04406 * SQRT( salinprac(ji,jj,jk) ) - 0.002847 * salinprac(ji,jj,jk)
                   zkox   = ( 10.** zkox ) * spd
                   zkox   = zkox * MAX( 1.e-6, zoxy) / ( chemo2(ji,jj,jk) + rtrn )
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 6966)
@@ -128,8 +128,8 @@
             ! DUMMY VARIABLES FOR DIC, H+, AND BORATE
             zfact = rhop(ji,jj,1) / 1000. + rtrn
-            zdic  = trb(ji,jj,1,jpdic) / zfact
+            zdic  = trb(ji,jj,1,jpdic) 
             zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
             ! CALCULATE [H2CO3]
-            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)*zfact
+            zh2co3(ji,jj) = zdic/(1. + ak13(ji,jj,1)/zph + ak13(ji,jj,1)*ak23(ji,jj,1)/zph**2)
          END DO
       END DO
@@ -146,5 +146,5 @@
 !CDIR NOVERRCHK
          DO ji = 1, jpi
-            ztc  = MIN( 35., tsn(ji,jj,1,jp_tem) )
+            ztc  = MIN( 35., tempis(ji,jj,1) )
             ztc2 = ztc * ztc
             ztc3 = ztc * ztc2 
@@ -169,6 +169,6 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ztkel = tsn(ji,jj,1,jp_tem) + 273.15
-            zsal  = tsn(ji,jj,1,jp_sal) + ( 1.- tmask(ji,jj,1) ) * 35.
+            ztkel = tempis(ji,jj,1) + 273.15
+            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
             zvapsw = exp(24.4543 - 67.4509*(100.0/ztkel) - 4.8489*log(ztkel/100) - 0.000544*zsal)
             zpco2atm(ji,jj) = satmco2(ji,jj) * ( patm(ji,jj) - zvapsw )
@@ -232,5 +232,5 @@
             trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1)
             trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1)
-            trc2d(:,:,jp_pcs0_2d + 3) = ( satmco2(:,:) * patm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
+            trc2d(:,:,jp_pcs0_2d + 3) = ( zpco2atm(:,:) - zh2co3(:,:) / ( chemc(:,:,1) + rtrn ) ) * tmask(:,:,1)
          ENDIF
       ENDIF
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 6966)
@@ -66,10 +66,10 @@
       REAL(wp) ::   zomegaca, zexcess, zexcess0
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss, zhinit, zhi
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss, zhinit, zhi, zco3sat
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_lys')
       !
-      CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi )
+      CALL wrk_alloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat )
       !
       zco3    (:,:,:) = 0.
@@ -85,8 +85,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zfact = rhop(ji,jj,jk) / 1000. + rtrn
                zco3(ji,jj,jk) = trb(ji,jj,jk,jpdic) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
                &                + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
-               hi(ji,jj,jk)   = zhi(ji,jj,jk) * zfact
+               hi(ji,jj,jk)   = zhi(ji,jj,jk) * rhop(ji,jj,jk) / 1000.
             END DO
          END DO
@@ -105,7 +104,8 @@
                ! DEVIATION OF [CO3--] FROM SATURATION VALUE
                ! Salinity dependance in zomegaca and divide by rhop/1000 to have good units
-               zcalcon  = calcon * ( tsn(ji,jj,jk,jp_sal) / 35._wp )
+               zcalcon  = calcon * ( salinprac(ji,jj,jk) / 35._wp )
                zfact    = rhop(ji,jj,jk) / 1000._wp
-               zomegaca = ( zcalcon * zco3(ji,jj,jk) * zfact ) / aksp(ji,jj,jk) 
+               zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zfact + rtrn )
+               zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zfact / ( zcalcon + rtrn )
 
                ! SET DEGREE OF UNDER-/SUPERSATURATION
@@ -135,6 +135,6 @@
          IF( iom_use( "PH"     ) ) CALL iom_put( "PH"    , -1. * LOG10( hi(:,:,:) )          * tmask(:,:,:) )
          IF( iom_use( "CO3"    ) ) CALL iom_put( "CO3"   , zco3(:,:,:) * 1.e+3               * tmask(:,:,:) )
-         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", aksp(:,:,:) * 1.e+3 / calcon      * tmask(:,:,:) )
-         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) )
+         IF( iom_use( "CO3sat" ) ) CALL iom_put( "CO3sat", zco3sat(:,:,:) * 1.e+3            * tmask(:,:,:) )
+         IF( iom_use( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r * tmask(:,:,:) )
 
       ELSE
@@ -142,5 +142,5 @@
             trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:)
             trc3d(:,:,:,jp_pcs0_3d + 1) = zco3(:,:,:)              * tmask(:,:,:)
-            trc3d(:,:,:,jp_pcs0_3d + 2) = aksp(:,:,:) / calcon     * tmask(:,:,:)
+            trc3d(:,:,:,jp_pcs0_3d + 2) = zco3sat(:,:,:)           * tmask(:,:,:)
          ENDIF
       ENDIF
@@ -153,5 +153,5 @@
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi )
+      CALL wrk_dealloc( jpi, jpj, jpk, zco3, zcaldiss, zhinit, zhi, zco3sat )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_lys')
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 6966)
@@ -125,9 +125,7 @@
       !                                        !  --------------------------------------
       IF( l_trcdm2dc ) THEN                     !  diurnal cycle
-         ! 1% of qsr to compute euphotic layer
-         zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr
          !
          zqsr_corr(:,:) = qsr_mean(:,:) / ( 1. - fr_i(:,:) + rtrn )
-         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pe100 = zqsr100 )
          !
          DO jk = 1, nksrp      
@@ -148,9 +146,7 @@
          !
       ELSE
-         ! 1% of qsr to compute euphotic layer
-         zqsr100(:,:) = 0.01 * qsr(:,:)
          !
          zqsr_corr(:,:) = qsr(:,:) / ( 1. - fr_i(:,:) + rtrn )
-         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3 )
+         CALL p4z_opt_par( kt, zqsr_corr, ze1, ze2, ze3, pe100 = zqsr100 )
          !
          DO jk = 1, nksrp      
@@ -184,5 +180,5 @@
          DO jj = 1, jpj
            DO ji = 1, jpi
-              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * zqsr100(ji,jj) )  THEN
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= zqsr100(ji,jj) )  THEN
                  neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
                  !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
@@ -190,6 +186,5 @@
               ENDIF
               IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.50 )  THEN
-                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
-                 heup_01(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
+                 heup_01(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth (light level definition)
               ENDIF
            END DO
@@ -273,5 +268,5 @@
    END SUBROUTINE p4z_opt
 
-   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 ) 
+   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0, pe100 ) 
       !!----------------------------------------------------------------------
       !!                  ***  routine p4z_opt_par  ***
@@ -285,5 +280,6 @@
       REAL(wp), DIMENSION(jpi,jpj)    , INTENT(in)              ::  pqsr          !   shortwave
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)           ::  pe1 , pe2 , pe3   !  PAR ( R-G-B)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0  
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL ::  pe0
+      REAL(wp), DIMENSION(jpi,jpj),     INTENT(inout), OPTIONAL ::  pe100
       !! * local variables
       INTEGER    ::   ji, jj, jk     ! dummy loop indices
@@ -294,4 +290,8 @@
       IF( ln_varpar ) THEN  ;  zqsr(:,:) = par_varsw(:,:) * pqsr(:,:)
       ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:)
+      ENDIF
+      !
+      IF( PRESENT( pe100 ) ) THEN
+         pe100(:,:) = 3. * 0.01 * zqsr(:,:)
       ENDIF
       !
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 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zpoc.F90	(revision 6966)
@@ -63,5 +63,5 @@
       !
       INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zremip, zremig, zdep, zorem2, zofer
+      REAL(wp) ::   zremip, zremig, zdep, zorem, zorem2, zofer
       REAL(wp) ::   zsizek, zsizek1, alphat, remint
       REAL(wp) ::   solgoc, zpoc
@@ -220,4 +220,5 @@
                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
@@ -389,5 +390,5 @@
                  ! --------------------------------------------------------
                  zremip          = zremipoc(ji,jj,jk) * zstep * tgfunc(ji,jj,jk)
-                 orem(ji,jj,jk)  = zremip * trb(ji,jj,jk,jppoc)
+                 zorem           = zremip * trb(ji,jj,jk,jppoc)
                  zofer           = zremip * trb(ji,jj,jk,jpsfe)
 #if defined key_kriest
@@ -398,7 +399,8 @@
                  ! -------------------------------------
 
-                 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + orem(ji,jj,jk)
+                 tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zorem
+                 orem(ji,jj,jk)      = orem(ji,jj,jk) + zorem
                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + zofer
-                 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - orem(ji,jj,jk)
+                 tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zorem
 #if defined key_kriest
                  tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zorem2
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 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 6966)
@@ -33,9 +33,9 @@
    !! * Shared module variables
    LOGICAL , PUBLIC ::  ln_newprod      !:
-   REAL(wp), PUBLIC ::  pislope         !:
-   REAL(wp), PUBLIC ::  pislope2        !:
+   REAL(wp), PUBLIC ::  pislopen         !:
+   REAL(wp), PUBLIC ::  pisloped        !:
    REAL(wp), PUBLIC ::  xadap           !:
-   REAL(wp), PUBLIC ::  excret          !:
-   REAL(wp), PUBLIC ::  excret2         !:
+   REAL(wp), PUBLIC ::  excretn          !:
+   REAL(wp), PUBLIC ::  excretd         !:
    REAL(wp), PUBLIC ::  bresp           !:
    REAL(wp), PUBLIC ::  chlcnm          !:
@@ -51,6 +51,6 @@
    
    REAL(wp) :: r1_rday                !: 1 / rday
-   REAL(wp) :: texcret                !: 1 - excret 
-   REAL(wp) :: texcret2               !: 1 - excret2        
+   REAL(wp) :: texcretn               !: 1 - excretn 
+   REAL(wp) :: texcretd               !: 1 - excretd        
 
 
@@ -78,14 +78,13 @@
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zsilfac, znanotot, zdiattot, zconctemp, zconctemp2
-      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap
-      REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zproreg, zproreg2
-      REAL(wp) ::   zmxltst, zmxlday, zmaxday, zdocprod
-      REAL(wp) ::   zpislopen  , zpislope2n
+      REAL(wp) ::   zratio, zmax, zsilim, ztn, zadap, zlim, zsilfac2, zsiborn
+      REAL(wp) ::   zprod, zproreg, zproreg2, zprochln, zprochld
+      REAL(wp) ::   zmaxday, zdocprod, zpislopen, zpisloped
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
       REAL(wp) ::   zfact
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zstrn, zw2d
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt, zw3d   
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopeadn, zpislopeadd, zprdia, zprbio, zprdch, zprnch, zysopt, zw3d   
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zprorcan, zprorcad, zprofed, zprofen, zpronewn, zpronewd
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxl_fac, zmxl_chl
       !!---------------------------------------------------------------------
@@ -95,21 +94,12 @@
       !  Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj,      zstrn )
-      CALL wrk_alloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zpislopeadn, zpislopeadd, zprdia, zprbio, zprdch, zprnch, zysopt ) 
       CALL wrk_alloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
-      CALL wrk_alloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd )
-      !
-      zprorca (:,:,:) = 0._wp
-      zprorcad(:,:,:) = 0._wp
-      zprofed (:,:,:) = 0._wp
-      zprofen (:,:,:) = 0._wp
-      zprochln(:,:,:) = 0._wp
-      zprochld(:,:,:) = 0._wp
-      zpronew (:,:,:) = 0._wp
-      zpronewd(:,:,:) = 0._wp
-      zprdia  (:,:,:) = 0._wp
-      zprbio  (:,:,:) = 0._wp
-      zprdch  (:,:,:) = 0._wp
-      zprnch  (:,:,:) = 0._wp
-      zysopt  (:,:,:) = 0._wp
+      CALL wrk_alloc( jpi, jpj, jpk, zprorcan, zprorcad, zprofed, zprofen, zpronewn, zpronewd )
+      !
+      zprorcan(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp ; zprofed (:,:,:) = 0._wp
+      zprofen (:,:,:) = 0._wp ; zysopt  (:,:,:) = 0._wp
+      zpronewn(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp ; zprdia  (:,:,:) = 0._wp
+      zprbio  (:,:,:) = 0._wp ; zprdch  (:,:,:) = 0._wp ; zprnch  (:,:,:) = 0._wp 
 
       ! Computation of the optimal production
@@ -144,6 +134,6 @@
                      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) ) ) * ( 1. - fr_i(ji,jj) )
-                  zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+                  zmxl_fac(ji,jj,jk) = ( 1. - exp( -0.2 * zmxl_fac(ji,jj,jk) ) )
+                  zmxl_chl(ji,jj,jk) = zmxl_chl(ji,jj,jk)
                ENDIF
             END DO
@@ -166,8 +156,8 @@
                   zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
                   !
-                  zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
+                  zpislopeadn(ji,jj,jk) = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
                   &                   * trb(ji,jj,jk,jpnch) /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
                   !
-                  zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn )   &
+                  zpislopeadd(ji,jj,jk) = (pislopen * zconctemp2 + pisloped * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn )   &
                   &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
                ENDIF
@@ -183,16 +173,16 @@
                       ! Computation of production function for Carbon
                       !  ---------------------------------------------
-                      zpislopen  = zpislopead (ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+                      zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
                       &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
-                      zpislope2n = zpislopead2(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+                      zpisloped = zpislopeadd(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
                       &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
-                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) )  )
-                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) )  )
+                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
+                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
                       !  Computation of production function for Chlorophyll
                       !--------------------------------------------------
-                      zpislopen  = zpislopead (ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-                      zpislope2n = zpislopead2(ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                      zpislopen = zpislopeadn(ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+                      zpisloped = zpislopeadd(ji,jj,jk) / ( prmax(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
+                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
                   ENDIF
                END DO
@@ -206,14 +196,14 @@
                       ! Computation of production function for Carbon
                       !  ---------------------------------------------
-                      zpislopen =  zpislopead(ji,jj,jk)  / ( zprbio(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
-                      zpislope2n = zpislopead2(ji,jj,jk) / ( zprdia(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
-                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
-                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                      zpislopen = zpislopeadn(ji,jj,jk)  / ( zprbio(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
+                      zpisloped = zpislopeadd(ji,jj,jk) / ( zprdia(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
+                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
+                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
                       !  Computation of production function for Chlorophyll
                       !--------------------------------------------------
-                      zpislopen = zpislopen   * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-                      zpislope2n = zpislope2n * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                      zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+                      zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) ) )
+                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) ) )
                   ENDIF
                END DO
@@ -221,5 +211,4 @@
          END DO
       ENDIF
-
 
       !  Computation of a proxy of the N/C ratio
@@ -264,4 +253,12 @@
       END DO
 
+      !  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
+
+
       ! Computation of the various production terms 
 !CDIR NOVERRCHK
@@ -272,23 +269,21 @@
             DO ji = 1, jpi
                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-                  !  production terms for nanophyto.
-                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
-                  zpronew(ji,jj,jk) = zprorca(ji,jj,jk) * xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
+                  !  production terms for nanophyto. (C)
+                  zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
+                  zpronewn(ji,jj,jk)  = zprorcan(ji,jj,jk)* xnanono3(ji,jj,jk) / ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) + rtrn )
                   !
-                  zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) + rtrn )
-                  zratio = zratio / fecnm 
+                  zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) * fecnm + rtrn )
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
-                  zprofen(ji,jj,jk) = fecnm * prmax(ji,jj,jk)  &
+                  zprofen(ji,jj,jk) = fecnm * prmax(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
                   &             * ( 4. - 4.5 * xlimnfe(ji,jj,jk) / ( xlimnfe(ji,jj,jk) + 0.5 ) )    &
                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concnfe(ji,jj,jk) )  &
                   &             * zmax * trb(ji,jj,jk,jpphy) * rfact2
-                  !  production terms for diatomees
+                  !  production terms for diatoms (C)
                   zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
                   zpronewd(ji,jj,jk) = zprorcad(ji,jj,jk) * xdiatno3(ji,jj,jk) / ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) + rtrn )
                   !
-                  zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn )
-                  zratio = zratio / fecdm 
+                  zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) * fecdm + rtrn )
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
-                  zprofed(ji,jj,jk) = fecdm * prmax(ji,jj,jk)  &
+                  zprofed(ji,jj,jk) = fecdm * prmax(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
                   &             * ( 4. - 4.5 * xlimdfe(ji,jj,jk) / ( xlimdfe(ji,jj,jk) + 0.5 ) )    &
                   &             * biron(ji,jj,jk) / ( biron(ji,jj,jk) + concdfe(ji,jj,jk) )  &
@@ -299,4 +294,5 @@
       END DO
 
+      ! Computation of the chlorophyll production terms
       DO jk = 1, jpkm1
          DO jj = 1, jpj
@@ -305,14 +301,17 @@
                   !  production terms for nanophyto. ( chlorophyll )
                   znanotot = enano(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
-                  zprod    = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * xlimphy(ji,jj,jk)
-                  zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
-                  zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 12. * zprod / &
-                                        & (  zpislopead(ji,jj,jk) * znanotot +rtrn)
-                  !  production terms for diatomees ( chlorophyll )
+                  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 / &
+                                        & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
+                  !  production terms for diatoms ( chlorophyll )
                   zdiattot = ediat(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
                   zprod    = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * xlimdia(ji,jj,jk)
-                  zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
-                  zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 12. * zprod / &
-                  & ( zpislopead2(ji,jj,jk) * zdiattot +rtrn )
+                  zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
+                  zprochld = zprochld + (chlcdm-chlcmin) * 12. * zprod / &
+                                        & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
+                  !   Update the arrays TRA which contain the Chla sources and sinks
+                  tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln * texcretn
+                  tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld * texcretd
                ENDIF
             END DO
@@ -324,30 +323,30 @@
          DO jj = 1, jpj
            DO ji =1 ,jpi
-              zproreg  = zprorca(ji,jj,jk) - zpronew(ji,jj,jk)
-              zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk)
-              zdocprod = excret2 * zprorcad(ji,jj,jk) + excret * zprorca(ji,jj,jk)
-              tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
-              tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronew(ji,jj,jk) - zpronewd(ji,jj,jk)
-              tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproreg - zproreg2
-              tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorca(ji,jj,jk) * texcret
-              tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln(ji,jj,jk) * texcret
-              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcret
-              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zdocprod
-              zfeup = texcret * zprofen(ji,jj,jk) + texcret2 * zprofed(ji,jj,jk)
+              IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+                 zproreg  = zprorcan(ji,jj,jk) - zpronewn(ji,jj,jk)
+                 zproreg2 = zprorcad(ji,jj,jk) - zpronewd(ji,jj,jk)
+                 zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
+                 tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk)
+                 tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) - zpronewn(ji,jj,jk) - zpronewd(ji,jj,jk)
+                 tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproreg - zproreg2
+                 tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorcan(ji,jj,jk) * texcretn
+                 tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcretn
+                 tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(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,jpdoc) = tra(ji,jj,jk,jpdoc) + zdocprod
+                 zfeup = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
 #if defined key_ligand
-              tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp  - zfeup * plig(ji,jj,jk) * lthet
+                 tra(ji,jj,jk,jplgw) = tra(ji,jj,jk,jplgw) + zdocprod * ldocp  - zfeup * plig(ji,jj,jk) * lthet
 #endif
-              tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
-                 &                + ( o2ut + o2nit ) * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) )
-              tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) - zfeup
-              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
-              tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) - zprorca(ji,jj,jk) - zprorcad(ji,jj,jk)
-              tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronew(ji,jj,jk) + zpronewd(ji,jj,jk) ) &
-                 &                                      - rno3 * ( zproreg + zproreg2 )
-          END DO
+                 tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproreg + zproreg2) &
+                 &                   + ( o2ut + o2nit ) * ( zpronewn(ji,jj,jk) + zpronewd(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,jpdic) = tra(ji,jj,jk,jpdic) - zprorcan(ji,jj,jk) - zprorcad(ji,jj,jk)
+                 tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3 * ( zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk) ) &
+                 &                                         - rno3 * ( zproreg + zproreg2 )
+              ENDIF
+           END DO
         END DO
      END DO
@@ -356,5 +355,5 @@
     ! Total primary production per year
     IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
-         & tpp = glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
+         & tpp = glob_sum( ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
 
     IF( lk_iomput ) THEN
@@ -365,5 +364,5 @@
           !
           IF( iom_use( "PPPHY" ) .OR. iom_use( "PPPHY2" ) )  THEN
-              zw3d(:,:,:) = zprorca (:,:,:) * zfact * tmask(:,:,:)  ! primary production by nanophyto
+              zw3d(:,:,:) = zprorcan(:,:,:) * zfact * tmask(:,:,:)  ! primary production by nanophyto
               CALL iom_put( "PPPHY"  , zw3d )
               !
@@ -372,5 +371,5 @@
           ENDIF
           IF( iom_use( "PPNEWN" ) .OR. iom_use( "PPNEWD" ) )  THEN
-              zw3d(:,:,:) = zpronew (:,:,:) * zfact * tmask(:,:,:)  ! new primary production by nanophyto
+              zw3d(:,:,:) = zpronewn(:,:,:) * zfact * tmask(:,:,:)  ! new primary production by nanophyto
               CALL iom_put( "PPNEWN"  , zw3d )
               !
@@ -408,9 +407,9 @@
           ENDIF
           IF( iom_use( "TPP" ) )  THEN
-              zw3d(:,:,:) = ( zprorca(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  ! total primary production
+              zw3d(:,:,:) = ( zprorcan(:,:,:) + zprorcad(:,:,:) ) * zfact * tmask(:,:,:)  ! total primary production
               CALL iom_put( "TPP"  , zw3d )
           ENDIF
           IF( iom_use( "TPNEW" ) )  THEN
-              zw3d(:,:,:) = ( zpronew(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ! total new production
+              zw3d(:,:,:) = ( zpronewn(:,:,:) + zpronewd(:,:,:) ) * zfact * tmask(:,:,:)  ! total new production
               CALL iom_put( "TPNEW"  , zw3d )
           ENDIF
@@ -422,5 +421,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-               zw2d(:,:) = zw2d(:,:) + zprorca (:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated  primary produc. by nano
+               zw2d(:,:) = zw2d(:,:) + zprorcan(:,:,jk) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated  primary produc. by nano
              ENDDO
              CALL iom_put( "INTPPPHY" , zw2d )
@@ -435,5 +434,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-                zw2d(:,:) = zw2d(:,:) + ( zprorca(:,:,jk) + zprorcad(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated pp
+                zw2d(:,:) = zw2d(:,:) + ( zprorcan(:,:,jk) + zprorcad(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk) ! vert. integrated pp
              ENDDO
              CALL iom_put( "INTPP" , zw2d )
@@ -442,5 +441,5 @@
              zw2d(:,:) = 0.
              DO jk = 1, jpkm1
-                zw2d(:,:) = zw2d(:,:) + ( zpronew(:,:,jk) + zpronewd(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated new prod
+                zw2d(:,:) = zw2d(:,:) + ( zpronewn(:,:,jk) + zpronewd(:,:,jk) ) * fse3t(:,:,jk) * zfact * tmask(:,:,jk)  ! vert. integrated new prod
              ENDDO
              CALL iom_put( "INTPNEW" , zw2d )
@@ -468,7 +467,7 @@
         IF( ln_diatrc ) THEN
            zfact = 1.e+3 * rfact2r
-           trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorca (:,:,:) * zfact * tmask(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 4)  = zprorcan(:,:,:) * zfact * tmask(:,:,:)
            trc3d(:,:,:,jp_pcs0_3d + 5)  = zprorcad(:,:,:) * zfact * tmask(:,:,:)
-           trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronew (:,:,:) * zfact * tmask(:,:,:)
+           trc3d(:,:,:,jp_pcs0_3d + 6)  = zpronewn(:,:,:) * zfact * tmask(:,:,:)
            trc3d(:,:,:,jp_pcs0_3d + 7)  = zpronewd(:,:,:) * zfact * tmask(:,:,:)
            trc3d(:,:,:,jp_pcs0_3d + 8)  = zprorcad(:,:,:) * zfact * tmask(:,:,:) * zysopt(:,:,:)
@@ -487,7 +486,7 @@
      !
      CALL wrk_dealloc( jpi, jpj,      zstrn )
-     CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopead2, zprdia, zprbio, zprdch, zprnch, zysopt ) 
+     CALL wrk_dealloc( jpi, jpj, jpk, zpislopeadn, zpislopeadd, zprdia, zprbio, zprdch, zprnch, zysopt ) 
      CALL wrk_dealloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
-     CALL wrk_dealloc( jpi, jpj, jpk, zprorca, zprorcad, zprofed, zprofen, zprochln, zprochld, zpronew, zpronewd )
+     CALL wrk_dealloc( jpi, jpj, jpk, zprorcan, zprorcad, zprofed, zprofen, zpronewn, zpronewd )
      !
      IF( nn_timing == 1 )  CALL timing_stop('p4z_prod')
@@ -508,5 +507,5 @@
       !!----------------------------------------------------------------------
       !
-      NAMELIST/nampisprod/ pislope, pislope2, xadap, ln_newprod, bresp, excret, excret2,  &
+      NAMELIST/nampisprod/ pislopen, pisloped, xadap, ln_newprod, bresp, excretn, excretd,  &
          &                 chlcnm, chlcdm, chlcmin, fecnm, fecdm, grosip
       INTEGER :: ios                 ! Local integer output status for namelist read
@@ -526,15 +525,15 @@
          WRITE(numout,*) ' Namelist parameters for phytoplankton growth, nampisprod'
          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '    Enable new parame. of production (T/F)   ln_newprod   =', ln_newprod
+         WRITE(numout,*) '    Enable new parame. of production (T/F)   ln_newprod    =', ln_newprod
          WRITE(numout,*) '    mean Si/C ratio                           grosip       =', grosip
-         WRITE(numout,*) '    P-I slope                                 pislope      =', pislope
-         WRITE(numout,*) '    Acclimation factor to low light           xadap       =', xadap
-         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excret       =', excret
-         WRITE(numout,*) '    excretion ratio of diatoms                excret2      =', excret2
+         WRITE(numout,*) '    P-I slope                                 pislopen     =', pislopen
+         WRITE(numout,*) '    Acclimation factor to low light           xadap        =', xadap
+         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excretn      =', excretn
+         WRITE(numout,*) '    excretion ratio of diatoms                excretd      =', excretd
          IF( ln_newprod )  THEN
             WRITE(numout,*) '    basal respiration in phytoplankton        bresp        =', bresp
             WRITE(numout,*) '    Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
          ENDIF
-         WRITE(numout,*) '    P-I slope  for diatoms                    pislope2     =', pislope2
+         WRITE(numout,*) '    P-I slope  for diatoms                    pisloped     =', pisloped
          WRITE(numout,*) '    Minimum Chl/C in nanophytoplankton        chlcnm       =', chlcnm
          WRITE(numout,*) '    Minimum Chl/C in diatoms                  chlcdm       =', chlcdm
@@ -544,6 +543,6 @@
       !
       r1_rday   = 1._wp / rday 
-      texcret   = 1._wp - excret
-      texcret2  = 1._wp - excret2
+      texcretn  = 1._wp - excretn
+      texcretd  = 1._wp - excretd
       tpp       = 0._wp
       !
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 6966)
@@ -19,4 +19,5 @@
    USE p4zbio          !  Biological model
    USE p4zche          !  Chemical model
+   USE p4zlys          !  Calcite saturation
    USE p4zflx          !  Gas exchange
    USE p4zsbc          !  External source of nutrients
@@ -132,4 +133,5 @@
          !
          CALL p4z_bio( kt, jnt )   ! Biology
+         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
          CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
          CALL p4z_sed( kt, jnt )   ! Sedimentation
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zbio.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zbio.F90	(revision 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zbio.F90	(revision 6966)
@@ -28,6 +28,4 @@
    USE p5zpoc          !  Remineralisation of organic particles
    USE p5zagg          !  Aggregation of particles
-   USE p4zlys          !  Dissolution of calcite
-   USE p4zfechem       !  Iron chemistry
    USE p4zligand       !  Remineralization of ligands
    USE prtctl_trc      !  print control for debugging
@@ -82,5 +80,4 @@
       CALL p4z_opt  ( kt, knt )     ! Optic: PAR in the water column
       CALL p5z_sink ( kt, knt )     ! vertical flux of particulate organic matter
-      CALL p4z_lys   (kt, knt )     ! Dissolution of calcite
       CALL p4z_fechem(kt, knt )     ! Iron chemistry/scavenging
       CALL p5z_lim  ( kt, knt )     ! co-limitations by the various nutrients
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 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zpoc.F90	(revision 6966)
@@ -67,5 +67,5 @@
       INTEGER  ::   ji, jj, jk, jn
       REAL(wp) ::   zremip, zremig, zdep, zstep
-      REAL(wp) ::   zopon, zopop, zopoc2, zopon2, zopop2, zofer
+      REAL(wp) ::   zopon, zopop, zopoc, zopoc2, zopon2, zopop2, zofer
       REAL(wp) ::   zsizek, zsizek1, alphat, remint
       REAL(wp) ::   solgoc, zpoc
@@ -211,4 +211,5 @@
                zremig = zremigoc(ji,jj,jk) * zstep * tgfunc(ji,jj,jk)
                zopoc2 = zremig  * trb(ji,jj,jk,jpgoc)
+               orem(ji,jj,jk) = zopoc2
                zopon2 = xremipn / xremipc * zremig * trb(ji,jj,jk,jpgon)
                zopop2 = xremipp / xremipc * zremig * trb(ji,jj,jk,jpgop)
@@ -396,6 +397,6 @@
                ! --------------------------------------------------------
                zremip = zremipoc(ji,jj,jk) * zstep * tgfunc(ji,jj,jk)
-
-               orem(ji,jj,jk)  = zremip * trb(ji,jj,jk,jppoc)
+               zopoc  = zremip * trb(ji,jj,jk,jppoc)
+               orem(ji,jj,jk)  = orem(ji,jj,jk) + zopoc
                zopon  = xremipn / xremipc * zremip * trb(ji,jj,jk,jppon)
                zopop  = xremipp / xremipc * zremip * trb(ji,jj,jk,jppop)
@@ -407,9 +408,9 @@
                ! Update the appropriate tracers trends
                ! -------------------------------------
-               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - orem(ji,jj,jk)
+               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) - zopoc
                tra(ji,jj,jk,jppon) = tra(ji,jj,jk,jppon) - zopon
                tra(ji,jj,jk,jppop) = tra(ji,jj,jk,jppop) - zopop
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) - zofer
-               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + orem(ji,jj,jk)
+               tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + zopoc
                tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + zopon 
                tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + zopop 
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 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zprod.F90	(revision 6966)
@@ -33,11 +33,11 @@
 
    !! * Shared module variables
-   REAL(wp), PUBLIC ::  pislope         !:
+   REAL(wp), PUBLIC ::  pislopen        !:
    REAL(wp), PUBLIC ::  pislopep        !:
-   REAL(wp), PUBLIC ::  pislope2        !:
+   REAL(wp), PUBLIC ::  pisloped        !:
    REAL(wp), PUBLIC ::  xadap           !:
-   REAL(wp), PUBLIC ::  excret          !:
+   REAL(wp), PUBLIC ::  excretn         !:
    REAL(wp), PUBLIC ::  excretp         !:
-   REAL(wp), PUBLIC ::  excret2         !:
+   REAL(wp), PUBLIC ::  excretd         !:
    REAL(wp), PUBLIC ::  bresp           !:
    REAL(wp), PUBLIC ::  thetanpm        !:
@@ -56,7 +56,7 @@
    
    REAL(wp) :: r1_rday                !: 1 / rday
-   REAL(wp) :: texcret                !: 1 - excret 
+   REAL(wp) :: texcretn               !: 1 - excret 
    REAL(wp) :: texcretp               !: 1 - excretp 
-   REAL(wp) :: texcret2               !: 1 - excret2        
+   REAL(wp) :: texcretd               !: 1 - excret2        
 
 
@@ -87,16 +87,15 @@
       REAL(wp) ::   zpronmax, zpropmax, zprofmax, zrat
       REAL(wp) ::   zlim, zsilfac2, zsiborn, zprod, zprontot, zproptot, zprodtot
-      REAL(wp) ::   zmxltst, zmxlday, zprnutmax, zdocprod
-      REAL(wp) ::   zpislopen, zpislopep, zpislope2n
+      REAL(wp) ::   zprnutmax, zdocprod, zprochln, zprochld, zprochlp
+      REAL(wp) ::   zpislopen, zpislopep, zpisloped
       REAL(wp) ::   zrum, zcodel, zargu, zval, zfeup
       REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zmixnano, zmixpico, zmixdiat, zstrn
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopead, zpislopeadp, zpislopead2
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpislopeadn, zpislopeadp, zpislopead
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zprbio, zprpic, zprdia, zysopt
       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(:,:,:) :: zprochln, zprochlp, zprochld
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zpronewn, zpronewp, zpronewd
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zproregn, zproregp, zproregd
@@ -105,4 +104,5 @@
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zrespn, zrespp, zrespd, zprnut
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zcroissn, zcroissp, zcroissd
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxl_fac, zmxl_chl
       !!---------------------------------------------------------------------
       !
@@ -111,7 +111,8 @@
       !  Allocate temporary workspace
       CALL wrk_alloc( jpi, jpj,      zmixnano, zmixpico, zmixdiat, zstrn )
-      CALL wrk_alloc( jpi, jpj, jpk, zpislopead, zpislopeadp, zpislopead2, zysopt ) 
+      CALL wrk_alloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
+      CALL wrk_alloc( jpi, jpj, jpk, zpislopeadn, zpislopeadp, zpislopead, zysopt ) 
       CALL wrk_alloc( jpi, jpj, jpk, zprdia, zprpic, zprbio, zprorcan, zprorcap, zprorcad )
-      CALL wrk_alloc( jpi, jpj, jpk, zprofed, zprofep, zprofen, zprochln, zprochlp, zprochld )
+      CALL wrk_alloc( jpi, jpj, jpk, zprofed, zprofep, zprofen )
       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 )
@@ -121,5 +122,4 @@
       zprorcan(:,:,:) = 0._wp ; zprorcap(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp
       zprofed (:,:,:) = 0._wp ; zprofep (:,:,:) = 0._wp ; zprofen (:,:,:) = 0._wp
-      zprochln(:,:,:) = 0._wp ; zprochlp(:,:,:) = 0._wp ; zprochld(:,:,:) = 0._wp
       zpronewn(:,:,:) = 0._wp ; zpronewp(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp
       zproregn(:,:,:) = 0._wp ; zproregp(:,:,:) = 0._wp ; zproregd(:,:,:) = 0._wp 
@@ -159,14 +159,24 @@
          DO jj = 1 ,jpj
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   zval = MAX( 1., zstrn(ji,jj) )
-                  zval = 1.5 * zval / (12. + zval) * (1. - fr_i(ji,jj))
-                  zprbio(ji,jj,jk) = prmaxn(ji,jj,jk) * zval
-                  zprpic(ji,jj,jk) = prmaxp(ji,jj,jk) * zval
-                  zprdia(ji,jj,jk) = prmaxd(ji,jj,jk) * zval
+                  zmxl_fac(ji,jj,jk) = zval
+                  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)
                ENDIF
             END DO
          END DO
       END DO
+
+      zprbio(:,:,:) = prmaxn(:,:,:) * zmxl_fac(:,:,:)
+      zprdia(:,:,:) = prmaxd(:,:,:) * zmxl_fac(:,:,:)
+      zprpic(:,:,:) = prmaxp(:,:,:) * zmxl_fac(:,:,:)
+
 
       ! Maximum light intensity
@@ -180,33 +190,34 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   ! 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       = xadap * ztn / ( 2.+ ztn )
-                  znanotot    = enano(ji,jj,jk) * zstrn(ji,jj)
-                  zpicotot    = epico(ji,jj,jk) * zstrn(ji,jj)
-                  zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
                   !
-                  zpislopead (ji,jj,jk) = pislope * trb(ji,jj,jk,jpnch)    &
+                  zpislopeadn(ji,jj,jk) = pislopen * trb(ji,jj,jk,jpnch)    &
                   &                       /( trb(ji,jj,jk,jpphy) * 12. + rtrn)
-                  zpislopeadp(ji,jj,jk) = pislopep * ( 1. + zadap * EXP( -0.5 * zpicotot ) )   &
+                  zpislopeadp(ji,jj,jk) = pislopep * ( 1. + zadap * EXP( -0.25 * epico(ji,jj,jk) ) )   &
                   &                       * trb(ji,jj,jk,jppch) /( trb(ji,jj,jk,jppic) * 12. + rtrn)
-                  zpislopead2(ji,jj,jk) = pislope2 * trb(ji,jj,jk,jpdch)    &
+                  zpislopeadd(ji,jj,jk) = pisloped * trb(ji,jj,jk,jpdch)    &
                      &                    /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
-                  zpislopen   = zpislopead (ji,jj,jk) / ( prmaxn(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
-                  zpislopep   = zpislopeadp(ji,jj,jk) / ( prmaxp(ji,jj,jk) * rday * xlimpic(ji,jj,jk) + rtrn )
-                  zpislope2n  = zpislopead2(ji,jj,jk) / ( prmaxd(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
+                  !
+                  zpislopen = zpislopeadn(ji,jj,jk) / ( prmaxn(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
+                  zpislopep = zpislopeadp(ji,jj,jk) / ( prmaxp(ji,jj,jk) * rday * xlimpic(ji,jj,jk) + rtrn )
+                  zpisloped = zpislopeadd(ji,jj,jk) / ( prmaxd(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
 
                   ! Computation of production function for Carbon
                   !  ---------------------------------------------
-                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot )  )
-                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep  * zpicotot )  )
-                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot )  )
+                  zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
+                  zprpic(ji,jj,jk) = zprpic(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) )  )
+                  zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
 
                   ! Computation of production function for Chlorophyll
                   !  -------------------------------------------------
-                  zprchln(ji,jj,jk) = prmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) )  )
-                  zprchlp(ji,jj,jk) = prmaxp(ji,jj,jk) * ( 1.- EXP( -zpislopep  * epico(ji,jj,jk) )  )
-                  zprchld(ji,jj,jk) = prmaxd(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) )  )
+                  zpislopen = zpislopen * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+                  zpisloped = zpisloped * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+                  zpislopep = zpislopep * zmxl_fac(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+                  zprchln(ji,jj,jk) = prmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
+                  zprchlp(ji,jj,jk) = prmaxp(ji,jj,jk) * ( 1.- EXP( -zpislopep * epico(ji,jj,jk) )  )
+                  zprchld(ji,jj,jk) = prmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
                ENDIF
             END DO
@@ -217,5 +228,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-              IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+
+                IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   !    Si/C of diatoms
                   !    ------------------------
@@ -238,27 +250,12 @@
       END DO
 
-      !  Computation of the limitation term due to a mixed layer deeper than the euphotic depth
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zmxltst = MAX( 0.e0, hmld(ji,jj) - heup_01(ji,jj) )
-            zmxlday = zmxltst * zmxltst * r1_rday
-            zmixnano(ji,jj) = 1. - zmxlday / ( zlimmxln + zmxlday )
-            zmixpico(ji,jj) = 1. - zmxlday / ( zlimmxlp + zmxlday )
-            zmixdiat(ji,jj) = 1. - zmxlday / ( zlimmxld + zmxlday )
-         END DO
-      END DO
- 
-      !  Mixed-layer effect on production                                                                               
-      DO jk = 1, jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                  zprbio(ji,jj,jk)  = zprbio(ji,jj,jk) * zmixnano(ji,jj)
-                  zprpic(ji,jj,jk)  = zprpic(ji,jj,jk) * zmixpico(ji,jj)
-                  zprdia(ji,jj,jk)  = zprdia(ji,jj,jk) * zmixdiat(ji,jj)
-                  zprchln(ji,jj,jk) = zprchln(ji,jj,jk) * zmixnano(ji,jj)
-                  zprchlp(ji,jj,jk) = zprchlp(ji,jj,jk) * zmixpico(ji,jj)
-                  zprchld(ji,jj,jk) = zprchld(ji,jj,jk) * zmixdiat(ji,jj)
-               ENDIF
+      !  Sea-ice effect on production                                                                               
+      DO jk = 1, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zprbio(ji,jj,jk)  = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+               zprpic(ji,jj,jk)  = zprpic(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
+               zprdia(ji,jj,jk)  = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) ) 
+               zprnut(ji,jj,jk)  = zprnut(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
             END DO
          END DO
@@ -269,5 +266,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   !  production terms for nanophyto.
                   zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trb(ji,jj,jk,jpphy) * rfact2
@@ -306,5 +303,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   !  production terms for picophyto.
                   zprorcap(ji,jj,jk) = zprpic(ji,jj,jk)  * xlimpic(ji,jj,jk) * trb(ji,jj,jk,jppic) * rfact2
@@ -343,5 +340,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                   !  production terms for diatomees
                   zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trb(ji,jj,jk,jpdia) * rfact2
@@ -381,17 +378,24 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+               IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                      !  production terms for nanophyto. ( chlorophyll )
+                  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(ji,jj,jk) = thetannm * zprod / ( zpislopead(ji,jj,jk) * enano(ji,jj,jk) +rtrn )
-                  zprochln(ji,jj,jk) = MAX(zprochln(ji,jj,jk), chlcmin * 12. * zprorcan (ji,jj,jk) )
+                  zprochln = thetannm * 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(ji,jj,jk) = thetanpm * zprod / ( zpislopeadp(ji,jj,jk) * epico(ji,jj,jk) +rtrn )
-                  zprochlp(ji,jj,jk) = MAX(zprochlp(ji,jj,jk), chlcmin * 12. * zprorcap(ji,jj,jk) )
+                  zprochlp = thetanpm * 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(ji,jj,jk) = thetandm * zprod / ( zpislopead2(ji,jj,jk) * ediat(ji,jj,jk) +rtrn )
-                  zprochld(ji,jj,jk) = MAX(zprochld(ji,jj,jk), chlcmin * 12. * zprorcad(ji,jj,jk) )
+                  zprochld = thetandm * 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
+                  tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln * texcretn
+                  tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld * texcretd
+                  tra(ji,jj,jk,jppch) = tra(ji,jj,jk,jppch) + zprochlp * texcretp
                ENDIF
             END DO
@@ -406,5 +410,5 @@
               zproptot = zpronewp(ji,jj,jk) + zproregp(ji,jj,jk)
               zprodtot = zpronewd(ji,jj,jk) + zproregd(ji,jj,jk)
-              zdocprod = excret2 * zprorcad(ji,jj,jk) + excret * zprorcan(ji,jj,jk)  &
+              zdocprod = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
               &          + excretp * zprorcap(ji,jj,jk)
               tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) - zpropo4n(ji,jj,jk) - zpropo4d(ji,jj,jk)  &
@@ -414,13 +418,12 @@
               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) - zproregn(ji,jj,jk) - zproregd(ji,jj,jk)  &
               &                     - zproregp(ji,jj,jk)
-              tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorcan(ji,jj,jk) * texcret     &
+              tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) + zprorcan(ji,jj,jk) * texcretn    &
                  &                  - zpsino3 * zpronewn(ji,jj,jk) - zpsinh4 * zproregn(ji,jj,jk)   &
                  &                  - zrespn(ji,jj,jk) 
               zcroissn(ji,jj,jk) = tra(ji,jj,jk,jpphy) / rfact2/ (trb(ji,jj,jk,jpphy) + rtrn)
-              tra(ji,jj,jk,jpnph) = tra(ji,jj,jk,jpnph) + zprontot * texcret
-              tra(ji,jj,jk,jppph) = tra(ji,jj,jk,jppph) + zpropo4n(ji,jj,jk) * texcret   &
-              &                     + zprodopn(ji,jj,jk) * texcret
-              tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) + zprochln(ji,jj,jk) * texcret
-              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcret
+              tra(ji,jj,jk,jpnph) = tra(ji,jj,jk,jpnph) + zprontot * texcretn
+              tra(ji,jj,jk,jppph) = tra(ji,jj,jk,jppph) + zpropo4n(ji,jj,jk) * texcretn   &
+              &                     + zprodopn(ji,jj,jk) * texcretn
+              tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) + zprofen(ji,jj,jk) * texcretn
               tra(ji,jj,jk,jppic) = tra(ji,jj,jk,jppic) + zprorcap(ji,jj,jk) * texcretp     &
                  &                  - zpsino3 * zpronewp(ji,jj,jk) - zpsinh4 * zproregp(ji,jj,jk)   &
@@ -430,22 +433,20 @@
               tra(ji,jj,jk,jpppi) = tra(ji,jj,jk,jpppi) + zpropo4p(ji,jj,jk) * texcretp   &
               &                     + zprodopp(ji,jj,jk) * texcretp
-              tra(ji,jj,jk,jppch) = tra(ji,jj,jk,jppch) + zprochlp(ji,jj,jk) * texcretp
               tra(ji,jj,jk,jppfe) = tra(ji,jj,jk,jppfe) + zprofep(ji,jj,jk) * texcretp
-              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcret2   &
+              tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) + zprorcad(ji,jj,jk) * texcretd   &
                  &                  - zpsino3 * zpronewd(ji,jj,jk) - zpsinh4 * zproregd(ji,jj,jk)   &
                  &                  - zrespd(ji,jj,jk) 
               zcroissd(ji,jj,jk) = tra(ji,jj,jk,jpdia) / rfact2 / (trb(ji,jj,jk,jpdia) + rtrn)
-              tra(ji,jj,jk,jpndi) = tra(ji,jj,jk,jpndi) + zprodtot * texcret2
-              tra(ji,jj,jk,jppdi) = tra(ji,jj,jk,jppdi) + zpropo4d(ji,jj,jk) * texcret2   &
-              &                     + zprodopd(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) + zprochld(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) + zprofed(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) + zprorcad(ji,jj,jk) * zysopt(ji,jj,jk) * texcret2
-              tra(ji,jj,jk,jpdoc) = tra(ji,jj,jk,jpdoc) + excret2 * zprorcad(ji,jj,jk) + excret * zprorcan(ji,jj,jk)  &
+              tra(ji,jj,jk,jpndi) = tra(ji,jj,jk,jpndi) + zprodtot * texcretd
+              tra(ji,jj,jk,jppdi) = tra(ji,jj,jk,jppdi) + zpropo4d(ji,jj,jk) * texcretd   &
+              &                     + 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,jpdoc) = tra(ji,jj,jk,jpdoc) + excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)  &
               &                     + excretp * zprorcap(ji,jj,jk)
-              tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + excret2 * zprodtot + excret * zprontot   &
+              tra(ji,jj,jk,jpdon) = tra(ji,jj,jk,jpdon) + excretd * zprodtot + excretn * zprontot   &
               &                     + excretp * zproptot
-              tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + excret2 * zpropo4d(ji,jj,jk) + excret * zpropo4n(ji,jj,jk)   &
-              &    - texcret * zprodopn(ji,jj,jk) - texcret2 * zprodopd(ji,jj,jk) + excretp * zpropo4p(ji,jj,jk)     &
+              tra(ji,jj,jk,jpdop) = tra(ji,jj,jk,jpdop) + excretd * zpropo4d(ji,jj,jk) + excretn * zpropo4n(ji,jj,jk)   &
+              &    - texcretn * zprodopn(ji,jj,jk) - texcretd * zprodopd(ji,jj,jk) + excretp * zpropo4p(ji,jj,jk)     &
               &    - texcretp * zprodopp(ji,jj,jk)
               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2ut * ( zproregn(ji,jj,jk) + zproregd(ji,jj,jk)   &
@@ -453,8 +454,8 @@
                  &                + zpronewd(ji,jj,jk) + zpronewp(ji,jj,jk) )   &
                  &                - o2ut * ( zrespn(ji,jj,jk) + zrespp(ji,jj,jk) + zrespd(ji,jj,jk) )
-              zfeup    = texcret * zprofen(ji,jj,jk) + texcret2 * zprofed(ji,jj,jk)  &
+              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) - texcret2 * zprorcad(ji,jj,jk) * zysopt(ji,jj,jk)
+              tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) - texcretd * zprorcad(ji,jj,jk) * zysopt(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)   &
@@ -522,7 +523,8 @@
       !
       CALL wrk_dealloc( jpi, jpj,      zmixnano, zmixpico, zmixdiat, zstrn )
-      CALL wrk_dealloc( jpi, jpj, jpk, zpislopead, zpislopeadp, zpislopead2, zysopt )                           
+      CALL wrk_dealloc( jpi, jpj, jpk, zmxl_fac, zmxl_chl )
+      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, zprochln, zprochlp, zprochld ) 
+      CALL wrk_dealloc( jpi, jpj, jpk, zprofed, zprofep, zprofen ) 
       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 )
@@ -547,5 +549,5 @@
       !!----------------------------------------------------------------------
       !
-      NAMELIST/nampisprod/ pislope, pislopep, pislope2, xadap, bresp, excret, excretp, excret2,  &
+      NAMELIST/nampisprod/ pislopen, pislopep, pisloped, xadap, bresp, excretn, excretp, excretd,  &
          &                 thetannm, thetanpm, thetandm, chlcmin, grosip, zlimmxln,           &
          &                 zlimmxlp, zlimmxld
@@ -568,12 +570,12 @@
          WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'
          WRITE(numout,*) '    mean Si/C ratio                           grosip       =', grosip
-         WRITE(numout,*) '    P-I slope                                 pislope      =', pislope
-         WRITE(numout,*) '    Acclimation factor to low light           xadap       =', xadap
-         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excret       =', excret
+         WRITE(numout,*) '    P-I slope                                 pislopen     =', pislopen
+         WRITE(numout,*) '    Acclimation factor to low light           xadap        =', xadap
+         WRITE(numout,*) '    excretion ratio of nanophytoplankton      excretn      =', excretn
          WRITE(numout,*) '    excretion ratio of picophytoplankton      excretp      =', excretp
-         WRITE(numout,*) '    excretion ratio of diatoms                excret2      =', excret2
+         WRITE(numout,*) '    excretion ratio of diatoms                excretd      =', excretd
          WRITE(numout,*) '    basal respiration in phytoplankton        bresp        =', bresp
          WRITE(numout,*) '    Maximum Chl/C in phytoplankton            chlcmin      =', chlcmin
-         WRITE(numout,*) '    P-I slope  for diatoms                    pislope2     =', pislope2
+         WRITE(numout,*) '    P-I slope  for diatoms                    pisloped     =', pisloped
          WRITE(numout,*) '    P-I slope  for picophytoplankton          pislopep     =', pislopep
          WRITE(numout,*) '    Minimum Chl/N in nanophytoplankton        thetannm     =', thetannm
@@ -586,7 +588,7 @@
       !
       r1_rday   = 1._wp / rday 
-      texcret   = 1._wp - excret
+      texcretn  = 1._wp - excretn
       texcretp  = 1._wp - excretp
-      texcret2  = 1._wp - excret2
+      texcretd  = 1._wp - excretd
       tpp       = 0._wp
       !
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90	(revision 6931)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/TOP_SRC/PISCES/P5Z/p5zsms.F90	(revision 6966)
@@ -20,4 +20,5 @@
    USE p5zbio          !  Biological model
    USE p4zche          !  Chemical model
+   USE p4zlys          !  Calcite saturation
    USE p4zflx          !  Gas exchange
    USE p4zsbc          !  External source of nutrients
@@ -132,4 +133,5 @@
          !
          CALL p5z_bio (kt, knt)    ! Biology
+         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
          CALL p4z_flx( kt, knt )   ! Compute surface fluxes
          CALL p5z_sed (kt, knt)    ! Sedimentation
