Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/NST/agrif_top_sponge.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/NST/agrif_top_sponge.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/NST/agrif_top_sponge.F90	(revision 13373)
@@ -67,14 +67,14 @@
       !
       INTEGER  ::   ji, jj, jk, jn   ! dummy loop indices
-      REAL(wp) ::   zabe1, zabe2, ztrelax
-      REAL(wp), DIMENSION(i1:i2,j1:j2)               ::   ztu, ztv
-      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,1:jptra) ::   trbdiff
+      INTEGER  ::   iku, ikv
+      REAL(wp) ::   ztra, zabe1, zabe2, zbtr
+      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk)  ::   ztu, ztv
+      REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::   trbdiff
       ! vertical interpolation:
-      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,1:jptra) ::tabres_child
+      REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
       REAL(wp), DIMENSION(k1:k2,1:jptra) :: tabin
       REAL(wp), DIMENSION(k1:k2) :: h_in
       REAL(wp), DIMENSION(1:jpk) :: h_out
       INTEGER :: N_in, N_out
-      REAL(wp) :: h_diff
       !!----------------------------------------------------------------------
       !
@@ -94,5 +94,5 @@
             DO jj=j1,j2
                DO ji=i1,i2
-                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 
+                  tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)
                END DO
             END DO
@@ -128,7 +128,7 @@
                DO jk=1,jpkm1
 # if defined key_vertical
-                  trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra)
+                  trbdiff(ji,jj,jk,1:jptra) = ( tr(ji,jj,jk,1:jptra,Kbb_a) - tabres_child(ji,jj,jk,1:jptra) ) * tmask(ji,jj,jk)
 # else
-                  trbdiff(ji,jj,jk,1:jptra) = tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra)
+                  trbdiff(ji,jj,jk,1:jptra) = ( tr(ji,jj,jk,1:jptra,Kbb_a) - tabres(ji,jj,jk,1:jptra)) * tmask(ji,jj,jk)
 # endif
                ENDDO
@@ -136,32 +136,48 @@
          ENDDO
 
-         !* set relaxation time scale
-         IF( l_1st_euler .AND. lk_agrif_fstep ) THEN   ;   ztrelax =   rn_trelax_tra  / (        rn_Dt )
-         ELSE                                          ;   ztrelax =   rn_trelax_tra  / (2._wp * rn_Dt )
-         ENDIF
-
          DO jn = 1, jptra
             DO jk = 1, jpkm1
-               DO jj = j1,j2-1
+               ztu(i1:i2,j1:j2,jk) = 0._wp
+               DO jj = j1,j2
                   DO ji = i1,i2-1
-                     zabe1 = rn_sponge_tra * fspu(ji,jj) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a) * umask(ji,jj,jk)
-                     zabe2 = rn_sponge_tra * fspv(ji,jj) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk)
-                     ztu(ji,jj) = zabe1 * ( trbdiff(ji+1,jj  ,jk,jn) - trbdiff(ji,jj,jk,jn) )
-                     ztv(ji,jj) = zabe2 * ( trbdiff(ji  ,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) )
+                     zabe1 = rn_sponge_tra * r1_Dt * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
+                     ztu(ji,jj,jk) = zabe1 * ( trbdiff(ji+1,jj,jk,jn) - trbdiff(ji,jj,jk,jn) )
+                  END DO
+               END DO
+               ztv(i1:i2,j1:j2,jk) = 0._wp
+               DO ji = i1,i2
+                  DO jj = j1,j2-1
+                     zabe2 = rn_sponge_tra * r1_Dt * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
+                     ztv(ji,jj,jk) = zabe2 * ( trbdiff(ji,jj+1,jk,jn) - trbdiff(ji,jj,jk,jn) )
                   END DO
                END DO
                !
+               IF( ln_zps ) THEN      ! set gradient at partial step level
+                  DO jj = j1,j2
+                     DO ji = i1,i2
+                        ! last level
+                        iku = mbku(ji,jj)
+                        ikv = mbkv(ji,jj)
+                        IF( iku == jk )   ztu(ji,jj,jk) = 0._wp
+                        IF( ikv == jk )   ztv(ji,jj,jk) = 0._wp
+                     END DO
+                  END DO
+               ENDIF
+            END DO
+            !
+            DO jk = 1, jpkm1
                DO jj = j1+1,j2-1
                   DO ji = i1+1,i2-1
-                     IF( .NOT. tabspongedone_trn(ji,jj) ) THEN 
-                        tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + (  ztu(ji,jj) - ztu(ji-1,jj  )     &
-                           &                                   + ztv(ji,jj) - ztv(ji  ,jj-1)  )  &
-                           &                                * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)  &
-                           &                                - ztrelax * fspt(ji,jj) * trbdiff(ji,jj,jk,jn)
+                     IF( .NOT. tabspongedone_trn(ji,jj) ) THEN
+                        zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
+                                                ! horizontal diffusive trends
+                        ztra = zbtr * (  ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk)  ) &
+                             &  - rn_trelax_tra * r1_Dt * fspt(ji,jj) * trbdiff(ji,jj,jk,jn)
+                        ! add it to the general tracer trends
+                        tr(ji,jj,jk,jn,Krhs_a) = tr(ji,jj,jk,jn,Krhs_a) + ztra
                      ENDIF
                   END DO
                END DO
             END DO
-            !
          END DO
          !
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/NST/agrif_user.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/NST/agrif_user.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/NST/agrif_user.F90	(revision 13373)
@@ -749,5 +749,5 @@
       CALL Agrif_Bc_variable(trn_sponge_id,calledweight=1.,procname=interptrn_sponge)
       ! reset tsa to zero
-      tra(:,:,:,:) = 0._wp
+      tr(:,:,:,:,Krhs_a) = 0._wp
 
       ! 3. Some controls
@@ -808,6 +808,6 @@
 # else
 ! LAURENT: STRANGE why (3,3) here ?
-      CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,jptra/),trn_id)
-      CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,jptra/),trn_sponge_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,jptra/),trn_id)
+      CALL agrif_declare_variable((/2,2,0,0/),(/ind2,ind3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/jpi,jpj,jpk,jptra/),trn_sponge_id)
 # endif
 
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/OCE/IOM/iom_def.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/OCE/IOM/iom_def.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/OCE/IOM/iom_def.F90	(revision 13373)
@@ -25,5 +25,5 @@
    INTEGER, PARAMETER, PUBLIC ::   jp_i1    = 204      !: write INTEGER(1)
 
-   INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 100  !: maximum number of simultaneously opened file
+   INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 200  !: maximum number of simultaneously opened file
    INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 1200 !: maximum number of variables in one file
    INTEGER, PARAMETER, PUBLIC ::   jpmax_dims   =  4   !: maximum number of dimensions for one variable
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zche.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zche.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zche.F90	(revision 13373)
@@ -188,259 +188,238 @@
       ! CHEMICAL CONSTANTS - SURFACE LAYER
       ! ----------------------------------
