Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P2Z/p2zopt.F90	(revision 5462)
@@ -89,5 +89,8 @@
 
       !                                          ! surface irradiance
-      zpar0m (:,:)   = qsr   (:,:) * 0.43        ! ------------------
+      !                                          ! ------------------
+      IF( ln_dm2dc ) THEN   ;   zpar0m(:,:) = qsr_mean(:,:) * 0.43
+      ELSE                  ;   zpar0m(:,:) = qsr     (:,:) * 0.43
+      ENDIF
       zpar100(:,:)   = zpar0m(:,:) * 0.01
       zparr  (:,:,1) = zpar0m(:,:) * 0.5
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zbio.F90	(revision 5462)
@@ -44,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_bio ( kt, jnt )
+   SUBROUTINE p4z_bio ( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_bio  ***
@@ -54,10 +54,6 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt, jnt
-      INTEGER  ::  ji, jj, jk, jn
-      REAL(wp) ::  ztra
-#if defined key_kriest
-      REAL(wp) ::  zcoef1, zcoef2
-#endif
+      INTEGER, INTENT(in) :: kt, knt
+      INTEGER             :: ji, jj, jk, jn
       CHARACTER (len=25) :: charout
 
@@ -80,54 +76,21 @@
 
           
-      CALL p4z_opt  ( kt, jnt )     ! Optic: PAR in the water column
-      CALL p4z_sink ( kt, jnt )     ! vertical flux of particulate organic matter
-      CALL p4z_fechem(kt, jnt )     ! Iron chemistry/scavenging
-      CALL p4z_lim  ( kt, jnt )     ! co-limitations by the various nutrients
-      CALL p4z_prod ( kt, jnt )     ! phytoplankton growth rate over the global ocean. 
+      CALL p4z_opt  ( kt, knt )     ! Optic: PAR in the water column
+      CALL p4z_sink ( kt, knt )     ! vertical flux of particulate organic matter
+      CALL p4z_fechem(kt, knt )     ! Iron chemistry/scavenging
+      CALL p4z_lim  ( kt, knt )     ! co-limitations by the various nutrients
+      CALL p4z_prod ( kt, knt )     ! phytoplankton growth rate over the global ocean. 
       !                             ! (for each element : C, Si, Fe, Chl )
       CALL p4z_mort ( kt      )     ! phytoplankton mortality
-      !                             ! zooplankton sources/sinks routines 
-      CALL p4z_micro( kt, jnt )           ! microzooplankton
-      CALL p4z_meso ( kt, jnt )           ! mesozooplankton
-      CALL p4z_rem  ( kt, jnt )     ! remineralization terms of organic matter+scavenging of Fe
+     !                             ! zooplankton sources/sinks routines 
+      CALL p4z_micro( kt, knt )           ! microzooplankton
+      CALL p4z_meso ( kt, knt )           ! mesozooplankton
+      CALL p4z_rem  ( kt, knt )     ! remineralization terms of organic matter+scavenging of Fe
       !                             ! test if tracers concentrations fall below 0.
-      xnegtr(:,:,:) = 1.e0
-      DO jn = jp_pcs0, jp_pcs1
-         DO jk = 1, jpk
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( ( trn(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
-                     ztra             = ABS( trn(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
-
-                     xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
-                  ENDIF
-              END DO
-            END DO
-         END DO
-      END DO
-      !                                ! where at least 1 tracer concentration becomes negative
-      !                                ! 
-      DO jn = jp_pcs0, jp_pcs1
-         trn(:,:,:,jn) = trn(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
-      END DO
-
-
-      tra(:,:,:,:) = 0.e0
-
-#if defined key_kriest
-      ! 
-      zcoef1 = 1.e0 / xkr_massp 
-      zcoef2 = 1.e0 / xkr_massp / 1.1
-      DO jk = 1,jpkm1
-         trn(:,:,jk,jpnum) = MAX(  trn(:,:,jk,jpnum), trn(:,:,jk,jppoc) * zcoef1 / xnumm(jk)  )
-         trn(:,:,jk,jpnum) = MIN(  trn(:,:,jk,jpnum), trn(:,:,jk,jppoc) * zcoef2              )
-      END DO
-#endif
-
-      !
+      !                                                             !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
          WRITE(charout, FMT="('bio ')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
       ENDIF
       !
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zfechem.F90	(revision 5462)
@@ -48,5 +48,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_fechem( kt, jnt )
+   SUBROUTINE p4z_fechem( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_fechem  ***
@@ -62,5 +62,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       !
       INTEGER  ::   ji, jj, jk, jic
@@ -101,5 +101,5 @@
       ! -------------------------------------------------
       IF( ln_ligvar ) THEN
-         ztotlig(:,:,:) =  0.09 * trn(:,:,:,jpdoc) * 1E6 + ligand * 1E9
+         ztotlig(:,:,:) =  0.09 * trb(:,:,:,jpdoc) * 1E6 + ligand * 1E9
          ztotlig(:,:,:) =  MIN( ztotlig(:,:,:), 10. )
       ELSE
@@ -127,5 +127,5 @@
                   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   = trn(ji,jj,jk,jpoxy) * ( rhop(ji,jj,jk) / 1.e3 )
+                  zoxy   = trb(ji,jj,jk,jpoxy) * ( rhop(ji,jj,jk) / 1.e3 )
                   ! 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 )  &
@@ -137,5 +137,5 @@
                   zkph1 = zkph2 / 5.
                   ! pass the dfe concentration from PISCES
-                  ztfe = trn(ji,jj,jk,jpfer) * 1e9
+                  ztfe = trb(ji,jj,jk,jpfer) * 1e9
                   ! ----------------------------------------------------------
                   ! ANALYTICAL SOLUTION OF ROOTS OF THE FE3+ EQUATION
@@ -204,5 +204,5 @@
                   zkeq           = fekeq(ji,jj,jk)
                   zfesatur       = zTL1(ji,jj,jk) * 1E-9
-                  ztfe           = trn(ji,jj,jk,jpfer) 
+                  ztfe           = trb(ji,jj,jk,jpfer) 
                   ! Fe' is the root of a 2nd order polynom
                   zFe3 (ji,jj,jk) = ( -( 1. + zfesatur * zkeq - zkeq * ztfe )               &
@@ -210,5 +210,5 @@
                      &               + 4. * ztfe * zkeq) ) / ( 2. * zkeq )
                   zFe3 (ji,jj,jk) = zFe3(ji,jj,jk) * 1E9
-                  zFeL1(ji,jj,jk) = MAX( 0., trn(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
+                  zFeL1(ji,jj,jk) = MAX( 0., trb(ji,jj,jk,jpfer) * 1E9 - zFe3(ji,jj,jk) )
               END DO
             END DO
@@ -240,7 +240,7 @@
                ENDIF
 #if defined key_kriest
-               ztrc   = ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpgsi) ) * 1.e6 
+               ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
 #else
-               ztrc   = ( trn(ji,jj,jk,jppoc) + trn(ji,jj,jk,jpgoc) + trn(ji,jj,jk,jpcal) + trn(ji,jj,jk,jpgsi) ) * 1.e6 
+               ztrc   = ( trb(ji,jj,jk,jppoc) + trb(ji,jj,jk,jpgoc) + trb(ji,jj,jk,jpcal) + trb(ji,jj,jk,jpgsi) ) * 1.e6 
 #endif
                IF( ln_dust )  zdust  = dust(ji,jj) / ( wdust / rday ) * tmask(ji,jj,jk) ! dust in kg/m2/s
@@ -251,7 +251,7 @@
                ! to later allocate scavenged iron to the different organic pools
                ! ---------------------------------------------------------
-               zdenom1 = xlam1 * trn(ji,jj,jk,jppoc) / zlam1b
+               zdenom1 = xlam1 * trb(ji,jj,jk,jppoc) / zlam1b
 #if ! defined key_kriest
-               zdenom2 = xlam1 * trn(ji,jj,jk,jpgoc) / zlam1b
+               zdenom2 = xlam1 * trb(ji,jj,jk,jpgoc) / zlam1b
 #endif
 
@@ -262,6 +262,6 @@
                zlamfac = MIN( 1.  , zlamfac )
                zdep    = MIN( 1., 1000. / fsdept(ji,jj,jk) )
-               zlam1b  = xlam1 * MAX( 0.e0, ( trn(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
-               zcoag   = zfeequi * zlam1b * zstep + 1E-4 * ( 1. - zlamfac ) * zdep * zstep * trn(ji,jj,jk,jpfer)
+               zlam1b  = xlam1 * MAX( 0.e0, ( trb(ji,jj,jk,jpfer) * 1.e9 - ztotlig(ji,jj,jk) ) )
+               zcoag   = zfeequi * zlam1b * zstep + 1E-4 * ( 1. - zlamfac ) * zdep * zstep * trb(ji,jj,jk,jpfer)
 
                !  Compute the coagulation of colloidal iron. This parameterization 
@@ -269,6 +269,6 @@
                !  It requires certainly some more work as it is very poorly constrained.
                !  ----------------------------------------------------------------
-               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) )
+               zlam1a  = ( 0.369  * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4  * trb(ji,jj,jk,jppoc) ) * xdiss(ji,jj,jk)    &
+                   &   + ( 114.   * 0.3 * trb(ji,jj,jk,jpdoc) + 5.09E3 * trb(ji,jj,jk,jppoc) )
                zaggdfea = zlam1a * zstep * zfecoll
 #if defined key_kriest
@@ -278,5 +278,5 @@
                tra(ji,jj,jk,jpsfe) = tra(ji,jj,jk,jpsfe) + zscave * zdenom1 + zaggdfea + zaggdfeb
 #else
-               zlam1b = 3.53E3 *   trn(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
+               zlam1b = 3.53E3 *   trb(ji,jj,jk,jpgoc) * xdiss(ji,jj,jk)
                zaggdfeb = zlam1b * zstep * zfecoll
                !
@@ -292,12 +292,12 @@
       !  ----------------------------------------
       IF( ln_fechem ) THEN
-          biron(:,:,:) = MAX( 0., trn(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
+          biron(:,:,:) = MAX( 0., trb(:,:,:,jpfer) - zFeP(:,:,:) * 1E-9 )
       ELSE
-          biron(:,:,:) = trn(:,:,:,jpfer) 
+          biron(:,:,:) = trb(:,:,:,jpfer) 
       ENDIF
 
       !  Output of some diagnostics variables
       !     ---------------------------------
-      IF( lk_iomput .AND. jnt == nrdttrc ) THEN
+      IF( lk_iomput .AND. knt == nrdttrc ) THEN
          IF( iom_use("Fe3")    )  CALL iom_put("Fe3"    , zFe3   (:,:,:)       * tmask(:,:,:) )   ! Fe3+
          IF( iom_use("FeL1")   )  CALL iom_put("FeL1"   , zFeL1  (:,:,:)       * tmask(:,:,:) )   ! FeL1
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zflx.F90	(revision 5462)
@@ -68,5 +68,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_flx ( kt )
+   SUBROUTINE p4z_flx ( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_flx  ***
@@ -81,5 +81,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt   !
+      INTEGER, INTENT(in) ::   kt, knt   !
       !
       INTEGER  ::   ji, jj, jm, iind, iindm1
@@ -101,5 +101,5 @@
       !     IS USED TO COMPUTE AIR-SEA FLUX OF CO2
 
-      IF( kt /= nit000 ) CALL p4z_patm( kt )    ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
+      IF( kt /= nit000 .AND. knt == 1 ) CALL p4z_patm( kt )    ! Get sea-level pressure (E&K [1981] climatology) for use in flux calcs
 
       IF( ln_co2int ) THEN 
@@ -130,7 +130,7 @@
                zbot  = borat(ji,jj,1)
                zfact = rhop(ji,jj,1) / 1000. + rtrn
-               zdic  = trn(ji,jj,1,jpdic) / zfact
+               zdic  = trb(ji,jj,1,jpdic) / zfact
                zph   = MAX( hi(ji,jj,1), 1.e-10 ) / zfact
-               zalka = trn(ji,jj,1,jptal) / zfact
+               zalka = trb(ji,jj,1,jptal) / zfact
 
                ! CALCULATE [ALK]([CO3--], [HCO3-])
@@ -184,13 +184,13 @@
             zfld = satmco2(ji,jj) * patm(ji,jj) * tmask(ji,jj,1) * chemc(ji,jj,1) * zkgco2(ji,jj)   ! (mol/L) * (m/s)
             zflu = zh2co3(ji,jj) * tmask(ji,jj,1) * zkgco2(ji,jj)                                   ! (mol/L) (m/s) ?
-            oce_co2(ji,jj) = ( zfld - zflu ) * rfact * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
+            oce_co2(ji,jj) = ( zfld - zflu ) * rfact2 * e1e2t(ji,jj) * tmask(ji,jj,1) * 1000.
             ! compute the trend
-            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) / fse3t(ji,jj,1)
+            tra(ji,jj,1,jpdic) = tra(ji,jj,1,jpdic) + ( zfld - zflu ) * rfact2 / fse3t(ji,jj,1)
 
             ! Compute O2 flux 
             zfld16 = atcox * patm(ji,jj) * chemc(ji,jj,2) * tmask(ji,jj,1) * zkgo2(ji,jj)          ! (mol/L) * (m/s)
-            zflu16 = trn(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
+            zflu16 = trb(ji,jj,1,jpoxy) * tmask(ji,jj,1) * zkgo2(ji,jj)
             zoflx(ji,jj) = zfld16 - zflu16
-            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) / fse3t(ji,jj,1)
+            tra(ji,jj,1,jpoxy) = tra(ji,jj,1,jpoxy) + zoflx(ji,jj) * rfact2 / fse3t(ji,jj,1)
          END DO
       END DO
@@ -207,8 +207,8 @@
       ENDIF
 
-      IF( lk_iomput ) THEN
+      IF( lk_iomput .AND. knt == nrdttrc ) THEN
          CALL wrk_alloc( jpi, jpj, zw2d )  
          IF( iom_use( "Cflx"  ) )  THEN
-            zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) / rfact
+            zw2d(:,:) = oce_co2(:,:) / e1e2t(:,:) * rfact2r
             CALL iom_put( "Cflx"     , zw2d ) 
          ENDIF
@@ -226,5 +226,5 @@
          ENDIF
          IF( iom_use( "Dpo2" ) )  THEN
-           zw2d(:,:) = ( atcox * patm(:,:) - trn(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1)
+           zw2d(:,:) = ( atcox * patm(:,:) - trb(:,:,1,jpoxy) / ( chemc(:,:,2) + rtrn ) ) * tmask(:,:,1)
            CALL iom_put( "Dpo2"  , zw2d )
          ENDIF
@@ -235,5 +235,5 @@
       ELSE
          IF( ln_diatrc ) THEN
-            trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) / rfact 
+            trc2d(:,:,jp_pcs0_2d    ) = oce_co2(:,:) / e1e2t(:,:) * rfact2r 
             trc2d(:,:,jp_pcs0_2d + 1) = zoflx(:,:) * 1000 * tmask(:,:,1) 
             trc2d(:,:,jp_pcs0_2d + 2) = zkgco2(:,:) * tmask(:,:,1) 
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zint.F90	(revision 5462)
@@ -56,5 +56,5 @@
       DO ji = 1, jpi
          DO jj = 1, jpj
-            zvar = trn(ji,jj,1,jpsil) * trn(ji,jj,1,jpsil)
+            zvar = trb(ji,jj,1,jpsil) * trb(ji,jj,1,jpsil)
             xksimax(ji,jj) = MAX( xksimax(ji,jj), ( 1.+ 7.* zvar / ( xksilim * xksilim + zvar ) ) * 1e-6 )
          END DO
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlim.F90	(revision 5462)
@@ -62,5 +62,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_lim( kt, jnt )
+   SUBROUTINE p4z_lim( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_lim  ***
@@ -72,10 +72,10 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in)  :: kt, jnt
+      INTEGER, INTENT(in)  :: kt, knt
       !
       INTEGER  ::   ji, jj, jk
       REAL(wp) ::   zlim1, zlim2, zlim3, zlim4, zno3, zferlim
       REAL(wp) ::   zconcd, zconcd2, zconcn, zconcn2
-      REAL(wp) ::   z1_trndia, z1_trnphy, ztem1, ztem2, zetot1, zetot2
+      REAL(wp) ::   z1_trbdia, z1_trbphy, ztem1, ztem2, zetot1, zetot2
       REAL(wp) ::   zdenom, zratio, zironmin
       REAL(wp) ::   zconc1d, zconc1dnh4, zconc0n, zconc0nnh4   
@@ -90,37 +90,37 @@
                ! Tuning of the iron concentration to a minimum level that is set to the detection limit
                !-------------------------------------
-               zno3    = trn(ji,jj,jk,jpno3) / 40.e-6
+               zno3    = trb(ji,jj,jk,jpno3) / 40.e-6
                zferlim = MAX( 3e-11 * zno3 * zno3, 5e-12 )
                zferlim = MIN( zferlim, 7e-11 )
-               trn(ji,jj,jk,jpfer) = MAX( trn(ji,jj,jk,jpfer), zferlim )
+               trb(ji,jj,jk,jpfer) = MAX( trb(ji,jj,jk,jpfer), zferlim )
 
                ! Computation of a variable Ks for iron on diatoms taking into account
                ! that increasing biomass is made of generally bigger cells
                !------------------------------------------------
-               zconcd   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - xsizedia )
-               zconcd2  = trn(ji,jj,jk,jpdia) - zconcd
-               zconcn   = MAX( 0.e0 , trn(ji,jj,jk,jpphy) - xsizephy )
-               zconcn2  = trn(ji,jj,jk,jpphy) - zconcn
-               z1_trnphy   = 1. / ( trn(ji,jj,jk,jpphy) + rtrn )
-               z1_trndia   = 1. / ( trn(ji,jj,jk,jpdia) + rtrn )
-
-               concdfe(ji,jj,jk) = MAX( concdfer, ( zconcd2 * concdfer + concdfer * xsizerd * zconcd ) * z1_trndia )
-               zconc1d           = MAX( concdno3, ( zconcd2 * concdno3 + concdno3 * xsizerd * zconcd ) * z1_trndia )
-               zconc1dnh4        = MAX( concdnh4, ( zconcd2 * concdnh4 + concdnh4 * xsizerd * zconcd ) * z1_trndia )
-
-               concnfe(ji,jj,jk) = MAX( concnfer, ( zconcn2 * concnfer + concnfer * xsizern * zconcn ) * z1_trnphy )
-               zconc0n           = MAX( concnno3, ( zconcn2 * concnno3 + concnno3 * xsizern * zconcn ) * z1_trnphy )
-               zconc0nnh4        = MAX( concnnh4, ( zconcn2 * concnnh4 + concnnh4 * xsizern * zconcn ) * z1_trnphy )
+               zconcd   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
+               zconcd2  = trb(ji,jj,jk,jpdia) - zconcd
+               zconcn   = MAX( 0.e0 , trb(ji,jj,jk,jpphy) - xsizephy )
+               zconcn2  = trb(ji,jj,jk,jpphy) - zconcn
+               z1_trbphy   = 1. / ( trb(ji,jj,jk,jpphy) + rtrn )
+               z1_trbdia   = 1. / ( trb(ji,jj,jk,jpdia) + rtrn )
+
+               concdfe(ji,jj,jk) = MAX( concdfer, ( zconcd2 * concdfer + concdfer * xsizerd * zconcd ) * z1_trbdia )
+               zconc1d           = MAX( concdno3, ( zconcd2 * concdno3 + concdno3 * xsizerd * zconcd ) * z1_trbdia )
+               zconc1dnh4        = MAX( concdnh4, ( zconcd2 * concdnh4 + concdnh4 * xsizerd * zconcd ) * z1_trbdia )
+
+               concnfe(ji,jj,jk) = MAX( concnfer, ( zconcn2 * concnfer + concnfer * xsizern * zconcn ) * z1_trbphy )
+               zconc0n           = MAX( concnno3, ( zconcn2 * concnno3 + concnno3 * xsizern * zconcn ) * z1_trbphy )
+               zconc0nnh4        = MAX( concnnh4, ( zconcn2 * concnnh4 + concnnh4 * xsizern * zconcn ) * z1_trbphy )
 
                ! 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
+               zdenom = 1. /  ( concbno3 * concbnh4 + concbnh4 * trb(ji,jj,jk,jpno3) + concbno3 * trb(ji,jj,jk,jpnh4) )
+               xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * concbnh4 * zdenom
+               xnanonh4(ji,jj,jk) = trb(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) )
+               zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concbnh4 )
+               zlim3    = trb(ji,jj,jk,jpfer) / ( concbfe + trb(ji,jj,jk,jpfer) )
+               zlim4    = trb(ji,jj,jk,jpdoc) / ( xkdoc   + trb(ji,jj,jk,jpdoc) )
                xlimbacl(ji,jj,jk) = MIN( zlim1, zlim2, zlim3 )
                xlimbac (ji,jj,jk) = MIN( zlim1, zlim2, zlim3 ) * zlim4
@@ -128,12 +128,12 @@
                ! Michaelis-Menten Limitation term for nutrients Small flagellates
                ! -----------------------------------------------
-               zdenom = 1. /  ( zconc0n * zconc0nnh4 + zconc0nnh4 * trn(ji,jj,jk,jpno3) + zconc0n * trn(ji,jj,jk,jpnh4) )
-               xnanono3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * zconc0nnh4 * zdenom
-               xnanonh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * zconc0n    * zdenom
+               zdenom = 1. /  ( zconc0n * zconc0nnh4 + zconc0nnh4 * trb(ji,jj,jk,jpno3) + zconc0n * trb(ji,jj,jk,jpnh4) )
+               xnanono3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc0nnh4 * zdenom
+               xnanonh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * zconc0n    * zdenom
                !
                zlim1    = xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk)
-               zlim2    = trn(ji,jj,jk,jppo4) / ( trn(ji,jj,jk,jppo4) + zconc0nnh4 )
-               zratio   = trn(ji,jj,jk,jpnfe) * z1_trnphy 
-               zironmin = xcoef1 * trn(ji,jj,jk,jpnch) * z1_trnphy + xcoef2 * zlim1 + xcoef3 * xnanono3(ji,jj,jk)
+               zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc0nnh4 )
+               zratio   = trb(ji,jj,jk,jpnfe) * z1_trbphy 
+               zironmin = xcoef1 * trb(ji,jj,jk,jpnch) * z1_trbphy + xcoef2 * zlim1 + xcoef3 * xnanono3(ji,jj,jk)
                zlim3    = MAX( 0.,( zratio - zironmin ) / qnfelim )
                xnanopo4(ji,jj,jk) = zlim2
@@ -143,13 +143,13 @@
                !   Michaelis-Menten Limitation term for nutrients Diatoms
                !   ----------------------------------------------
-               zdenom   = 1. / ( zconc1d * zconc1dnh4 + zconc1dnh4 * trn(ji,jj,jk,jpno3) + zconc1d * trn(ji,jj,jk,jpnh4) )
-               xdiatno3(ji,jj,jk) = trn(ji,jj,jk,jpno3) * zconc1dnh4 * zdenom
-               xdiatnh4(ji,jj,jk) = trn(ji,jj,jk,jpnh4) * zconc1d    * zdenom
+               zdenom   = 1. / ( zconc1d * zconc1dnh4 + zconc1dnh4 * trb(ji,jj,jk,jpno3) + zconc1d * trb(ji,jj,jk,jpnh4) )
+               xdiatno3(ji,jj,jk) = trb(ji,jj,jk,jpno3) * zconc1dnh4 * zdenom
+               xdiatnh4(ji,jj,jk) = trb(ji,jj,jk,jpnh4) * zconc1d    * zdenom
                !
                zlim1    = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
-               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) * z1_trndia
-               zironmin = xcoef1 * trn(ji,jj,jk,jpdch) * z1_trndia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
+               zlim2    = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + zconc1dnh4  )
+               zlim3    = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi(ji,jj) )
+               zratio   = trb(ji,jj,jk,jpdfe) * z1_trbdia
+               zironmin = xcoef1 * trb(ji,jj,jk,jpdch) * z1_trbdia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
                zlim4    = MAX( 0., ( zratio - zironmin ) / qdfelim )
                xdiatpo4(ji,jj,jk) = zlim2
@@ -166,16 +166,16 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zlim1 =  ( trn(ji,jj,jk,jpno3) * concnnh4 + trn(ji,jj,jk,jpnh4) * concnno3 )    &
-                  &   / ( 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) +  5.E-11   )
+               zlim1 =  ( trb(ji,jj,jk,jpno3) * concnnh4 + trb(ji,jj,jk,jpnh4) * concnno3 )    &
+                  &   / ( concnno3 * concnnh4 + concnnh4 * trb(ji,jj,jk,jpno3) + concnno3 * trb(ji,jj,jk,jpnh4) ) 
+               zlim2  = trb(ji,jj,jk,jppo4) / ( trb(ji,jj,jk,jppo4) + concnnh4 )
+               zlim3  = trb(ji,jj,jk,jpfer) / ( trb(ji,jj,jk,jpfer) +  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 = 30. / ( 30. + etot(ji,jj,jk) ) 
+               zetot1 = MAX( 0., etot_ndcy(ji,jj,jk) - 1.) / ( 4. + etot_ndcy(ji,jj,jk) ) 
+               zetot2 = 30. / ( 30. + etot_ndcy(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. )  &
+                  &                       * MAX( 1., trb(ji,jj,jk,jpphy) * 1.e6 / 2. )  &
                   &                       * zetot1 * zetot2               &
                   &                       * ( 1. + EXP(-ztem2 * ztem2 / 25. ) )         &
@@ -188,5 +188,5 @@
       !
       !
-      IF( lk_iomput .AND. jnt == nrdttrc ) THEN        ! save output diagnostics
+      IF( lk_iomput .AND. knt == nrdttrc ) THEN        ! save output diagnostics
         IF( iom_use( "xfracal" ) ) CALL iom_put( "xfracal", xfracal(:,:,:) * tmask(:,:,:) )  ! euphotic layer deptht
         IF( iom_use( "LNnut"   ) ) CALL iom_put( "LNnut"  , xlimphy(:,:,:) * tmask(:,:,:) )  ! Nutrient limitation term
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zlys.F90	(revision 5462)
@@ -48,5 +48,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_lys( kt )
+   SUBROUTINE p4z_lys( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_lys  ***
@@ -59,10 +59,9 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       INTEGER  ::   ji, jj, jk, jn
       REAL(wp) ::   zalk, zdic, zph, zah2
       REAL(wp) ::   zdispot, zfact, zcalcon, zalka, zaldi
       REAL(wp) ::   zomegaca, zexcess, zexcess0
-      REAL(wp) ::   zrfact2
       CHARACTER (len=25) :: charout
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zco3, zcaldiss   
@@ -89,6 +88,6 @@
                   zfact = rhop(ji,jj,jk) / 1000. + rtrn
                   zph  = hi(ji,jj,jk) * tmask(ji,jj,jk) / zfact + ( 1.-tmask(ji,jj,jk) ) * 1.e-9 ! [H+]
-                  zdic  = trn(ji,jj,jk,jpdic) / zfact
-                  zalka = trn(ji,jj,jk,jptal) / zfact
+                  zdic  = trb(ji,jj,jk,jpdic) / zfact
+                  zalka = trb(ji,jj,jk,jptal) / zfact
                   ! CALCULATE [ALK]([CO3--], [HCO3-])
                   zalk  = zalka - ( akw3(ji,jj,jk) / zph - zph + borat(ji,jj,jk) / ( 1. + zph / akb3(ji,jj,jk) ) )
@@ -130,5 +129,5 @@
                !       (ACCORDING TO THIS FORMULATION ALSO SOME PARTICULATE
                !       CACO3 GETS DISSOLVED EVEN IN THE CASE OF OVERSATURATION)
-               zdispot = kdca * zexcess * trn(ji,jj,jk,jpcal)
+               zdispot = kdca * zexcess * trb(ji,jj,jk,jpcal)
 # if defined key_degrad
                zdispot = zdispot * facvol(ji,jj,jk)
@@ -136,6 +135,6 @@
               !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
               !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
-              zcaldiss(ji,jj,jk)  = zdispot / rmtss  ! calcite dissolution
-              zco3(ji,jj,jk)      = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk) * rfact
+              zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
+              zco3(ji,jj,jk)      = zco3(ji,jj,jk) + zcaldiss(ji,jj,jk)
               !
               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + 2. * zcaldiss(ji,jj,jk)
@@ -147,9 +146,9 @@
       !
 
-      IF( lk_iomput ) THEN
+      IF( lk_iomput .AND. knt == nrdttrc ) THEN
          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( "DCAL"   ) ) CALL iom_put( "DCAL"  , zcaldiss(:,:,:) * 1.e+3 * rfact2r   * tmask(:,:,:) )
       ELSE
          trc3d(:,:,:,jp_pcs0_3d    ) = -1. * LOG10( hi(:,:,:) ) * tmask(:,:,:)
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmeso.F90	(revision 5462)
@@ -60,5 +60,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_meso( kt, jnt )
+   SUBROUTINE p4z_meso( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_meso  ***
@@ -68,5 +68,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       INTEGER  :: ji, jj, jk
       REAL(wp) :: zcompadi, zcompaph, zcompapoc, zcompaz, zcompam
@@ -97,5 +97,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zcompam   = MAX( ( trn(ji,jj,jk,jpmes) - 1.e-9 ), 0.e0 )
+               zcompam   = MAX( ( trb(ji,jj,jk,jpmes) - 1.e-9 ), 0.e0 )
 # if defined key_degrad
                zstep     = xstep * facvol(ji,jj,jk)
@@ -107,5 +107,5 @@
                !  Respiration rates of both zooplankton
                !  -------------------------------------
-               zrespz2   = resrat2 * zfact * trn(ji,jj,jk,jpmes) / ( xkmort + trn(ji,jj,jk,jpmes) )  &
+               zrespz2   = resrat2 * zfact * trb(ji,jj,jk,jpmes) / ( xkmort + trb(ji,jj,jk,jpmes) )  &
                   &      + resrat2 * zfact * 3. * nitrfac(ji,jj,jk)
 
@@ -113,14 +113,14 @@
                !  no real reason except that it seems to be more stable and may mimic predation
                !  ---------------------------------------------------------------
-               ztortz2   = mzrat2 * 1.e6 * zfact * trn(ji,jj,jk,jpmes)
+               ztortz2   = mzrat2 * 1.e6 * zfact * trb(ji,jj,jk,jpmes)
                !
-               zcompadi  = MAX( ( trn(ji,jj,jk,jpdia) - xthresh2dia ), 0.e0 )
-               zcompaz   = MAX( ( trn(ji,jj,jk,jpzoo) - xthresh2zoo ), 0.e0 )
+               zcompadi  = MAX( ( trb(ji,jj,jk,jpdia) - xthresh2dia ), 0.e0 )
+               zcompaz   = MAX( ( trb(ji,jj,jk,jpzoo) - xthresh2zoo ), 0.e0 )
                ! Size effect of nanophytoplankton on grazing : the smaller it is, the less prone
                ! it is to predation by mesozooplankton
                ! -------------------------------------------------------------------------------
-               zcompaph  = MAX( ( trn(ji,jj,jk,jpphy) - xthresh2phy ), 0.e0 ) &
+               zcompaph  = MAX( ( trb(ji,jj,jk,jpphy) - xthresh2phy ), 0.e0 ) &
                   &      * MIN(1., MAX( 0., ( quotan(ji,jj,jk) - 0.2) / 0.3 ) )
-               zcompapoc = MAX( ( trn(ji,jj,jk,jppoc) - xthresh2poc ), 0.e0 )
+               zcompapoc = MAX( ( trb(ji,jj,jk,jppoc) - xthresh2poc ), 0.e0 )
 
                zfood     = xprefc * zcompadi + xprefz * zcompaz + xprefp * zcompaph + xprefpoc * zcompapoc 
@@ -128,5 +128,5 @@
                zdenom    = zfoodlim / ( xkgraz2 + zfoodlim )
                zdenom2   = zdenom / ( zfood + rtrn )
-               zgraze2   = grazrat2 * zstep * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpmes) 
+               zgraze2   = grazrat2 * zstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpmes) 
 
                zgrazd    = zgraze2  * xprefc   * zcompadi  * zdenom2 
@@ -135,7 +135,7 @@
                zgrazpoc  = zgraze2  * xprefpoc * zcompapoc * zdenom2 
 
-               zgraznf   = zgrazn   * trn(ji,jj,jk,jpnfe) / ( trn(ji,jj,jk,jpphy) + rtrn)
-               zgrazf    = zgrazd   * trn(ji,jj,jk,jpdfe) / ( trn(ji,jj,jk,jpdia) + rtrn)
-               zgrazpof  = zgrazpoc * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn)
+               zgraznf   = zgrazn   * trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) + rtrn)
+               zgrazf    = zgrazd   * trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn)
+               zgrazpof  = zgrazpoc * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn)
 
                !  Mesozooplankton flux feeding on GOC
