diff --git a/cfgs/SHARED/field_def_nemo-pisces.xml b/cfgs/SHARED/field_def_nemo-pisces.xml
index 5ed5b7043bf43d7f5971c1f101fe2b2503cb1ea5..6b63d355ed943ae8fe3412009efeca44863d1e9d 100644
--- a/cfgs/SHARED/field_def_nemo-pisces.xml
+++ b/cfgs/SHARED/field_def_nemo-pisces.xml
@@ -173,23 +173,23 @@
   <!-- PISCES additional diagnostics on T grid  -->
 
   <field_group id="diad_T" grid_ref="grid_T_2D">
-    <field id="PH"          long_name="PH"                                      unit="1"          grid_ref="grid_T_3D" />
-    <field id="CO3"         long_name="Bicarbonates"                            unit="mol/m3"     grid_ref="grid_T_3D" />
-    <field id="CO3sat"      long_name="CO3 saturation"                          unit="mol/m3"     grid_ref="grid_T_3D" />
+    <field id="PH"          long_name="PH"                                      unit="1"          grid_ref="grid_T_3D_inner" />
+    <field id="CO3"         long_name="Bicarbonates"                            unit="mol/m3"     grid_ref="grid_T_3D_inner" />
+    <field id="CO3sat"      long_name="CO3 saturation"                          unit="mol/m3"     grid_ref="grid_T_3D_inner" />
     <field id="PAR"         long_name="Photosynthetically Available Radiation"  unit="W/m2"       grid_ref="grid_T_3D" />
-    <field id="PPPHYN"      long_name="Primary production of nanophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PPPHYP"      long_name="Primary production of picophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PPPHYD"      long_name="Primary production of diatoms"           unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PPNEWN"      long_name="New Primary production of nanophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D" />
+    <field id="PPPHYN"      long_name="Primary production of nanophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="PPPHYP"      long_name="Primary production of picophyto"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="PPPHYD"      long_name="Primary production of diatoms"           unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="PPNEWN"      long_name="New Primary production of nanophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="PPNEWP"      long_name="New Primary production of picophyto"     unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PPNEWD"      long_name="New Primary production of diatoms"       unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PBSi"        long_name="Primary production of Si diatoms"        unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PFeN"        long_name="Primary production of nano iron"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PFeP"        long_name="Primary production of pico iron"         unit="molC/m3/s"  grid_ref="grid_T_3D" />
-    <field id="PFeD"        long_name="Primary production of diatoms iron"      unit="mol/m3/s"   grid_ref="grid_T_3D" />
+    <field id="PPNEWD"      long_name="New Primary production of diatoms"       unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="PBSi"        long_name="Primary production of Si diatoms"        unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="PFeN"        long_name="Primary production of nano iron"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="PFeP"        long_name="Primary production of pico iron"         unit="molC/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="PFeD"        long_name="Primary production of diatoms iron"      unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <field id="xfracal"     long_name="Calcifying fraction"                     unit="1"          grid_ref="grid_T_3D" />
     <field id="PCAL"        long_name="Calcite production"                      unit="mol/m3/s"   grid_ref="grid_T_3D" />
-    <field id="DCAL"        long_name="Calcite dissolution"                     unit="mol/m3/s"   grid_ref="grid_T_3D" />
+    <field id="DCAL"        long_name="Calcite dissolution"                     unit="mol/m3/s"   grid_ref="grid_T_3D_inner" />
     <field id="GRAZ1"       long_name="Grazing by microzooplankton"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="GRAZ2"       long_name="Grazing by mesozooplankton"              unit="mol/m3/s"   grid_ref="grid_T_3D" />
     <field id="REMIN"       long_name="Oxic remineralization of OM"             unit="mol/m3/s"   grid_ref="grid_T_3D" />
@@ -197,22 +197,22 @@
     <field id="REMINP"       long_name="Oxic remineralization rate of POC"      unit="d-1"        grid_ref="grid_T_3D" />
     <field id="REMING"       long_name="Oxic remineralization rate of GOC"      unit="d-1"        grid_ref="grid_T_3D" />
     <field id="Nfix"        long_name="Nitrogen fixation"                       unit="mol/m3/s"   grid_ref="grid_T_3D" />
-    <field id="Mumax"       long_name="Maximum growth rate"                     unit="s-1"        grid_ref="grid_T_3D" />
-    <field id="MuN"         long_name="Realized growth rate for nanophyto"      unit="s-1"        grid_ref="grid_T_3D" />
-    <field id="MuP"         long_name="Realized growth rate for picophyto"      unit="s-1"        grid_ref="grid_T_3D" />
-    <field id="MuD"         long_name="Realized growth rate for diatomes"       unit="s-1"        grid_ref="grid_T_3D" />
+    <field id="Mumax"       long_name="Maximum growth rate"                     unit="s-1"        grid_ref="grid_T_3D_inner" />
+    <field id="MuN"         long_name="Realized growth rate for nanophyto"      unit="s-1"        grid_ref="grid_T_3D_inner" />
+    <field id="MuP"         long_name="Realized growth rate for picophyto"      unit="s-1"        grid_ref="grid_T_3D_inner" />
+    <field id="MuD"         long_name="Realized growth rate for diatomes"       unit="s-1"        grid_ref="grid_T_3D_inner" />
     <field id="MunetN"      long_name="Net growth rate for nanophyto"           unit="s-1"        grid_ref="grid_T_3D" />
     <field id="MunetP"      long_name="Net growth rate for picophyto"           unit="s-1"        grid_ref="grid_T_3D" />
     <field id="MunetD"      long_name="Net growth rate for diatomes"            unit="s-1"        grid_ref="grid_T_3D" />
-    <field id="LNnut"       long_name="Nutrient limitation term in Nanophyto"   unit=""           grid_ref="grid_T_3D" />
-    <field id="LPnut"       long_name="Nutrient limitation term in Picophyto"   unit="-"          grid_ref="grid_T_3D" />
-    <field id="LDnut"       long_name="Nutrient limitation term in Diatoms"     unit=""           grid_ref="grid_T_3D" />
+    <field id="LNnut"       long_name="Nutrient limitation term in Nanophyto"   unit=""           grid_ref="grid_T_3D_inner" />
+    <field id="LPnut"       long_name="Nutrient limitation term in Picophyto"   unit="-"          grid_ref="grid_T_3D_inner" />
+    <field id="LDnut"       long_name="Nutrient limitation term in Diatoms"     unit=""           grid_ref="grid_T_3D_inner" />
     <field id="LNFe"        long_name="Iron limitation term in Nanophyto"       unit=""           grid_ref="grid_T_3D" />
     <field id="LPFe"        long_name="Iron limitation term in Picophyto"       unit="-"          grid_ref="grid_T_3D" />
     <field id="LDFe"        long_name="Iron limitation term in Diatoms"         unit=""           grid_ref="grid_T_3D" />