-!CDIR NOVERRCHK
-      DO jj = 1, jpj
-!CDIR NOVERRCHK
-         DO ji = 1, jpi
-            !                             ! SET ABSOLUTE TEMPERATURE
-            ztkel = tempis(ji,jj,1) + 273.15
-            zt    = ztkel * 0.01
-            zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
-            !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
-            !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
-            zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
-            &       + 0.0047036e-4*ztkel**2)
-            chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm)
-            chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
-            chemc(ji,jj,3) = 57.7 - 0.118*ztkel
-            !
-         END DO
-      END DO
+      DO_2D( 1, 1, 1, 1 )
+         !                             ! SET ABSOLUTE TEMPERATURE
+         ztkel = tempis(ji,jj,1) + 273.15
+         zt    = ztkel * 0.01
+         zsal  = salinprac(ji,jj,1) + ( 1.- tmask(ji,jj,1) ) * 35.
+         !                             ! LN(K0) OF SOLUBILITY OF CO2 (EQ. 12, WEISS, 1980)
+         !                             !     AND FOR THE ATMOSPHERE FOR NON IDEAL GAS
+         zcek1 = 9345.17/ztkel - 60.2409 + 23.3585 * LOG(zt) + zsal*(0.023517 - 0.00023656*ztkel    &
+         &       + 0.0047036e-4*ztkel**2)
+         chemc(ji,jj,1) = EXP( zcek1 ) * 1E-6 * rhop(ji,jj,1) / 1000. ! mol/(L atm)
+         chemc(ji,jj,2) = -1636.75 + 12.0408*ztkel - 0.0327957*ztkel**2 + 0.0000316528*ztkel**3
+         chemc(ji,jj,3) = 57.7 - 0.118*ztkel
+      END_2D
 
       ! OXYGEN SOLUBILITY - DEEP OCEAN
       ! -------------------------------