@@ -144,10 +144,10 @@
 # if ! defined key_kriest
                zgrazffeg = grazflux  * zstep * wsbio4(ji,jj,jk)      &
-               &           * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpmes)
-               zgrazfffg = zgrazffeg * trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)
+               &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpgoc) * trb(ji,jj,jk,jpmes)
+               zgrazfffg = zgrazffeg * trb(ji,jj,jk,jpbfe) / (trb(ji,jj,jk,jpgoc) + rtrn)
 # endif
                zgrazffep = grazflux  * zstep *  wsbio3(ji,jj,jk)     &
-               &           * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpmes)
-               zgrazfffp = zgrazffep * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
+               &           * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpmes)
+               zgrazfffp = zgrazffep * trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
               !
 # if ! defined key_kriest
@@ -158,10 +158,10 @@
               ! diatoms based aggregates are more prone to fractionation
               ! since they are more porous (marine snow instead of fecal pellets)
-              zratio    = trn(ji,jj,jk,jpgsi) / ( trn(ji,jj,jk,jpgoc) + rtrn )
+              zratio    = trb(ji,jj,jk,jpgsi) / ( trb(ji,jj,jk,jpgoc) + rtrn )
               zratio2   = zratio * zratio
               zfrac     = zproport * grazflux  * zstep * wsbio4(ji,jj,jk)      &
-               &          * trn(ji,jj,jk,jpgoc) * trn(ji,jj,jk,jpmes)          &
+               &          * trb(ji,jj,jk,jpgoc) * trb(ji,jj,jk,jpmes)          &
                &          * ( 0.2 + 3.8 * zratio2 / ( 1.**2 + zratio2 ) )
-              zfracfe   = zfrac * trn(ji,jj,jk,jpbfe) / (trn(ji,jj,jk,jpgoc) + rtrn)
+              zfracfe   = zfrac * trb(ji,jj,jk,jpbfe) / (trb(ji,jj,jk,jpgoc) + rtrn)
 
               zgrazffep = zproport * zgrazffep
@@ -215,8 +215,8 @@
                tra(ji,jj,jk,jpzoo) = tra(ji,jj,jk,jpzoo) - zgrazz
                tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazn
-               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazn * trn(ji,jj,jk,jpnch) / ( trn(ji,jj,jk,jpphy) + rtrn )
-               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazd * trn(ji,jj,jk,jpdch) / ( trn(ji,jj,jk,jpdia) + rtrn )
-               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zgrazd * trn(ji,jj,jk,jpdsi) / ( trn(ji,jj,jk,jpdia) + rtrn )
-               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zgrazd * trn(ji,jj,jk,jpdsi) / ( trn(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazn * trb(ji,jj,jk,jpnch) / ( trb(ji,jj,jk,jpphy) + rtrn )
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazd * trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zgrazd * trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn )
+               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zgrazd * trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn )
                tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgraznf
                tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazf
@@ -231,5 +231,5 @@
                tra(ji,jj,jk,jpcal) = tra(ji,jj,jk,jpcal) + zprcaca
 #if defined key_kriest
-              znumpoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
+              znumpoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
               tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zmortzgoc - zgrazpoc - zgrazffep + zgrapoc2
               tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) - zgrazpoc * znumpoc + zgrapoc2 * xkr_dmeso      &
@@ -248,5 +248,5 @@
       END DO
       !