-    <field id="LNlight"     long_name="Light limitation term in Nanophyto"      unit=""           grid_ref="grid_T_3D" />
-    <field id="LPlight"     long_name="Light limitation term in Picophyto"      unit="-"          grid_ref="grid_T_3D" />
-    <field id="LDlight"     long_name="Light limitation term in Diatoms"        unit=""           grid_ref="grid_T_3D" />
+    <field id="LNlight"     long_name="Light limitation term in Nanophyto"      unit=""           grid_ref="grid_T_3D_inner" />
+    <field id="LPlight"     long_name="Light limitation term in Picophyto"      unit="-"          grid_ref="grid_T_3D_inner" />
+    <field id="LDlight"     long_name="Light limitation term in Diatoms"        unit=""           grid_ref="grid_T_3D_inner" />
     <field id="SIZEN"       long_name="Mean relative size of nanophyto."        unit="-"          grid_ref="grid_T_3D" />
     <field id="SIZEP"       long_name="Mean relative size of picophyto."        unit="-"          grid_ref="grid_T_3D" />
     <field id="SIZED"       long_name="Mean relative size of diatoms"           unit="-"          grid_ref="grid_T_3D" />
@@ -276,9 +276,9 @@
 
 
     <!-- dbio_T on T grid : variables available with diaar5 -->
-    <field id="TPP"         long_name="Total Primary production of phyto"                   unit="mol/m3/s"  grid_ref="grid_T_3D" />
-    <field id="TPNEW"       long_name="New Primary production of phyto"                     unit="mol/m3/s"  grid_ref="grid_T_3D" />
-    <field id="TPBFE"       long_name="Total biogenic iron production"                      unit="mol/m3/s"  grid_ref="grid_T_3D" />
+    <field id="TPP"         long_name="Total Primary production of phyto"                   unit="mol/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="TPNEW"       long_name="New Primary production of phyto"                     unit="mol/m3/s"  grid_ref="grid_T_3D_inner" />
+    <field id="TPBFE"       long_name="Total biogenic iron production"                      unit="mol/m3/s"  grid_ref="grid_T_3D_inner" />
     <field id="INTDIC"      long_name="DIC content"                                         unit="kg/m2"                          />
     <field id="O2MIN"       long_name="Oxygen minimum concentration"                        unit="mol/m3"                          />
     <field id="ZO2MIN"      long_name="Depth of oxygen minimum concentration"               unit="m"                              />
diff --git a/src/OCE/DOM/domtile.F90 b/src/OCE/DOM/domtile.F90
index dd391b709044048c4aade3bd090ddb86c87d713b..3c1f99b5944e36b343479b30ac6e96b6d633b0fe 100644
--- a/src/OCE/DOM/domtile.F90
+++ b/src/OCE/DOM/domtile.F90
@@ -52,16 +52,6 @@ CONTAINS
       !!----------------------------------------------------------------------
       IF( ln_tile .AND. nn_hls /= 2 ) CALL ctl_stop('dom_tile_init: Tiling is only supported for nn_hls = 2')
 
-      ntile = 0                     ! Initialise to full domain
-      nijtile = 1
-      ntsi = Nis0
-      ntsj = Njs0
-      ntei = Nie0
-      ntej = Nje0
-      nthl = 0
-      nthr = 0
-      nthb = 0
-      ntht = 0
       l_istiled = .FALSE.
 
       IF( ln_tile ) THEN            ! Calculate tile domain indices
diff --git a/src/OCE/LBC/mppini.F90 b/src/OCE/LBC/mppini.F90
index da174b39072a4294fd5303db17f984660faa4851..6697cc4151a1ccbf0ba68e1dc61383141882e16d 100644
--- a/src/OCE/LBC/mppini.F90
+++ b/src/OCE/LBC/mppini.F90
@@ -1415,6 +1415,17 @@ ENDIF
       !
       jpkm1 = jpk-1                             !   "           "
       !
+      ntile = 0                     ! Initialise to full domain
+      nijtile = 1
+      ntsi = Nis0
+      ntsj = Njs0
+      ntei = Nie0
+      ntej = Nje0
+      nthl = 0
+      nthr = 0
+      nthb = 0
+      ntht = 0
+      !
    END SUBROUTINE init_doloop
 
    
diff --git a/src/TOP/AGE/trcsms_age.F90 b/src/TOP/AGE/trcsms_age.F90
index a8283edf6df534771cbac7aaa40fa40946b6f4e1..f3352a5c41c2cac32bddf15065c31e071cdb196b 100644
--- a/src/TOP/AGE/trcsms_age.F90
+++ b/src/TOP/AGE/trcsms_age.F90
@@ -46,7 +46,7 @@ CONTAINS
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt              ! ocean time-step index
       INTEGER, INTENT(in) ::   Kbb, Kmm, Krhs  ! ocean time level
-      INTEGER ::   jn, jk   ! dummy loop index
+      INTEGER ::   jk   ! dummy loop index
       !!----------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('trc_sms_age')
diff --git a/src/TOP/AGE/trcwri_age.F90 b/src/TOP/AGE/trcwri_age.F90
index 98b7c6863ed6ba0d6eac44a892864fc0a6980dc3..3df7b3e7e99e766d7d7f4bd5f7f57555bbfb9998 100644
--- a/src/TOP/AGE/trcwri_age.F90
+++ b/src/TOP/AGE/trcwri_age.F90
@@ -28,7 +28,6 @@ CONTAINS
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in)  :: Kmm  ! time level indices
       CHARACTER (len=20)   :: cltra
-      INTEGER              :: jn
       !!---------------------------------------------------------------------
 
       ! write the tracer concentrations in the file