-!CDIR NOVERRCHK
-      DO jk = 1, jpk
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-              ztkel = tempis(ji,jj,jk) + 273.15
-              zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
-              zsal2 = zsal * zsal
-              ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
-              ztgg2 = ztgg  * ztgg
-              ztgg3 = ztgg2 * ztgg
-              ztgg4 = ztgg3 * ztgg
-              ztgg5 = ztgg4 * ztgg
-
-              zoxy  = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3    &
-              &       + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3   &
-              &       - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 )   &
-              &       - 3.11680e-7 * zsal2
-              chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox     ! mol/(L atm)
-            END DO
-          END DO
-        END DO
+      DO_3D( 1, 1, 1, 1, 1, jpk )
+         ztkel = tempis(ji,jj,jk) + 273.15
+         zsal  = salinprac(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 35.
+         zsal2 = zsal * zsal
+         ztgg  = LOG( ( 298.15 - tempis(ji,jj,jk) ) / ztkel )  ! Set the GORDON & GARCIA scaled temperature
+         ztgg2 = ztgg  * ztgg
+         ztgg3 = ztgg2 * ztgg
+         ztgg4 = ztgg3 * ztgg
+         ztgg5 = ztgg4 * ztgg
+
+         zoxy  = 2.00856 + 3.22400 * ztgg + 3.99063 * ztgg2 + 4.80299 * ztgg3    &
+         &       + 9.78188e-1 * ztgg4 + 1.71069 * ztgg5 + zsal * ( -6.24097e-3   &
+         &       - 6.93498e-3 * ztgg - 6.90358e-3 * ztgg2 - 4.29155e-3 * ztgg3 )   &
+         &       - 3.11680e-7 * zsal2
+         chemo2(ji,jj,jk) = ( EXP( zoxy ) * o2atm ) * oxyco * atcox     ! mol/(L atm)
+      END_3D
 
       ! CHEMICAL CONSTANTS - DEEP OCEAN
       ! -------------------------------
-!CDIR NOVERRCHK
-      DO jk = 1, jpk
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-
-               ! SET PRESSION ACCORDING TO SAUNDER (1980)
-               zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
-               zc1 = 5.92E-3 + zplat**2 * 5.25E-3
-               zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6
-               zpres = zpres / 10.0
-
-               ! SET ABSOLUTE TEMPERATURE
-               ztkel   = tempis(ji,jj,jk) + 273.15
-               zsal    = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
-               zsqrt  = SQRT( zsal )
-               zsal15  = zsqrt * zsal
-               zlogt  = LOG( ztkel )
-               ztr    = 1. / ztkel
-               zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
-               zis2   = zis * zis
-               zisqrt = SQRT( zis )
-               ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
-
-               ! CHLORINITY (WOOSTER ET AL., 1969)
-               zcl     = zsal / 1.80655
-
-               ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
-               zst     = 0.14 * zcl /96.062
-
-               ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
-               zft     = 0.000067 * zcl /18.9984
-
-               ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
-               zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         &
-               &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
-               &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    &
-               &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
-               &         + LOG(1.0 - 0.001005 * zsal))
-
-               ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
-               zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   &
-               &         + LOG(1.0d0 - 0.001005d0*zsal)            &
-               &         + LOG(1.0d0 + zst/zcks))
-
-               ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
-               zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        &
-               &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
-               &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
-               &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      & 
-               &      * zlogt + 0.053105*zsqrt*ztkel
-
-               ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
-               ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
-               zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
-                  - 0.011555*zsal + 0.0001152*zsal*zsal)
-               zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
-                  - 0.01781*zsal + 0.0001122*zsal*zsal)
-
-               ! PKW (H2O) (MILLERO, 1995) from composite data
-               zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
-                         - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
-
-               ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
-              zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
-              &          + (-106.736*ztr + 0.69171) * zsqrt       &
-              &          + (-0.65643*ztr - 0.01844) * zsal
-
-              zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
-              &          + (-160.340*ztr + 1.3566)*zsqrt          &
-              &          + (0.37335*ztr - 0.05778)*zsal
-
-              zck3p    = -3070.75*ztr - 18.126                    &
-              &          + (17.27039*ztr + 2.81197) * zsqrt       &
-              &          + (-44.99486*ztr - 0.09984) * zsal 
-
-              ! CONSTANT FOR SILICATE, MILLERO (1995)
-              zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
-              &          + (-458.79*ztr + 3.5913) * zisqrt       &
-              &          + (188.74*ztr - 1.5998) * zis           &
-              &          + (-12.1652*ztr + 0.07871) * zis2       &
-              &          + LOG(1.0 - 0.001005*zsal)
-
-               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
-               !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
-               zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
-                  &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
-                  &      - 0.07711*zsal + 0.0041249*zsal15
-
-               ! CONVERT FROM DIFFERENT PH SCALES
-               total2free  = 1.0/(1.0 + zst/zcks)
-               free2SWS    = 1. + zst/zcks + zft/(zckf*total2free)
-               total2SWS   = total2free * free2SWS
-               SWS2total   = 1.0 / total2SWS
-
-               ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
-               zak1    = 10**(zck1) * total2SWS
-               zak2    = 10**(zck2) * total2SWS
-               zakb    = EXP( zckb ) * total2SWS
-               zakw    = EXP( zckw )
-               zaksp1  = 10**(zaksp0)
-               zak1p   = exp( zck1p )
-               zak2p   = exp( zck2p )
-               zak3p   = exp( zck3p )
-               zaksi   = exp( zcksi )
-               zckf    = zckf * total2SWS
-
-               ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
-               !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
-               !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
-               !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
-               !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
-               !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
-               !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
-               !        & GIESKES (1970), P. 1285-1286 (THE SMALL
-               !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
-               !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
-               zcpexp  = zpres / (rgas*ztkel)
-               zcpexp2 = zpres * zcpexp
-
-               ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
-               !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
-               !        (CF. BROECKER ET AL., 1982)
-
-               zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
-               ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
-               ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
-               akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
-               akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
-               aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
-               akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
-               ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
-               ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk19 + devk29 * ztc + devk39 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk49 + devk59 * ztc )
-               ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
-               aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               ! CONVERT FROM DIFFERENT PH SCALES
-               total2free  = 1.0/(1.0 + zst/aks3(ji,jj,jk))
-               free2SWS    = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk)
-               total2SWS   = total2free * free2SWS
-               SWS2total   = 1.0 / total2SWS
-
-               ! Convert to total scale
-               ak13(ji,jj,jk)  = ak13(ji,jj,jk)  * SWS2total
-               ak23(ji,jj,jk)  = ak23(ji,jj,jk)  * SWS2total
-               akb3(ji,jj,jk)  = akb3(ji,jj,jk)  * SWS2total
-               akw3(ji,jj,jk)  = akw3(ji,jj,jk)  * SWS2total
-               ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total
-               ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total
-               ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total
-               aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total
-               akf3(ji,jj,jk)  = akf3(ji,jj,jk)  / total2free
-
-               ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
-               !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
-               !        (P. 1285) AND BERNER (1976)
-               zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
-               zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
-               aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
-
-               ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
-               borat(ji,jj,jk) = 0.0002414 * zcl / 10.811
-               sulfat(ji,jj,jk) = zst
-               fluorid(ji,jj,jk) = zft 
-
-               ! Iron and SIO3 saturation concentration from ...
-               sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
-               fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel )
-
-               ! Liu and Millero (1999) only valid 5 - 50 degC
-               ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16
-               fesol(ji,jj,jk,1) = 10**(-13.486 - 0.1856* zis**0.5 + 0.3073*zis + 5254.0/ztkel1)
-               fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139 * zis - 1320.0/ztkel1 )
-               fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 )
-               fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 )
-               fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**0.5 - 7980.0/ztkel1 )
-            END DO
-         END DO
-      END DO
+      DO_3D( 1, 1, 1, 1, 1, jpk )
+          ! SET PRESSION ACCORDING TO SAUNDER (1980)
+          zplat   = SIN ( ABS(gphit(ji,jj)*3.141592654/180.) )
+          zc1 = 5.92E-3 + zplat**2 * 5.25E-3
+          zpres = ((1-zc1)-SQRT(((1-zc1)**2)-(8.84E-6*gdept(ji,jj,jk,Kmm)))) / 4.42E-6
+          zpres = zpres / 10.0
+
+          ! SET ABSOLUTE TEMPERATURE
+          ztkel   = tempis(ji,jj,jk) + 273.15
+          zsal    = salinprac(ji,jj,jk) + ( 1.-tmask(ji,jj,jk) ) * 35.
+          zsqrt  = SQRT( zsal )
+          zsal15  = zsqrt * zsal
+          zlogt  = LOG( ztkel )
+          ztr    = 1. / ztkel
+          zis    = 19.924 * zsal / ( 1000.- 1.005 * zsal )
+          zis2   = zis * zis
+          zisqrt = SQRT( zis )
+          ztc     = tempis(ji,jj,jk) + ( 1.- tmask(ji,jj,jk) ) * 20.
+
+          ! CHLORINITY (WOOSTER ET AL., 1969)
+          zcl     = zsal / 1.80655
+
+          ! TOTAL SULFATE CONCENTR. [MOLES/kg soln]
+          zst     = 0.14 * zcl /96.062
+
+          ! TOTAL FLUORIDE CONCENTR. [MOLES/kg soln]
+          zft     = 0.000067 * zcl /18.9984
+
+          ! DISSOCIATION CONSTANT FOR SULFATES on free H scale (Dickson 1990)
+          zcks    = EXP(-4276.1 * ztr + 141.328 - 23.093 * zlogt         &
+          &         + (-13856. * ztr + 324.57 - 47.986 * zlogt) * zisqrt &
+          &         + (35474. * ztr - 771.54 + 114.723 * zlogt) * zis    &
+          &         - 2698. * ztr * zis**1.5 + 1776.* ztr * zis2         &
+          &         + LOG(1.0 - 0.001005 * zsal))
+
+          ! DISSOCIATION CONSTANT FOR FLUORIDES on free H scale (Dickson and Riley 79)
+          zckf    = EXP( 1590.2*ztr - 12.641 + 1.525*zisqrt   &
+          &         + LOG(1.0d0 - 0.001005d0*zsal)            &
+          &         + LOG(1.0d0 + zst/zcks))
+
+          ! DISSOCIATION CONSTANT FOR CARBONATE AND BORATE
+          zckb=  (-8966.90 - 2890.53*zsqrt - 77.942*zsal        &
+          &      + 1.728*zsal15 - 0.0996*zsal*zsal)*ztr         &
+          &      + (148.0248 + 137.1942*zsqrt + 1.62142*zsal)   &
+          &      + (-24.4344 - 25.085*zsqrt - 0.2474*zsal)      & 
+          &      * zlogt + 0.053105*zsqrt*ztkel
+
+          ! DISSOCIATION COEFFICIENT FOR CARBONATE ACCORDING TO 
+          ! MEHRBACH (1973) REFIT BY MILLERO (1995), seawater scale
+          zck1    = -1.0*(3633.86*ztr - 61.2172 + 9.6777*zlogt  &
+             - 0.011555*zsal + 0.0001152*zsal*zsal)
+          zck2    = -1.0*(471.78*ztr + 25.9290 - 3.16967*zlogt      &
+             - 0.01781*zsal + 0.0001122*zsal*zsal)
+
+          ! PKW (H2O) (MILLERO, 1995) from composite data
+          zckw    = -13847.26 * ztr + 148.9652 - 23.6521 * zlogt + ( 118.67 * ztr    &
+                    - 5.977 + 1.0495 * zlogt ) * zsqrt - 0.01615 * zsal
+
+          ! CONSTANTS FOR PHOSPHATE (MILLERO, 1995)
+         zck1p    = -4576.752*ztr + 115.540 - 18.453*zlogt   &
+         &          + (-106.736*ztr + 0.69171) * zsqrt       &
+         &          + (-0.65643*ztr - 0.01844) * zsal
+
+         zck2p    = -8814.715*ztr + 172.1033 - 27.927*zlogt  &
+         &          + (-160.340*ztr + 1.3566)*zsqrt          &
+         &          + (0.37335*ztr - 0.05778)*zsal
+
+         zck3p    = -3070.75*ztr - 18.126                    &
+         &          + (17.27039*ztr + 2.81197) * zsqrt       &
+         &          + (-44.99486*ztr - 0.09984) * zsal 
+
+         ! CONSTANT FOR SILICATE, MILLERO (1995)
+         zcksi    = -8904.2*ztr  + 117.400 - 19.334*zlogt   &
+         &          + (-458.79*ztr + 3.5913) * zisqrt       &
+         &          + (188.74*ztr - 1.5998) * zis           &
+         &          + (-12.1652*ztr + 0.07871) * zis2       &
+         &          + LOG(1.0 - 0.001005*zsal)
+
+          ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE IN SEAWATER
+          !       (S=27-43, T=2-25 DEG C) at pres =0 (atmos. pressure) (MUCCI 1983)
+          zaksp0  = -171.9065 -0.077993*ztkel + 2839.319*ztr + 71.595*LOG10( ztkel )   &
+             &      + (-0.77712 + 0.00284263*ztkel + 178.34*ztr) * zsqrt  &
+             &      - 0.07711*zsal + 0.0041249*zsal15
+
+          ! CONVERT FROM DIFFERENT PH SCALES
+          total2free  = 1.0/(1.0 + zst/zcks)
+          free2SWS    = 1. + zst/zcks + zft/(zckf*total2free)
+          total2SWS   = total2free * free2SWS
+          SWS2total   = 1.0 / total2SWS
+
+          ! K1, K2 OF CARBONIC ACID, KB OF BORIC ACID, KW (H2O) (LIT.?)
+          zak1    = 10**(zck1) * total2SWS
+          zak2    = 10**(zck2) * total2SWS
+          zakb    = EXP( zckb ) * total2SWS
+          zakw    = EXP( zckw )
+          zaksp1  = 10**(zaksp0)
+          zak1p   = exp( zck1p )
+          zak2p   = exp( zck2p )
+          zak3p   = exp( zck3p )
+          zaksi   = exp( zcksi )
+          zckf    = zckf * total2SWS
+
+          ! FORMULA FOR CPEXP AFTER EDMOND & GIESKES (1970)
+          !        (REFERENCE TO CULBERSON & PYTKOQICZ (1968) AS MADE
+          !        IN BROECKER ET AL. (1982) IS INCORRECT; HERE RGAS IS
+          !        TAKEN TENFOLD TO CORRECT FOR THE NOTATION OF pres  IN
+          !        DBAR INSTEAD OF BAR AND THE EXPRESSION FOR CPEXP IS
+          !        MULTIPLIED BY LN(10.) TO ALLOW USE OF EXP-FUNCTION
+          !        WITH BASIS E IN THE FORMULA FOR AKSPP (CF. EDMOND
+          !        & GIESKES (1970), P. 1285-1286 (THE SMALL
+          !        FORMULA ON P. 1286 IS RIGHT AND CONSISTENT WITH THE
+          !        SIGN IN PARTIAL MOLAR VOLUME CHANGE AS SHOWN ON P. 1285))
+          zcpexp  = zpres / (rgas*ztkel)
+          zcpexp2 = zpres * zcpexp
+
+          ! KB OF BORIC ACID, K1,K2 OF CARBONIC ACID PRESSURE
+          !        CORRECTION AFTER CULBERSON AND PYTKOWICZ (1968)
+          !        (CF. BROECKER ET AL., 1982)
+
+          zbuf1  = -     ( devk10 + devk20 * ztc + devk30 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk40 + devk50 * ztc )
+          ak13(ji,jj,jk) = zak1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk11 + devk21 * ztc + devk31 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk41 + devk51 * ztc )
+          ak23(ji,jj,jk) = zak2 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk12 + devk22 * ztc + devk32 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk42 + devk52 * ztc )
+          akb3(ji,jj,jk) = zakb * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk13 + devk23 * ztc + devk33 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk43 + devk53 * ztc )
+          akw3(ji,jj,jk) = zakw * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk14 + devk24 * ztc + devk34 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk44 + devk54 * ztc )
+          aks3(ji,jj,jk) = zcks * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk15 + devk25 * ztc + devk35 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk45 + devk55 * ztc )
+          akf3(ji,jj,jk) = zckf * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk17 + devk27 * ztc + devk37 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk47 + devk57 * ztc )
+          ak1p3(ji,jj,jk) = zak1p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk18 + devk28 * ztc + devk38 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk48 + devk58 * ztc )
+          ak2p3(ji,jj,jk) = zak2p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk19 + devk29 * ztc + devk39 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk49 + devk59 * ztc )
+          ak3p3(ji,jj,jk) = zak3p * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          zbuf1  =     - ( devk110 + devk210 * ztc + devk310 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk410 + devk510 * ztc )
+          aksi3(ji,jj,jk) = zaksi * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          ! CONVERT FROM DIFFERENT PH SCALES
+          total2free  = 1.0/(1.0 + zst/aks3(ji,jj,jk))
+          free2SWS    = 1. + zst/aks3(ji,jj,jk) + zft/akf3(ji,jj,jk)
+          total2SWS   = total2free * free2SWS
+          SWS2total   = 1.0 / total2SWS
+
+          ! Convert to total scale
+          ak13(ji,jj,jk)  = ak13(ji,jj,jk)  * SWS2total
+          ak23(ji,jj,jk)  = ak23(ji,jj,jk)  * SWS2total
+          akb3(ji,jj,jk)  = akb3(ji,jj,jk)  * SWS2total
+          akw3(ji,jj,jk)  = akw3(ji,jj,jk)  * SWS2total
+          ak1p3(ji,jj,jk) = ak1p3(ji,jj,jk) * SWS2total
+          ak2p3(ji,jj,jk) = ak2p3(ji,jj,jk) * SWS2total
+          ak3p3(ji,jj,jk) = ak3p3(ji,jj,jk) * SWS2total
+          aksi3(ji,jj,jk) = aksi3(ji,jj,jk) * SWS2total
+          akf3(ji,jj,jk)  = akf3(ji,jj,jk)  / total2free
+
+          ! APPARENT SOLUBILITY PRODUCT K'SP OF CALCITE 
+          !        AS FUNCTION OF PRESSURE FOLLOWING MILLERO
+          !        (P. 1285) AND BERNER (1976)
+          zbuf1  =     - ( devk16 + devk26 * ztc + devk36 * ztc * ztc )
+          zbuf2  = 0.5 * ( devk46 + devk56 * ztc )
+          aksp(ji,jj,jk) = zaksp1 * EXP( zbuf1 * zcpexp + zbuf2 * zcpexp2 )
+
+          ! TOTAL F, S, and BORATE CONCENTR. [MOLES/L]
+          borat(ji,jj,jk) = 0.0002414 * zcl / 10.811
+          sulfat(ji,jj,jk) = zst
+          fluorid(ji,jj,jk) = zft 
+
+          ! Iron and SIO3 saturation concentration from ...
+          sio3eq(ji,jj,jk) = EXP(  LOG( 10.) * ( 6.44 - 968. / ztkel )  ) * 1.e-6
+          fekeq (ji,jj,jk) = 10**( 17.27 - 1565.7 / ztkel ) 
+          ! Liu and Millero (1999) only valid 5 - 50 degC
+          ztkel1 = MAX( 5. , tempis(ji,jj,jk) ) + 273.16
+          fesol(ji,jj,jk,1) = 10**(-13.486 - 0.1856* zis**0.5 + 0.3073*zis + 5254.0/ztkel1)
+          fesol(ji,jj,jk,2) = 10**(2.517 - 0.8885*zis**0.5 + 0.2139 * zis - 1320.0/ztkel1 )
+          fesol(ji,jj,jk,3) = 10**(0.4511 - 0.3305*zis**0.5 - 1996.0/ztkel1 )
+          fesol(ji,jj,jk,4) = 10**(-0.2965 - 0.7881*zis**0.5 - 4086.0/ztkel1 )
+          fesol(ji,jj,jk,5) = 10**(4.4466 - 0.8505*zis**0.5 - 7980.0/ztkel1 )
+      END_3D
       !
       IF( ln_timing )  CALL timing_stop('p4z_che')
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zint.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zint.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zint.F90	(revision 13373)
@@ -19,4 +19,5 @@
    REAL(wp) ::   xksilim = 16.5e-6_wp   ! Half-saturation constant for the Si half-saturation constant computation
 
