Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/C14/trcatm_c14.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/C14/trcatm_c14.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/C14/trcatm_c14.F90	(revision 13571)
@@ -120,5 +120,5 @@
             IF( ierr3 /= 0 )   CALL ctl_stop( 'STOP', 'trc_atm_c14_ini: unable to allocate fareaz' )
       !
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )                 ! from C14b package
               IF( gphit(ji,jj) >= yn40 ) THEN
                  fareaz(ji,jj,1) = 0.
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/CFC/trcsms_cfc.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/CFC/trcsms_cfc.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/CFC/trcsms_cfc.F90	(revision 13571)
@@ -126,6 +126,6 @@
          
          !                                                         !------------!
-         DO_2D( 1, 1, 1, 1 )
- 
+         DO_2D( 1, 1, 1, 1 )                                       !  i-j loop  !
+            !                                                      !------------!
             ! space interpolation
             zpp_cfc  =       xphem(ji,jj)   * zpatm(1,jl)   &
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P2Z/p2zopt.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P2Z/p2zopt.F90	(revision 13571)
@@ -95,5 +95,5 @@
       !                                          ! Photosynthetically Available Radiation (PAR)
       zcoef = 12 * redf / rcchl / rpig           ! --------------------------------------
-      DO_3D( 1, 1, 1, 1, 2, jpk )
+      DO_3D( 1, 1, 1, 1, 2, jpk )                     ! local par at w-levels
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk-1,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
@@ -102,5 +102,5 @@
          zparg(ji,jj,jk) = zparg(ji,jj,jk-1) * EXP( -zkg * e3t(ji,jj,jk-1,Kmm) )
       END_3D
-      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )                   ! mean par at t-levels
          zpig = LOG(  MAX( TINY(0.), tr(ji,jj,jk,jpphy,Kmm) ) * zcoef  )
          zkr  = xkr0 + xkrp * EXP( xlr * zpig )
@@ -114,5 +114,5 @@
       !                                          ! --------------
       neln(:,:) = 1                                   ! euphotic layer level
-      DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+      DO_3D( 1, 1, 1, 1, 1, jpkm1 )                   ! (i.e. 1rst T-level strictly below EL bottom)
         IF( etot(ji,jj,jk) >= zpar100(ji,jj) )   neln(ji,jj) = jk + 1 
       END_3D
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zfechem.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zfechem.F90	(revision 13571)
@@ -118,8 +118,4 @@
          !
          zfeequi = zFe3(ji,jj,jk) * 1E-9
-         zhplus  = max( rtrn, hi(ji,jj,jk) )
-         fe3sol  = fesol(ji,jj,jk,1) * ( zhplus**3 + fesol(ji,jj,jk,2) * zhplus**2  &
-            &         + fesol(ji,jj,jk,3) * zhplus + fesol(ji,jj,jk,4)     &
-            &         + fesol(ji,jj,jk,5) / zhplus )
          zfecoll = 0.5 * zFeL1(ji,jj,jk) * 1E-9
          ! precipitation of Fe3+, creation of nanoparticles
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zlim.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zlim.F90	(revision 13571)
@@ -161,5 +161,5 @@
          zlim1    = xdiatno3(ji,jj,jk) + xdiatnh4(ji,jj,jk)
          zlim2    = tr(ji,jj,jk,jppo4,Kbb) / ( tr(ji,jj,jk,jppo4,Kbb) + zconc1dnh4  )
-         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) )
+         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) + rtrn )
          zratio   = tr(ji,jj,jk,jpdfe,Kbb) * z1_trbdia
          zironmin = xcoef1 * tr(ji,jj,jk,jpdch,Kbb) * z1_trbdia + xcoef2 * zlim1 + xcoef3 * xdiatno3(ji,jj,jk)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13571)
@@ -37,8 +37,4 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   par_varsw      ! PAR fraction of shortwave
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ekb, ekg, ekr  ! wavelength (Red-Green-Blue)
-
-   INTEGER  ::   nksrp   ! levels below which the light cannot penetrate ( depth larger than 391 m)
-
-   REAL(wp), DIMENSION(3,61) ::   xkrgb   ! tabulated attenuation coefficients for RGB absorption
    
    !! * Substitutions