diff --git a/src/TOP/PISCES/P4Z/p4zlys.F90 b/src/TOP/PISCES/P4Z/p4zlys.F90
index cc0f80ff519095c8e27752d39de1e53ed3c00b73..729daab9dfa7e0e306bfb81915ddf6cc4a2fefee 100644
--- a/src/TOP/PISCES/P4Z/p4zlys.F90
+++ b/src/TOP/PISCES/P4Z/p4zlys.F90
@@ -64,21 +64,15 @@ CONTAINS
       !
       INTEGER  ::   ji, jj, jk, jn
       REAL(wp) ::   zdispot, zrhd, zcalcon
-      REAL(wp) ::   zomegaca, zexcess, zexcess0, zkd, zco3, ztra
+      REAL(wp) ::   zomegaca, zexcess, zexcess0, zkd
       CHARACTER (len=25) ::   charout
-      REAL(wp), DIMENSION(A2D(0),jpk) :: zhinit, zhi
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  zcaldiss, zbicarb, zco3sat
+      REAL(wp), DIMENSION(A2D(0),jpk) ::   zco3, zcaldiss, zhinit, zhi, zco3sat
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )  CALL timing_start('p4z_lys')
       !
-      IF( iom_use( "CO3"    ) ) ALLOCATE( zbicarb (A2D(0),jpk) )  
-      IF( iom_use( "CO3sat" ) ) ALLOCATE( zco3sat (A2D(0),jpk) ) 
-      IF( iom_use( "DCAL"   ) ) ALLOCATE( zcaldiss(A2D(0),jpk) ) 
 
-      DO_3D( 0, 0, 0, 0, 1, jpkm1)
-        zhinit(ji,jj,jk) = hi(ji,jj,jk) / ( rhd(ji,jj,jk) + 1._wp )
-      END_3D
+      zhinit  (:,:,:) = hi(:,:,:) / ( rhd(:,:,:) + 1._wp )
       !
       !     -------------------------------------------
       !     COMPUTE [CO3--] and [H+] CONCENTRATIONS
@@ -87,20 +81,25 @@ CONTAINS
       CALL solve_at_general( zhinit, zhi, Kbb )
 
       DO_3D( 0, 0, 0, 0, 1, jpkm1)
-         !
-         zco3           = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
+         zco3(ji,jj,jk) = tr(ji,jj,jk,jpdic,Kbb) * ak13(ji,jj,jk) * ak23(ji,jj,jk) / (zhi(ji,jj,jk)**2   &
             &             + ak13(ji,jj,jk) * zhi(ji,jj,jk) + ak13(ji,jj,jk) * ak23(ji,jj,jk) + rtrn )
          hi  (ji,jj,jk) = zhi(ji,jj,jk) * ( rhd(ji,jj,jk) + 1._wp )
+      END_3D
 
-         ! CALCULATE DEGREE OF CACO3 SATURATION AND CORRESPONDING
-         ! DISSOLOUTION AND PRECIPITATION OF CACO3 (BE AWARE OF MGCO3)
+      !     ---------------------------------------------------------
+      !        CALCULATE DEGREE OF CACO3 SATURATION AND CORRESPONDING
+      !        DISSOLOUTION AND PRECIPITATION OF CACO3 (BE AWARE OF
+      !        MGCO3)
+      !     ---------------------------------------------------------
 
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
 
          ! DEVIATION OF [CO3--] FROM SATURATION VALUE
          ! Salinity dependance in zomegaca and divide by rhd to have good units
          zcalcon  = calcon * ( salinprac(ji,jj,jk) / 35._wp )
          zrhd    = rhd(ji,jj,jk) + 1._wp
-         zomegaca = ( zcalcon * zco3 ) / ( aksp(ji,jj,jk) * zrhd + rtrn )
+         zomegaca = ( zcalcon * zco3(ji,jj,jk) ) / ( aksp(ji,jj,jk) * zrhd + rtrn )
+         zco3sat(ji,jj,jk) = aksp(ji,jj,jk) * zrhd / ( zcalcon + rtrn )
 
          ! SET DEGREE OF UNDER-/SUPERSATURATION
          excess(ji,jj,jk) = 1._wp - zomegaca
@@ -118,36 +117,44 @@ CONTAINS
 
         !  CHANGE OF [CO3--] , [ALK], PARTICULATE [CACO3],
         !       AND [SUM(CO2)] DUE TO CACO3 DISSOLUTION/PRECIPITATION
-        ztra =  zdispot * rfact2 / rmtss ! calcite dissolution
-        !
-        tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + 2. * ztra
-        tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) -      ztra
-        tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) +      ztra
-        ! 
-        IF( iom_use( "CO3"    ) ) zbicarb (ji,jj,jk) = zco3                       !  bicarbonate
-        IF( iom_use( "CO3sat" ) ) zco3sat (ji,jj,jk) = zco3 / ( zomegaca + rtrn )  ! calcite saturation
-        IF( iom_use( "DCAL"   ) ) zcaldiss(ji,jj,jk) = ztra                        ! calcite dissolution
+        zcaldiss(ji,jj,jk)  = zdispot * rfact2 / rmtss ! calcite dissolution
         !
+        tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + 2. * zcaldiss(ji,jj,jk)
+        tr(ji,jj,jk,jpcal,Krhs) = tr(ji,jj,jk,jpcal,Krhs) -      zcaldiss(ji,jj,jk)
+        tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) +      zcaldiss(ji,jj,jk)
       END_3D
       !
 
       IF( lk_iomput .AND. knt == nrdttrc ) THEN
-         IF( iom_use( "PH" ) )  CALL iom_put( "PH" , -1. * LOG10( MAX( hi(A2D(0),1:jpk), rtrn ) ) * tmask(A2D(0),1:jpk) )
-         IF( iom_use( "CO3"    ) ) THEN  ! bicarbonate 
-            zbicarb(A2D(0),jpk)  = 0. 
-            CALL iom_put( "CO3"   , zbicarb(A2D(0),1:jpk) * 1.e+3 * tmask(A2D(0),1:jpk) ) 
-            DEALLOCATE( zbicarb )
+          IF( l_dia ) THEN
+             ALLOCATE( zw3d(A2D(0),jpk) ) 
+             zw3d(A2D(0),jpk) = 0._wp
+          ENDIF
+          IF( iom_use( "PH" ) ) THEN
+             DO_3D( 0, 0, 0, 0, 1, jpkm1)
+                zw3d(ji,jj,jk) = -1. * LOG10( MAX( hi(ji,jj,jk), rtrn ) ) * tmask(ji,jj,jk) 
+             END_3D
+             CALL iom_put( "PH" , zw3d )
+          ENDIF
+          IF( iom_use( "CO3" ) ) THEN  ! bicarbonate 
+             DO_3D( 0, 0, 0, 0, 1, jpkm1)
+                zw3d(ji,jj,jk) = zco3(ji,jj,jk) * 1.e+3 * tmask(ji,jj,jk) 
+             END_3D
+             CALL iom_put( "CO3", zw3d ) 
          ENDIF