+#  include "do_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/TOP 4.0 , NEMO Consortium (2018)
@@ -49,10 +50,8 @@
       ! Computation of the silicon dependant half saturation  constant for silica uptake
       ! ---------------------------------------------------
-      DO ji = 1, jpi
-         DO jj = 1, jpj
-            zvar = tr(ji,jj,1,jpsil,Kbb) * tr(ji,jj,1,jpsil,Kbb)
-            xksimax(ji,jj) = MAX( xksimax(ji,jj), ( 1.+ 7.* zvar / ( xksilim * xksilim + zvar ) ) * 1e-6 )
-         END DO
-      END DO
+      DO_2D( 1, 1, 1, 1 )
+         zvar = tr(ji,jj,1,jpsil,Kbb) * tr(ji,jj,1,jpsil,Kbb)
+         xksimax(ji,jj) = MAX( xksimax(ji,jj), ( 1.+ 7.* zvar / ( xksilim * xksilim + zvar ) ) * 1e-6 )
+      END_2D
       !
       IF( nday_year == nyear_len(1) ) THEN
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zopt.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zopt.F90	(revision 13373)
@@ -32,7 +32,7 @@
    REAL(wp) ::   xsi0r       ! 1. /rn_si0
 
+  INTEGER  ::   nksr1        ! levels below which the light cannot penetrate ( depth larger than 391 m)
+
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_par      ! structure of input par
-   INTEGER , PARAMETER :: nbtimes = 366  !: maximum number of times record in a file
-   INTEGER  :: ntimes_par                ! number of time steps in a file
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   par_varsw      ! PAR fraction of shortwave
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ekb, ekg, ekr  ! wavelength (Red-Green-Blue)
@@ -279,16 +279,10 @@
          pe3(:,:,1) = zqsr(:,:)
          !