@@ -94,7 +90,7 @@
          irgb = NINT( 41 + 20.* LOG10( zchl ) + rtrn )
          !                                                         
-         ekb(ji,jj,jk) = xkrgb(1,irgb) * e3t(ji,jj,jk,Kmm)
-         ekg(ji,jj,jk) = xkrgb(2,irgb) * e3t(ji,jj,jk,Kmm)
-         ekr(ji,jj,jk) = xkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
+         ekb(ji,jj,jk) = rkrgb(1,irgb) * e3t(ji,jj,jk,Kmm)
+         ekg(ji,jj,jk) = rkrgb(2,irgb) * e3t(ji,jj,jk,Kmm)
+         ekr(ji,jj,jk) = rkrgb(3,irgb) * e3t(ji,jj,jk,Kmm)
       END_3D
       !                                        !* Photosynthetically Available Radiation (PAR)
@@ -106,5 +102,5 @@
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100 ) 
          !
-         DO jk = 1, nksrp      
+         DO jk = 1, nksr      
             etot_ndcy(:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
             enano    (:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
@@ -112,5 +108,5 @@
          END DO
          IF( ln_p5z ) THEN
-            DO jk = 1, nksrp      
+            DO jk = 1, nksr      
               epico  (:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
             END DO
@@ -121,5 +117,5 @@
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3 ) 
          !
-         DO jk = 1, nksrp      
+         DO jk = 1, nksr      
             etot(:,:,jk) =  ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk)
          END DO
@@ -131,5 +127,5 @@
          CALL p4z_opt_par( kt, Kmm, zqsr_corr, ze1, ze2, ze3, pqsr100 = zqsr100  ) 
          !
-         DO jk = 1, nksrp      
+         DO jk = 1, nksr      
             etot (:,:,jk) =        ze1(:,:,jk) +        ze2(:,:,jk) +       ze3(:,:,jk)
             enano(:,:,jk) =  1.85 * ze1(:,:,jk) + 0.69 * ze2(:,:,jk) + 0.46 * ze3(:,:,jk)
@@ -137,5 +133,5 @@
          END DO
          IF( ln_p5z ) THEN
-            DO jk = 1, nksrp      
+            DO jk = 1, nksr      
               epico(:,:,jk) =  1.94 * ze1(:,:,jk) + 0.66 * ze2(:,:,jk) + 0.4 * ze3(:,:,jk)
             END DO
@@ -150,5 +146,5 @@
          !
          etot3(:,:,1) =  qsr(:,:) * tmask(:,:,1)
-         DO jk = 2, nksrp + 1
+         DO jk = 2, nksr + 1
             etot3(:,:,jk) =  ( ze0(:,:,jk) + ze1(:,:,jk) + ze2(:,:,jk) + ze3(:,:,jk) ) * tmask(:,:,jk)
          END DO
@@ -160,5 +156,5 @@
       heup_01(:,:) = gdepw(:,:,2,Kmm)
 
-      DO_3D( 1, 1, 1, 1, 2, nksrp )
+      DO_3D( 1, 1, 1, 1, 2, nksr )
         IF( etot_ndcy(ji,jj,jk) * tmask(ji,jj,jk) >=  zqsr100(ji,jj) )  THEN
            neln(ji,jj) = jk+1                    ! Euphotic level : 1rst T-level strictly below Euphotic layer
@@ -178,5 +174,5 @@
       zetmp2 (:,:)   = 0.e0
 
-      DO_3D( 1, 1, 1, 1, 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             zetmp1 (ji,jj) = zetmp1 (ji,jj) + etot     (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! remineralisation
@@ -189,5 +185,5 @@
       zpar(:,:,:) = etot_ndcy(:,:,:)  ! diagnostic : PAR with no diurnal cycle 
       !
-      DO_3D( 1, 1, 1, 1, 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
@@ -201,5 +197,5 @@
       zetmp4 (:,:)   = 0.e0
       !
-      DO_3D( 1, 1, 1, 1, 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
             zetmp3 (ji,jj) = zetmp3 (ji,jj) + enano    (ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
@@ -211,5 +207,5 @@
       ediatm(:,:,:) = ediat(:,:,:)
       !
-      DO_3D( 1, 1, 1, 1, 1, nksrp )
+      DO_3D( 1, 1, 1, 1, 1, nksr )
          IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
             z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
@@ -221,5 +217,5 @@
       IF( ln_p5z ) THEN
          ALLOCATE( zetmp5(jpi,jpj) )  ;   zetmp5 (:,:) = 0.e0
-         DO_3D( 1, 1, 1, 1, 1, nksrp )
+         DO_3D( 1, 1, 1, 1, 1, nksr )
             IF( gdepw(ji,jj,jk+1,Kmm) <= MIN(hmld(ji,jj), heup_01(ji,jj)) ) THEN
                zetmp5(ji,jj)  = zetmp5 (ji,jj) + epico(ji,jj,jk) * e3t(ji,jj,jk,Kmm) ! production
@@ -229,5 +225,5 @@
          epicom(:,:,:) = epico(:,:,:)
          !
-         DO_3D( 1, 1, 1, 1, 1, nksrp )
+         DO_3D( 1, 1, 1, 1, 1, nksr )
             IF( gdepw(ji,jj,jk+1,Kmm) <= hmld(ji,jj) ) THEN
                z1_dep = 1. / ( zdepmoy(ji,jj) + rtrn )
@@ -283,5 +279,5 @@
          pe3(:,:,1) = zqsr(:,:)
          !
-         DO jk = 2, nksrp + 1
+         DO jk = 2, nksr + 1
             DO jj = 1, jpj
                DO ji = 1, jpi
@@ -302,5 +298,5 @@
         pe3(:,:,1) = zqsr(:,:) * EXP( -0.5 * ekr(:,:,1) )
         !
-        DO_3D( 1, 1, 1, 1, 2, nksrp )
+        DO_3D( 1, 1, 1, 1, 2, nksr )
            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) ) )
@@ -400,9 +396,4 @@
          ntimes_par = iom_getszuld( numpar )   ! get number of record in file
       ENDIF
-      !
-      CALL trc_oce_rgb( xkrgb )                  ! tabulated attenuation coefficients
-      nksrp = trc_oce_ext_lev( r_si2, 0.33e2_wp )     ! max level of light extinction (Blue Chl=0.01)
-      !
-      IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksrp, ' ref depth = ', gdepw_1d(nksrp+1), ' m'
       !
                          ekr      (:,:,:) = 0._wp
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zsed.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zsed.F90	(revision 13571)
@@ -313,5 +313,5 @@
       ENDIF
       !
-      IF(sn_cfctl%l_prttrc) THEN  ! print mean trends (USEd for debugging)
+      IF(sn_cfctl%l_prttrc) THEN  ! print mean trneds (USEd for debugging)
          WRITE(charout, fmt="('sed ')")
          CALL prt_ctl_info( charout, cdcomp = 'top' )
@@ -366,4 +366,6 @@
       lk_sed = ln_sediment .AND. ln_sed_2way 
       !
+      nitrpot(:,:,jpk) = 0._wp   ! define last level for iom_put
+      !
    END SUBROUTINE p4z_sed_init
 
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13571)
@@ -69,5 +69,5 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:    ) :: zw2d
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:  ) :: zw3d
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrdt   ! 4D workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jp_pisces) :: ztrbbio
 
       !!---------------------------------------------------------------------
