Index: NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/TRP/trcsbc.F90
===================================================================
--- NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/TRP/trcsbc.F90	(revision 15379)
+++ NEMO/branches/2021/dev_r14318_RK3_stage1/src/TOP/TRP/trcsbc.F90	(revision 15380)
@@ -232,16 +232,12 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) ::   ptr       ! passive tracers and RHS of tracer equation
       !
-      INTEGER  ::   ji, jj, jn                      ! dummy loop indices
-      REAL(wp) ::   zse3t, zrtrn, zfact     ! local scalars
-      REAL(wp) ::   zftra, zdtra, ztfx, ztra   !   -      -
+      INTEGER  ::   ji, jj, jn           ! dummy loop indices
+      REAL(wp) ::   z1_rho0_e3t          ! local scalars
+      REAL(wp) ::   zftra, zdtra, ztfx   !   -      -
       CHARACTER (len=22) :: charout
-      REAL(wp), DIMENSION(jpi,jpj)   ::   zmfx
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrtrd
       !!---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('trc_sbc_RK3')
-      !
-      ! Allocate temporary workspace
-      IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) )
       !
       IF( kt == nittrc000 ) THEN
@@ -250,98 +246,109 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
-!!st: not sure of the way to deal about this
+      !
+!!st note that trc_sbc can be removed only re-use in atf (not relevant for RK3)
       SELECT CASE( kstg )
-      !
+         !
       CASE( 1 , 2 )                       !=  stage 1 and 2  =!   only in non linear ssh
          !
          IF( .NOT.ln_linssh ) THEN           !* only passive tracer fluxes associated with mass fluxes
+            !                                        ! no passive tracer concentration modification due to ssh variation
+!!st emp includes fmm see iceupdate.F90
+!!not sure about trc_i case... (1)
+            DO jn = 1, jptra
+               DO_2D( 0, 0, 0, 1 )              !!st WHY 1 : exterior here ? 
+                  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) * ptr(ji,jj,1,jn,Kmm) * z1_rho0_e3t
+               END_2D
+            END DO
+            !
+         ENDIF
+         !
+      CASE( 3 )
+         !
+         ! Allocate temporary workspace
+         IF( l_trdtrc )  ALLOCATE( ztrtrd(jpi,jpj,jpk) )
+         !
+         DO jn = 1, jptra
+            IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends
+         END DO
+         !
+         IF( ln_linssh ) THEN                !* linear free surface (add concentration/dilution effect artificially since no volume variation)
             !
             SELECT CASE ( nn_ice_tr )
                !
-            CASE ( -1 , 0 , 1 )                      ! no passive tracer concentration modification due to ssh variation
-               !
-!!st emp includes fmm see iceupdate.F90
+            CASE ( -1 ) ! No tracers in sea ice (null concentration in sea ice)
+               !
                DO jn = 1, jptra
                   DO_2D( 0, 0, 0, 1 )
-                     sbc_trc(ji,jj,jn) =  - emp(ji,jj) * ptr(ji,jj,1,jn,Kmm) * r1_rho0
+                     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) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
                   END_2D
                END DO
                !
-            END SELECT
-            !
-         ENDIF
-         !
-      CASE( 3 )
-!!st copy of existing code
-         !
-         IF( .NOT.ln_linssh ) THEN    !!st concentration/dilution effect due to volume variation 
-            SELECT CASE ( nn_ice_tr )
-            !
-            CASE ( 0 )  ! Same concentration in sea ice and in the ocean
+            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 )
-                     sbc_trc(ji,jj,jn) =  - fmmflx(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
+                     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
                END DO
                !
-            CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice 
+            CASE ( 1 )  ! Specific treatment of sea ice fluxes with an imposed concentration in sea ice !!st TODO : check Christian new implementation
                !
                DO jn = 1, jptra
                   DO_2D( 0, 0, 0, 1 )
-                  ! tracer flux at the ice/ocean interface (tracer/m2/s)
-                  zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
-                  !                                         ! only used in the levitating sea ice case
-                  ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
-                  ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
-                  ztfx  = zftra                        ! net tracer flux
-                  !
-                  zdtra = r1_rho0 * ( ztfx -  fmmflx(ji,jj) * ptr(ji,jj,1,jn,Kmm) ) 
-                  IF ( zdtra < 0. ) THEN
-                     zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise
-                  ENDIF
-                  sbc_trc(ji,jj,jn) =  zdtra 
+                     z1_rho0_e3t = r1_rho0  / e3t(ji,jj,1,Kmm)
+                     ! tracer flux at the ice/ocean interface (tracer/m2/s)
+                     zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
+                     !                                         ! only used in the levitating sea ice case
+                     ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
+                     ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
+                     ztfx  = zftra                        ! net tracer flux
+                     !
+                     zdtra = r1_rho0 * ( ztfx +  ( emp(ji,jj) - fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) 
+                     IF ( zdtra < 0. ) THEN
+                        zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise
+                     ENDIF
+                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + zdtra
                   END_2D
                END DO
                !
             END SELECT
-         ELSE                         !linear ssh !!st need to mimic concentration/dilution effect since no volume variation 
+            !
+         ELSE                                !* non linear free surface (concentration/dilution effect due to volume variation)
+            !
             SELECT CASE ( nn_ice_tr )
-            !
-            CASE ( -1 ) ! No tracers in sea ice (null concentration in sea ice)
+            ! CASE ( -1 ) natural concentration/dilution effect due to volume variation : nothing to do
+            !
+            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 )
-                  sbc_trc(ji,jj,jn) = emp(ji,jj) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
+                     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
                END DO
                !
-            CASE ( 0 )  ! Same concentration in sea ice and in the ocean
+            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 )
-                     sbc_trc(ji,jj,jn) = ( emp(ji,jj) - fmmflx(ji,jj) ) * r1_rho0 * ptr(ji,jj,1,jn,Kmm)
+                     ! tracer flux at the ice/ocean interface (tracer/m2/s)
+                     zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
+                     !                                         ! only used in the levitating sea ice case
+                     ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
+                     ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
+                     ztfx  = zftra                        ! net tracer flux
+                     !
+                     zdtra = r1_rho0 * ( ztfx -  fmmflx(ji,jj) * ptr(ji,jj,1,jn,Kmm) ) 
+                     IF ( zdtra < 0. ) THEN
+                        zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise
+                     ENDIF
+                     ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) + zdtra
                   END_2D
                END DO
                !