-         IF( iom_use( "CO3sat" ) ) THEN   ! calcite saturation
-            zco3sat(A2D(0),jpk)  = 0. 
-            CALL iom_put( "CO3sat", zco3sat(A2D(0),1:jpk) * 1.e+3 * tmask(A2D(0),1:jpk) ) 
-            DEALLOCATE( zco3sat )
+         IF( iom_use( "CO3sat" ) ) THEN  ! calcite saturation
+             DO_3D( 0, 0, 0, 0, 1, jpkm1)
+                zw3d(ji,jj,jk) = zco3sat(ji,jj,jk) * 1.e+3 * tmask(ji,jj,jk) 
+             END_3D
+             CALL iom_put( "CO3sat", zw3d ) 
          ENDIF
          IF( iom_use( "DCAL" ) ) THEN  ! calcite dissolution
-            zcaldiss(A2D(0),jpk)  = 0. 
-            CALL iom_put( "DCAL", zcaldiss(A2D(0),1:jpk) * 1.e+3 * rfact2r * tmask(A2D(0),1:jpk) ) 
-            DEALLOCATE( zcaldiss )
+             DO_3D( 0, 0, 0, 0, 1, jpkm1)
+                zw3d(ji,jj,jk) = zcaldiss(ji,jj,jk) * 1.e+3 * rfact2r * tmask(ji,jj,jk) 
+             END_3D
+             CALL iom_put( "DCAL", zw3d ) 
          ENDIF              
+         IF( l_dia ) DEALLOCATE( zw3d )
       ENDIF
       !
       IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
diff --git a/src/TOP/PISCES/P4Z/p4zprod.F90 b/src/TOP/PISCES/P4Z/p4zprod.F90
index 357611278a1c4d7e49cead02c90f1a584c39989a..00bbcfe8c63b53a28727230aff851fad6a46992a 100644
--- a/src/TOP/PISCES/P4Z/p4zprod.F90
+++ b/src/TOP/PISCES/P4Z/p4zprod.F90
@@ -75,28 +75,48 @@ CONTAINS
       REAL(wp) ::   zproddoc, zprodsil, zprodfer, zprodlig
       REAL(wp) ::   zpislopen, zpisloped, zfact
       REAL(wp) ::   zratiosi, zmaxsi, zlimfac, zsizetmp, zfecnm, zfecdm
-      REAL(wp) ::   zprod, zval, zmxl_fac, zmxl_chl
-      REAL(wp) ::   zprmax, zpislopeadn,zpislopeadd,zprchld, zprchln   
+      REAL(wp) ::   zprod, zval
       CHARACTER (len=25) :: charout
-      REAL(wp), DIMENSION(A2D(0),jpk) :: zprdia, zprbio, zysopt, zmxl    
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprmaxn,zprmaxd
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zpislopeadn, zpislopeadd, zysopt  
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zprdia, zprbio, zprchld, zprchln   
       REAL(wp), DIMENSION(A2D(0),jpk) :: zprorcan, zprorcad, zprofed, zprofen
       REAL(wp), DIMENSION(A2D(0),jpk) :: zpronewn, zpronewd
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
+      REAL(wp), DIMENSION(A2D(0),jpk) :: zmxl_fac, zmxl_chl
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zpligprod
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('p4z_prod')
       !
+      !  Allocate temporary workspace
+      !
+      zprorcan(:,:,:) = 0._wp ; zprorcad(:,:,:) = 0._wp ; zprofed(:,:,:) = 0._wp
+      zprofen (:,:,:) = 0._wp ; zysopt  (:,:,:) = 0._wp
+      zpronewn(:,:,:) = 0._wp ; zpronewd(:,:,:) = 0._wp ; zprdia(:,:,:) = 0._wp
+      zprbio  (:,:,:) = 0._wp ; zprchld (:,:,:) = 0._wp ; zprchln(:,:,:) = 0._wp 
+      zmxl_fac(:,:,:) = 0._wp ; zmxl_chl(:,:,:) = 0._wp 
+
+      ! Computation of the maximimum production. Based on a Q10 description
+      ! of the thermal dependency
+      ! Parameters are taken from Bissinger et al. (2008)
+      zprmaxn(:,:,:) = 0.65_wp * r1_rday * tgfunc(:,:,:)
+      zprmaxd(:,:,:) = zprmaxn(:,:,:)
 
       ! Intermittency is supposed to have a similar effect on production as 
-      ! day length (Shatwell et al., 2012).  
+      ! day length (Shatwell et al., 2012). The correcting factor is zmxl_fac. 
+      ! zmxl_chl is the fractional day length and is used to compute the mean
+      ! PAR during daytime. The effect of mixing is computed using the 
+      ! absolute light level definition of the euphotic zone
       ! ------------------------------------------------------------------------- 
       IF ( ln_p4z_dcyc ) THEN    ! Diurnal cycle in PISCES
+
          DO_3D( 0, 0, 0, 0, 1, jpkm1)
             IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
                IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
                   zval = MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
                ENDIF
-               zmxl(ji,jj,jk) = zval
+               zmxl_chl(ji,jj,jk) = zval / 24.
+               zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval )
             ENDIF
          END_3D
  
@@ -107,32 +127,21 @@ CONTAINS
                IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
                   zval = zval * MIN(1., heup_01(ji,jj) / ( hmld(ji,jj) + rtrn ))
                ENDIF
-               zmxl(ji,jj,jk) = zval
+               zmxl_chl(ji,jj,jk) = zval / 24.
+               zmxl_fac(ji,jj,jk) = 1.0 - exp( -0.26 * zval )
             ENDIF
          END_3D
 
       ENDIF
 