-      IF( lk_iomput .AND. jnt == nrdttrc ) THEN
+      IF( lk_iomput .AND. knt == nrdttrc ) THEN
          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
          IF( iom_use( "GRAZ2" ) ) THEN
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmicro.F90	(revision 5462)
@@ -59,5 +59,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_micro( kt, jnt )
+   SUBROUTINE p4z_micro( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_micro  ***
@@ -68,5 +68,5 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::  kt  ! ocean time step
-      INTEGER, INTENT(in) ::  jnt 
+      INTEGER, INTENT(in) ::  knt 
       !
       INTEGER  :: ji, jj, jk
@@ -90,5 +90,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zcompaz = MAX( ( trn(ji,jj,jk,jpzoo) - 1.e-9 ), 0.e0 )
+               zcompaz = MAX( ( trb(ji,jj,jk,jpzoo) - 1.e-9 ), 0.e0 )
                zstep   = xstep
 # if defined key_degrad
@@ -99,5 +99,5 @@
                !  Respiration rates of both zooplankton
                !  -------------------------------------
-               zrespz = resrat * zfact * trn(ji,jj,jk,jpzoo) / ( xkmort + trn(ji,jj,jk,jpzoo) )  &
+               zrespz = resrat * zfact * trb(ji,jj,jk,jpzoo) / ( xkmort + trb(ji,jj,jk,jpzoo) )  &
                   &   + resrat * zfact * 3. * nitrfac(ji,jj,jk)
 
@@ -105,9 +105,9 @@
                !  no real reason except that it seems to be more stable and may mimic predation.
                !  ---------------------------------------------------------------
-               ztortz = mzrat * 1.e6 * zfact * trn(ji,jj,jk,jpzoo)
-
-               zcompadi  = MIN( MAX( ( trn(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 ), xsizedia )
-               zcompaph  = MAX( ( trn(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
-               zcompapoc = MAX( ( trn(ji,jj,jk,jppoc) - xthreshpoc ), 0.e0 )
+               ztortz = mzrat * 1.e6 * zfact * trb(ji,jj,jk,jpzoo)
+
+               zcompadi  = MIN( MAX( ( trb(ji,jj,jk,jpdia) - xthreshdia ), 0.e0 ), xsizedia )
+               zcompaph  = MAX( ( trb(ji,jj,jk,jpphy) - xthreshphy ), 0.e0 )
+               zcompapoc = MAX( ( trb(ji,jj,jk,jppoc) - xthreshpoc ), 0.e0 )
                
                !     Microzooplankton grazing
@@ -117,5 +117,5 @@
                zdenom    = zfoodlim / ( xkgraz + zfoodlim )
                zdenom2   = zdenom / ( zfood + rtrn )
-               zgraze    = grazrat * zstep * tgfunc2(ji,jj,jk) * trn(ji,jj,jk,jpzoo) 
+               zgraze    = grazrat * zstep * tgfunc2(ji,jj,jk) * trb(ji,jj,jk,jpzoo) 
 
                zgrazp    = zgraze  * xpref2p * zcompaph  * zdenom2 
@@ -123,7 +123,7 @@
                zgrazsd   = zgraze  * xpref2d * zcompadi  * zdenom2 
 
-               zgrazpf   = zgrazp  * trn(ji,jj,jk,jpnfe) / (trn(ji,jj,jk,jpphy) + rtrn)
-               zgrazmf   = zgrazm  * trn(ji,jj,jk,jpsfe) / (trn(ji,jj,jk,jppoc) + rtrn)
-               zgrazsf   = zgrazsd * trn(ji,jj,jk,jpdfe) / (trn(ji,jj,jk,jpdia) + rtrn)
+               zgrazpf   = zgrazp  * trb(ji,jj,jk,jpnfe) / (trb(ji,jj,jk,jpphy) + rtrn)
+               zgrazmf   = zgrazm  * trb(ji,jj,jk,jpsfe) / (trb(ji,jj,jk,jppoc) + rtrn)
+               zgrazsf   = zgrazsd * trb(ji,jj,jk,jpdfe) / (trb(ji,jj,jk,jpdia) + rtrn)
                !
                zgraztot  = zgrazp  + zgrazm  + zgrazsd 
@@ -165,8 +165,8 @@
                tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zgrazp
                tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zgrazsd
-               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazp  * trn(ji,jj,jk,jpnch)/(trn(ji,jj,jk,jpphy)+rtrn)
-               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazsd * trn(ji,jj,jk,jpdch)/(trn(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zgrazsd * trn(ji,jj,jk,jpdsi)/(trn(ji,jj,jk,jpdia)+rtrn)
-               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zgrazsd * trn(ji,jj,jk,jpdsi)/(trn(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zgrazp  * trb(ji,jj,jk,jpnch)/(trb(ji,jj,jk,jpphy)+rtrn)
+               tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zgrazsd * trb(ji,jj,jk,jpdch)/(trb(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpdsi) = tra(ji,jj,jk,jpdsi) - zgrazsd * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
+               tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) + zgrazsd * trb(ji,jj,jk,jpdsi)/(trb(ji,jj,jk,jpdia)+rtrn)
                tra(ji,jj,jk,jpnfe) = tra(ji,jj,jk,jpnfe) - zgrazpf
                tra(ji,jj,jk,jpdfe) = tra(ji,jj,jk,jpdfe) - zgrazsf
@@ -184,5 +184,5 @@
 #if defined key_kriest
                tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zmortz * xkr_dmicro &
-                                                         - zgrazm * trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
+                                                         - zgrazm * trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
 #endif
             END DO
@@ -190,5 +190,5 @@
       END DO
       !
-      IF( lk_iomput .AND. jnt == nrdttrc ) THEN
+      IF( lk_iomput .AND. knt == nrdttrc ) THEN
          CALL wrk_alloc( jpi, jpj, jpk, zw3d )
          IF( iom_use( "GRAZ1" ) ) THEN
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zmort.F90	(revision 5462)
@@ -85,5 +85,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zcompaph = MAX( ( trn(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
+               zcompaph = MAX( ( trb(ji,jj,jk,jpphy) - 1e-8 ), 0.e0 )
                zstep    = xstep
 # if defined key_degrad
@@ -94,5 +94,5 @@
                !     due to turbulence is negligible. Mortality is also set
                !     to 0
-               zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * trn(ji,jj,jk,jpphy)
+               zsizerat = MIN(1., MAX( 0., (quotan(ji,jj,jk) - 0.2) / 0.3) ) * trb(ji,jj,jk,jpphy)
                !     Squared mortality of Phyto similar to a sedimentation term during
                !     blooms (Doney et al. 1996)
@@ -102,5 +102,5 @@
                !     increased when nutrients are limiting phytoplankton growth
                !     as observed for instance in case of iron limitation.
-               ztortp = mprat * xstep * zcompaph / ( xkmort + trn(ji,jj,jk,jpphy) ) * zsizerat
+               ztortp = mprat * xstep * zcompaph / ( xkmort + trb(ji,jj,jk,jpphy) ) * zsizerat
 
                zmortp = zrespp + ztortp
@@ -108,6 +108,6 @@
                !   Update the arrays TRA which contains the biological sources and sinks
 
-               zfactfe = trn(ji,jj,jk,jpnfe)/(trn(ji,jj,jk,jpphy)+rtrn)
-               zfactch = trn(ji,jj,jk,jpnch)/(trn(ji,jj,jk,jpphy)+rtrn)
+               zfactfe = trb(ji,jj,jk,jpnfe)/(trb(ji,jj,jk,jpphy)+rtrn)
+               zfactch = trb(ji,jj,jk,jpnch)/(trb(ji,jj,jk,jpphy)+rtrn)
                tra(ji,jj,jk,jpphy) = tra(ji,jj,jk,jpphy) - zmortp
                tra(ji,jj,jk,jpnch) = tra(ji,jj,jk,jpnch) - zmortp * zfactch
@@ -172,5 +172,5 @@
             DO ji = 1, jpi
 
-               zcompadi = MAX( ( trn(ji,jj,jk,jpdia) - 1e-9), 0. )
+               zcompadi = MAX( ( trb(ji,jj,jk,jpdia) - 1e-9), 0. )
 
                !    Aggregation term for diatoms is increased in case of nutrient
@@ -186,9 +186,9 @@
                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)
+               zrespp2 = 1.e6 * zstep * (  wchld + wchldm * zlim1 ) * xdiss(ji,jj,jk) * zcompadi * trb(ji,jj,jk,jpdia)
 
                !     Phytoplankton mortality. 
                !     ------------------------
-               ztortp2 = mprat2 * zstep * trn(ji,jj,jk,jpdia)  / ( xkmort + trn(ji,jj,jk,jpdia) ) * zcompadi 
+               ztortp2 = mprat2 * zstep * trb(ji,jj,jk,jpdia)  / ( xkmort + trb(ji,jj,jk,jpdia) ) * zcompadi 
 
                zmortp2 = zrespp2 + ztortp2
@@ -196,7 +196,7 @@
                !   Update the arrays tra which contains the biological sources and sinks
                !   ---------------------------------------------------------------------
-               zfactch = trn(ji,jj,jk,jpdch) / ( trn(ji,jj,jk,jpdia) + rtrn )
-               zfactfe = trn(ji,jj,jk,jpdfe) / ( trn(ji,jj,jk,jpdia) + rtrn )
-               zfactsi = trn(ji,jj,jk,jpdsi) / ( trn(ji,jj,jk,jpdia) + rtrn )
+               zfactch = trb(ji,jj,jk,jpdch) / ( trb(ji,jj,jk,jpdia) + rtrn )
+               zfactfe = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn )
+               zfactsi = trb(ji,jj,jk,jpdsi) / ( trb(ji,jj,jk,jpdia) + rtrn )
                tra(ji,jj,jk,jpdia) = tra(ji,jj,jk,jpdia) - zmortp2 
                tra(ji,jj,jk,jpdch) = tra(ji,jj,jk,jpdch) - zmortp2 * zfactch
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90	(revision 5462)
@@ -35,4 +35,5 @@
    REAL(wp) :: parlux      !: Fraction of shortwave as PAR
    REAL(wp) :: xparsw                 !: parlux/3
+   REAL(wp) :: xsi0r                 !:  1. /rn_si0
 
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
@@ -42,5 +43,7 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: enano, ediat   !: PAR for phyto, nano and diat 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etot_ndcy      !: PAR over 24h in case of diurnal cycle
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emoy           !: averaged PAR in the mixed layer
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ekb, ekg, ekr  !: wavelength (Red-Green-Blue)
 
    INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
@@ -57,5 +60,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_opt( kt, jnt )
+   SUBROUTINE p4z_opt( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_opt  ***
@@ -67,12 +70,12 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt, jnt   ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt   ! ocean time step
       !
       INTEGER  ::   ji, jj, jk
       INTEGER  ::   irgb
-      REAL(wp) ::   zchl, zxsi0r
+      REAL(wp) ::   zchl
       REAL(wp) ::   zc0 , zc1 , zc2, zc3, z1_dep
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp, zetmp1, zetmp2
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zekg, zekr, zekb, ze0, ze1, ze2, ze3
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4, zqsr100
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpar, ze0, ze1, ze2, ze3
       !!---------------------------------------------------------------------
       !
@@ -80,15 +83,14 @@
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 ) 
-      CALL wrk_alloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
-
-      IF( jnt == 1 .AND. ln_varpar ) CALL p4z_optsbc( kt )
+      CALL wrk_alloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_alloc( jpi, jpj, jpk, zpar, ze0, ze1, ze2, ze3 )
+
+      IF( knt == 1 .AND. ln_varpar ) CALL p4z_opt_sbc( kt )
 
       !     Initialisation of variables used to compute PAR
       !     -----------------------------------------------
-      ze1(:,:,jpk) = 0._wp
-      ze2(:,:,jpk) = 0._wp
-      ze3(:,:,jpk) = 0._wp
-
+      ze1(:,:,:) = 0._wp
+      ze2(:,:,:) = 0._wp
+      ze3(:,:,:) = 0._wp
       !                                        !* attenuation coef. function of Chlorophyll and wavelength (Red-Green-Blue)
       DO jk = 1, jpkm1                         !  --------------------------------------------------------
@@ -97,103 +99,59 @@
 !CDIR NOVERRCHK
             DO ji = 1, jpi
-               zchl = ( trn(ji,jj,jk,jpnch) + trn(ji,jj,jk,jpdch) + rtrn ) * 1.e6
+               zchl = ( trb(ji,jj,jk,jpnch) + trb(ji,jj,jk,jpdch) + rtrn ) * 1.e6
                zchl = MIN(  10. , MAX( 0.05, zchl )  )
                irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
                !                                                         
-               zekb(ji,jj,jk) = xkrgb(1,irgb) * fse3t(ji,jj,jk)
-               zekg(ji,jj,jk) = xkrgb(2,irgb) * fse3t(ji,jj,jk)
-               zekr(ji,jj,jk) = xkrgb(3,irgb) * fse3t(ji,jj,jk)
+               ekb(ji,jj,jk) = xkrgb(1,irgb) * fse3t(ji,jj,jk)
+               ekg(ji,jj,jk) = xkrgb(2,irgb) * fse3t(ji,jj,jk)
+               ekr(ji,jj,jk) = xkrgb(3,irgb) * fse3t(ji,jj,jk)
             END DO
          END DO
       END DO
-
-
       !                                        !* Photosynthetically Available Radiation (PAR)
       !                                        !  --------------------------------------
-
-      IF( ln_varpar ) THEN
-         ze1(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
-         ze2(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
-         ze3(:,:,1) = par_varsw(:,:) * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
+      IF( l_trcdm2dc ) THEN                     !  diurnal cycle
+         ! 1% of qsr to compute euphotic layer
+         zqsr100(:,:) = 0.01 * qsr_mean(:,:)     !  daily mean qsr
+         !
+         CALL p4z_opt_par( kt, qsr_mean, ze1, ze2, ze3 ) 
+         !
+         DO jk = 1, nksrp      
+            etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
+            enano    (:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+            ediat    (:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk)
+         END DO
+         !
+         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         !
+         DO jk = 1, nksrp      
+            etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
+         END DO
+         !
       ELSE
-         ze1(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekb(:,:,1) )
-         ze2(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekg(:,:,1) )
-         ze3(:,:,1) = xparsw         * qsr(:,:) * EXP( -0.5 * zekr(:,:,1) )
-      ENDIF
-
-!CDIR NOVERRCHK
-      DO jj = 1, jpj
-!CDIR NOVERRCHK
-         DO ji = 1, jpi
-            zc1 = ze1(ji,jj,1)
-            zc2 = ze2(ji,jj,1) 
-            zc3 = ze3(ji,jj,1)
-            etot (ji,jj,1) = (       zc1 +        zc2 +       zc3 )
-            enano(ji,jj,1) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
-            ediat(ji,jj,1) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
-         END DO
-      END DO
-
-    
-      DO jk = 2, nksrp      
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-               zc1 = ze1(ji,jj,jk-1) * EXP( -0.5 * ( zekb(ji,jj,jk-1) + zekb(ji,jj,jk) ) )
-               zc2 = ze2(ji,jj,jk-1) * EXP( -0.5 * ( zekg(ji,jj,jk-1) + zekg(ji,jj,jk) ) )
-               zc3 = ze3(ji,jj,jk-1) * EXP( -0.5 * ( zekr(ji,jj,jk-1) + zekr(ji,jj,jk) ) )
-               ze1  (ji,jj,jk) = zc1
-               ze2  (ji,jj,jk) = zc2
-               ze3  (ji,jj,jk) = zc3
-               etot (ji,jj,jk) = (       zc1 +        zc2 +       zc3 )
-               enano(ji,jj,jk) = ( 2.1 * zc1 + 0.42 * zc2 + 0.4 * zc3 )
-               ediat(ji,jj,jk) = ( 1.6 * zc1 + 0.69 * zc2 + 0.7 * zc3 )
-            END DO
-         END DO
-      END DO
+         ! 1% of qsr to compute euphotic layer
+         zqsr100(:,:) = 0.01 * qsr(:,:)
+         !
+         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3 ) 
+         !
+         DO jk = 1, nksrp      
+            etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
+            enano(:,:,jk) =  2.1 * ze1(:,:,jk) + 0.42 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
+            ediat(:,:,jk) =  1.6 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.7 * ze3(:,:,jk)
+         END DO
+         etot_ndcy(:,:,:) =  etot(:,:,:) 
+      ENDIF
+
 
       IF( ln_qsr_bio ) THEN                    !* heat flux accros w-level (used in the dynamics)
          !                                     !  ------------------------
-         zxsi0r = 1.e0 / rn_si0
-         !
-         ze0(:,:,1) = rn_abs * qsr(:,:)
-         !                                                    ! surface value : separation in R-G-B + near surface
-         IF( ln_varpar ) THEN
-            ze0(:,:,1) = ( 1. - 3. * par_varsw(:,:) ) * qsr(:,:)
-            ze1(:,:,1) = par_varsw(:,:)               * qsr(:,:)         
-            ze2(:,:,1) = par_varsw(:,:)               * qsr(:,:)
-            ze3(:,:,1) = par_varsw(:,:)               * qsr(:,:)
-         ELSE
-            ze0(:,:,1) = ( 1. - 3. * xparsw )  * qsr(:,:)
-            ze1(:,:,1) = xparsw                * qsr(:,:)         
-            ze2(:,:,1) = xparsw                * qsr(:,:)
-            ze3(:,:,1) = xparsw                * qsr(:,:)
-         ENDIF
+         CALL p4z_opt_par( kt, qsr, ze1, ze2, ze3, pe0=ze0 )
+         !
          etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
-         !
-         !
          DO jk = 2, nksrp + 1
-!CDIR NOVERRCHK
-            DO jj = 1, jpj
-!CDIR NOVERRCHK
-               DO ji = 1, jpi
-                  zc0 = ze0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * zxsi0r )
-                  zc1 = ze1(ji,jj,jk-1) * EXP( -zekb(ji,jj,jk-1 ) )
-                  zc2 = ze2(ji,jj,jk-1) * EXP( -zekg(ji,jj,jk-1 ) )
-                  zc3 = ze3(ji,jj,jk-1) * EXP( -zekr(ji,jj,jk-1 ) )
-                  ze0(ji,jj,jk) = zc0
-                  ze1(ji,jj,jk) = zc1
-                  ze2(ji,jj,jk) = zc2
-                  ze3(ji,jj,jk) = zc3
-                  etot3(ji,jj,jk) = ( zc0 + zc1 + zc2 + zc3 ) * tmask(ji,jj,jk)
-              END DO
-              !
-            END DO
-            !
-        END DO
-        !
-      ENDIF
-
+            etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
+         END DO
+         !                                     !  ------------------------
+      ENDIF
       !                                        !* Euphotic depth and level
       neln(:,:) = 1                            !  ------------------------
@@ -203,7 +161,7 @@
          DO jj = 1, jpj
            DO ji = 1, jpi
-              IF( etot(ji,jj,jk) * tmask(ji,jj,jk) >= 0.0043 * qsr(ji,jj) )  THEN
+              IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >= 0.43 * 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_mxl_trc_zint
+                 !                                     ! nb: ensure the compatibility with nmld_trc definition in trd_mld_trc_zint
                  heup(ji,jj) = fsdepw(ji,jj,jk+1)      ! Euphotic layer depth
               ENDIF
@@ -211,12 +169,12 @@
         END DO
       END DO
- 
+      !
       heup(:,:) = MIN( 300., heup(:,:) )
-
       !                                        !* mean light over the mixed layer
       zdepmoy(:,:)   = 0.e0                    !  -------------------------------
-      zetmp  (:,:)   = 0.e0
       zetmp1 (:,:)   = 0.e0
       zetmp2 (:,:)   = 0.e0
+      zetmp3 (:,:)   = 0.e0
+      zetmp4 (:,:)   = 0.e0
 
       DO jk = 1, nksrp
@@ -226,7 +184,8 @@
             DO ji = 1, jpi
                IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
-                  zetmp  (ji,jj) = zetmp  (ji,jj) + etot (ji,jj,jk) * fse3t(ji,jj,jk)
-                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + enano(ji,jj,jk) * fse3t(ji,jj,jk)
-                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + ediat(ji,jj,jk) * fse3t(ji,jj,jk)
+                  zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * fse3t(ji,jj,jk) ! remineralisation
+                  zetmp2 (ji,jj) = zetmp2 (ji,jj) + etot_ndcy(ji,jj,jk) * fse3t(ji,jj,jk) ! production
+                  zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
+                  zetmp4 (ji,jj) = zetmp4 (ji,jj) + ediat    (ji,jj,jk) * fse3t(ji,jj,jk) ! production
                   zdepmoy(ji,jj) = zdepmoy(ji,jj) + fse3t(ji,jj,jk)
                ENDIF
@@ -235,5 +194,6 @@
       END DO
       !
-      emoy(:,:,:) = etot(:,:,:)
+      emoy(:,:,:) = etot(:,:,:)       ! remineralisation
+      zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
       DO jk = 1, nksrp
@@ -244,26 +204,28 @@
                IF( fsdepw(ji,jj,jk+1) <= hmld(ji,jj) ) THEN
                   z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
-                  emoy (ji,jj,jk) = zetmp (ji,jj) * z1_dep
-                  enano(ji,jj,jk) = zetmp1(ji,jj) * z1_dep
-                  ediat(ji,jj,jk) = zetmp2(ji,jj) * z1_dep
+                  emoy (ji,jj,jk) = zetmp1(ji,jj) * z1_dep
+                  zpar (ji,jj,jk) = zetmp2(ji,jj) * z1_dep
+                  enano(ji,jj,jk) = zetmp3(ji,jj) * z1_dep
+                  ediat(ji,jj,jk) = zetmp4(ji,jj) * z1_dep
                ENDIF
             END DO
          END DO
       END DO
-
+      !
       IF( lk_iomput ) THEN
-        IF( jnt == nrdttrc  ) THEN
-           IF( iom_use( "Heup" ) ) CALL iom_put( "Heup", heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
-           IF( iom_use( "PAR"  ) ) CALL iom_put( "PAR" , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
+        IF( knt == nrdttrc ) THEN
+           IF( iom_use( "Heup"  ) ) CALL iom_put( "Heup" , heup(:,:  ) * tmask(:,:,1) )  ! euphotic layer deptht
+           IF( iom_use( "PARDM" ) ) CALL iom_put( "PARDM", zpar(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
+           IF( iom_use( "PAR"   ) ) CALL iom_put( "PAR"  , emoy(:,:,:) * tmask(:,:,:) )  ! Photosynthetically Available Radiation
         ENDIF
       ELSE
          IF( ln_diatrc ) THEN        ! save output diagnostics
-            trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)  
+            trc2d(:,:,  jp_pcs0_2d + 10) = heup(:,:  ) * tmask(:,:,1)
             trc3d(:,:,:,jp_pcs0_3d + 3)  = etot(:,:,:) * tmask(:,:,:)
          ENDIF
       ENDIF
       !
-      CALL wrk_dealloc( jpi, jpj,      zdepmoy, zetmp, zetmp1, zetmp2 )
-      CALL wrk_dealloc( jpi, jpj, jpk, zekg, zekr, zekb, ze0, ze1, ze2, ze3 )
+      CALL wrk_dealloc( jpi, jpj,      zqsr100, zdepmoy, zetmp1, zetmp2, zetmp3, zetmp4 )
+      CALL wrk_dealloc( jpi, jpj, jpk, zpar,  ze0, ze1, ze2, ze3 )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_opt')
@@ -271,7 +233,75 @@
    END SUBROUTINE p4z_opt
 
-   SUBROUTINE p4z_optsbc( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  routine p4z_optsbc  ***
+   SUBROUTINE p4z_opt_par( kt, pqsr, pe1, pe2, pe3, pe0 ) 
+      !!----------------------------------------------------------------------
+      !!                  ***  routine p4z_opt_par  ***
+      !!
+      !! ** purpose :   compute PAR of each wavelength (Red-Green-Blue)
+      !!                for a given shortwave radiation
+      !!
+      !!----------------------------------------------------------------------
+      !! * arguments
+      INTEGER, INTENT(in)                                       ::  kt            !   ocean time-step
+      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  
+      !! * local variables
+      INTEGER    ::   ji, jj, jk     ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj)     ::  zqsr          !   shortwave
+      !!----------------------------------------------------------------------
+
+      !  Real shortwave
+      IF( ln_varpar ) THEN  ;  zqsr(:,:) = par_varsw(:,:) * pqsr(:,:)
+      ELSE                  ;  zqsr(:,:) = xparsw         * pqsr(:,:)
+      ENDIF
+      !
+      IF( PRESENT( pe0 ) ) THEN     !  W-level
+         !
+         pe0(:,:,1) = pqsr(:,:) - 3. * zqsr(:,:)    !   ( 1 - 3 * alpha ) * q
+         pe1(:,:,1) = zqsr(:,:)         
+         pe2(:,:,1) = zqsr(:,:)
+         pe3(:,:,1) = zqsr(:,:)
+         !
+         DO jk = 2, nksrp + 1
+!CDIR NOVERRCHK
+            DO jj = 1, jpj
+!CDIR NOVERRCHK
+               DO ji = 1, jpi
+                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -fse3t(ji,jj,jk-1) * xsi0r )
+                  pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -ekb(ji,jj,jk-1 ) )
+                  pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -ekg(ji,jj,jk-1 ) )
+                  pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -ekr(ji,jj,jk-1 ) )
+               END DO
+              !
+            END DO
+            !
+         END DO
+        !
+      ELSE   ! T- level
+        !
+        pe1(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekb(:,:,1) )
+        pe2(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekg(:,:,1) )
+        pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
+        !
+        DO jk = 2, nksrp      
+!CDIR NOVERRCHK
+           DO jj = 1, jpj
+!CDIR NOVERRCHK
+              DO ji = 1, jpi
+                 pe1(ji,jj,jk) = pe1(ji,jj,jk-1) * EXP( -0.5 * ( ekb(ji,jj,jk-1) + ekb(ji,jj,jk) ) )
+                 pe2(ji,jj,jk) = pe2(ji,jj,jk-1) * EXP( -0.5 * ( ekg(ji,jj,jk-1) + ekg(ji,jj,jk) ) )
+                 pe3(ji,jj,jk) = pe3(ji,jj,jk-1) * EXP( -0.5 * ( ekr(ji,jj,jk-1) + ekr(ji,jj,jk) ) )
+              END DO
+           END DO
+        END DO    
+        !
+      ENDIF
+      ! 
+   END SUBROUTINE p4z_opt_par
+
+
+   SUBROUTINE p4z_opt_sbc( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine p4z_opt_sbc  ***
       !!
       !! ** purpose :   read and interpolate the variable PAR fraction
@@ -284,5 +314,5 @@
       !!----------------------------------------------------------------------
       !! * arguments
-      INTEGER, INTENT( in  ) ::   kt   ! ocean time step
+      INTEGER ,                INTENT(in) ::   kt     ! ocean time step
 
       !! * local declarations
@@ -297,5 +327,5 @@
          IF( kt == nit000 .OR. ( kt /= nit000 .AND. ntimes_par > 1 ) ) THEN
             CALL fld_read( kt, 1, sf_par )
-            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) )/3.0
+            par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
          ENDIF
       ENDIF
@@ -303,5 +333,5 @@
       IF( nn_timing == 1 )  CALL timing_stop('p4z_optsbc')
       !
-   END SUBROUTINE p4z_optsbc
+   END SUBROUTINE p4z_opt_sbc
 
    SUBROUTINE p4z_opt_init
@@ -347,4 +377,5 @@
       !
       xparsw = parlux / 3.0
+      xsi0r  = 1.e0 / rn_si0
       !
       ! Variable PAR at the surface of the ocean
@@ -372,8 +403,12 @@
       IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
       !
-                         etot (:,:,:) = 0._wp
-                         enano(:,:,:) = 0._wp
-                         ediat(:,:,:) = 0._wp
-      IF( ln_qsr_bio )   etot3(:,:,:) = 0._wp
+                         ekr      (:,:,:) = 0._wp
+                         ekb      (:,:,:) = 0._wp
+                         ekg      (:,:,:) = 0._wp
+                         etot     (:,:,:) = 0._wp
+                         etot_ndcy(:,:,:) = 0._wp
+                         enano    (:,:,:) = 0._wp
+                         ediat    (:,:,:) = 0._wp
+      IF( ln_qsr_bio )   etot3    (:,:,:) = 0._wp
       ! 
       IF( nn_timing == 1 )  CALL timing_stop('p4z_opt_init')
@@ -386,5 +421,7 @@
       !!                     ***  ROUTINE p4z_opt_alloc  ***
       !!----------------------------------------------------------------------
-      ALLOCATE( enano(jpi,jpj,jpk), ediat(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 
+      ALLOCATE( ekb(jpi,jpj,jpk)      , ekr(jpi,jpj,jpk), ekg(jpi,jpj,jpk),   &
+        &       enano(jpi,jpj,jpk)    , ediat(jpi,jpj,jpk), &
+        &       etot_ndcy(jpi,jpj,jpk), emoy (jpi,jpj,jpk), STAT=p4z_opt_alloc ) 
          !
       IF( p4z_opt_alloc /= 0 ) CALL ctl_warn('p4z_opt_alloc : failed to allocate arrays.')
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zprod.F90	(revision 5462)
@@ -64,5 +64,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_prod( kt , jnt )
+   SUBROUTINE p4z_prod( kt , knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_prod  ***
@@ -74,5 +74,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER, INTENT(in) :: kt, knt
       !
       INTEGER  ::   ji, jj, jk
@@ -129,19 +129,17 @@
       END DO
 
-      IF( ln_newprod ) THEN
-         ! Impact of the day duration on phytoplankton growth
-         DO jk = 1, jpkm1
-            DO jj = 1 ,jpj
-               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)
-                  ENDIF
-               END DO
-            END DO
-         END DO
-      ENDIF
+      ! Impact of the day duration on phytoplankton growth
+      DO jk = 1, jpkm1
+         DO jj = 1 ,jpj
+            DO ji = 1, jpi
+               IF( etot_ndcy(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)
+               ENDIF
+            END DO
+         END DO
+      END DO
 
       ! Maximum light intensity
@@ -157,17 +155,17 @@
                DO ji = 1, jpi
                   ! Computation of the P-I slope for nanos and diatoms
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                       ztn         = MAX( 0., tsn(ji,jj,jk,jp_tem) - 15. )
                       zadap       = xadap * ztn / ( 2.+ ztn )
-                      zconctemp   = MAX( 0.e0 , trn(ji,jj,jk,jpdia) - xsizedia )
-                      zconctemp2  = trn(ji,jj,jk,jpdia) - zconctemp
+                      zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
+                      zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
                       znanotot    = enano(ji,jj,jk) * zstrn(ji,jj)
                       zdiattot    = ediat(ji,jj,jk) * zstrn(ji,jj)
                       !
                       zpislopead (ji,jj,jk) = pislope * ( 1.+ zadap  * EXP( -znanotot ) )  &
-                         &                   * trn(ji,jj,jk,jpnch) /( trn(ji,jj,jk,jpphy) * 12. + rtrn)
+                         &                   * trb(ji,jj,jk,jpnch) /( trb(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)
+                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp) / ( trb(ji,jj,jk,jpdia) + rtrn )   &
+                         &                   * trb(ji,jj,jk,jpdch) /( trb(ji,jj,jk,jpdia) * 12. + rtrn)
 
                       ! Computation of production function for Carbon
@@ -196,30 +194,32 @@
 
                   ! Computation of the P-I slope for nanos and diatoms
-                  IF( etot(ji,jj,jk) > 1.E-3 ) THEN
+                  IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                       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
+                      zconctemp   = MAX( 0.e0 , trb(ji,jj,jk,jpdia) - xsizedia )
+                      zconctemp2  = trb(ji,jj,jk,jpdia) - zconctemp
+                      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 * enano(ji,jj,jk) ) )
-                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trn(ji,jj,jk,jpdia) + rtrn )
-
-                      zpislopen =  zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch)                &
-                        &          / ( trn(ji,jj,jk,jpphy) * 12.                  + rtrn )   &
+                      zpislopead (ji,jj,jk) = pislope  * ( 1.+ zadap  * EXP( -znanotot ) )
+                      zpislopead2(ji,jj,jk) = (pislope * zconctemp2 + pislope2 * zconctemp)  / ( trb(ji,jj,jk,jpdia) + rtrn )
+
+                      zpislopen =  zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch)                &
+                        &          / ( trb(ji,jj,jk,jpphy) * 12.                  + rtrn )   &
                         &          / ( prmax(ji,jj,jk) * rday * xlimphy(ji,jj,jk) + rtrn )
 
-                      zpislope2n = zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch)                &
-                        &          / ( trn(ji,jj,jk,jpdia) * 12.                  + rtrn )   &
+                      zpislope2n = zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch)                &
+                        &          / ( trb(ji,jj,jk,jpdia) * 12.                  + rtrn )   &
                         &          / ( prmax(ji,jj,jk) * rday * xlimdia(ji,jj,jk) + rtrn )
 
                       ! Computation of production function for Carbon
                       !  ---------------------------------------------
-                      zprbio(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) ) )
-                      zprdia(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) ) )
+                      zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen  * znanotot ) )
+                      zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpislope2n * zdiattot ) )
 
                       !  Computation of production function for Chlorophyll
                       !--------------------------------------------------