@@ -93,11 +93,4 @@
       rfact = rDt_trc
       !
-      ! trends computation initialisation
-      IF( l_trdtrc )  THEN
-         ALLOCATE( ztrdt(jpi,jpj,jpk,jp_pisces) )  !* store now fields before applying the Asselin filter
-         ztrdt(:,:,:,:)  = tr(:,:,:,:,Kmm)
-      ENDIF
-      !
-
       IF( ( ln_top_euler .AND. kt == nittrc000 )  .OR. ( .NOT.ln_top_euler .AND. kt <= nittrc000 + 1 ) ) THEN
          rfactr  = 1. / rfact
@@ -117,4 +110,9 @@
          END DO
       ENDIF
+
+      DO jn = jp_pcs0, jp_pcs1              !   Store the tracer concentrations before entering PISCES
+         ztrbbio(:,:,:,jn) = tr(:,:,:,jn,Kbb)
+      END DO
+
       !
       IF( ll_bc )    CALL p4z_bc( kt, Kbb, Kmm, Krhs )   ! external sources of nutrients 
@@ -198,32 +196,26 @@
          END DO
          !
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               tr(:,:,:,jn,Kmm) = tr(:,:,:,jn,Kbb)
-            END DO
-         ENDIF
+      END DO
+      !
+#endif
+      !
+      IF( ln_sediment ) THEN 
+         !
+         CALL sed_model( kt, Kbb, Kmm, Krhs )     !  Main program of Sediment model
+         !
+      ENDIF
+      !
+      DO jn = jp_pcs0, jp_pcs1
+         tr(:,:,:,jn,Krhs) = ( tr(:,:,:,jn,Kbb) - ztrbbio(:,:,:,jn) ) * rfactr
+         tr(:,:,:,jn,Kbb ) = ztrbbio(:,:,:,jn)
+         ztrbbio(:,:,:,jn) = 0._wp
       END DO
       !
       IF( l_trdtrc ) THEN
          DO jn = jp_pcs0, jp_pcs1