+      zprbio(:,:,:) = zprmaxn(:,:,:) * zmxl_fac(:,:,:)
+      zprdia(:,:,:) = zprmaxd(:,:,:) * zmxl_fac(:,:,:)
+
       ! The formulation proposed by Geider et al. (1997) has been modified 
       ! to exclude the effect of nutrient limitation and temperature in the PI
       ! curve following Vichi et al. (2007)
       ! -----------------------------------------------------------------------
       DO_3D( 0, 0, 0, 0, 1, jpkm1)
-         !
-         ! Computation of the maximimum production. Based on a Q10 description
-         ! of the thermal dependency. Parameters are taken from Bissinger et al. (2008)
-         zprmax = 0.65_wp * r1_rday * tgfunc(ji,jj,jk)
-         !
-         ! The correcting factor of Intermittency is zmxl_fac. 
-         ! zmxl_chl is the fractional day length and is used to compute the mean
-         ! PAR during daytime. The effect of mixing is computed using the 
-         ! absolute light level definition of the euphotic zone
-         zmxl_fac = 1.0 - EXP( -0.26 * zmxl(ji,jj,jk) )
-         zmxl_chl = zmxl(ji,jj,jk) / 24.
-         !
-         zprbio(ji,jj,jk) = zprmax * zmxl_fac 
-         zprdia(ji,jj,jk) = zprmax * zmxl_fac
-         !
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             ztn         = MAX( 0., ts(ji,jj,jk,jp_tem,Kmm) - 15. )
             zadap       = xadap * ztn / ( 2.+ ztn )
@@ -145,59 +154,66 @@ CONTAINS
             ! improved or removed in future versions of the model
 
             ! Nanophytoplankton
-            zpislopeadn = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
+            zpislopeadn(ji,jj,jk) = pislopen * ( 1.+ zadap  * EXP( -0.25 * enano(ji,jj,jk) ) )  &
             &                   * tr(ji,jj,jk,jpnch,Kbb) /( tr(ji,jj,jk,jpphy,Kbb) * 12. + rtrn)
 
             ! Diatoms
-            zpislopeadd = (pislopen * zconctemp2 + pisloped * zconctemp) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )   &
+            zpislopeadd(ji,jj,jk) = (pislopen * zconctemp2 + pisloped * zconctemp) / ( tr(ji,jj,jk,jpdia,Kbb) + rtrn )   &
             &                   * tr(ji,jj,jk,jpdch,Kbb) /( tr(ji,jj,jk,jpdia,Kbb) * 12. + rtrn)
-            !
-            ! Computation of production function for Carbon ; Actual light levels are used here 
-            zfact = 1._wp / ( ( r1_rday + bresp * r1_rday ) * zmxl_fac * rday + rtrn)
-            zpislopen = zpislopeadn * zfact
-            zpisloped = zpislopeadd * zfact
-            !
-            zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
-            zprdia(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
          ENDIF
-         !  Computation of a proxy of the N/C quota from nutrient limitation 
-         !  and light limitation. Steady state is assumed to allow the computation
-         !  ----------------------------------------------------------------------
-         zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
-         &      * zprmax / ( zprbio(ji,jj,jk) + rtrn )
-         quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
-         !
-         zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   &
-         &      * zprmax / ( zprdia(ji,jj,jk) + rtrn )
-         quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
-
-         ! Sea-ice effect on production : No production is assumed below sea ice
-         zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
-         zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+      END_3D
 
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
-            !  Computation of production function for Chlorophyll
-            !  Mean light level in the mixed layer (when appropriate)
-            !  is used here (acclimation is in general slower than 
-            !  the characteristic time scales of vertical mixing)
-            !  ------------------------------------------------------
-            zfact = 1._wp / ( zprmax * zmxl_chl * rday + rtrn )
-            zpislopen = zpislopeadn * zfact 
-            zpisloped = zpislopeadd * zfact
-            zprchln = zprmax * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
-            zprchld = zprmax * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
-            ! Si/C of diatoms
-            ! ------------------------
-            ! Si/C increases with iron stress and silicate availability
-            ! Si/C is arbitrariliy increased for very high Si concentrations
-            ! to mimic the very high ratios observed in the Southern Ocean (zsilfac)
-            ! A parameterization derived from Flynn (2003) is used for the control
-            ! when Si is not limiting which is similar to the parameterisation
-            ! proposed by Gurney and Davidson (1999).
-            ! -----------------------------------------------------------------------
+             ! Computation of production function for Carbon
+             ! Actual light levels are used here 
+             ! ----------------------------------------------
+             zpislopen = zpislopeadn(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+             &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
+             zpisloped = zpislopeadd(ji,jj,jk) / ( ( r1_rday + bresp * r1_rday ) &
+             &            * zmxl_fac(ji,jj,jk) * rday + rtrn)
+             zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1.- EXP( -zpislopen * enano(ji,jj,jk) )  )
+             zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediat(ji,jj,jk) )  )
+
+             !  Computation of production function for Chlorophyll
+             !  Mean light level in the mixed layer (when appropriate)
+             !  is used here (acclimation is in general slower than 
+             !  the characteristic time scales of vertical mixing)
+             !  ------------------------------------------------------
+             zpislopen = zpislopeadn(ji,jj,jk) / ( zprmaxn(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+             zpisloped = zpislopeadd(ji,jj,jk) / ( zprmaxd(ji,jj,jk) * zmxl_chl(ji,jj,jk) * rday + rtrn )
+             zprchln(ji,jj,jk) = zprmaxn(ji,jj,jk) * ( 1.- EXP( -zpislopen * enanom(ji,jj,jk) ) )
+             zprchld(ji,jj,jk) = zprmaxd(ji,jj,jk) * ( 1.- EXP( -zpisloped * ediatm(ji,jj,jk) ) )
+         ENDIF
+      END_3D
+
+      !  Computation of a proxy of the N/C quota from nutrient limitation 
+      !  and light limitation. Steady state is assumed to allow the computation
+      !  ----------------------------------------------------------------------
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+          zval = MIN( xnanopo4(ji,jj,jk), ( xnanonh4(ji,jj,jk) + xnanono3(ji,jj,jk) ) )   &
+          &      * zprmaxn(ji,jj,jk) / ( zprbio(ji,jj,jk) + rtrn )
+          quotan(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
+          zval = MIN( xdiatpo4(ji,jj,jk), ( xdiatnh4(ji,jj,jk) + xdiatno3(ji,jj,jk) ) )   &
+          &      * zprmaxd(ji,jj,jk) / ( zprdia(ji,jj,jk) + rtrn )
+          quotad(ji,jj,jk) = MIN( 1., 0.3 + 0.7 * zval )
+      END_3D
+
+
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+
+          IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
+             ! Si/C of diatoms
+             ! ------------------------
+             ! Si/C increases with iron stress and silicate availability
+             ! Si/C is arbitrariliy increased for very high Si concentrations
+             ! to mimic the very high ratios observed in the Southern Ocean (zsilfac)
+             ! A parameterization derived from Flynn (2003) is used for the control
+             ! when Si is not limiting which is similar to the parameterisation
+             ! proposed by Gurney and Davidson (1999).
+             ! -----------------------------------------------------------------------
             zlim  = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi1 )
-            !
-            zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmax + rtrn )
+            zsilim = xlimdia(ji,jj,jk) * zprdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
             zsiborn = tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb) * tr(ji,jj,jk,jpsil,Kbb)
             IF (gphit(ji,jj) < -30 ) THEN
               zsilfac = 1. + 2. * zsiborn / ( zsiborn + xksi2**3 )