-         DO jk = 2, nksr + 1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * 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
+         DO_3D( 1, 1, 1, 1, 2, nksr1 )
+            pe0(ji,jj,jk) = pe0(ji,jj,jk-1) * EXP( -e3t(ji,jj,jk-1,Kmm) * 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_3D
         !
       ELSE   ! T- level
@@ -331,8 +325,6 @@
       ! Compute par_varsw at nit000 or only if there is more than 1 time record in par coefficient file
       IF( ln_varpar ) THEN
-         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
-         ENDIF
+         CALL fld_read( kt, 1, sf_par )
+         par_varsw(:,:) = ( sf_par(1)%fnow(:,:,1) ) / 3.0
       ENDIF
       !
@@ -377,4 +369,5 @@
       xparsw = parlux / 3.0
       xsi0r  = 1.e0 / rn_si0
+      nksr1  = nksr + 1
       !
       ! Variable PAR at the surface of the ocean
@@ -392,7 +385,4 @@
                                    ALLOCATE( sf_par(1)%fnow(jpi,jpj,1)   )
          IF( sn_par%ln_tint )      ALLOCATE( sf_par(1)%fdta(jpi,jpj,1,2) )
-
-         CALL iom_open (  TRIM( sn_par%clname ) , numpar )
-         ntimes_par = iom_getszuld( numpar )   ! get number of record in file
       ENDIF
       !
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 13373)
@@ -11,5 +11,4 @@
    USE oce_trc         ! shared variables between ocean and passive tracers
    USE trc             ! passive tracers common variables 
-   USE trcdta          ! 
    USE sms_pisces      ! PISCES Source Minus Sink variables
    USE p4zbio          ! Biological model
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/oce_sed.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/oce_sed.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/oce_sed.F90	(revision 13373)
@@ -8,6 +8,7 @@
    !!----------------------------------------------------------------------
    USE par_sed
+   USE par_trc , ONLY : rtrn  => rtrn
+   USE par_pisces
    USE timing
-   USE par_trc
 
    USE dom_oce , ONLY :   glamt     =>   glamt          !: longitude of t-point (degre)
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/par_sed.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/par_sed.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/par_sed.F90	(revision 13373)
@@ -20,13 +20,12 @@
       jp_sal   =>   jp_sal     !: indice of salintity
 
-   INTEGER, PARAMETER :: jpdta = 17
+   INTEGER, PUBLIC, PARAMETER :: jpdta = 17
 
    ! Vertical sediment geometry
-   INTEGER, PUBLIC   ::      &
-      jpksed   = 11 ,        &
-      jpksedm1 = 10
+   INTEGER, PUBLIC  :: jpksed  = 11 
+   INTEGER, PUBLIC  :: jpksedm1 
 
    ! sediment tracer species
-   INTEGER, PARAMETER ::    &
+   INTEGER, PUBLIC, PARAMETER ::    &
       jpsol =  8,           &  !: number of solid component
       jpwat = 10,           &   !: number of pore water component
@@ -36,5 +35,5 @@
    
    ! pore water components       