-                      zprnch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislopen  * enano(ji,jj,jk) * zstrn(ji,jj) ) )
-                      zprdch(ji,jj,jk) = prmax(ji,jj,jk) * ( 1.- EXP( -zpislope2n * ediat(ji,jj,jk) * zstrn(ji,jj) ) )
+                      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) ) )
                   ENDIF
                END DO
@@ -252,5 +252,5 @@
             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
                    !    ------------------------
@@ -258,8 +258,8 @@
                    !    Si/C is arbitrariliy increased for very high Si concentrations
                    !    to mimic the very high ratios observed in the Southern Ocean (silpot2)
-                  zlim  = trn(ji,jj,jk,jpsil) / ( trn(ji,jj,jk,jpsil) + xksi1 )
+                  zlim  = trb(ji,jj,jk,jpsil) / ( trb(ji,jj,jk,jpsil) + xksi1 )
                   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 = trn(ji,jj,jk,jpsil) * trn(ji,jj,jk,jpsil) * trn(ji,jj,jk,jpsil)
+                  zsiborn = trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil) * trb(ji,jj,jk,jpsil)
                   IF (gphit(ji,jj) < -30 ) THEN
                     zsilfac2 = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
@@ -302,10 +302,10 @@
 !CDIR NOVERRCHK
             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.
-                  zprorca(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * trn(ji,jj,jk,jpphy) * rfact2
+                  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 )
                   !
-                  zratio = trn(ji,jj,jk,jpnfe) / ( trn(ji,jj,jk,jpphy) + rtrn )
+                  zratio = trb(ji,jj,jk,jpnfe) / ( trb(ji,jj,jk,jpphy) + rtrn )
                   zratio = zratio / fecnm 
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
@@ -313,10 +313,10 @@
                   &             * ( 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 * trn(ji,jj,jk,jpphy) * rfact2
+                  &             * zmax * trb(ji,jj,jk,jpphy) * rfact2
                   !  production terms for diatomees
-                  zprorcad(ji,jj,jk) = zprdia(ji,jj,jk) * xlimdia(ji,jj,jk) * trn(ji,jj,jk,jpdia) * rfact2
+                  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 = trn(ji,jj,jk,jpdfe) / ( trn(ji,jj,jk,jpdia) + rtrn )
+                  zratio = trb(ji,jj,jk,jpdfe) / ( trb(ji,jj,jk,jpdia) + rtrn )
                   zratio = zratio / fecdm 
                   zmax   = MAX( 0., ( 1. - zratio ) / ABS( 1.05 - zratio ) ) 
@@ -324,5 +324,5 @@
                   &             * ( 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) )  &
-                  &             * zmax * trn(ji,jj,jk,jpdia) * rfact2
+                  &             * zmax * trb(ji,jj,jk,jpdia) * rfact2
                ENDIF
             END DO
@@ -341,5 +341,5 @@
                      zprdch(ji,jj,jk) = zprdch(ji,jj,jk) * zmixdiat(ji,jj)
                   ENDIF
-                  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) * zstrn(ji,jj)
@@ -365,17 +365,17 @@
 !CDIR NOVERRCHK
                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) * zstrn(ji,jj)
-                     zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trn(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk)
+                     znanotot = enano(ji,jj,jk)
+                     zprod = rday * zprorca(ji,jj,jk) * zprnch(ji,jj,jk) * trb(ji,jj,jk,jpphy) * xlimphy(ji,jj,jk)
                      zprochln(ji,jj,jk) = chlcmin * 12. * zprorca (ji,jj,jk)
                      zprochln(ji,jj,jk) = zprochln(ji,jj,jk) + (chlcnm-chlcmin) * 144. * zprod            &
-                     &                    / ( zpislopead(ji,jj,jk) * trn(ji,jj,jk,jpnch) * znanotot +rtrn )
+                     &                    / ( zpislopead(ji,jj,jk) * trb(ji,jj,jk,jpnch) * znanotot +rtrn )
                      !  production terms for diatomees ( chlorophyll )
-                     zdiattot = ediat(ji,jj,jk) * zstrn(ji,jj)
-                     zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trn(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk)
+                     zdiattot = ediat(ji,jj,jk)
+                     zprod = rday * zprorcad(ji,jj,jk) * zprdch(ji,jj,jk) * trb(ji,jj,jk,jpdia) * xlimdia(ji,jj,jk)
                      zprochld(ji,jj,jk) = chlcmin * 12. * zprorcad(ji,jj,jk)
                      zprochld(ji,jj,jk) = zprochld(ji,jj,jk) + (chlcdm-chlcmin) * 144. * zprod             &
-                     &                    / ( zpislopead2(ji,jj,jk) * trn(ji,jj,jk,jpdch) * zdiattot +rtrn )
+                     &                    / ( zpislopead2(ji,jj,jk) * trb(ji,jj,jk,jpdch) * zdiattot +rtrn )
                   ENDIF
                END DO
@@ -414,9 +414,9 @@
 
     ! Total primary production per year
-    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. jnt == nrdttrc )  )  &
+    IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
          & tpp = glob_sum( ( zprorca(:,:,:) + zprorcad(:,:,:) ) * cvol(:,:,:) )
 
     IF( lk_iomput ) THEN
-       IF( jnt == nrdttrc ) THEN
+       IF( knt == nrdttrc ) THEN
           CALL wrk_alloc( jpi, jpj,      zw2d )
           CALL wrk_alloc( jpi, jpj, jpk, zw3d )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zrem.F90	(revision 5462)
@@ -59,5 +59,5 @@
 CONTAINS
 
-   SUBROUTINE p4z_rem( kt, jnt )
+   SUBROUTINE p4z_rem( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_rem  ***
@@ -68,5 +68,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       !
       INTEGER  ::   ji, jj, jk
@@ -104,5 +104,5 @@
                zdep = MAX( hmld(ji,jj), heup(ji,jj) )
                IF( fsdept(ji,jj,jk) < zdep ) THEN
-                  zdepbac(ji,jj,jk) = MIN( 0.7 * ( trn(ji,jj,jk,jpzoo) + 2.* trn(ji,jj,jk,jpmes) ), 4.e-6 )
+                  zdepbac(ji,jj,jk) = MIN( 0.7 * ( trb(ji,jj,jk,jpzoo) + 2.* trb(ji,jj,jk,jpmes) ), 4.e-6 )
                   ztempbac(ji,jj)   = zdepbac(ji,jj,jk)
                ELSE
@@ -119,6 +119,6 @@
             DO ji = 1, jpi
                ! denitrification factor computed from O2 levels
-               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trn(ji,jj,jk,jpoxy) )    &
-                  &                                / ( oxymin + trn(ji,jj,jk,jpoxy) )  )
+               nitrfac(ji,jj,jk) = MAX(  0.e0, 0.4 * ( 6.e-6  - trb(ji,jj,jk,jpoxy) )    &
+                  &                                / ( oxymin + trb(ji,jj,jk,jpoxy) )  )
                nitrfac(ji,jj,jk) = MIN( 1., nitrfac(ji,jj,jk) )
             END DO
@@ -140,10 +140,10 @@
                ! Ammonification in oxic waters with oxygen consumption
                ! -----------------------------------------------------
-               zolimit = zremik * ( 1.- nitrfac(ji,jj,jk) ) * trn(ji,jj,jk,jpdoc) 
-               zolimi(ji,jj,jk) = MIN( ( trn(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimit ) 
+               zolimit = zremik * ( 1.- nitrfac(ji,jj,jk) ) * trb(ji,jj,jk,jpdoc) 
+               zolimi(ji,jj,jk) = MIN( ( trb(ji,jj,jk,jpoxy) - rtrn ) / o2ut, zolimit ) 
                ! Ammonification in suboxic waters with denitrification
                ! -------------------------------------------------------
-               denitr(ji,jj,jk)  = MIN(  ( trn(ji,jj,jk,jpno3) - rtrn ) / rdenit,   &
-                  &                     zremik * nitrfac(ji,jj,jk) * trn(ji,jj,jk,jpdoc)  )
+               denitr(ji,jj,jk)  = MIN(  ( trb(ji,jj,jk,jpno3) - rtrn ) / rdenit,   &
+                  &                     zremik * nitrfac(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  )
                !
                zolimi (ji,jj,jk) = MAX( 0.e0, zolimi (ji,jj,jk) )
@@ -165,6 +165,6 @@
                ! below 2 umol/L. Inhibited at strong light 
                ! ----------------------------------------------------------
-               zonitr  =nitrif * zstep * trn(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) ) * ( 1.- nitrfac(ji,jj,jk) ) 
-               denitnh4(ji,jj,jk) = nitrif * zstep * trn(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 
+               zonitr  =nitrif * zstep * trb(ji,jj,jk,jpnh4) / ( 1.+ emoy(ji,jj,jk) ) * ( 1.- nitrfac(ji,jj,jk) ) 
+               denitnh4(ji,jj,jk) = nitrif * zstep * trb(ji,jj,jk,jpnh4) * nitrfac(ji,jj,jk) 
                ! Update of the tracers trends
                ! ----------------------------
@@ -192,5 +192,5 @@
                ! ----------------------------------------------------------
                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) )    &
+                  &              * trb(ji,jj,jk,jpfer) / ( 2.5E-10 + trb(ji,jj,jk,jpfer) )    &
                   &              * zdepprod(ji,jj,jk) * zdepbac(ji,jj,jk)
 #if defined key_kriest
@@ -228,11 +228,11 @@
                ! means a disaggregation constant about 0.5 the value in oxic zones
                ! -----------------------------------------------------------------
-               zorem  = zremip * trn(ji,jj,jk,jppoc)
-               zofer  = zremip * trn(ji,jj,jk,jpsfe)
+               zorem  = zremip * trb(ji,jj,jk,jppoc)
+               zofer  = zremip * trb(ji,jj,jk,jpsfe)
 #if ! defined key_kriest
-               zorem2 = zremip * trn(ji,jj,jk,jpgoc)
-               zofer2 = zremip * trn(ji,jj,jk,jpbfe)
+               zorem2 = zremip * trb(ji,jj,jk,jpgoc)
+               zofer2 = zremip * trb(ji,jj,jk,jpbfe)
 #else
-               zorem2 = zremip * trn(ji,jj,jk,jpnum)
+               zorem2 = zremip * trb(ji,jj,jk,jpnum)
 #endif
 
@@ -272,5 +272,5 @@
                ! Remineralization rate of BSi depedant on T and saturation
                ! ---------------------------------------------------------
-               zsatur   = ( sio3eq(ji,jj,jk) - trn(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
+               zsatur   = ( sio3eq(ji,jj,jk) - trb(ji,jj,jk,jpsil) ) / ( sio3eq(ji,jj,jk) + rtrn )
                zsatur   = MAX( rtrn, zsatur )
                zsatur2  = ( 1. + tsn(ji,jj,jk,jp_tem) / 400.)**37
@@ -287,5 +287,5 @@
                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)
+               zosil    = zsiremin * trb(ji,jj,jk,jpgsi)
                !
                tra(ji,jj,jk,jpgsi) = tra(ji,jj,jk,jpgsi) - zosil
@@ -315,5 +315,5 @@
       END DO
 
-      IF( jnt == nrdttrc ) THEN
+      IF( knt == nrdttrc ) THEN
           CALL wrk_alloc( jpi, jpj, jpk, zw3d )
           zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsbc.F90	(revision 5462)
@@ -85,5 +85,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Header:$ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -136,5 +136,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  zcoef = ryyss * cvol(ji,jj,1) 
+                  zcoef = ryyss * e1e2t(ji,jj) * h_rnf(ji,jj) 
                   rivalk(ji,jj) =   sf_river(jr_dic)%fnow(ji,jj,1)                                    &
                      &              * 1.E3        / ( 12. * zcoef + rtrn )
@@ -187,4 +187,6 @@
       INTEGER  :: ierr, ierr1, ierr2, ierr3
       INTEGER  :: ios                 ! Local integer output status for namelist read
+      INTEGER  :: ik50                !  last level where depth less than 50 m
+      INTEGER  :: isrow             ! index for ORCA1 starting row
       REAL(wp) :: zexpide, zdenitide, zmaskt
       REAL(wp) :: ztimes_dust, ztimes_riv, ztimes_ndep 
@@ -216,4 +218,13 @@
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nampissbc in configuration namelist', lwp )
       IF(lwm) WRITE ( numonp, nampissbc )
+
+      IF ( ( nn_ice_tr >= 0 ) .AND. ln_ironice ) THEN
+         IF(lwp) THEN
+            WRITE(numout,*) ' ln_ironice incompatible with nn_ice_tr = ', nn_ice_tr
+            WRITE(numout,*) ' Specify your sea ice iron concentration in nampisice instead '
+            WRITE(numout,*) ' ln_ironice is forced to .FALSE. '
+            ln_ironice = .FALSE.
+         ENDIF
+      ENDIF
 
       IF(lwp) THEN
@@ -247,4 +258,11 @@
       ENDIF
 
+      ! set the number of level over which river runoffs are applied 
+      ! online configuration : computed in sbcrnf
+      IF( lk_offline ) THEN
+        nk_rnf(:,:) = 1
+        h_rnf (:,:) = fsdept(:,:,1)
+      ENDIF
+
       ! dust input from the atmosphere
       ! ------------------------------
@@ -358,5 +376,4 @@
          rivalkinput = 0._wp
       END IF 
-
       ! nutrient input from dust
       ! ------------------------
@@ -410,5 +427,12 @@
          CALL iom_close( numiron )
          !
-         DO jk = 1, 5
+         ik50 = 5        !  last level where depth less than 50 m
+         DO jk = jpkm1, 1, -1
+            IF( gdept_1d(jk) > 50. )  ik50 = jk - 1
+         END DO
+         IF (lwp) WRITE(numout,*)
+         IF (lwp) WRITE(numout,*) ' Level corresponding to 50m depth ',  ik50,' ', gdept_1d(ik50+1)
+         IF (lwp) WRITE(numout,*)
+         DO jk = 1, ik50
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -421,30 +445,7 @@
             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
             DO jj = 1, jpj
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsed.F90	(revision 5462)
@@ -21,5 +21,4 @@
    USE p4zopt          !  optical model
    USE p4zlim          !  Co-limitations of differents nutrients
-   USE p4zrem          !  Remineralisation of organic matter
    USE p4zsbc          !  External source of nutrients 
    USE p4zint          !  interpolation and computation of various fields
@@ -30,11 +29,12 @@
    PRIVATE
 
-   PUBLIC   p4z_sed   
+   PUBLIC   p4z_sed  
+   PUBLIC   p4z_sed_alloc
+ 
 
    !! * Module variables
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: nitrpot    !: Nitrogen fixation 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:  ) :: sdenit     !: Nitrate reduction in the sediments
    REAL(wp) :: r1_rday                  !: inverse of rday
-
-   INTEGER ::  numnit  
-
 
    !!* Substitution
@@ -42,10 +42,10 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Header:$ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE p4z_sed( kt, jnt )
+   SUBROUTINE p4z_sed( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sed  ***
@@ -58,5 +58,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) ::   kt, jnt ! ocean time step
+      INTEGER, INTENT(in) ::   kt, knt ! ocean time step
       INTEGER  ::   ji, jj, jk, ikt
 #if ! defined key_sed
@@ -69,25 +69,20 @@
       REAL(wp) ::  zsiloss, zcaloss, zws3, zws4, zwsc, zdep, zwstpoc
       REAL(wp) ::  ztrfer, ztrpo4, zwdust, zlight
-      REAL(wp) ::  zrdenittot, zsdenittot, znitrpottot
       !
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3, zwork4
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zpdep, zsidep, zwork1, zwork2, zwork3
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zdenit2d, zironice, zbureff
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zwsbio3, zwsbio4, zwscal
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: znitrpot, zirondep, zsoufer
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zirondep, zsoufer
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_sed')
       !
-      IF( kt == nittrc000 .AND. jnt == 1 )  THEN
-         r1_rday  = 1. / rday
-         IF( ln_check_mass .AND. lwp)  &
-           &  CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
-      ENDIF
+      IF( kt == nittrc000 .AND. knt == 1 )   r1_rday  = 1. / rday
       !
       ! Allocate temporary workspace
-      CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zwork4, zbureff )
+      CALL wrk_alloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff )
       CALL wrk_alloc( jpi, jpj, zwsbio3, zwsbio4, zwscal )