@@ -212,9 +228,21 @@ CONTAINS
             ELSE
                zysopt(ji,jj,jk) = zlim * zsilfac * grosip * 1.0 * zsilim**0.7 * zmaxsi
             ENDIF
+        ENDIF
+      END_3D
 
-            ! Computation of the various production  and nutrient uptake terms
-            ! ---------------------------------------------------------------
+      ! Sea-ice effect on production
+      ! No production is assumed below sea ice
+      ! -------------------------------------- 
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         zprbio(ji,jj,jk) = zprbio(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+         zprdia(ji,jj,jk) = zprdia(ji,jj,jk) * ( 1. - fr_i(ji,jj) )
+      END_3D
+
+      ! Computation of the various production  and nutrient uptake terms
+      ! ---------------------------------------------------------------
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             !  production terms for nanophyto. (C)
             zprorcan(ji,jj,jk) = zprbio(ji,jj,jk)  * xlimphy(ji,jj,jk) * tr(ji,jj,jk,jpphy,Kbb) * rfact2
 
@@ -227,7 +255,7 @@ CONTAINS
             ! size at time step t+1 and is thus updated at the end of the 
             ! current time step
             ! --------------------------------------------------------------------
-            zlimfac = xlimphy(ji,jj,jk) * zprchln / ( zprmax + rtrn )
+            zlimfac = xlimphy(ji,jj,jk) * zprchln(ji,jj,jk) / ( zprmaxn(ji,jj,jk) + rtrn )
             zsizetmp = 1.0 + 1.3 * ( xsizern - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
             sizena(ji,jj,jk) = min(xsizern, max( sizena(ji,jj,jk), zsizetmp ) )
 
@@ -238,7 +266,7 @@ CONTAINS
             zfecnm = xqfuncfecn(ji,jj,jk) + ( fecnm - xqfuncfecn(ji,jj,jk) ) * ( xnanono3(ji,jj,jk) + xnanonh4(ji,jj,jk) )
             zratio = 1.0 - MIN(1.0,tr(ji,jj,jk,jpnfe,Kbb) / ( tr(ji,jj,jk,jpphy,Kbb) * zfecnm + rtrn ) )
             zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
-            zprofen(ji,jj,jk) = zfecnm * zprmax * ( 1.0 - fr_i(ji,jj) )  &
+            zprofen(ji,jj,jk) = zfecnm * zprmaxn(ji,jj,jk) * ( 1.0 - fr_i(ji,jj) )  &
             &          * (1. + 0.8 * xnanono3(ji,jj,jk) / ( rtrn + xnanono3(ji,jj,jk)  &
             &          + xnanonh4(ji,jj,jk) ) * (1. - xnanofer(ji,jj,jk) ) )   &
             &          * xnanofer(ji,jj,jk) * zmax * tr(ji,jj,jk,jpphy,Kbb) * rfact2
@@ -254,7 +282,7 @@ CONTAINS
             ! size at time step t+1 and is thus updated at the end of the 
             ! current time step. 
             ! --------------------------------------------------------------------
-            zlimfac = zprchld * xlimdia(ji,jj,jk) / ( zprmax + rtrn )
+            zlimfac = zprchld(ji,jj,jk) * xlimdia(ji,jj,jk) / ( zprmaxd(ji,jj,jk) + rtrn )
             zsizetmp = 1.0 + 1.3 * ( xsizerd - 1.0 ) * zlimfac**3/(0.3 + zlimfac**3)
             sizeda(ji,jj,jk) = min(xsizerd, max( sizeda(ji,jj,jk), zsizetmp ) )
 
@@ -265,37 +293,45 @@ CONTAINS
             zfecdm = xqfuncfecd(ji,jj,jk) + ( fecdm - xqfuncfecd(ji,jj,jk) ) * ( xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk) )
             zratio = 1.0 - MIN(1.0, tr(ji,jj,jk,jpdfe,Kbb) / ( tr(ji,jj,jk,jpdia,Kbb) * zfecdm + rtrn ) )
             zmax   = MAX( 0., MIN( 1.0, zratio**2/ (0.05**2+zratio**2) ) ) 
-            zprofed(ji,jj,jk) = zfecdm * zprmax * (1.0 - fr_i(ji,jj) )  &
+            zprofed(ji,jj,jk) = zfecdm * zprmaxd(ji,jj,jk) * (1.0 - fr_i(ji,jj) )  &
             &          * (1. + 0.8 * xdiatno3(ji,jj,jk) / ( rtrn + xdiatno3(ji,jj,jk)  &
             &          + xdiatnh4(ji,jj,jk) ) * (1. - xdiatfer(ji,jj,jk) ) )   &
             &          * xdiatfer(ji,jj,jk) * zmax * tr(ji,jj,jk,jpdia,Kbb) * rfact2
+         ENDIF
+      END_3D
 
-            ! Computation of the chlorophyll production terms
-            ! The parameterization is taken from Geider et al. (1997)
-
+      ! Computation of the chlorophyll production terms
+      ! The parameterization is taken from Geider et al. (1997)
+      ! -------------------------------------------------------
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+         IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
             !  production terms for nanophyto. ( chlorophyll )
-            znanotot = enanom(ji,jj,jk) / ( zmxl_chl + rtrn )
-            zprod    = rday * zprorcan(ji,jj,jk) * zprchln * xlimphy(ji,jj,jk)
+            znanotot = enanom(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+            zprod    = rday * zprorcan(ji,jj,jk) * zprchln(ji,jj,jk) * xlimphy(ji,jj,jk)
             zprochln = chlcmin * 12. * zprorcan (ji,jj,jk)
             zprochln = zprochln + (chlcnm - chlcmin) * 12. * zprod / &
-                                  & (  zpislopeadn * znanotot +rtrn)
+                                  & (  zpislopeadn(ji,jj,jk) * znanotot +rtrn)
 
             !  production terms for diatoms ( chlorophyll )
-            zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl + rtrn )
-            zprod    = rday * zprorcad(ji,jj,jk) * zprchld * xlimdia(ji,jj,jk)
+            zdiattot = ediatm(ji,jj,jk) / ( zmxl_chl(ji,jj,jk) + rtrn )
+            zprod    = rday * zprorcad(ji,jj,jk) * zprchld(ji,jj,jk) * xlimdia(ji,jj,jk)
             zprochld = chlcmin * 12. * zprorcad(ji,jj,jk)
             zprochld = zprochld + (chlcdm - chlcmin) * 12. * zprod / &
-                                  & ( zpislopeadd * zdiattot +rtrn )
+                                  & ( zpislopeadd(ji,jj,jk) * zdiattot +rtrn )
 
             !   Update the arrays TRA which contain the Chla sources and sinks
             tr(ji,jj,jk,jpnch,Krhs) = tr(ji,jj,jk,jpnch,Krhs) + zprochln * texcretn
             tr(ji,jj,jk,jpdch,Krhs) = tr(ji,jj,jk,jpdch,Krhs) + zprochld * texcretd
+         ENDIF
+      END_3D
 
-            !   Update the arrays TRA which contain the biological sources and sinks
+      !   Update the arrays TRA which contain the biological sources and sinks
+      DO_3D( 0, 0, 0, 0, 1, jpkm1)
+        IF( etot_ndcy(ji,jj,jk) > 1.E-3 ) THEN
            zpptot   = zprorcan(ji,jj,jk) + zprorcad(ji,jj,jk)
            zpnewtot = zpronewn(ji,jj,jk) + zpronewd(ji,jj,jk)
            zpregtot = zpptot - zpnewtot
-           zprodsil  = zprmax * zysopt(ji,jj,jk) * rfact2 * tr(ji,jj,jk,jpdia,Kbb)
+           zprodsil  = zprmaxd(ji,jj,jk) * zysopt(ji,jj,jk) * rfact2 * tr(ji,jj,jk,jpdia,Kbb)
            zproddoc  = excretd * zprorcad(ji,jj,jk) + excretn * zprorcan(ji,jj,jk)
            zprodfer  = texcretn * zprofen(ji,jj,jk) + texcretd * zprofed(ji,jj,jk)
            !
@@ -317,7 +353,7 @@ CONTAINS
            !
            tr(ji,jj,jk,jpdic,Krhs) = tr(ji,jj,jk,jpdic,Krhs) - zpptot
            tr(ji,jj,jk,jptal,Krhs) = tr(ji,jj,jk,jptal,Krhs) + rno3 * ( zpnewtot - zpregtot )
-         ENDIF
+        ENDIF
       END_3D
 
      ! Production and uptake of ligands by phytoplankton. This part is activated 
@@ -337,6 +373,7 @@ CONTAINS
          END_3D
      ENDIF
 
+
     ! Output of the diagnostics
     ! Total primary production per year
     IF( iom_use( "tintpp" ) .OR. ( ln_check_mass .AND. kt == nitend .AND. knt == nrdttrc )  )  THEN
@@ -345,6 +382,10 @@ CONTAINS
        tpp = glob_sum( 'p4zprod', zw3d )
        DEALLOCATE ( zw3d ) 
     ENDIF
+    
+    IF( lk_iomput .AND.  knt == nrdttrc ) THEN
+       zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
+       !
 
     IF( lk_iomput .AND.  knt == nrdttrc ) THEN
        zfact = 1.e+3 * rfact2r  !  conversion from mol/l/kt to  mol/m3/s
@@ -398,7 +439,7 @@ CONTAINS
        ENDIF
        !
        IF( iom_use ( "Mumax" ) ) THEN    ! Maximum growth rate
-         zw3d(A2D(0),1:jpkm1) = 0.65_wp * r1_rday * tgfunc(A2D(0),1:jpkm1)  * tmask(A2D(0),1:jpkm1) 
+         zw3d(A2D(0),1:jpkm1) = zprmaxn(A2D(0),1:jpkm1)  * tmask(A2D(0),1:jpkm1) 
          CALL iom_put( "Mumax", zw3d )  
        ENDIF
        !
@@ -413,12 +454,12 @@ CONTAINS
        ENDIF
        !
        IF( iom_use ( "LNlight" ) ) THEN    ! light limitation term for nano
-         zw3d(A2D(0),1:jpkm1) = zprbio(A2D(0),1:jpkm1) / ( 0.65_wp * r1_rday * tgfunc(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1) 
+         zw3d(A2D(0),1:jpkm1) = zprbio(A2D(0),1:jpkm1) / ( zprmaxn(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1) 
          CALL iom_put( "LNlight", zw3d )  
        ENDIF
        !
        IF( iom_use ( "LDlight" ) ) THEN    ! light limitation term for diatomes
-         zw3d(A2D(0),1:jpkm1) = zprdia(A2D(0),1:jpkm1) / ( 0.65_wp * r1_rday * tgfunc(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1) 
+         zw3d(A2D(0),1:jpkm1) = zprdia(A2D(0),1:jpkm1) / ( zprmaxd(A2D(0),1:jpkm1) + rtrn ) * tmask(A2D(0),1:jpkm1) 
          CALL iom_put( "LDlight", zw3d )  
        ENDIF
        !
@@ -440,6 +481,7 @@ CONTAINS
        IF( iom_use( "tintpp") )  CALL iom_put( "tintpp"  , tpp * zfact )  !  global total integrated primary production molC/s
        !
        DEALLOCATE( zw3d )
+
      ENDIF
 
      IF(sn_cfctl%l_prttrc)   THEN  ! print mean trends (used for debugging)
diff --git a/src/TOP/TRP/trcadv.F90 b/src/TOP/TRP/trcadv.F90
index 86776d8784ecf5cd192ec8e6917b8d0642f096ad..e98590e1693612ef16395d623624c477273d89d9 100644
--- a/src/TOP/TRP/trcadv.F90
+++ b/src/TOP/TRP/trcadv.F90
@@ -123,19 +123,19 @@ CONTAINS
          !
          IF( ln_wave .AND. ln_sdw )  THEN
             DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )                            ! eulerian transport + Stokes Drift
-               zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * ( zptu(ji,jj,jk) + usd(ji,jj,jk) )
-               zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * ( zptv(ji,jj,jk) + vsd(ji,jj,jk) )
+               zuu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * ( zptu(ji,jj,jk) + usd(ji,jj,jk) )
+               zvv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * ( zptv(ji,jj,jk) + vsd(ji,jj,jk) )
             END_3D
             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
-               zww(ji,jj,jk) = e1e2t(ji,jj)                     * ( zptw(ji,jj,jk) + wsd(ji,jj,jk) )
+               zww(ji,jj,jk) = e1e2t(ji,jj)                   * ( zptw(ji,jj,jk) + wsd(ji,jj,jk) )
             END_3D
          ELSE
             DO_3D( nn_hls, nn_hls-1, nn_hls, nn_hls-1, 1, jpkm1 )
-               zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * zptu(ji,jj,jk)           ! eulerian transport
-               zvv(ji,jj,jk) = e1v  (ji,jj) * e3v(ji,jj,jk,Kmm) * zptv(ji,jj,jk)
+               zuu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zptu(ji,jj,jk)           ! eulerian transport
+               zvv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zptv(ji,jj,jk)
             END_3D
             DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
-               zww(ji,jj,jk) = e1e2t(ji,jj)                     * zptw(ji,jj,jk)
+               zww(ji,jj,jk) = e1e2t(ji,jj)                   * zptw(ji,jj,jk)
             END_3D
          ENDIF
          !
diff --git a/src/TOP/TRP/trcsbc.F90 b/src/TOP/TRP/trcsbc.F90
index f817b677b011d92edd8a7515e009f0f22f2a0272..88fd617f1efbeeb03ff099a56348c5dd3f6c1c30 100644
--- a/src/TOP/TRP/trcsbc.F90
+++ b/src/TOP/TRP/trcsbc.F90
@@ -115,14 +115,14 @@ CONTAINS
       CASE ( -1 ) ! ! No tracers in sea ice ( trc_i = 0 )
          !
          DO jn = 1, jptra
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 0, 0, 0 )
                sbc_trc(ji,jj,jn) = 0._wp
             END_2D
          END DO
          !
          IF( ln_linssh ) THEN  !* linear free surface  
             DO jn = 1, jptra
-               DO_2D( 0, 0, 0, 1 )
+               DO_2D( 0, 0, 0, 0 )
                   sbc_trc(ji,jj,jn) = sbc_trc(ji,jj,jn) + r1_rho0 * emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) !==>> add concentration/dilution effect due to constant volume cell
                END_2D
             END DO
@@ -131,14 +131,14 @@ CONTAINS
       CASE ( 0 )  ! Same concentration in sea ice and in the ocean ( trc_i = ptr(...,Kmm)  )
          !
          DO jn = 1, jptra
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 0, 0, 0 )
                sbc_trc(ji,jj,jn) = - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
             END_2D
          END DO
          !
          IF( ln_linssh ) THEN  !* linear free surface  
             DO jn = 1, jptra