-   INTEGER, PARAMETER :: &
+   INTEGER, PUBLIC, PARAMETER :: &
       jwsil  = 1,        & !: silic acid
       jwoxy  = 2,        & !: oxygen
@@ -49,5 +48,5 @@
 
    ! solid components       
-   INTEGER, PARAMETER ::  &
+   INTEGER, PUBLIC, PARAMETER ::  &
       jsopal  = 1,        & !: opal sediment
       jsclay  = 2,        & !: clay
@@ -59,5 +58,5 @@
       jsfes   = 8           !: FeS
 
-   INTEGER, PARAMETER ::  &
+   INTEGER, PUBLIC, PARAMETER ::  &
       jptrased   = jpsol + jpwat , &
       jpdia3dsed = 2             , &
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sed.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sed.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sed.F90	(revision 13373)
@@ -7,5 +7,4 @@
    !!        !  06-12  (C. Ethe)  Orignal
    !!----------------------------------------------------------------------
-   USE par_sed
    USE oce_sed
    USE in_out_manager
@@ -63,5 +62,4 @@
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  solcp      !: solid sediment data at given time-step
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  solcp0     !: solid sediment data at initial time
-   REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  trc_dta
    REAL(wp), PUBLIC, DIMENSION(:,:,:), ALLOCATABLE ::  diff
 
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedarr.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedarr.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedarr.F90	(revision 13373)
@@ -10,7 +10,8 @@
    !!----------------------------------------------------------------------
    !! * Modules used
+   USE par_oce
    USE par_sed
-   USE dom_oce
-   USE sed
+   USE in_out_manager, ONLY : ln_timing
+   USE timing
 
    IMPLICIT NONE
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedchem.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedchem.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedchem.F90	(revision 13373)
@@ -6,4 +6,5 @@
    !!======================================================================
    !!   modules used
+   USE par_sed
    USE sed     ! sediment global variable
    USE sedarr
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedini.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedini.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedini.F90	(revision 13373)
@@ -9,4 +9,6 @@
    !!----------------------------------------------------------------------
    !! * Modules used
+   USE par_trc        ! need jptra, number of passive tracers
+   USE par_sed
    USE sed     ! sediment global variable
    USE sed_oce
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedwri.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedwri.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/sedwri.F90	(revision 13373)
@@ -4,4 +4,5 @@
    !!         Sediment diagnostics :  write sediment output files
    !!======================================================================
+   USE par_sed
    USE sed
    USE sedarr
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/trcdmp_sed.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/SED/trcdmp_sed.F90	(revision 13373)
@@ -91,5 +91,5 @@
                !
                jl = n_trc_index(jn) 
-               CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
+               CALL trc_dta( kt, jl, ztrcdta )   ! read tracer data at nit000
                !
                DO_2D( 1, 1, 1, 1 )
@@ -108,5 +108,5 @@
          WRITE(charout, FMT="('dmp ')")
          CALL prt_ctl_info( charout, cdcomp = 'top' )
-         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
+         CALL prt_ctl( tab4d_1=tr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm,clinfo3='trd' )
       ENDIF
       !
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/trcini_pisces.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/trcini_pisces.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/PISCES/trcini_pisces.F90	(revision 13373)
@@ -79,15 +79,9 @@
       !
       INTEGER, INTENT(in)  ::  Kmm      ! time level indices
-      REAL(wp), SAVE ::   sco2   =  2.312e-3_wp
-      REAL(wp), SAVE ::   alka0  =  2.426e-3_wp
-      REAL(wp), SAVE ::   oxyg0  =  177.6e-6_wp 
-      REAL(wp), SAVE ::   po4    =  2.165e-6_wp 
-      REAL(wp), SAVE ::   bioma0 =  1.000e-8_wp  
-      REAL(wp), SAVE ::   silic1 =  91.51e-6_wp  
-      REAL(wp), SAVE ::   no3    =  30.9e-6_wp * 7.625_wp
       !
       INTEGER  ::  ji, jj, jk, jn, ierr
       REAL(wp) ::  zcaralk, zbicarb, zco3
       REAL(wp) ::  ztmas, ztmas1
+      REAL(wp) ::  sco2, alka0, oxyg0, po4, bioma0, silic1, no3   
       CHARACTER(len = 20)  ::  cltra
       !!----------------------------------------------------------------------
@@ -103,4 +97,12 @@
          ENDIF
       ENDIF
+      !
+      sco2   =  2.312e-3_wp
+      alka0  =  2.426e-3_wp
+      oxyg0  =  177.6e-6_wp 
+      po4    =  2.165e-6_wp 
+      bioma0 =  1.000e-8_wp  
+      silic1 =  91.51e-6_wp  
+      no3    =  30.9e-6_wp * 7.625_wp
       !
       ! Allocate PISCES arrays
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcadv.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcadv.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcadv.F90	(revision 13373)
@@ -16,4 +16,5 @@
    !!   trc_adv_ini   : control the different options of advection scheme
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc        ! ocean dynamics and active tracers
    USE trc            ! ocean passive tracers variables
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcatf.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcatf.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcatf.F90	(revision 13373)
@@ -27,4 +27,5 @@
    !!   trc_atf     : time stepping on passive tracers
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc         ! ocean dynamics and tracers variables
    USE trc             ! ocean passive tracers variables
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcbbl.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcbbl.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcbbl.F90	(revision 13373)
@@ -20,4 +20,5 @@
    !!    trc_bbl      : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc        ! ocean dynamics and passive tracers variables
    USE trc            ! ocean passive tracers variables
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcdmp.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcdmp.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcdmp.F90	(revision 13373)
@@ -16,4 +16,5 @@
    !!   trc_dmp_init : initialization, namlist read, parameters control
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc         ! ocean dynamics and tracers variables
    USE trc             ! ocean passive tracers variables
@@ -108,5 +109,5 @@
                !
                jl = n_trc_index(jn) 
-               CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
+               CALL trc_dta( kt, jl, ztrcdta )   ! read tracer data at nit000
                !
                SELECT CASE ( nn_zdmp_tr )
@@ -350,5 +351,5 @@
             IF( ln_trc_ini(jn) ) THEN      ! update passive tracers arrays with input data read from file
                 jl = n_trc_index(jn)
-                CALL trc_dta( kt, Kmm, sf_trcdta(jl), rf_trfac(jl), ztrcdta )   ! read tracer data at nit000
+                CALL trc_dta( kt, jl, ztrcdta )   ! read tracer data at nit000
                 DO jc = 1, npncts
                    DO jk = 1, jpkm1
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcldf.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcldf.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcldf.F90	(revision 13373)
@@ -15,4 +15,5 @@
    !!   trc_ldf_ini   : initialization, namelist read, and parameters control
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE trc            ! ocean passive tracers variables
    USE oce_trc        ! ocean dynamics and active tracers
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcsbc.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcsbc.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trcsbc.F90	(revision 13373)
@@ -16,4 +16,5 @@
    !!   trc_sbc      : update the tracer trend at ocean surface
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc         ! ocean dynamics and active tracers variables
    USE trc             ! ocean  passive tracers variables
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trctrp.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trctrp.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trctrp.F90	(revision 13373)
@@ -13,4 +13,5 @@
    !!   trc_trp        : passive tracer transport
    !!----------------------------------------------------------------------