-           ztrdt(:,:,:,jn) = ( tr(:,:,:,jn,Kbb) - ztrdt(:,:,:,jn) ) * rfactr 
            CALL trd_trc( tr(:,:,:,jn,Krhs), jn, jptra_sms, kt, Kmm )   ! save trends
          END DO
-         DEALLOCATE( ztrdt ) 
       END IF
-#endif
-      !
-      IF( ln_sediment ) THEN 
-         !
-         CALL sed_model( kt, Kbb, Kmm, Krhs )     !  Main program of Sediment model
-         !
-         IF( ln_top_euler ) THEN
-            DO jn = jp_pcs0, jp_pcs1
-               tr(:,:,:,jn,Kmm) = tr(:,:,:,jn,Kbb)
-            END DO
-         ENDIF
-         !
-      ENDIF
-      !
+      !  
       IF( lrst_trc )  CALL p4z_rst( kt, Kbb, Kmm,  'WRITE' )           !* Write PISCES informations in restart file 
       !
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p5zlim.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/PISCES/P4Z/p5zlim.F90	(revision 13571)
@@ -306,5 +306,5 @@
          &          / (xqndmax(ji,jj,jk) - 2. * xqndmin(ji,jj,jk) ) )   &
          &          * xqndmax(ji,jj,jk) / (zration + rtrn)
-         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) )
+         zlim3    = tr(ji,jj,jk,jpsil,Kbb) / ( tr(ji,jj,jk,jpsil,Kbb) + xksi(ji,jj) + rtrn )
          zlim4    = MAX( 0., ( zratiof - zqfemd ) / qfdopt )
          xlimdfe(ji,jj,jk) = MIN( 1., zlim4 )
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trcrad.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trcrad.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trcrad.F90	(revision 13571)
@@ -168,5 +168,5 @@
               IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,itime)                       ! save input tr(:,:,:,:,Kbb) for trend computation           
               !
-              DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+              DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
                  IF( ztrneg(ji,jj,jn) /= 0. ) THEN                                 ! if negative values over the 3x3 box
                     !
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trdmxl_trc.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trdmxl_trc.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/TRP/trdmxl_trc.F90	(revision 13571)
@@ -148,5 +148,5 @@
          ! ... Weights for vertical averaging
          wkx_trc(:,:,:) = 0.e0