-               DO_2D( 0, 0, 0, 1 )
+               DO_2D( 0, 0, 0, 0 )
                   sbc_trc(ji,jj,jn) = sbc_trc(ji,jj,jn) + r1_rho0 * emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) !==>> add concentration/dilution effect due to constant volume cell
                END_2D
             END DO
@@ -147,21 +147,21 @@ CONTAINS
       CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice 
          !
          DO jn = 1, jptra
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 0, 0, 0 )
                sbc_trc(ji,jj,jn) = - fmmflx(ji,jj) * r1_rho0 * trc_i(ji,jj,jn)
             END_2D
          END DO
          !
          IF( ln_linssh ) THEN  !* linear free surface  
             DO jn = 1, jptra
-               DO_2D( 0, 0, 0, 1 )
+               DO_2D( 0, 0, 0, 0 )
                   sbc_trc(ji,jj,jn) = sbc_trc(ji,jj,jn) + r1_rho0 * emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) !==>> add concentration/dilution effect due to constant volume cell
                END_2D
             END DO
          ENDIF
          !
          DO jn = 1, jptra
-            DO_2D( 0, 0, 0, 1 )
+            DO_2D( 0, 0, 0, 0 )
                zse3t = rDt_trc / e3t(ji,jj,1,Kmm)
                zdtra = ptr(ji,jj,1,jn,Kmm) + sbc_trc(ji,jj,jn) * zse3t 
                IF( zdtra < 0. ) sbc_trc(ji,jj,jn) = MAX( zdtra, -ptr(ji,jj,1,jn,Kmm) / zse3t  ) ! avoid negative concentration that can occurs if trc_i > ptr 
@@ -176,7 +176,7 @@ CONTAINS
          !
          IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends
          !
-         DO_2D( 0, 0, 0, 1 )
+         DO_2D( 0, 0, 0, 0 )
             zse3t = zfact / e3t(ji,jj,1,Kmm)
             ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( sbc_trc_b(ji,jj,jn) + sbc_trc(ji,jj,jn) ) * zse3t
          END_2D
@@ -295,7 +295,7 @@ CONTAINS
             CASE ( 0 )  ! Same concentration in sea ice and in the ocean fmm contribution to concentration/dilution effect has to be removed
                !
                DO jn = 1, jptra
-                  DO_2D( 0, 0, 0, 1 )
+                  DO_2D( 0, 0, 0, 0 )
                      z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm)
                      ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + ( emp(ji,jj) - fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
                   END_2D
@@ -331,7 +331,7 @@ CONTAINS
             CASE ( 0 )  ! Same concentration in sea ice and in the ocean : correct concentration/dilution effect due to "freezing - melting"
                !
                DO jn = 1, jptra
-                  DO_2D( 0, 0, 0, 1 )
+                  DO_2D( 0, 0, 0, 0 )
                      z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm)
                      ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
                   END_2D