+   USE par_trc         ! need jptra, number of passive tracers
    USE oce_trc         ! ocean dynamics and active tracers variables
    USE trc             ! ocean passive tracers variables 
@@ -27,4 +28,5 @@
    USE bdy_oce   , ONLY: ln_bdy
    USE trcbdy          ! BDY open boundaries
+   USE in_out_manager
 
 #if defined key_agrif
@@ -62,4 +64,11 @@
       IF( .NOT. lk_c1d ) THEN
          !
+         !                                                         ! Partial top/bottom cell: GRADh( trb )  
+         IF( ln_zps ) THEN
+            IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, Kmm, jptra, tr(:,:,:,:,Kbb), pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
+            ELSE                 ; CALL zps_hde    ( kt, Kmm, jptra, tr(:,:,:,:,Kbb), gtru, gtrv )                                      !  only bottom
+            ENDIF
+         ENDIF
+         !
                                 CALL trc_sbc    ( kt,      Kmm, tr, Krhs )      ! surface boundary condition
          IF( ln_trcbc .AND. lltrcbc .AND. kt /= nit000 )  &
@@ -68,16 +77,9 @@
          IF( ln_trcdmp )        CALL trc_dmp    ( kt, Kbb, Kmm, tr, Krhs )      ! internal damping trends
          IF( ln_bdy )           CALL trc_bdy_dmp( kt, Kbb,      Krhs )      ! BDY damping trends
-                                CALL trc_adv    ( kt, Kbb, Kmm, tr, Krhs )      ! horizontal & vertical advection 
-         !                                                         ! Partial top/bottom cell: GRADh( trb )  
-         IF( ln_zps ) THEN
-           IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, Kmm, jptra, tr(:,:,:,:,Kbb), pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom
-           ELSE                 ; CALL zps_hde    ( kt, Kmm, jptra, tr(:,:,:,:,Kbb), gtru, gtrv )                                      !  only bottom
-           ENDIF
-         ENDIF
-         !                                                      
-                                CALL trc_ldf    ( kt, Kbb, Kmm,       tr, Krhs )  ! lateral mixing
 #if defined key_agrif
          IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc       ! tracers sponge
 #endif
+                                CALL trc_adv    ( kt, Kbb, Kmm, tr, Krhs )      ! horizontal & vertical advection 
+                                CALL trc_ldf    ( kt, Kbb, Kmm,       tr, Krhs )  ! lateral mixing
                                 CALL trc_zdf    ( kt, Kbb, Kmm, Krhs, tr, Kaa  )  ! vert. mixing & after tracer	==> after
                                 CALL trc_atf    ( kt, Kbb, Kmm, Kaa , tr )        ! time filtering of "now" tracer fields    
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trczdf.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trczdf.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/TRP/trczdf.F90	(revision 13373)
@@ -14,11 +14,9 @@
    !!   trc_zdf      : update the tracer trend with the vertical diffusion
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE trc           ! ocean passive tracers variables
    USE oce_trc       ! ocean dynamics and active tracers
    USE trd_oce       ! trends: ocean variables
    USE trazdf        ! tracer: vertical diffusion
-!!gm do we really need this ?
-   USE trcldf        ! passive tracers: lateral diffusion
-!!gm
    USE trdtra        ! trends manager: tracers 
    USE prtctl        ! Print control
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcdta.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcdta.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcdta.F90	(revision 13373)
@@ -31,11 +31,9 @@
    PUBLIC   trc_dta_ini     ! called in trcini.F90 
 
-   INTEGER  , SAVE, PUBLIC                             :: nb_trcdta   ! number of tracers to be initialised with data
-   INTEGER  , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_index ! indice of tracer which is initialised with data
-   INTEGER  , SAVE, PUBLIC                             :: ntra        ! MAX( 1, nb_trcdta ) to avoid compilation error with bounds checking
-   REAL(wp) , SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: rf_trfac    ! multiplicative factor for tracer values
-!$AGRIF_DO_NOT_TREAT
-   TYPE(FLD), SAVE, PUBLIC, ALLOCATABLE, DIMENSION(:)  :: sf_trcdta   ! structure of input SST (file informations, fields read)
-!$AGRIF_END_DO_NOT_TREAT
+   INTEGER  , PUBLIC                             :: nb_trcdta   ! number of tracers to be initialised with data
+   INTEGER  , PUBLIC, ALLOCATABLE, DIMENSION(:)  :: n_trc_index ! indice of tracer which is initialised with data
+   INTEGER  , PUBLIC                             :: ntra        ! MAX( 1, nb_trcdta ) to avoid compilation error with bounds checking
+   REAL(wp) ,         ALLOCATABLE, DIMENSION(:)  :: rf_trfac    ! multiplicative factor for tracer values
+   TYPE(FLD),         ALLOCATABLE, DIMENSION(:)  :: sf_trcdta   ! structure of input SST (file informations, fields read)
 
    !! Substitutions
@@ -155,5 +153,5 @@
 
 
-   SUBROUTINE trc_dta( kt, Kmm, sf_trcdta, ptrcfac, ptrcdta)
+   SUBROUTINE trc_dta( kt, kjl, ptrcdta)
       !!----------------------------------------------------------------------
       !!                   ***  ROUTINE trc_dta  ***
@@ -168,7 +166,5 @@
       !!----------------------------------------------------------------------
       INTEGER                          , INTENT(in   )   ::   kt         ! ocean time-step
-      INTEGER                          , INTENT(in   )   ::   Kmm        ! time level index
-      TYPE(FLD), DIMENSION(1)          , INTENT(inout)   ::   sf_trcdta  ! array of information on the field to read
-      REAL(wp)                         , INTENT(in   )   ::   ptrcfac    ! multiplication factor
+      INTEGER                          , INTENT(in   )   ::   kjl        ! tracer index
       REAL(wp),  DIMENSION(jpi,jpj,jpk), INTENT(inout  ) ::   ptrcdta    ! 3D data array
       !
@@ -191,5 +187,5 @@
          ! read data at kt time step
          CALL fld_read( kt, 1, sf_trcdta )
-         ptrcdta(:,:,:) = sf_trcdta(1)%fnow(:,:,:) * tmask(:,:,:)
+         ptrcdta(:,:,:) = sf_trcdta(kjl)%fnow(:,:,:) * tmask(:,:,:)
          ! 
          IF( ln_sco ) THEN                !== s- or mixed s-zps-coordinate  ==!