-            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 )
-                  ! tracer flux at the ice/ocean interface (tracer/m2/s)
-                  zftra = - trc_i(ji,jj,jn) * fmmflx(ji,jj) ! uptake of tracer in the sea ice
-                  !                                         ! only used in the levitating sea ice case
-                  ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
-                  ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
-                  ztfx  = zftra                        ! net tracer flux
-                  !
-                  zdtra = r1_rho0 * ( ztfx +  ( emp(ji,jj) - fmmflx(ji,jj) ) * ptr(ji,jj,1,jn,Kmm) ) 
-                  IF ( zdtra < 0. ) THEN
-                     zdtra  = MAX(zdtra, -ptr(ji,jj,1,jn,Kmm) * e3t(ji,jj,1,Kmm) / rDt_trc )   ! avoid negative concentrations to arise
-                  ENDIF
-                  sbc_trc(ji,jj,jn) =  zdtra 
-                  END_2D
-               END DO
-               !
             END SELECT
             !
@@ -349,29 +356,22 @@
          !
          !
+!!st useless trc_sbc only in the interior even in MLF case         CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1.0_wp )
+         !                                       Concentration dilution effect on tracers due to evaporation & precipitation 
+         DO jn = 1, jptra
+            !
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) 'trc_sbc_RK3 : Runge Kutta 3rd order at stage ', kstg, jn
+            IF(lwp) WRITE(numout,*)
+            !
+            IF( l_trdtrc ) THEN
+               ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) - ztrtrd(:,:,:)
+               CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_nsr, ztrtrd )
+            END IF
+            !
+         END DO
+         !
+         IF( l_trdtrc )  DEALLOCATE( ztrtrd )
          !
       END SELECT
-      !
-      CALL lbc_lnk( 'trcsbc', sbc_trc(:,:,:), 'T', 1.0_wp )
-      !                                       Concentration dilution effect on tracers due to evaporation & precipitation 
-      DO jn = 1, jptra
-         !
-         IF( l_trdtrc )   ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs)  ! save trends
-         !
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'trc_sbc_RK3 : Runge Kutta 3rd order at stage ', kstg, jn
-         IF(lwp) WRITE(numout,*) 'emp', MAXVAL(emp(:,:))
-         IF(lwp) WRITE(numout,*) 'sbc_trc', MAXVAL(sbc_trc(:,:,jn))
-         IF(lwp) WRITE(numout,*)
-         DO_2D( 0, 0, 0, 1 )
-            zse3t = 1._wp / e3t(ji,jj,1,Kmm)
-            ptr(ji,jj,1,jn,Krhs) = ptr(ji,jj,1,jn,Krhs) +  sbc_trc(ji,jj,jn) * zse3t
-         END_2D
-         !
-         IF( l_trdtrc ) THEN
-            ztrtrd(:,:,:) = ptr(:,:,:,jn,Krhs) - ztrtrd(:,:,:)
-            CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_nsr, ztrtrd )
-         END IF
-         !
-      END DO
       !
       IF( sn_cfctl%l_prttrc )   THEN
@@ -379,6 +379,4 @@
                                            CALL prt_ctl( tab4d_1=ptr(:,:,:,:,Krhs), mask1=tmask, clinfo=ctrcnm, clinfo3='trd' )
       ENDIF
-      !
-      IF( l_trdtrc )  DEALLOCATE( ztrtrd )
       !
       IF( ln_timing )   CALL timing_stop('trc_sbc_RK3')