-      CALL wrk_alloc( jpi, jpj, jpk, znitrpot, zsoufer )
+      CALL wrk_alloc( jpi, jpj, jpk, zsoufer )
 
       zdenit2d(:,:) = 0.e0
@@ -96,5 +91,4 @@
       zwork2  (:,:) = 0.e0
       zwork3  (:,:) = 0.e0
-      zwork4  (:,:) = 0.e0
 
       ! Iron input/uptake due to sea ice : Crude parameterization based on Lancelot et al.
@@ -108,5 +102,5 @@
                zdep    = rfact2 / fse3t(ji,jj,1)
                zwflux  = fmmflx(ji,jj) / 1000._wp
-               zfminus = MIN( 0._wp, -zwflux ) * trn(ji,jj,1,jpfer) * zdep
+               zfminus = MIN( 0._wp, -zwflux ) * trb(ji,jj,1,jpfer) * zdep
                zfplus  = MAX( 0._wp, -zwflux ) * icefeinput * zdep
                zironice(ji,jj) =  zfplus + zfminus
@@ -114,7 +108,7 @@
          END DO
          !
-         trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + zironice(:,:) 
+         tra(:,:,1,jpfer) = tra(:,:,1,jpfer) + zironice(:,:) 
          ! 
-         IF( lk_iomput .AND. jnt == nrdttrc .AND. iom_use( "Ironice" ) )   &
+         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironice" ) )   &
             &   CALL iom_put( "Ironice", zironice(:,:) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! iron flux from ice
          !
@@ -144,10 +138,10 @@
          END DO
          !                                              ! Iron solubilization of particles in the water column
-         trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + zpdep   (:,:)
-         trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + zsidep  (:,:)
-         trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + zirondep(:,:,:) 
+         tra(:,:,1,jppo4) = tra(:,:,1,jppo4) + zpdep   (:,:)
+         tra(:,:,1,jpsil) = tra(:,:,1,jpsil) + zsidep  (:,:)
+         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + zirondep(:,:,:) 
          ! 
          IF( lk_iomput ) THEN
-            IF( jnt == nrdttrc ) THEN
+            IF( knt == nrdttrc ) THEN
                 IF( iom_use( "Irondep" ) )   &
                 &  CALL iom_put( "Irondep", zirondep(:,:,1) * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1) ) ! surface downward dust depo of iron
@@ -167,10 +161,16 @@
       ! ----------------------------------------------------------
       IF( ln_river ) THEN
-         trn(:,:,1,jppo4) = trn(:,:,1,jppo4) + rivdip(:,:) * rfact2
-         trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + rivdin(:,:) * rfact2
-         trn(:,:,1,jpfer) = trn(:,:,1,jpfer) + rivdic(:,:) * 5.e-5 * rfact2
-         trn(:,:,1,jpsil) = trn(:,:,1,jpsil) + rivdsi(:,:) * rfact2
-         trn(:,:,1,jpdic) = trn(:,:,1,jpdic) + rivdic(:,:) * rfact2
-         trn(:,:,1,jptal) = trn(:,:,1,jptal) + ( rivalk(:,:) - rno3 * rivdin(:,:) ) * rfact2
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               DO jk = 1, nk_rnf(ji,jj)
+                  tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) +  rivdip(ji,jj) * rfact2
+                  tra(ji,jj,jk,jpno3) = tra(ji,jj,jk,jpno3) +  rivdin(ji,jj) * rfact2
+                  tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) +  rivdic(ji,jj) * 5.e-5 * rfact2
+                  tra(ji,jj,jk,jpsil) = tra(ji,jj,jk,jpsil) +  rivdsi(ji,jj) * rfact2
+                  tra(ji,jj,jk,jpdic) = tra(ji,jj,jk,jpdic) +  rivdic(ji,jj) * rfact2
+                  tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) +  ( rivalk(ji,jj) - rno3 * rivdin(ji,jj) ) * rfact2
+               ENDDO
+            ENDDO
+         ENDDO
       ENDIF
       
@@ -178,6 +178,6 @@
       ! ----------------------------------------------------------
       IF( ln_ndepo ) THEN
-         trn(:,:,1,jpno3) = trn(:,:,1,jpno3) + nitdep(:,:) * rfact2
-         trn(:,:,1,jptal) = trn(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2
+         tra(:,:,1,jpno3) = tra(:,:,1,jpno3) + nitdep(:,:) * rfact2
+         tra(:,:,1,jptal) = tra(:,:,1,jptal) - rno3 * nitdep(:,:) * rfact2
       ENDIF
 
@@ -185,7 +185,7 @@
       ! ------------------------------------------------------
       IF( ln_ironsed ) THEN
-         trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
+         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + ironsed(:,:,:) * rfact2
          !
-         IF( lk_iomput .AND. jnt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
+         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "Ironsed" ) )   &
             &   CALL iom_put( "Ironsed", ironsed(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! iron inputs from sediments
       ENDIF
@@ -194,7 +194,7 @@
       ! ------------------------------------------------------
       IF( ln_hydrofe ) THEN
-         trn(:,:,:,jpfer) = trn(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
+         tra(:,:,:,jpfer) = tra(:,:,:,jpfer) + hydrofe(:,:,:) * rfact2
          !
-         IF( lk_iomput .AND. jnt == nrdttrc .AND. iom_use( "HYDR" ) )   &
+         IF( lk_iomput .AND. knt == nrdttrc .AND. iom_use( "HYDR" ) )   &
             &   CALL iom_put( "HYDR", hydrofe(:,:,:) * 1.e+3 * tmask(:,:,:) ) ! hydrothermal iron input
       ENDIF
@@ -222,12 +222,12 @@
               ikt = mbkt(ji,jj)
 # if defined key_kriest
-              zflx =    trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)    * 1E3 * 1E6 / 1E4
+              zflx =    trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)    * 1E3 * 1E6 / 1E4
 # else
-              zflx = (  trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
-                &     + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4
+              zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
+                &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) )  * 1E3 * 1E6 / 1E4
 #endif
               zflx  = LOG10( MAX( 1E-3, zflx ) )
-              zo2   = LOG10( MAX( 10. , trn(ji,jj,ikt,jpoxy) * 1E6 ) )
-              zno3  = LOG10( MAX( 1.  , trn(ji,jj,ikt,jpno3) * 1E6 * rno3 ) )
+              zo2   = LOG10( MAX( 10. , trb(ji,jj,ikt,jpoxy) * 1E6 ) )
+              zno3  = LOG10( MAX( 1.  , trb(ji,jj,ikt,jpno3) * 1E6 * rno3 ) )
               zdep  = LOG10( fsdepw(ji,jj,ikt+1) )
               zdenit2d(ji,jj) = -2.2567 - 1.185 * zflx - 0.221 * zflx**2 - 0.3995 * zno3 * zo2 + 1.25 * zno3    &
@@ -235,6 +235,6 @@
               zdenit2d(ji,jj) = 10.0**( zdenit2d(ji,jj) )
               !
-              zflx = (  trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
-                &     + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6
+              zflx = (  trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj)   &
+                &     + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) ) * 1E6
               zbureff(ji,jj) = 0.013 + 0.53 * zflx**2 / ( 7.0 + zflx )**2
            ENDIF
@@ -251,14 +251,14 @@
                ikt = mbkt(ji,jj) 
 # if defined key_kriest
-               zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * zwscal (ji,jj)
-               zwork2(ji,jj) = trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)
+               zwork1(ji,jj) = trb(ji,jj,ikt,jpgsi) * zwscal (ji,jj)
+               zwork2(ji,jj) = trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj)
 # else
-               zwork1(ji,jj) = trn(ji,jj,ikt,jpgsi) * zwsbio4(ji,jj)
-               zwork2(ji,jj) = trn(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj) + trn(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) 
+               zwork1(ji,jj) = trb(ji,jj,ikt,jpgsi) * zwsbio4(ji,jj)
+               zwork2(ji,jj) = trb(ji,jj,ikt,jpgoc) * zwsbio4(ji,jj) + trb(ji,jj,ikt,jppoc) * zwsbio3(ji,jj) 
 # endif
                ! For calcite, burial efficiency is made a function of saturation
                zfactcal      = MIN( excess(ji,jj,ikt), 0.2 )
                zfactcal      = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
-               zwork3(ji,jj) = trn(ji,jj,ikt,jpcal) * zwscal(ji,jj) * 2.e0 * zfactcal
+               zwork3(ji,jj) = trb(ji,jj,ikt,jpcal) * zwscal(ji,jj) * 2.e0 * zfactcal
             ENDIF
          END DO
@@ -279,23 +279,23 @@
          DO ji = 1, jpi
             ikt  = mbkt(ji,jj)
-            zdep = xstep / fse3t(ji,jj,ikt)
+            zdep = xstep / fse3t(ji,jj,ikt) 
             zws4 = zwsbio4(ji,jj) * zdep
             zwsc = zwscal (ji,jj) * zdep
 # if defined key_kriest
-            zsiloss = trn(ji,jj,ikt,jpgsi) * zws4
+            zsiloss = trb(ji,jj,ikt,jpgsi) * zws4
 # else
-            zsiloss = trn(ji,jj,ikt,jpgsi) * zwsc
+            zsiloss = trb(ji,jj,ikt,jpgsi) * zwsc
 # endif
-            zcaloss = trn(ji,jj,ikt,jpcal) * zwsc
+            zcaloss = trb(ji,jj,ikt,jpcal) * zwsc
             !
-            trn(ji,jj,ikt,jpgsi) = trn(ji,jj,ikt,jpgsi) - zsiloss
-            trn(ji,jj,ikt,jpcal) = trn(ji,jj,ikt,jpcal) - zcaloss
+            tra(ji,jj,ikt,jpgsi) = tra(ji,jj,ikt,jpgsi) - zsiloss
+            tra(ji,jj,ikt,jpcal) = tra(ji,jj,ikt,jpcal) - zcaloss
 #if ! defined key_sed