@@ -201,5 +197,5 @@
             DO_2D( 1, 1, 1, 1 )
                DO jk = 1, jpk                        ! determines the intepolated T-S profiles at each (i,j) points
-                  zl = gdept(ji,jj,jk,Kmm)
+                  zl = gdept_0(ji,jj,jk)
                   IF(     zl < gdept_1d(1  ) ) THEN         ! above the first level of data
                      ztp(jk) = ptrcdta(ji,jj,1)
@@ -223,25 +219,22 @@
          ELSE                                !==   z- or zps- coordinate   ==!
             ! zps-coordinate (partial steps) interpolation at the last ocean level
-!            IF( ln_zps ) THEN
-!               DO jj = 1, jpj
-!                  DO ji = 1, jpi
-!                     ik = mbkt(ji,jj) 
-!                     IF( ik > 1 ) THEN
-!                        zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
-!                        ptrcdta(ji,jj,ik) = (1.-zl) * ptrcdta(ji,jj,ik) + zl * ptrcdta(ji,jj,ik-1)
-!                     ENDIF
-!                     ik = mikt(ji,jj)
-!                     IF( ik > 1 ) THEN
-!                        zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
-!                        ptrcdta(ji,jj,ik) = (1.-zl) * ptrcdta(ji,jj,ik) + zl * ptrcdta(ji,jj,ik+1)
-!                     ENDIF
-!                  END DO
-!              END DO
-!            ENDIF
+            IF( ln_zps ) THEN
+                DO_2D( 1, 1, 1, 1 )
+                   ik = mbkt(ji,jj)
+                   IF( ik > 1 ) THEN
+                      zl = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
+                      ptrcdta(ji,jj,ik) = (1.-zl) * ptrcdta(ji,jj,ik) + zl * ptrcdta(ji,jj,ik-1)
+                   ENDIF
+                   ik = mikt(ji,jj)
+                   IF( ik > 1 ) THEN
+                      zl = ( gdept_0(ji,jj,ik) - gdept_1d(ik) ) / ( gdept_1d(ik+1) - gdept_1d(ik) )
+                      ptrcdta(ji,jj,ik) = (1.-zl) * ptrcdta(ji,jj,ik) + zl * ptrcdta(ji,jj,ik+1)
+                   ENDIF
+                 END_2D
+            ENDIF
             !
          ENDIF
-         !
          ! Scale by multiplicative factor
-         ptrcdta(:,:,:) = ptrcdta(:,:,:) * ptrcfac
+         ptrcdta(:,:,:) = ptrcdta(:,:,:) * rf_trfac(kjl)
          !
       ENDIF
@@ -256,5 +249,5 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE trc_dta( kt, sf_trcdta, ptrcfac, ptrcdta)        ! Empty routine
+   SUBROUTINE trc_dta( kt, kjl, ptrcdta)        ! Empty routine
       WRITE(*,*) 'trc_dta: You should not have seen this print! error?', kt
    END SUBROUTINE trc_dta
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcice.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcice.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcice.F90	(revision 13373)
@@ -12,4 +12,5 @@
    !!   trc_ice       :  Call the appropriate sea ice tracer subroutine
    !!----------------------------------------------------------------------
+   USE par_trc         ! need jptra, number of passive tracers
    USE oce_trc        ! shared variables between ocean and passive tracers
    USE trc            ! passive tracers common variables
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcini.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcini.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcini.F90	(revision 13373)
@@ -16,4 +16,5 @@
    !!   top_alloc :   allocate the TOP arrays
    !!----------------------------------------------------------------------
+   USE par_trc         ! need jptra, number of passive tracers
    USE oce_trc         ! shared variables between ocean and passive tracers
    USE trc             ! passive tracers common variables
@@ -95,5 +96,4 @@
       CHARACTER (len=25) :: charout
       REAL(wp), DIMENSION(jpi,jpj,jpk,jptra) :: zzmsk
-      CHARACTER (len=25), DIMENSION(jptra) :: clseb  
       !!----------------------------------------------------------------------
       !
@@ -131,9 +131,4 @@
          CALL prt_ctl_info( charout, cdcomp = 'top' )
          CALL prt_ctl( tab4d_1=tr(:,:,:,:,Kmm), mask1=tmask, clinfo=ctrcnm )
-         DO jn = 1, jptra
-            zzmsk(:,:,:,jn) = tmask(:,:,:)
-            WRITE(clseb(jn),'(a,i2.2)') 'seb ', jn
-         END DO
-         CALL prt_ctl( tab4d_1=zzmsk, mask1=tmask, clinfo=clseb )
       ENDIF
 9000  FORMAT('      tracer nb : ',i2,'      name :',a10,'      initial content :',e18.10)
@@ -254,5 +249,4 @@
         !
       ELSE                             ! Initialisation of tracer from a file that may also be used for damping
-!!gm BUG ?   if damping and restart, what's happening ?
         IF( ln_trcdta .AND. nb_trcdta > 0 ) THEN
             ! update passive tracers arrays with input data read from file
@@ -260,13 +254,5 @@
                IF( ln_trc_ini(jn) ) THEN
                   jl = n_trc_index(jn) 
-                  CALL trc_dta( nit000, Kmm, sf_trcdta(jl), rf_trfac(jl), tr(:,:,:,jn,Kmm) )
-                  !
-                  ! deallocate data structure if data are not used for damping
-                  IF( .NOT.ln_trcdmp .AND. .NOT.ln_trcdmp_clo ) THEN
-                     IF(lwp) WRITE(numout,*) 'trc_ini_state: deallocate data arrays as they are only used to initialize the run'
-                                                  DEALLOCATE( sf_trcdta(jl)%fnow )
-                     IF( sf_trcdta(jl)%ln_tint )  DEALLOCATE( sf_trcdta(jl)%fdta )
-                     !
-                  ENDIF
+                  CALL trc_dta( nit000, jl, tr(:,:,:,jn,Kmm) )
                ENDIF
             END DO
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcnam.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcnam.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcnam.F90	(revision 13373)
@@ -18,4 +18,5 @@
    !!   trc_nam    :  Read and print options for the passive tracer run (namelist)
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc     ! shared variables between ocean and passive tracers
    USE trc         ! passive tracers common variables
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcrst.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcrst.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcrst.F90	(revision 13373)
@@ -19,4 +19,5 @@
    !!   trc_rst_wri    : write restart file
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc
    USE trc
Index: /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcstp.F90
===================================================================
--- /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcstp.F90	(revision 13372)
+++ /NEMO/branches/2020/dev_r13333_TOP-05_Ethe_Agrif/src/TOP/trcstp.F90	(revision 13373)
@@ -11,4 +11,5 @@
    !!   trc_stp       : passive tracer system time-stepping
    !!----------------------------------------------------------------------
+   USE par_trc        ! need jptra, number of passive tracers
    USE oce_trc        ! ocean dynamics and active tracers variables
    USE sbc_oce