-         DO_3D( 1, 1, 1, 1, 1, jpktrd_trc )
+         DO_3D( 1, 1, 1, 1, 1, jpktrd_trc )                       ! initialize wkx_trc with vertical scale factor in mixed-layer
             IF( jk - nmld_trc(ji,jj) < 0 )   wkx_trc(ji,jj,jk) = e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk)
          END_3D
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/oce_trc.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/oce_trc.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/oce_trc.F90	(revision 13571)
@@ -85,4 +85,6 @@
    USE traqsr  , ONLY :   rn_abs     =>    rn_abs     !: fraction absorbed in the very near surface
    USE traqsr  , ONLY :   rn_si0     =>    rn_si0     !: very near surface depth of extinction
+   USE traqsr  , ONLY :   nksr       =>    nksr       !: levels below which the light cannot penetrate (depth larger than 391 m)
+   USE traqsr  , ONLY :   rkrgb      =>    rkrgb      !: tabulated attenuation coefficients for RGB absorption
    USE traqsr  , ONLY :   ln_qsr_bio =>    ln_qsr_bio !: flag to use or not the biological fluxes for light
    USE sbcrnf  , ONLY :   rnfmsk     =>    rnfmsk     !: mixed adv scheme in runoffs vicinity (hori.) 
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trc.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trc.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trc.F90	(revision 13571)
@@ -21,5 +21,5 @@
    INTEGER, PUBLIC ::   numonr     = -1   !: reference passive tracer namelist output output.namelist.top
    INTEGER, PUBLIC ::   numstr            !: tracer statistics
-   INTEGER, PUBLIC ::   numrtr            !: trc restart (read )
+   INTEGER, PUBLIC ::   numrtr     = -1   !: trc restart (read )
    INTEGER, PUBLIC ::   numrtw            !: trc restart ( write )
    CHARACTER(:), ALLOCATABLE, PUBLIC ::   numnat_ref   !: character buffer for reference passive tracer namelist_top_ref
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcbdy.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcbdy.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcbdy.F90	(revision 13571)
@@ -49,5 +49,4 @@
       INTEGER                           :: ib_bdy ,ir, jn ,igrd ! Loop indices
       REAL(wp), POINTER, DIMENSION(:,:) ::  ztrc
-      REAL(wp), POINTER                 ::  zfac
       LOGICAL                           :: llrim0               ! indicate if rim 0 is treated
       LOGICAL, DIMENSION(4)             :: llsend1, llrecv1     ! indicate how communications are to be carried out
@@ -61,30 +60,38 @@
          IF( ir == 0 ) THEN   ;   llrim0 = .TRUE.
          ELSE                 ;   llrim0 = .FALSE.
-         END IF
+         ENDIF
          DO ib_bdy=1, nb_bdy
+            !
             DO jn = 1, jptra
                !
-               ztrc => trcdta_bdy(jn,ib_bdy)%trc 
-               zfac => trcdta_bdy(jn,ib_bdy)%rn_fac
+               IF( ASSOCIATED(trcdta_bdy(jn,ib_bdy)%trc) .AND. trcdta_bdy(jn,ib_bdy)%cn_obc /= 'neumann' ) THEN
+                  IF( .NOT. ASSOCIATED(ztrc) )   ALLOCATE( ztrc(idx_bdy(ib_bdy)%nblen(igrd),jpk) )
+                  ztrc(:,:) = trcdta_bdy(jn,ib_bdy)%trc(:,:) * trcdta_bdy(jn,ib_bdy)%rn_fac
+               ENDIF
                !
-               SELECT CASE( TRIM(trcdta_bdy(jn,ib_bdy)%cn_obc) )
+               SELECT CASE( trcdta_bdy(jn,ib_bdy)%cn_obc )
                CASE('none'        )   ;   CYCLE
                CASE('frs'         )   ! treat the whole boundary at once
-                  IF( ir == 0 ) CALL bdy_frs( idx_bdy(ib_bdy),                tr(:,:,:,jn,Krhs), ztrc*zfac )
+                  IF( ir == 0 )           CALL bdy_frs( idx_bdy(ib_bdy),                   tr(:,:,:,jn,Krhs), ztrc )
                CASE('specified'   )   ! treat the whole rim      at once