-            trn(ji,jj,ikt,jpsil) = trn(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
+            tra(ji,jj,ikt,jpsil) = tra(ji,jj,ikt,jpsil) + zsiloss * zrivsil 
             zfactcal = MIN( excess(ji,jj,ikt), 0.2 )
             zfactcal = MIN( 1., 1.3 * ( 0.2 - zfactcal ) / ( 0.4 - zfactcal ) )
             zrivalk  =  1._wp - ( rivalkinput * r1_ryyss ) * zfactcal / ( zsumsedcal + rtrn )
-            trn(ji,jj,ikt,jptal) =  trn(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
-            trn(ji,jj,ikt,jpdic) =  trn(ji,jj,ikt,jpdic) + zcaloss * zrivalk
+            tra(ji,jj,ikt,jptal) =  tra(ji,jj,ikt,jptal) + zcaloss * zrivalk * 2.0
+            tra(ji,jj,ikt,jpdic) =  tra(ji,jj,ikt,jpdic) + zcaloss * zrivalk
 #endif
          END DO
@@ -304,20 +304,20 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            ikt     = mbkt(ji,jj)
-            zdep    = xstep / fse3t(ji,jj,ikt)
+            ikt  = mbkt(ji,jj)
+            zdep = xstep / fse3t(ji,jj,ikt) 
             zws4 = zwsbio4(ji,jj) * zdep
             zws3 = zwsbio3(ji,jj) * zdep
             zrivno3 = 1. - zbureff(ji,jj)
 # if ! defined key_kriest
-            trn(ji,jj,ikt,jpgoc) = trn(ji,jj,ikt,jpgoc) - trn(ji,jj,ikt,jpgoc) * zws4
-            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zws3
-            trn(ji,jj,ikt,jpbfe) = trn(ji,jj,ikt,jpbfe) - trn(ji,jj,ikt,jpbfe) * zws4
-            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zws3
-            zwstpoc              =  trn(ji,jj,ikt,jpgoc) * zws4 + trn(ji,jj,ikt,jppoc) * zws3 
+            tra(ji,jj,ikt,jpgoc) = tra(ji,jj,ikt,jpgoc) - trb(ji,jj,ikt,jpgoc) * zws4 
+            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3
+            tra(ji,jj,ikt,jpbfe) = tra(ji,jj,ikt,jpbfe) - trb(ji,jj,ikt,jpbfe) * zws4
+            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3
+            zwstpoc              = trb(ji,jj,ikt,jpgoc) * zws4 + trb(ji,jj,ikt,jppoc) * zws3
 # else
-            trn(ji,jj,ikt,jpnum) = trn(ji,jj,ikt,jpnum) - trn(ji,jj,ikt,jpnum) * zws4
-            trn(ji,jj,ikt,jppoc) = trn(ji,jj,ikt,jppoc) - trn(ji,jj,ikt,jppoc) * zws3
-            trn(ji,jj,ikt,jpsfe) = trn(ji,jj,ikt,jpsfe) - trn(ji,jj,ikt,jpsfe) * zws3
-            zwstpoc = trn(ji,jj,ikt,jppoc) * zws3 
+            tra(ji,jj,ikt,jpnum) = tra(ji,jj,ikt,jpnum) - trb(ji,jj,ikt,jpnum) * zws4 
+            tra(ji,jj,ikt,jppoc) = tra(ji,jj,ikt,jppoc) - trb(ji,jj,ikt,jppoc) * zws3
+            tra(ji,jj,ikt,jpsfe) = tra(ji,jj,ikt,jpsfe) - trb(ji,jj,ikt,jpsfe) * zws3
+            zwstpoc = trb(ji,jj,ikt,jppoc) * zws3 
 # endif
 
@@ -325,16 +325,16 @@
             ! The 0.5 factor in zpdenit and zdenitt is to avoid negative NO3 concentration after both denitrification
             ! in the sediments and just above the sediments. Not very clever, but simpliest option.
-            zpdenit  = MIN( 0.5 * ( trn(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
+            zpdenit  = MIN( 0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, zdenit2d(ji,jj) * zwstpoc * zrivno3 )
             z1pdenit = zwstpoc * zrivno3 - zpdenit
-            zolimit = MIN( ( trn(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
-            zdenitt = MIN(  0.5 * ( trn(ji,jj,ikt,jpno3) - rtrn ) / rdenit, z1pdenit * nitrfac(ji,jj,ikt) )
-            trn(ji,jj,ikt,jpdoc) = trn(ji,jj,ikt,jpdoc) + z1pdenit - zolimit - zdenitt
-            trn(ji,jj,ikt,jppo4) = trn(ji,jj,ikt,jppo4) + zpdenit + zolimit + zdenitt
-            trn(ji,jj,ikt,jpnh4) = trn(ji,jj,ikt,jpnh4) + zpdenit + zolimit + zdenitt
-            trn(ji,jj,ikt,jpno3) = trn(ji,jj,ikt,jpno3) - rdenit * (zpdenit + zdenitt)
-            trn(ji,jj,ikt,jpoxy) = trn(ji,jj,ikt,jpoxy) - zolimit * o2ut
-            trn(ji,jj,ikt,jptal) = trn(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) )
-            trn(ji,jj,ikt,jpdic) = trn(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt
-            zwork4(ji,jj) = rdenit * zpdenit * fse3t(ji,jj,ikt)
+            zolimit = MIN( ( trb(ji,jj,ikt,jpoxy) - rtrn ) / o2ut, z1pdenit * ( 1.- nitrfac(ji,jj,ikt) ) )
+            zdenitt = MIN(  0.5 * ( trb(ji,jj,ikt,jpno3) - rtrn ) / rdenit, z1pdenit * nitrfac(ji,jj,ikt) )
+            tra(ji,jj,ikt,jpdoc) = tra(ji,jj,ikt,jpdoc) + z1pdenit - zolimit - zdenitt
+            tra(ji,jj,ikt,jppo4) = tra(ji,jj,ikt,jppo4) + zpdenit + zolimit + zdenitt
+            tra(ji,jj,ikt,jpnh4) = tra(ji,jj,ikt,jpnh4) + zpdenit + zolimit + zdenitt
+            tra(ji,jj,ikt,jpno3) = tra(ji,jj,ikt,jpno3) - rdenit * (zpdenit + zdenitt)
+            tra(ji,jj,ikt,jpoxy) = tra(ji,jj,ikt,jpoxy) - zolimit * o2ut
+            tra(ji,jj,ikt,jptal) = tra(ji,jj,ikt,jptal) + rno3 * (zolimit + (1.+rdenit) * (zpdenit + zdenitt) )
+            tra(ji,jj,ikt,jpdic) = tra(ji,jj,ikt,jpdic) + zpdenit + zolimit + zdenitt
+            sdenit(ji,jj) = rdenit * zpdenit * fse3t(ji,jj,ikt)
 #endif
          END DO
@@ -356,7 +356,7 @@
 #endif
                ztrfer = biron(ji,jj,jk)       / ( concfediaz + biron(ji,jj,jk)       )
-               ztrpo4 = trn  (ji,jj,jk,jppo4) / ( concnnh4   + trn  (ji,jj,jk,jppo4) ) 
-               zlight =  ( 1.- EXP( -etot(ji,jj,jk) / diazolight ) ) 
-               znitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   &
+               ztrpo4 = trb  (ji,jj,jk,jppo4) / ( concnnh4   + trb  (ji,jj,jk,jppo4) ) 
+               zlight =  ( 1.- EXP( -etot_ndcy(ji,jj,jk) / diazolight ) ) 
+               nitrpot(ji,jj,jk) =  MAX( 0.e0, ( 0.6 * tgfunc(ji,jj,jk) - 2.15 ) * r1_rday )   &
                  &         *  zfact * MIN( ztrfer, ztrpo4 ) * zlight
                zsoufer(ji,jj,jk) = zlight * 2E-11 / (2E-11 + biron(ji,jj,jk))
@@ -370,39 +370,23 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zfact = znitrpot(ji,jj,jk) * nitrfix
-               trn(ji,jj,jk,jpnh4) = trn(ji,jj,jk,jpnh4) +             zfact
-               trn(ji,jj,jk,jptal) = trn(ji,jj,jk,jptal) + rno3      * zfact
-               trn(ji,jj,jk,jpoxy) = trn(ji,jj,jk,jpoxy) + o2nit     * zfact 
-               trn(ji,jj,jk,jppo4) = trn(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trn(ji,jj,jk,jppo4) ) &
-               &                     * 0.002 * trn(ji,jj,jk,jpdoc) * rfact2 / rday
-               trn(ji,jj,jk,jpfer) = trn(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * rfact2 / rday
+               zfact = nitrpot(ji,jj,jk) * nitrfix
+               tra(ji,jj,jk,jpnh4) = tra(ji,jj,jk,jpnh4) +             zfact
+               tra(ji,jj,jk,jptal) = tra(ji,jj,jk,jptal) + rno3      * zfact
+               tra(ji,jj,jk,jpoxy) = tra(ji,jj,jk,jpoxy) + o2nit     * zfact 
+               tra(ji,jj,jk,jppo4) = tra(ji,jj,jk,jppo4) + concdnh4 / ( concdnh4 + trb(ji,jj,jk,jppo4) ) &
+               &                     * 0.002 * trb(ji,jj,jk,jpdoc) * xstep
+               tra(ji,jj,jk,jpfer) = tra(ji,jj,jk,jpfer) + 0.002 * 4E-10 * zsoufer(ji,jj,jk) * xstep
            END DO
          END DO 
       END DO
 
-      ! Global budget of N SMS : denitrification in the water column and in the sediment
-      !                          nitrogen fixation by the diazotrophs
-      ! --------------------------------------------------------------------------------
-      zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
-      zsdenittot   = glob_sum ( zwork4(:,:)   * e1e2t(:,:) )
-      znitrpottot  = glob_sum ( znitrpot(:,:,:) * nitrfix              * cvol(:,:,:) )
-      zfact = 1.e+3 * rfact2r * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/kt ----> TgN/m3/year
-      !
-      IF( ln_check_mass .AND. ( kt == nitend .AND. jnt == nrdttrc ) .AND. ( lwp )  )  &
-         &  WRITE(numnit,9100) ndastp, znitrpottot * zfact  , &
-         &                             zrdenittot  * zfact  , &
-         &                             zsdenittot  * zfact
-      !
       IF( lk_iomput ) THEN
-         IF( jnt == nrdttrc ) THEN
+         IF( knt == nrdttrc ) THEN
             zfact = 1.e+3 * rfact2r * rno3  !  conversion from molC/l/kt  to molN/m3/s
-            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix"  , znitrpot(:,:,:) * nitrfix * zfact * tmask(:,:,:) )  ! nitrogen fixation 
-            IF( iom_use("Sdenit" ) ) CALL iom_put( "Sdenit", zwork4(:,:) * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments
-            IF( iom_use("tnfix"  ) ) CALL iom_put( "tnfix"  , znitrpottot * zfact  )               ! Global  nitrogen fixation
-            IF( iom_use("tdenit" ) ) CALL iom_put( "tdenit" , zrdenittot  * zfact  )               ! Total denitrification
+            IF( iom_use("Nfix"   ) ) CALL iom_put( "Nfix", nitrpot(:,:,:) * nitrfix * zfact * tmask(:,:,:) )  ! nitrogen fixation 
             IF( iom_use("INTNFIX") ) THEN   ! nitrogen fixation rate in ocean ( vertically integrated )
                zwork1(:,:) = 0.
                DO jk = 1, jpkm1
-                 zwork1(:,:) = zwork1(:,:) + znitrpot(:,:,jk) * nitrfix * zfact * fse3t(:,:,jk) * tmask(:,:,jk)
+                 zwork1(:,:) = zwork1(:,:) + nitrpot(:,:,jk) * nitrfix * zfact * fse3t(:,:,jk) * tmask(:,:,jk)
                ENDDO
                CALL iom_put( "INTNFIX" , zwork1 ) 
@@ -411,5 +395,5 @@
       ELSE
          IF( ln_diatrc )  &
-            &  trc2d(:,:,jp_pcs0_2d + 12) = znitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
+            &  trc2d(:,:,jp_pcs0_2d + 12) = nitrpot(:,:,1) * nitrfix * rno3 * 1.e+3 * rfact2r * fse3t(:,:,1) * tmask(:,:,1)
       ENDIF
       !
@@ -417,10 +401,10 @@
          WRITE(charout, fmt="('sed ')")
          CALL prt_ctl_trc_info(charout)
-         CALL prt_ctl_trc(tab4d=trn, mask=tmask, clinfo=ctrcnm)
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zwork4, zbureff )
+         CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm)
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, zdenit2d, zwork1, zwork2, zwork3, zbureff )
       CALL wrk_dealloc( jpi, jpj, zwsbio3, zwsbio4, zwscal )
-      CALL wrk_dealloc( jpi, jpj, jpk, znitrpot, zsoufer )
+      CALL wrk_dealloc( jpi, jpj, jpk, zsoufer )
       !
       IF( nn_timing == 1 )  CALL timing_stop('p4z_sed')
@@ -429,4 +413,16 @@
       !
    END SUBROUTINE p4z_sed
+
+
+   INTEGER FUNCTION p4z_sed_alloc()
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE p4z_sed_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( nitrpot(jpi,jpj,jpk), sdenit(jpi,jpj), STAT=p4z_sed_alloc )
+      !
+      IF( p4z_sed_alloc /= 0 )   CALL ctl_warn('p4z_sed_alloc: failed to allocate arrays')
+      !
+   END FUNCTION p4z_sed_alloc
+
 
 #else
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsink.F90	(revision 5462)
@@ -79,5 +79,5 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE p4z_sink ( kt, jnt )
+   SUBROUTINE p4z_sink ( kt, knt )
       !!---------------------------------------------------------------------
       !!                     ***  ROUTINE p4z_sink  ***
@@ -88,5 +88,5 @@
       !! ** Method  : - ???
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER, INTENT(in) :: kt, knt
       INTEGER  ::   ji, jj, jk, jit
       INTEGER  ::   iiter1, iiter2
@@ -199,15 +199,15 @@
                zfact = zstep * xdiss(ji,jj,jk)
                !  Part I : Coagulation dependent on turbulence
-               zagg1 = 25.9  * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jppoc)
-               zagg2 = 4452. * zfact * trn(ji,jj,jk,jppoc) * trn(ji,jj,jk,jpgoc)
+               zagg1 = 25.9  * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
+               zagg2 = 4452. * zfact * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
 
                ! Part II : Differential settling
 
                !  Aggregation of small into large particles
-               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)
+               zagg3 =  47.1 * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jpgoc)
+               zagg4 =  3.3  * zstep * trb(ji,jj,jk,jppoc) * trb(ji,jj,jk,jppoc)
 
                zagg   = zagg1 + zagg2 + zagg3 + zagg4
-               zaggfe = zagg * trn(ji,jj,jk,jpsfe) / ( trn(ji,jj,jk,jppoc) + rtrn )
+               zaggfe = zagg * trb(ji,jj,jk,jpsfe) / ( trb(ji,jj,jk,jppoc) + rtrn )
 
                ! Aggregation of DOC to POC : 
@@ -215,12 +215,12 @@
                ! 2nd term is shear aggregation of DOC-POC
                ! 3rd term is differential settling of DOC-POC
-               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.369 * 0.3 * trb(ji,jj,jk,jpdoc) + 102.4 * trb(ji,jj,jk,jppoc) ) * zfact       &
+               &            + 2.4 * zstep * trb(ji,jj,jk,jppoc) ) * 0.3 * trb(ji,jj,jk,jpdoc)
                ! transfer of DOC to GOC : 
                ! 1st term is shear aggregation
                ! 2nd term is differential settling 
-               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trn(ji,jj,jk,jpgoc) * 0.3 * trn(ji,jj,jk,jpdoc)
+               zaggdoc2 = ( 3.53E3 * zfact + 0.1 * zstep ) * trb(ji,jj,jk,jpgoc) * 0.3 * trb(ji,jj,jk,jpdoc)
                ! tranfer of DOC to POC due to brownian motion
-               zaggdoc3 =  ( 5095. * trn(ji,jj,jk,jppoc) + 114. * 0.3 * trn(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trn(ji,jj,jk,jpdoc)
+               zaggdoc3 =  ( 5095. * trb(ji,jj,jk,jppoc) + 114. * 0.3 * trb(ji,jj,jk,jpdoc) ) *zstep * 0.3 * trb(ji,jj,jk,jpdoc)
 
                !  Update the trends
@@ -237,9 +237,9 @@
 
      ! Total carbon export per year
-     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. jnt == nrdttrc )  )  &
+     IF( iom_use( "tcexp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  &
         &   t_oce_co2_exp = glob_sum( ( sinking(:,:,ik100) + sinking2(:,:,ik100) ) * e1e2t(:,:) * tmask(:,:,1) )
      !
      IF( lk_iomput ) THEN
-       IF( jnt == nrdttrc ) THEN
+       IF( knt == nrdttrc ) THEN
           CALL wrk_alloc( jpi, jpj,      zw2d )
           CALL wrk_alloc( jpi, jpj, jpk, zw3d )
@@ -328,5 +328,5 @@
    !!----------------------------------------------------------------------
 
-   SUBROUTINE p4z_sink ( kt, jnt )
+   SUBROUTINE p4z_sink ( kt, knt )
       !!---------------------------------------------------------------------
       !!                ***  ROUTINE p4z_sink  ***
@@ -338,5 +338,5 @@
       !!---------------------------------------------------------------------
       !
-      INTEGER, INTENT(in) :: kt, jnt
+      INTEGER, INTENT(in) :: kt, knt
       !
       INTEGER  :: ji, jj, jk, jit, niter1, niter2
@@ -373,5 +373,5 @@
             DO ji = 1, jpi
                IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
-                  znum = trn(ji,jj,jk,jppoc) / ( trn(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
+                  znum = trb(ji,jj,jk,jppoc) / ( trb(ji,jj,jk,jpnum) + rtrn ) / xkr_massp
                   ! -------------- To avoid sinking speed over 50 m/day -------
                   znum  = MIN( xnumm(jk), znum )
@@ -435,5 +435,5 @@
                IF( tmask(ji,jj,jk) /= 0.e0 ) THEN
 
-                  znum = trn(ji,jj,jk,jppoc)/(trn(ji,jj,jk,jpnum)+rtrn) / xkr_massp
+                  znum = trb(ji,jj,jk,jppoc)/(trb(ji,jj,jk,jpnum)+rtrn) / xkr_massp
                   !-------------- To avoid sinking speed over 50 m/day -------
                   znum  = min(xnumm(jk),znum)
@@ -453,10 +453,10 @@
                   !    ----------------------------------------------
 
-                  zagg1 =  0.163 * trn(ji,jj,jk,jpnum)**2               &
+                  zagg1 =  0.163 * trb(ji,jj,jk,jpnum)**2               &
                      &            * 2.*( (zfm-1.)*(zfm*xkr_mass_max**3-xkr_mass_min**3)    &
                      &            * (zeps-1)/zdiv1 + 3.*(zfm*xkr_mass_max-xkr_mass_min)    &
                      &            * (zfm*xkr_mass_max**2-xkr_mass_min**2)                  &
                      &            * (zeps-1.)**2/(zdiv2*zdiv3)) 
-                  zagg2 =  2*0.163*trn(ji,jj,jk,jpnum)**2*zfm*                       &
+                  zagg2 =  2*0.163*trb(ji,jj,jk,jpnum)**2*zfm*                       &
                      &                   ((xkr_mass_max**3+3.*(xkr_mass_max**2          &
                      &                    *xkr_mass_min*(zeps-1.)/zdiv2                 &
@@ -466,5 +466,5 @@
                      &                    (zeps-2.)+(zeps-1.)/zdiv3)+(zeps-1.)/zdiv1))    
 
-                  zagg3 =  0.163*trn(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3  
+                  zagg3 =  0.163*trb(ji,jj,jk,jpnum)**2*zfm**2*8. * xkr_mass_max**3  
                   
                  !    Aggregation of small into large particles
@@ -472,5 +472,5 @@
                  !    ----------------------------------------------
 
-                  zagg4 =  2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2*                       &
+                  zagg4 =  2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2*                       &
                      &                 xkr_wsbio_min*(zeps-1.)**2                         &
                      &                 *(xkr_mass_min**2*((1.-zsm*zfm)/(zdiv3*zdiv4)      &
@@ -479,5 +479,5 @@
                      &                 *xkr_eta)/(zdiv*zdiv3*zdiv5) )   
 
-                  zagg5 =   2.*3.141*0.125*trn(ji,jj,jk,jpnum)**2                         &
+                  zagg5 =   2.*3.141*0.125*trb(ji,jj,jk,jpnum)**2                         &
                      &                 *(zeps-1.)*zfm*xkr_wsbio_min                        &
                      &                 *(zsm*(xkr_mass_min**2-zfm*xkr_mass_max**2)         &
@@ -489,5 +489,5 @@
                   !     ------------------------------------
 
-                  zfract = 2.*3.141*0.125*trn(ji,jj,jk,jpmes)*12./0.12/0.06**3*trn(ji,jj,jk,jpnum)  &
+                  zfract = 2.*3.141*0.125*trb(ji,jj,jk,jpmes)*12./0.12/0.06**3*trb(ji,jj,jk,jpnum)  &
                     &      * (0.01/xkr_mass_min)**(1.-zeps)*0.1**2  &
                     &      * 10000.*xstep
@@ -496,8 +496,8 @@
                   !     --------------------------------------
 
-                  zaggdoc = 0.83 * trn(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)   &
-                     &        + 0.005 * 231. * trn(ji,jj,jk,jpdoc) * xstep * trn(ji,jj,jk,jpdoc)
-                  zaggdoc1 = 271. * trn(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trn(ji,jj,jk,jpdoc)  &
-                     &  + 0.02 * 16706. * trn(ji,jj,jk,jppoc) * xstep * trn(ji,jj,jk,jpdoc)
+                  zaggdoc = 0.83 * trb(ji,jj,jk,jpdoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)   &
+                     &        + 0.005 * 231. * trb(ji,jj,jk,jpdoc) * xstep * trb(ji,jj,jk,jpdoc)
+                  zaggdoc1 = 271. * trb(ji,jj,jk,jppoc) * xstep * xdiss(ji,jj,jk) * trb(ji,jj,jk,jpdoc)  &
+                     &  + 0.02 * 16706. * trb(ji,jj,jk,jppoc) * xstep * trb(ji,jj,jk,jpdoc)
 
 # if defined key_degrad
@@ -514,5 +514,5 @@
                   zagg = 0.5 * xkr_stick * ( zaggsh + zaggsi )
                   !
-                  znumdoc = trn(ji,jj,jk,jpnum) / ( trn(ji,jj,jk,jppoc) + rtrn )
+                  znumdoc = trb(ji,jj,jk,jpnum) / ( trb(ji,jj,jk,jppoc) + rtrn )
                   tra(ji,jj,jk,jppoc) = tra(ji,jj,jk,jppoc) + zaggdoc + zaggdoc1
                   tra(ji,jj,jk,jpnum) = tra(ji,jj,jk,jpnum) + zfract + zaggdoc / xkr_massp - zagg
@@ -528,5 +528,5 @@
      !
      IF( lk_iomput ) THEN
-        IF( jnt == nrdttrc ) THEN
+        IF( knt == nrdttrc ) THEN
           CALL wrk_alloc( jpi, jpj,      zw2d )
           CALL wrk_alloc( jpi, jpj, jpk, zw3d )
@@ -800,5 +800,5 @@
       ztraz(:,:,:) = 0.e0
       zakz (:,:,:) = 0.e0
-      ztrb (:,:,:) = trn(:,:,:,jp_tra)
+      ztrb (:,:,:) = trb(:,:,:,jp_tra)
 
       DO jk = 1, jpkm1
@@ -815,5 +815,5 @@
          !  first guess of the slopes interior values
          DO jk = 2, jpkm1
-            ztraz(:,:,jk) = ( trn(:,:,jk-1,jp_tra) - trn(:,:,jk,jp_tra) ) * tmask(:,:,jk)
+            ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)
          END DO
          ztraz(:,:,1  ) = 0.0
@@ -846,5 +846,5 @@
                   zigma = zwsink2(ji,jj,jk+1) * zstep / fse3w(ji,jj,jk+1)
                   zew   = zwsink2(ji,jj,jk+1)
-                  psinkflx(ji,jj,jk+1) = -zew * ( trn(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
+                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
                END DO
             END DO
@@ -859,5 +859,5 @@
                DO ji = 1, jpi
                   zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / fse3t(ji,jj,jk)
-                  trn(ji,jj,jk,jp_tra) = trn(ji,jj,jk,jp_tra) + zflx
+                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
                END DO
             END DO
@@ -875,5 +875,5 @@
       END DO
 
-      trn(:,:,:,jp_tra) = ztrb(:,:,:)
+      trb(:,:,:,jp_tra) = ztrb(:,:,:)
       psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
       !
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zsms.F90	(revision 5462)
@@ -24,4 +24,5 @@
    USE p4zsed          !  Sedimentation
    USE p4zint          !  time interpolation
+   USE p4zrem          !  remineralisation
    USE iom             !  I/O manager
    USE trd_oce         !  Ocean trends variables
@@ -36,7 +37,14 @@
    PUBLIC   p4z_sms        ! called in p4zsms.F90
 
-   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget
-   INTEGER ::  numco2, numnut  !: logical unit for co2 budget
-
+   REAL(wp) :: alkbudget, no3budget, silbudget, ferbudget, po4budget
+   REAL(wp) :: xfact1, xfact2
+   INTEGER ::  numco2, numnut, numnit  !: logical unit for co2 budget
+
+   !!* Array used to indicate negative tracer values
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     !: ???
+
+
+   !! * Substitutions
+#  include "top_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
@@ -61,20 +69,17 @@
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index      
       !!
-      INTEGER ::   jnt, jn, jl
+      INTEGER ::   ji, jj, jk, jnt, jn, jl
+      REAL(wp) ::  ztra
+#if defined key_kriest
+      REAL(wp) ::  zcoef1, zcoef2
+#endif
       CHARACTER (len=25) :: charout
-      REAL(wp), POINTER, DIMENSION(:,:,:,:)  :: ztrdpis
       !!---------------------------------------------------------------------
       !
       IF( nn_timing == 1 )  CALL timing_start('p4z_sms')
       !
-      IF( l_trdtrc )  THEN
-         CALL wrk_alloc( jpi, jpj, jpk, jp_pisces, ztrdpis ) 
-         DO jn = 1, jp_pisces
-            jl = jn + jp_pcs0 - 1
-            ztrdpis(:,:,:,jn) = trn(:,:,:,jl)
-         ENDDO
-      ENDIF
-      !
       IF( kt == nittrc000 ) THEN
+        !
+        ALLOCATE( xnegtr(jpi,jpj,jpk) )
         !
         CALL p4z_che                              ! initialize the chemical constants
@@ -88,4 +93,26 @@
       IF( ln_pisdmp .AND. MOD( kt - nn_dttrc, nn_pisdmp ) == 0 )   CALL p4z_dmp( kt )      ! Relaxation of some tracers
       !
+      !                                                                    !   set time step size (Euler/Leapfrog)
+      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN   ;    rfact = rdttrc(1)     !  at nittrc000
+      ELSEIF( kt <= nittrc000 + nn_dttrc )                          THEN   ;    rfact = 2. * rdttrc(1)   ! at nittrc000 or nittrc000+nn_dttrc (Leapfrog)
+      ENDIF
+      !
+      IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + nn_dttrc ) ) THEN
+         rfactr  = 1. / rfact
+         rfact2  = rfact / FLOAT( nrdttrc )
+         rfact2r = 1. / rfact2
+         xstep = rfact2 / rday         ! Time step duration for biology
+         IF(lwp) WRITE(numout,*) 
+         IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdttra(1)
+         IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2
+         IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      IF( ( neuler == 0 .AND. kt == nittrc000 ) .OR. ln_top_euler ) THEN
+         DO jn = jp_pcs0, jp_pcs1              !   SMS on tracer without Asselin time-filter
+            trb(:,:,:,jn) = trn(:,:,:,jn)
+         END DO
+      ENDIF
+      !
       IF( ndayflxtr /= nday_year ) THEN      ! New days
          !
@@ -105,28 +132,57 @@
       DO jnt = 1, nrdttrc          ! Potential time splitting if requested
          !
-         CALL p4z_bio (kt, jnt)    ! Biology
-         CALL p4z_sed (kt, jnt)    ! Sedimentation
-         !
+         CALL p4z_bio( kt, jnt )   ! Biology
+         CALL p4z_sed( kt, jnt )   ! Sedimentation
+         CALL p4z_lys( kt, jnt )   ! Compute CaCO3 saturation
+         CALL p4z_flx( kt, jnt )   ! Compute surface fluxes
+         !
+         xnegtr(:,:,:) = 1.e0
          DO jn = jp_pcs0, jp_pcs1
-            trb(:,:,:,jn) = trn(:,:,:,jn)
-         ENDDO
-         !
+            DO jk = 1, jpk
+               DO jj = 1, jpj
+                  DO ji = 1, jpi
+                     IF( ( trb(ji,jj,jk,jn) + tra(ji,jj,jk,jn) ) < 0.e0 ) THEN
+                        ztra             = ABS( trb(ji,jj,jk,jn) ) / ( ABS( tra(ji,jj,jk,jn) ) + rtrn )
+                        xnegtr(ji,jj,jk) = MIN( xnegtr(ji,jj,jk),  ztra )
+                     ENDIF
+                 END DO
+               END DO
+            END DO
+         END DO
+         !                                ! where at least 1 tracer concentration becomes negative
+         !                                ! 
+         DO jn = jp_pcs0, jp_pcs1
+           trb(:,:,:,jn) = trb(:,:,:,jn) + xnegtr(:,:,:) * tra(:,:,:,jn)
+         END DO
+        !
+         DO jn = jp_pcs0, jp_pcs1
+            tra(:,:,:,jn) = 0._wp
+         END DO
+         !
+         IF( ln_top_euler ) THEN
+            DO jn = jp_pcs0, jp_pcs1
+               trn(:,:,:,jn) = trb(:,:,:,jn)
+            END DO
+         ENDIF
       END DO
 
-      IF( l_trdtrc )  THEN
-         DO jn = 1, jp_pisces
-            jl = jn + jp_pcs0 - 1
-            ztrdpis(:,:,:,jn) = ( ztrdpis(:,:,:,jn) - trn(:,:,:,jl) ) * rfact2r
-         ENDDO
-      ENDIF
-      CALL p4z_lys( kt )             ! Compute CaCO3 saturation
-      CALL p4z_flx( kt )             ! Compute surface fluxes
-
-      DO jn = jp_pcs0, jp_pcs1
-        CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
-        CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
-        CALL lbc_lnk( tra(:,:,:,jn), 'T', 1. )
+#if defined key_kriest
+      ! 
+      zcoef1 = 1.e0 / xkr_massp 
+      zcoef2 = 1.e0 / xkr_massp / 1.1
+      DO jk = 1,jpkm1
+         trb(:,:,jk,jpnum) = MAX(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef1 / xnumm(jk)  )
+         trb(:,:,jk,jpnum) = MIN(  trb(:,:,jk,jpnum), trb(:,:,jk,jppoc) * zcoef2              )
       END DO
       !
+#endif
+      !
+      !
+      IF( l_trdtrc ) THEN
+         DO jn = jp_pcs0, jp_pcs1
+           CALL trd_trc( tra(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
+         END DO
+      END IF
+      !
       IF( lk_sed ) THEN 
          !
@@ -134,5 +190,5 @@
          !
          DO jn = jp_pcs0, jp_pcs1
-           CALL lbc_lnk( trn(:,:,:,jn), 'T', 1. )
+           CALL lbc_lnk( trb(:,:,:,jn), 'T', 1. )
          END DO
          !
@@ -141,13 +197,5 @@
       IF( lrst_trc )  CALL p4z_rst( kt, 'WRITE' )  !* Write PISCES informations in restart file 
       !
-      IF( l_trdtrc ) THEN
-         DO jn = 1, jp_pisces
-            jl = jn + jp_pcs0 - 1
-             ztrdpis(:,:,:,jn) = ztrdpis(:,:,:,jn) + tra(:,:,:,jl)
-             CALL trd_trc( ztrdpis(:,:,:,jn), jn, jptra_sms, kt )   ! save trends
-          END DO
-          CALL wrk_dealloc( jpi, jpj, jpk, jp_pisces, ztrdpis ) 
-      END IF
-      !
+
       IF( lk_iomput .OR. ln_check_mass )  CALL p4z_chk_mass( kt ) ! Mass conservation checking
 
@@ -280,7 +328,7 @@
                ztmas   = tmask(ji,jj,jk)
                ztmas1  = 1. - tmask(ji,jj,jk)
-               zcaralk = trn(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
-               zco3    = ( zcaralk - trn(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
-               zbicarb = ( 2. * trn(ji,jj,jk,jpdic) - zcaralk )
+               zcaralk = trb(ji,jj,jk,jptal) - borat(ji,jj,jk) / (  1. + 1.E-8 / ( rtrn + akb3(ji,jj,jk) )  )
+               zco3    = ( zcaralk - trb(ji,jj,jk,jpdic) ) * ztmas + 0.5e-3 * ztmas1
+               zbicarb = ( 2. * trb(ji,jj,jk,jpdic) - zcaralk )
                hi(ji,jj,jk) = ( ak23(ji,jj,jk) * zbicarb / zco3 ) * ztmas + 1.e-9 * ztmas1
             END DO
@@ -361,5 +409,6 @@
       REAL(wp) ::  silmean = 91.51     ! mean value of silicate
       !
-      REAL(wp) :: zarea, zalksum, zpo4sum, zno3sum, zsilsum
+      REAL(wp) :: zarea, zalksumn, zpo4sumn, zno3sumn, zsilsumn
+      REAL(wp) :: zalksumb, zpo4sumb, zno3sumb, zsilsumb
       !!---------------------------------------------------------------------
 
@@ -374,23 +423,44 @@
          zarea          = 1._wp / glob_sum( cvol(:,:,:) ) * 1e6              
 
-         zalksum = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
-         zpo4sum = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
-         zno3sum = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
-         zsilsum = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
+         zalksumn = glob_sum( trn(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
+         zpo4sumn = glob_sum( trn(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
+         zno3sumn = glob_sum( trn(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
+         zsilsumn = glob_sum( trn(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
  
-         IF(lwp) WRITE(numout,*) '       TALK mean : ', zalksum
-         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksum
-
-         IF(lwp) WRITE(numout,*) '       PO4  mean : ', zpo4sum
-         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sum
-
-         IF(lwp) WRITE(numout,*) '       NO3  mean : ', zno3sum
-         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sum
-
-         IF(lwp) WRITE(numout,*) '       SiO3 mean : ', zsilsum
-         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsum )
-         !
-      ENDIF
-
+         IF(lwp) WRITE(numout,*) '       TALKN mean : ', zalksumn
+         trn(:,:,:,jptal) = trn(:,:,:,jptal) * alkmean / zalksumn
+
+         IF(lwp) WRITE(numout,*) '       PO4N  mean : ', zpo4sumn
+         trn(:,:,:,jppo4) = trn(:,:,:,jppo4) * po4mean / zpo4sumn
+
+         IF(lwp) WRITE(numout,*) '       NO3N  mean : ', zno3sumn
+         trn(:,:,:,jpno3) = trn(:,:,:,jpno3) * no3mean / zno3sumn
+
+         IF(lwp) WRITE(numout,*) '       SiO3N mean : ', zsilsumn
+         trn(:,:,:,jpsil) = MIN( 400.e-6,trn(:,:,:,jpsil) * silmean / zsilsumn )
+         !
+         !
+         IF( .NOT. ln_top_euler ) THEN
+            zalksumb = glob_sum( trb(:,:,:,jptal) * cvol(:,:,:)  ) * zarea
+            zpo4sumb = glob_sum( trb(:,:,:,jppo4) * cvol(:,:,:)  ) * zarea * po4r
+            zno3sumb = glob_sum( trb(:,:,:,jpno3) * cvol(:,:,:)  ) * zarea * rno3
+            zsilsumb = glob_sum( trb(:,:,:,jpsil) * cvol(:,:,:)  ) * zarea
+ 
+            IF(lwp) WRITE(numout,*) ' '
+            IF(lwp) WRITE(numout,*) '       TALKB mean : ', zalksumb
+            trb(:,:,:,jptal) = trb(:,:,:,jptal) * alkmean / zalksumb
+
+            IF(lwp) WRITE(numout,*) '       PO4B  mean : ', zpo4sumb
+            trb(:,:,:,jppo4) = trb(:,:,:,jppo4) * po4mean / zpo4sumb
+
+            IF(lwp) WRITE(numout,*) '       NO3B  mean : ', zno3sumb
+            trb(:,:,:,jpno3) = trb(:,:,:,jpno3) * no3mean / zno3sumb
+
+            IF(lwp) WRITE(numout,*) '       SiO3B mean : ', zsilsumb
+            trb(:,:,:,jpsil) = MIN( 400.e-6,trb(:,:,:,jpsil) * silmean / zsilsumb )
+        ENDIF
+        !
+      ENDIF
+        !
    END SUBROUTINE p4z_dmp
 
@@ -406,5 +476,11 @@
       INTEGER , INTENT( in ) ::   kt      ! ocean time-step index      
       REAL(wp)               ::  zfact       
-      !!
+      REAL(wp) ::  zrdenittot, zsdenittot, znitrpottot
+      CHARACTER(LEN=100)   ::   cltxt
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zvol
+      INTEGER :: jk
+      !!----------------------------------------------------------------------
+
+      !
       !!---------------------------------------------------------------------
 
@@ -413,7 +489,14 @@
             CALL ctl_opn( numco2, 'carbon.budget'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
             CALL ctl_opn( numnut, 'nutrient.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
+            CALL ctl_opn( numnit, 'nitrogen.budget', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, 6, .FALSE., narea )
+            xfact1 = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/yr
+            xfact2 = 1.e+3 * rno3 * 14. / 1.e12 * ryyss   ! conversion molC/l/s ----> TgN/m3/yr
+            cltxt='time-step   Alkalinity        Nitrate        Phosphorus         Silicate           Iron'
+            IF( lwp ) WRITE(numnut,*)  TRIM(cltxt)
+            IF( lwp ) WRITE(numnut,*) 
          ENDIF
       ENDIF
 
+      !
       IF( iom_use( "pno3tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          !   Compute the budget of NO3, ALK, Si, Fer
@@ -431,5 +514,18 @@
       ENDIF
       !
-      IF( iom_use( "psiltot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+      IF( iom_use( "ppo4tot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
+         po4budget = glob_sum( (   trn(:,:,:,jppo4)                     &
+            &                    + trn(:,:,:,jpphy) + trn(:,:,:,jpdia)  &
+            &                    + trn(:,:,:,jpzoo) + trn(:,:,:,jpmes)  &
+            &                    + trn(:,:,:,jppoc)                     &
+#if ! defined key_kriest
+            &                    + trn(:,:,:,jpgoc)                     &
+#endif
+            &                    + trn(:,:,:,jpdoc)                     ) * cvol(:,:,:)  )
+         po4budget = po4budget / areatot
+         CALL iom_put( "ppo4tot", po4budget )
+      ENDIF
+      !
+      IF( iom_use( "psiltot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          silbudget = glob_sum( (   trn(:,:,:,jpsil) + trn(:,:,:,jpgsi)  &
             &                    + trn(:,:,:,jpdsi)                     ) * cvol(:,:,:)  )
@@ -439,5 +535,5 @@
       ENDIF
       !
-      IF( iom_use( "palktot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+      IF( iom_use( "palktot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          alkbudget = glob_sum( (   trn(:,:,:,jpno3) * rno3              &
             &                    + trn(:,:,:,jptal)                     &
@@ -448,5 +544,5 @@
       ENDIF
       !
-      IF( iom_use( "pfertot" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+      IF( iom_use( "pfertot" ) .OR. ( ln_check_mass .AND. kt == nitend )  ) THEN
          ferbudget = glob_sum( (   trn(:,:,:,jpfer) + trn(:,:,:,jpnfe)  &
             &                    + trn(:,:,:,jpdfe)                     &
@@ -462,19 +558,45 @@
       ENDIF
       !
+
+      ! Global budget of N SMS : denitrification in the water column and in the sediment
+      !                          nitrogen fixation by the diazotrophs
+      ! --------------------------------------------------------------------------------
+      IF( iom_use( "tnfix" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+         znitrpottot  = glob_sum ( nitrpot(:,:,:) * nitrfix * cvol(:,:,:) )
+         CALL iom_put( "tnfix"  , znitrpottot * 1.e+3 * rno3 )  ! Global  nitrogen fixation molC/l  to molN/m3 
+      ENDIF
+      !
+      IF( iom_use( "tdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+         zrdenittot   = glob_sum ( denitr(:,:,:) * rdenit * xnegtr(:,:,:) * cvol(:,:,:) )
+         CALL iom_put( "tdenit"  , zrdenittot * 1.e+3 * rno3 )  ! Total denitrification molC/l to molN/m3 
+      ENDIF
+      !
+      IF( iom_use( "Sdenit" ) .OR.  ( ln_check_mass .AND. kt == nitend )  ) THEN
+         zsdenittot   = glob_sum ( sdenit(:,:) * e1e2t(:,:) )
+         CALL iom_put( "Sdenit", sdenit(:,:) * zfact * tmask(:,:,1) )  ! Nitrate reduction in the sediments
+      ENDIF
+
       IF( ln_check_mass .AND. kt == nitend ) THEN   ! Compute the budget of NO3, ALK, Si, Fer
-         zfact = rfact2r * 12. / 1.e15 * ryyss    ! conversion molC/kt --> PgC/year
          t_atm_co2_flx  = t_atm_co2_flx / glob_sum( e1e2t(:,:) )
-         t_oce_co2_flx  = t_oce_co2_flx         * zfact * (-1 )
-         tpp            = tpp           * 1000. * zfact
-         t_oce_co2_exp  = t_oce_co2_exp * 1000. * zfact
+         t_oce_co2_flx  = t_oce_co2_flx         * xfact1 * (-1 )
+         tpp            = tpp           * 1000. * xfact1
+         t_oce_co2_exp  = t_oce_co2_exp * 1000. * xfact1
          IF( lwp ) WRITE(numco2,9000) ndastp, t_atm_co2_flx, t_oce_co2_flx, tpp, t_oce_co2_exp
-         IF( lwp ) WRITE(numnut,9500) ndastp, alkbudget        * 1.e+06, &
+         IF( lwp ) WRITE(numnut,9100) ndastp, alkbudget        * 1.e+06, &
              &                                no3budget * rno3 * 1.e+06, &
+             &                                po4budget * po4r * 1.e+06, &
              &                                silbudget        * 1.e+06, &
              &                                ferbudget        * 1.e+09
+         !
+         IF( lwp ) WRITE(numnit,9200) ndastp, znitrpottot * xfact2  , &
+         &                             zrdenittot  * xfact2  , &
+         &                             zsdenittot  * xfact2
+
       ENDIF
       !
  9000  FORMAT(i8,f10.5,e18.10,f10.5,f10.5)
- 9500  FORMAT(i8,4e18.10)
+ 9100  FORMAT(i8,5e18.10)
+ 9200  FORMAT(i8,3f10.5)
+
        !
    END SUBROUTINE p4z_chk_mass
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/par_sed.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/par_sed.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/par_sed.F90	(revision 5462)
@@ -7,4 +7,5 @@
    !!        !  06-12  (C. Ethe)  Orignal
    !!----------------------------------------------------------------------
+   !! $Id$
 #if defined key_sed
    !! Domain characteristics
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sed.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sed.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sed.F90	(revision 5462)
@@ -160,4 +160,5 @@
    INTEGER, PUBLIC ::  numsed = 27    ! units
 
+   !! $Id$
 CONTAINS
 
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedadv.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedadv.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedadv.F90	(revision 5462)
@@ -23,4 +23,5 @@
    REAL(wp) :: eps = 1.e-13
 
+   !! $Id$
 CONTAINS
 
@@ -438,4 +439,5 @@
    !! MODULE sedbtb  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_adv( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedarr.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedarr.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedarr.F90	(revision 5462)
@@ -29,5 +29,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/TOP 3.3 , NEMO Consortium (2010)
-   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limtab.F90,v 1.2 2005/03/27 18:34:42 opalod Exp $ 
+   !! $Id$
    !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedbtb.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedbtb.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedbtb.F90	(revision 5462)
@@ -12,4 +12,5 @@
 
 
+   !! $Id$
 CONTAINS
    
@@ -77,4 +78,5 @@
    !! MODULE sedbtb  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_btb( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedchem.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedchem.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedchem.F90	(revision 5462)
@@ -163,4 +163,5 @@
    DATA Ddsw / 999.842594 , 6.793952E-2 , -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9/
 
+   !! $Id$
 CONTAINS
 
@@ -559,4 +560,5 @@
    !! MODULE sedchem  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_chem( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedco3.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedco3.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedco3.F90	(revision 5462)
@@ -23,4 +23,5 @@
    !!----------------------------------------------------------------------
 
+   !! $Id$
 CONTAINS
 
@@ -188,4 +189,5 @@
    !! MODULE sedco3  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_co3( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/seddsr.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/seddsr.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/seddsr.F90	(revision 5462)
@@ -20,4 +20,5 @@
    REAL(wp), DIMENSION(:), ALLOCATABLE, PUBLIC ::  dens_mol_wgt  ! molecular density 
 
+   !! $Id$
 CONTAINS
    
@@ -530,4 +531,5 @@
    !! MODULE seddsr  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_dsr ( kt )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/seddta.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/seddta.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/seddta.F90	(revision 5462)
@@ -28,4 +28,5 @@
 #endif
 
+   !! $Id$
 CONTAINS
 
@@ -268,4 +269,5 @@
    !! MODULE seddta  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_dta ( kt )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedini.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedini.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedini.F90	(revision 5462)
@@ -55,4 +55,5 @@
    PUBLIC sed_init          ! routine called by opa.F90
 
+   !! $Id$
 CONTAINS
 
@@ -856,4 +857,5 @@
    !!   Dummy module :                      NO Sediment model
    !!----------------------------------------------------------------------
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_ini              ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmat.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmat.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmat.F90	(revision 5462)
@@ -22,4 +22,5 @@
 
 
+   !! $Id$
  CONTAINS
 
@@ -257,4 +258,5 @@
    !! MODULE sedmat  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_mat         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmbc.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmbc.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmbc.F90	(revision 5462)
@@ -36,4 +36,5 @@
    REAL(wp)  :: src13ca  
 
+   !! $Id$
 CONTAINS
 
@@ -311,4 +312,5 @@
    !! MODULE sedmbc :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_mbc( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmodel.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmodel.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedmodel.F90	(revision 5462)
@@ -17,4 +17,5 @@
    LOGICAL, PUBLIC, PARAMETER ::   lk_sed = .TRUE.     !: sediment flag
 
+   !! $Id$
 CONTAINS
 
@@ -47,4 +48,5 @@
    !!======================================================================
    LOGICAL, PUBLIC, PARAMETER ::   lk_sed = .FALSE.     !: sediment flag
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_model( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedrst.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedrst.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedrst.F90	(revision 5462)
@@ -25,4 +25,5 @@
    
    
+   !! $Id$
 CONTAINS
 
@@ -270,4 +271,5 @@
    !! MODULE sedrst :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_rst_read                      ! Empty routines
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedsfc.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedsfc.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedsfc.F90	(revision 5462)
@@ -12,4 +12,5 @@
    PUBLIC sed_sfc
 
+   !! $Id$
 CONTAINS
 
@@ -67,4 +68,5 @@
    !! MODULE sedsfc  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_sfc ( kt )
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedstp.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedstp.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedstp.F90	(revision 5462)
@@ -23,4 +23,5 @@
    PUBLIC sed_stp  ! called by step.F90
 
+   !! $Id$
 CONTAINS
 
@@ -69,4 +70,5 @@
    !! MODULE sedstp  :   Dummy module 
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_stp( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedwri.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedwri.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/SED/sedwri.F90	(revision 5462)
@@ -25,4 +25,5 @@
    INTEGER, ALLOCATABLE, SAVE, DIMENSION(:) :: ndext51
 
+   !! $Id$
 CONTAINS
 
@@ -264,4 +265,5 @@
    !! MODULE sedwri  :   Dummy module
    !!======================================================================
+   !! $Id$
 CONTAINS
    SUBROUTINE sed_wri( kt )         ! Empty routine
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/par_pisces.F90	(revision 5462)
@@ -63,10 +63,10 @@
    INTEGER, PUBLIC, PARAMETER ::   jpdia = 11    !: Diatoms Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpmes = 12    !: Mesozooplankton Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: (big) Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: Diatoms Silicate Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpfer = 14    !: Iron Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpnum = 15    !: Big iron particles Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpsfe = 16    !: number of particulate organic phosphate concentration
    INTEGER, PUBLIC, PARAMETER ::   jpdfe = 17    !: Diatoms iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 18    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 18    !: (big) Silicate Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpnfe = 19    !: Nano iron Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpnch = 20    !: Nano Chlorophyll Concentration
@@ -102,5 +102,5 @@
    INTEGER, PUBLIC, PARAMETER ::   jpdia = 11    !: Diatoms Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpmes = 12    !: Mesozooplankton Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: (big) Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpdsi = 13    !: Diatoms Silicate Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpfer = 14    !: Iron Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpbfe = 15    !: Big iron particles Concentration
@@ -108,5 +108,5 @@
    INTEGER, PUBLIC, PARAMETER ::   jpsfe = 17    !: Small iron particles Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpdfe = 18    !: Diatoms iron Concentration
-   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 19    !: Diatoms Silicate Concentration
+   INTEGER, PUBLIC, PARAMETER ::   jpgsi = 19    !: (big) Silicate Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpnfe = 20    !: Nano iron Concentration
    INTEGER, PUBLIC, PARAMETER ::   jpnch = 21    !: Nano Chlorophyll Concentration
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/sms_pisces.F90	(revision 5462)
@@ -106,7 +106,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tgfunc2   !: Temp. dependancy of mesozooplankton rates
 
-   !!* Array used to indicate negative tracer values
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   xnegtr     !: ???
-
 #if defined key_kriest
    !!*  Kriest parameter for aggregation
@@ -131,5 +128,5 @@
       !!----------------------------------------------------------------------
       USE lib_mpp , ONLY: ctl_warn
-      INTEGER ::   ierr(6)        ! Local variables
+      INTEGER ::   ierr(5)        ! Local variables
       !!----------------------------------------------------------------------
       ierr(:) = 0
@@ -162,6 +159,4 @@
       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(6) )
 #endif
       !
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/trcice_pisces.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/trcice_pisces.F90	(revision 5462)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/trcice_pisces.F90	(revision 5462)
@@ -0,0 +1,266 @@
+MODULE trcice_pisces
+   !!======================================================================
+   !!                         ***  MODULE trcice_pisces  ***
+   !! TOP :   initialisation of the PISCES biochemical model
+   !!======================================================================
+   !! History :  3.5  ! 2013    (M. Vancoppenolle, O. Aumont, G. Madec), original code
+   !! Comment ! probably not properly done when the second particle export
+   !! scheme (kriest) is used
+   !!----------------------------------------------------------------------
+#if defined key_pisces || defined key_pisces_reduced
+   !!----------------------------------------------------------------------
+   !!   'key_pisces'                                       PISCES bio-model
+   !!----------------------------------------------------------------------
+   !! trc_ice_pisces   : PISCES fake sea ice model setting
+   !!----------------------------------------------------------------------
+   USE par_trc         ! TOP parameters
+   USE par_pisces      ! PISCES parameters
+   USE oce_trc         ! Shared variables between ocean and passive tracers
+   USE trc             ! Passive tracers common variables 
+   USE phycst          ! Ocean physics parameters
+   USE sms_pisces      ! PISCES Source Minus Sink variables
+   USE in_out_manager
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_ice_ini_pisces ! called by trcini.F90 module
+
+CONTAINS
+
+   SUBROUTINE trc_ice_ini_pisces
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE trc_ice_ini_pisces ***
+      !!
+      !! ** Purpose :   PISCES fake sea ice model setting
+      !!    Method  :   Assign prescribe values to tracer concentrations in sea ice
+      !!
+      !! For levitating sea ice, constant ocean tracer concentrations also have to be defined. 
+      !! This is done specifically for Global, Arctic, Antarctic and Baltic regions
+      !! 
+      !! Sea ice concentrations are by default prescribed as follows
+      !!  trc_i = zratio * trc_o
+      !!
+      !! This formulation is modulated by the namelist parameter trc_ice_ratio
+      !!
+      !! trc_ice_ratio  * betw 0 and 1: prescribed ice/ocean tracer concentration ratio
+      !!                * -1 => the ice-ocean tracer concentration ratio follows the
+      !!                         ice-ocean salinity ratio
+      !!                * -2 => no ice-ocean tracer concentration is used
+      !!                        instead, the tracer concentration in sea ice 
+      !!                        is prescribed to trc_ice_prescr
+      !! 
+      !! cn_trc_o  specifies which disinctions are made for prescribed tracer concentration
+      !!                * 'GL' use global ocean values making distinction for Baltic Sea only
+      !!                * 'AA' use Arctic/Antarctic contrasted values, + Baltic
+      !!
+      !!----------------------------------------------------------------------
+
+                                        !--- Dummy variables
+      REAL(wp), DIMENSION(jptra,2) &
+               ::  zratio            ! effective ice-ocean tracer cc ratio
+      REAL(wp), DIMENSION(2) :: zrs  ! ice-ocean salinity ratio, 1 - global, 2- Baltic
+      REAL(wp) :: zsice_bal          ! prescribed ice salinity in the Baltic
+      REAL(wp) :: zsoce_bal          ! prescribed ocean salinity in the Baltic
+      REAL(wp) :: zfeoce_glo         ! prescribed iron concentration in the global ocean
+      REAL(wp) :: zfeoce_bal         ! prescribed iron concentration in the global ocean
+      INTEGER  :: jn                 ! dummy loop index
+
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) ' trc_ice_ini_pisces: Prescribed sea ice biogeochemistry '
+      IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~'
+
+      !--------------------------------------------
+      ! Initialize ocean prescribed concentrations
+      !--------------------------------------------
+      ! values taken from a 500 yr equilibrium run
+      ! used only in the levitating sea ice case with virtual salt / tracer
+      ! fluxes
+
+      !--- Global case 
+      IF ( cn_trc_o(jpdic) == 'GL ' ) trc_o(:,:,jpdic) =  1.99e-3_wp 
+      IF ( cn_trc_o(jpdoc) == 'GL ' ) trc_o(:,:,jpdoc) =  2.04e-5_wp 
+      IF ( cn_trc_o(jptal) == 'GL ' ) trc_o(:,:,jptal) =  2.31e-3_wp 
+      IF ( cn_trc_o(jpoxy) == 'GL ' ) trc_o(:,:,jpoxy) =  2.47e-4_wp
+      IF ( cn_trc_o(jpcal) == 'GL ' ) trc_o(:,:,jpcal) =  1.04e-8_wp
+      IF ( cn_trc_o(jppo4) == 'GL ' ) trc_o(:,:,jppo4) =  5.77e-7_wp / po4r 
+      IF ( cn_trc_o(jppoc) == 'GL ' ) trc_o(:,:,jppoc) =  1.27e-6_wp  
+#  if ! defined key_kriest
+      IF ( cn_trc_o(jpgoc) == 'GL ' ) trc_o(:,:,jpgoc) =  5.23e-8_wp  
+      IF ( cn_trc_o(jpbfe) == 'GL ' ) trc_o(:,:,jpbfe) =  9.84e-13_wp 
+#  else
+      IF ( cn_trc_o(jpnum) == 'GL ' ) trc_o(:,:,jpnum) = 0. ! could not get this value since did not use it
+#  endif
+      IF ( cn_trc_o(jpsil) == 'GL ' ) trc_o(:,:,jpsil) =  7.36e-6_wp  
+      IF ( cn_trc_o(jpdsi) == 'GL ' ) trc_o(:,:,jpdsi) =  1.07e-7_wp 
+      IF ( cn_trc_o(jpgsi) == 'GL ' ) trc_o(:,:,jpgsi) =  1.53e-8_wp
+      IF ( cn_trc_o(jpphy) == 'GL ' ) trc_o(:,:,jpphy) =  9.57e-8_wp
+      IF ( cn_trc_o(jpdia) == 'GL ' ) trc_o(:,:,jpdia) =  4.24e-7_wp
+      IF ( cn_trc_o(jpzoo) == 'GL ' ) trc_o(:,:,jpzoo) =  6.07e-7_wp
+      IF ( cn_trc_o(jpmes) == 'GL ' ) trc_o(:,:,jpmes) =  3.44e-7_wp
+      IF ( cn_trc_o(jpfer) == 'GL ' ) trc_o(:,:,jpfer) =  4.06e-10_wp
+      IF ( cn_trc_o(jpsfe) == 'GL ' ) trc_o(:,:,jpsfe) =  2.51e-11_wp
+      IF ( cn_trc_o(jpdfe) == 'GL ' ) trc_o(:,:,jpdfe) =  6.57e-12_wp
+      IF ( cn_trc_o(jpnfe) == 'GL ' ) trc_o(:,:,jpnfe) =  1.76e-11_wp
+      IF ( cn_trc_o(jpnch) == 'GL ' ) trc_o(:,:,jpnch) =  1.67e-7_wp
+      IF ( cn_trc_o(jpdch) == 'GL ' ) trc_o(:,:,jpdch) =  1.02e-7_wp
+      IF ( cn_trc_o(jpno3) == 'GL ' ) trc_o(:,:,jpno3) =  5.79e-6_wp / rno3 
+      IF ( cn_trc_o(jpnh4) == 'GL ' ) trc_o(:,:,jpnh4) =  3.22e-7_wp / rno3
+
+      !--- Arctic specificities (dissolved inorganic & DOM)
+      IF ( cn_trc_o(jpdic) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpdic) =  1.98e-3_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdoc) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpdoc) =  6.00e-6_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jptal) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jptal) =  2.13e-3_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpoxy) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpoxy) =  3.65e-4_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpcal) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpcal) =  1.50e-9_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jppo4) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jppo4) =  4.09e-7_wp / po4r ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jppoc) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jppoc) =  4.05e-7_wp  ; END WHERE ; ENDIF
+#  if ! defined key_kriest
+      IF ( cn_trc_o(jpgoc) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpgoc) =  2.84e-8_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpbfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpbfe) =  7.03e-13_wp ; END WHERE ; ENDIF
+#  else
+      IF ( cn_trc_o(jpnum) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpnum) =  0.00e-00_wp ; END WHERE ; ENDIF
+#  endif
+      IF ( cn_trc_o(jpsil) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpsil) =  6.87e-6_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdsi) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpdsi) =  1.73e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpgsi) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpgsi) =  7.93e-9_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpphy) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpphy) =  5.25e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdia) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpdia) =  7.75e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpzoo) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpzoo) =  3.34e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpmes) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpmes) =  2.49e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpfer) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpfer) =  1.43e-9_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpsfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpsfe) =  2.21e-11_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpdfe) =  2.04e-11_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpnfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpnfe) =  1.75e-11_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpnch) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpnch) =  1.46e-07_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdch) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpdch) =  2.36e-07_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpno3) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpno3) =  3.51e-06_wp / rno3 ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpnh4) == 'AA ' ) THEN ; WHERE( gphit(:,:) >= 00._wp ) ; trc_o(:,:,jpnh4) =  6.15e-08_wp / rno3 ; END WHERE ; ENDIF
+
+      !--- Antarctic specificities (dissolved inorganic & DOM)
+      IF ( cn_trc_o(jpdic) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpdic) =  2.20e-3_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdoc) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpdoc) =  7.02e-6_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jptal) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jptal) =  2.37e-3_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpoxy) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpoxy) =  3.42e-4_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpcal) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpcal) =  3.17e-9_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jppo4) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jppo4) =  1.88e-6_wp / po4r  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jppoc) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jppoc) =  1.13e-6_wp  ; END WHERE ; ENDIF
+#  if ! defined key_kriest
+      IF ( cn_trc_o(jpgoc) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpgoc) =  2.89e-8_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpbfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpbfe) =  5.63e-13_wp ; END WHERE ; ENDIF
+#  else
+      IF ( cn_trc_o(jpnum) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpnum) =  0.00e-00_wp ; END WHERE ; ENDIF
+#  endif
+      IF ( cn_trc_o(jpsil) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpsil) =  4.96e-5_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdsi) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpdsi) =  5.63e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpgsi) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpgsi) =  5.35e-8_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpphy) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpphy) =  8.10e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdia) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpdia) =  5.77e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpzoo) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpzoo) =  6.68e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpmes) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpmes) =  3.55e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpfer) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpfer) =  1.62e-10_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpsfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpsfe) =  2.29e-11_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpdfe) =  8.75e-12_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpnfe) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpnfe) =  1.48e-11_wp ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpnch) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpnch) =  2.02e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpdch) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpdch) =  1.60e-7_wp  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpno3) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpno3) =  2.64e-5_wp / rno3  ; END WHERE ; ENDIF
+      IF ( cn_trc_o(jpnh4) == 'AA ' ) THEN ; WHERE( gphit(:,:) <  00._wp ) ; trc_o(:,:,jpnh4) =  3.39e-7_wp / rno3  ; END WHERE ; ENDIF
+
+      !--- Baltic Sea particular case for ORCA configurations
+      IF( cp_cfg == "orca" ) THEN            ! Baltic mask
+         WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.    &
+                54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
+         trc_o(:,:,jpdic) = 1.14e-3_wp
+         trc_o(:,:,jpdoc) = 1.06e-5_wp
+         trc_o(:,:,jptal) = 1.16e-3_wp
+         trc_o(:,:,jpoxy) = 3.71e-4_wp
+         trc_o(:,:,jpcal) = 1.51e-9_wp
+         trc_o(:,:,jppo4) = 2.85e-9_wp / po4r
+         trc_o(:,:,jppoc) = 4.84e-7_wp
+#  if ! defined key_kriest
+         trc_o(:,:,jpgoc) = 1.05e-8_wp
+         trc_o(:,:,jpbfe) = 4.97e-13_wp
+#  else
+         trc_o(:,:,jpnum) = 0. ! could not get this value
+#  endif
+         trc_o(:,:,jpsil) = 4.91e-5_wp
+         trc_o(:,:,jpdsi) = 3.25e-7_wp
+         trc_o(:,:,jpgsi) = 1.93e-8_wp
+         trc_o(:,:,jpphy) = 6.64e-7_wp
+         trc_o(:,:,jpdia) = 3.41e-7_wp
+         trc_o(:,:,jpzoo) = 3.83e-7_wp
+         trc_o(:,:,jpmes) = 0.225e-6_wp
+         trc_o(:,:,jpfer) = 2.45e-9_wp
+         trc_o(:,:,jpsfe) = 3.89e-11_wp
+         trc_o(:,:,jpdfe) = 1.33e-11_wp
+         trc_o(:,:,jpnfe) = 2.62e-11_wp
+         trc_o(:,:,jpnch) = 1.17e-7_wp
+         trc_o(:,:,jpdch) = 9.69e-8_wp
+         trc_o(:,:,jpno3) = 5.36e-5_wp / rno3
+         trc_o(:,:,jpnh4) = 7.18e-7_wp / rno3
+         END WHERE
+      ENDIF ! cfg
+
+      !-----------------------------
+      ! Assign ice-ocean cc ratios 
+      !-----------------------------
+      ! 0 means zero concentration in sea ice
+      ! 1 means same concentration in the sea ice as in the ocean
+
+      ! Ice ocean salinity ratio
+      zsoce_bal   = 4. ; zsice_bal   = 2. !! Baltic ocean and sea ice salinities
+      zrs(1) = sice / soce                !! ice-ocean salinity ratio, global case
+      zrs(2) = zsice_bal / zsoce_bal      !! ice-ocean salinity ratio, Baltic case
+
+      DO jn = jp_pcs0, jp_pcs1
+         IF ( trc_ice_ratio(jn) >= 0._wp )  zratio(jn,:) = trc_ice_ratio(jn)
+         IF ( trc_ice_ratio(jn) == -1._wp ) zratio(jn,:) = zrs(:)
+         IF ( trc_ice_ratio(jn) == -2._wp ) zratio(jn,:) = -9999.99_wp
+      END DO
+
+      !-------------------------------
+      ! Sea ice tracer concentrations
+      !-------------------------------
+      DO jn = jp_pcs0, jp_pcs1
+         !-- Everywhere but in the Baltic
+         IF ( trc_ice_ratio(jn) >= -1._wp ) THEN !! no prescribed concentration
+                                              !! (typically everything but iron) 
+            trc_i(:,:,jn) = zratio(jn,1) * trc_o(:,:,jn) 
+         ELSE                                 !! prescribed concentration
+            trc_i(:,:,jn) = trc_ice_prescr(jn)
+         ENDIF
+       
+         !-- Baltic
+         IF( cp_cfg == "orca" ) THEN !! Baltic treated seperately for ORCA configs
+            IF ( trc_ice_ratio(jn) >= - 1._wp ) THEN !! no prescribed concentration
+                                                 !! (typically everything but iron) 
+               WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.    &
+                      54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
+                     trc_i(:,:,jn) = zratio(jn,2) * trc_o(:,:,jn) 
+               END WHERE
+            ELSE                                 !! prescribed tracer concentration in ice
+               WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.    &
+                   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp )
+                     trc_i(:,:,jn) = trc_ice_prescr(jn)
+               END WHERE
+            ENDIF ! trc_ice_ratio
+         ENDIF
+      !
+      END DO ! jn
+
+   END SUBROUTINE trc_ice_ini_pisces
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module                            No PISCES biochemical model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_ice_ini_pisces         ! Empty routine
+   END SUBROUTINE trc_ice_ini_pisces
+#endif
+
+   !!======================================================================
+END MODULE trcice_pisces
Index: branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90
===================================================================
--- branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 5420)
+++ branches/UKMO/restart_datestamp/NEMOGCM/NEMO/TOP_SRC/PISCES/trcini_pisces.F90	(revision 5462)
@@ -71,12 +71,13 @@
       USE p4zmort         !  Mortality terms for phytoplankton
       USE p4zlys          !  Calcite saturation