-                  IF( ir == 0 ) CALL bdy_spe( idx_bdy(ib_bdy),                tr(:,:,:,jn,Krhs), ztrc*zfac )
-               CASE('neumann'     )   ;   CALL bdy_nmn( idx_bdy(ib_bdy), igrd         , tr(:,:,:,jn,Krhs) )   ! tra masked
-               CASE('orlanski'    )   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc*zfac, ll_npo=.false. )
-               CASE('orlanski_npo')   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc*zfac, ll_npo=.true. )
+                  IF( ir == 0 )           CALL bdy_spe( idx_bdy(ib_bdy),                   tr(:,:,:,jn,Krhs), ztrc )
+               CASE('neumann'     )   ;   CALL bdy_nmn( idx_bdy(ib_bdy), igrd            , tr(:,:,:,jn,Krhs),       llrim0 )   ! tra masked
+               CASE('orlanski'    )   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc, llrim0,   &
+                  &                                     ll_npo=.FALSE. )
+               CASE('orlanski_npo')   ;   CALL bdy_orl( idx_bdy(ib_bdy), tr(:,:,:,jn,Kbb), tr(:,:,:,jn,Krhs), ztrc, llrim0,   &
+                  &                                     ll_npo=.TRUE.  )
                CASE DEFAULT           ;   CALL ctl_stop( 'trc_bdy : unrecognised option for open boundaries for passive tracers' )
                END SELECT
                !
             END DO
+            !
+            IF( ASSOCIATED(ztrc) )   DEALLOCATE(ztrc)
+            !
          END DO
          !
          IF( nn_hls > 1 .AND. ir == 1 ) CYCLE   ! at least 2 halos will be corrected -> no need to correct rim 1 before rim 0
-         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   END IF
+         IF( nn_hls == 1 ) THEN   ;   llsend1(:) = .false.   ;   llrecv1(:) = .false.   ;   ENDIF
          DO ib_bdy=1, nb_bdy
-            SELECT CASE( TRIM(cn_tra(ib_bdy)) )
+            SELECT CASE( cn_tra(ib_bdy) )
             CASE('neumann')
                llsend1(:) = llsend1(:) .OR. lsend_bdyint(ib_bdy,1,:,ir)   ! possibly every direction, T points
@@ -97,5 +104,5 @@
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
             CALL lbc_lnk( 'trcbdy', tr(:,:,:,:,Krhs), 'T',  1.0_wp, kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
-         END IF
+         ENDIF
          !
       END DO   ! ir
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcdta.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcdta.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcdta.F90	(revision 13571)
@@ -199,5 +199,5 @@
                WRITE(numout,*) 'trc_dta: interpolates passive tracer data onto the s- or mixed s-z-coordinate mesh'
             ENDIF
-            DO_2D( 1, 1, 1, 1 )
+            DO_2D( 1, 1, 1, 1 )                 ! vertical interpolation of T & S
                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
                   zl = gdept(ji,jj,jk,Kmm)
Index: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcrst.F90
===================================================================
--- NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcrst.F90	(revision 13302)
+++ NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/TOP/trcrst.F90	(revision 13571)
@@ -237,5 +237,5 @@
                ! calculate start time in hours and minutes
                zdayfrac=adatrj-INT(adatrj)
-               ksecs = NINT(zdayfrac*86400)            ! Nearest second to catch rounding errors in adatrj              
+               ksecs = NINT(zdayfrac*86400)            ! Nearest second to catch rounding errors in adatrj
                ihour = INT(ksecs/3600)
                iminute = ksecs/60-ihour*60
@@ -258,7 +258,8 @@
                adatrj = INT(adatrj)                    ! adatrj set to integer as nn_time0 updated            
              ELSE
+               ndt05 = NINT( 0.5 * rn_Dt  )   !  --- WARNING --- not defined yet are we did not go through day_init
                ! parameters corresponding to nit000 - 1 (as we start the step
                ! loop with a call to day)
-               ndastp = ndate0 - 1       ! ndate0 read in the namelist in dom_nam
+               ndastp = ndate0        ! ndate0 read in the namelist in dom_nam
                nhour   =   nn_time0 / 100
                nminute = ( nn_time0 - nhour * 100 )