+      USE p4zsed          !  Sedimentation & burial
       !
       REAL(wp), SAVE :: sco2   =  2.312e-3_wp
-      REAL(wp), SAVE :: alka0  =  2.423e-3_wp
+      REAL(wp), SAVE :: alka0  =  2.426e-3_wp
       REAL(wp), SAVE :: oxyg0  =  177.6e-6_wp 
-      REAL(wp), SAVE :: po4    =  2.174e-6_wp 
+      REAL(wp), SAVE :: po4    =  2.165e-6_wp 
       REAL(wp), SAVE :: bioma0 =  1.000e-8_wp  
-      REAL(wp), SAVE :: silic1 =  91.65e-6_wp  
-      REAL(wp), SAVE :: no3    =  31.04e-6_wp * 7.625_wp
+      REAL(wp), SAVE :: silic1 =  91.51e-6_wp  
+      REAL(wp), SAVE :: no3    =  30.9e-6_wp * 7.625_wp
       !
       INTEGER  ::  ji, jj, jk, ierr
@@ -97,4 +98,5 @@
       ierr = ierr +  p4z_rem_alloc()
       ierr = ierr +  p4z_flx_alloc()
+      ierr = ierr +  p4z_sed_alloc()
       !
       IF( lk_mpp    )   CALL mpp_sum( ierr )
@@ -107,13 +109,4 @@
       CALL p4z_sms_init       !  Maint routine
       !                                            ! Time-step
-      rfact   = rdttrc(1)                          ! ---------
-      rfactr  = 1. / rfact
-      rfact2  = rfact / FLOAT( nrdttrc )
-      rfact2r = 1. / rfact2
-
-      IF(lwp) WRITE(numout,*) '    Passive Tracer  time step    rfact  = ', rfact, ' rdt = ', rdttra(1)
-      IF(lwp) write(numout,*) '    PISCES  Biology time step    rfact2 = ', rfact2
-
-
 
       ! Set biological ratios
@@ -165,6 +158,4 @@
       END IF
 
-      ! Time step duration for biology
-      xstep = rfact2 / rday
 
       CALL p4z_sink_init      !  vertical flux of particulate organic matter
