Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv_crs.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv_crs.F90	(revision 5105)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcadv_crs.F90	(revision 5105)
@@ -0,0 +1,228 @@
+MODULE trcadv_crs
+   !!==============================================================================
+   !!                       ***  MODULE  trcadv  ***
+   !! Ocean passive tracers:  advection trend 
+   !!==============================================================================
+   !! History :  2.0  !  05-11  (G. Madec)  Original code
+   !!            3.0  !  10-06  (C. Ethe)   Adapted to passive tracers
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!   trc_adv      : compute ocean tracer advection trend
+   !!   trc_adv_ctl  : control the different options of advection scheme
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! ocean dynamics and active tracers
+   !???USE oce_trc, ONLY: un,vn,wn
+   USE trc             ! ocean passive tracers variables
+   USE trcnam_trp      ! passive tracers transport namelist variables
+   USE traadv_cen2     ! 2nd order centered scheme (tra_adv_cen2   routine)
+   USE traadv_tvd_crs  ! TVD      scheme           (tra_adv_tvd    routine)
+   USE traadv_muscl    ! MUSCL    scheme           (tra_adv_muscl  routine)
+   USE traadv_muscl2   ! MUSCL2   scheme           (tra_adv_muscl2 routine)
+   USE traadv_ubs      ! UBS      scheme           (tra_adv_ubs    routine)
+   USE traadv_qck      ! QUICKEST scheme           (tra_adv_qck    routine)
+   USE traadv_eiv      ! eddy induced velocity     (tra_adv_eiv    routine)
+   USE ldftra_oce      ! lateral diffusion coefficient on tracers
+   USE prtctl_trc      ! Print control
+   USE crs , ONLY : e2e3u_msk , e1e3v_msk , e1e2w_msk,jpi_crs,jpj_crs
+   USE timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_adv_crs          ! routine called by step module
+   PUBLIC   trc_adv_alloc_crs   ! routine called by nemogcm module
+
+   INTEGER ::   nadv   ! choice of the type of advection scheme
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::   r2dt  ! vertical profile time-step, = 2 rdttra
+   !                                                    ! except at nitrrc000 (=rdttra) if neuler=0
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcadv.F90 3294 2012-01-28 16:44:18Z rblod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION trc_adv_alloc_crs()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_adv_alloc  ***
+      !!----------------------------------------------------------------------
+
+      ALLOCATE( r2dt(jpk), STAT=trc_adv_alloc_crs )
+
+      IF( trc_adv_alloc_crs /= 0 ) CALL ctl_warn('trc_adv_alloc : failed to allocate array.')
+
+   END FUNCTION trc_adv_alloc_crs
+
+
+   SUBROUTINE trc_adv_crs( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_adv  ***
+      !!
+      !! ** Purpose :   compute the ocean tracer advection trend.
+      !!
+      !! ** Method  : - Update the tracer with the advection term following nadv
+      !!----------------------------------------------------------------------
+      !!
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step index
+      !
+      INTEGER ::   jk 
+      INTEGER ::   ji,jj
+      CHARACTER (len=22) ::   charout
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zun, zvn, zwn  ! effective velocity
+      !!----------------------------------------------------------------------
+      !
+      
+      IF( nn_timing == 1 )  CALL timing_start('trc_adv')
+      !
+      CALL wrk_alloc( jpi, jpj, jpk, zun, zvn, zwn )
+      !
+
+      IF( kt == nittrc000 )   CALL trc_adv_ctl_crs          ! initialisation & control of options
+
+#if ! defined key_pisces
+      IF( neuler == 0 .AND. kt == nittrc000 ) THEN     ! at nittrc000
+         r2dt(:) =  rdttrc(:)           ! = rdttrc (restarting with Euler time stepping)
+      ELSEIF( kt <= nittrc000 + 1 ) THEN          ! at nittrc000 or nittrc000+1
+         r2dt(:) = 2. * rdttrc(:)       ! = 2 rdttrc (leapfrog)
+      ENDIF
+#else
+      r2dt(:) =  rdttrc(:)              ! = rdttrc (for PISCES use Euler time stepping)
+#endif
+
+  !    IF(lwp) WRITE(numout,*) 'TEST', e1e2t
+      !                                                   ! effective transport
+!         IF(lwp) WRITE(numout,*) 'un', maxval(un(:,:,:))
+!         IF(lwp) WRITE(numout,*) 'un', minval(un(:,:,:))
+!         IF(lwp) WRITE(numout,*) 'vn', maxval(vn(:,:,:))
+!         IF(lwp) WRITE(numout,*) 'vn', minval(vn(:,:,:))
+!         IF(lwp) WRITE(numout,*) 'wn', maxval(wn(:,:,:))
+!         IF(lwp) WRITE(numout,*) 'wn', minval(wn(:,:,:))
+      DO jk = 1, jpkm1
+         !                                                ! eulerian transport only
+         zun(:,:,jk) = e2e3u_msk(:,:,jk) * un(:,:,jk)
+         zvn(:,:,jk) = e1e3v_msk(:,:,jk) * vn(:,:,jk)
+         zwn(:,:,jk) = e1e2w_msk(:,:,jk) * wn(:,:,jk)
+         !
+      END DO
+
+         IF(lwp)WRITE(numout,*)"jpi_crs jpj_crs jpk ",jpi_crs,jpj_crs,jpk
+         DO jk=1,jpk
+           DO jj = 1, jpj_crs
+               DO ji = 1, jpi_crs
+                  IF( zwn(ji,jj,jk) .NE. zwn(ji,jj,jk) )WRITE(narea+200,*)"trcadv_zwn",zwn(ji,jj,jk) ; call flush(narea+200)
+               END DO
+            END DO
+         END DO
+
+
+      zwn(:,:,jpk) = 0.e0                                 ! no transport trough the bottom
+
+      IF( lk_traldf_eiv .AND. .NOT. ln_traldf_grif )   &  ! add the eiv transport (if necessary)
+         &              CALL tra_adv_eiv( kt, nittrc000, zun, zvn, zwn, 'TRC' )
+      !
+      SELECT CASE ( nadv )                            !==  compute advection trend and add it to general trend  ==!
+!cbr      CASE ( 1 )   ;    CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )   !  2nd order centered
+      CASE ( 2 )   ;    CALL tra_adv_tvd_crs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  TVD 
+!cbr      CASE ( 3 )   ;    CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra )   !  MUSCL 
+!cbr      CASE ( 4 )   ;    CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  MUSCL2 
+!cbr      CASE ( 5 )   ;    CALL tra_adv_ubs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  UBS 
+!cbr      CASE ( 6 )   ;    CALL tra_adv_qck   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )   !  QUICKEST 
+      !
+      CASE (-1 )                                      !==  esopa: test all possibility with control print  ==!
+!         CALL tra_adv_cen2  ( kt, nittrc000, 'TRC',       zun, zvn, zwn, trb, trn, tra, jptra )          
+!         WRITE(charout, FMT="('adv1')")  ; CALL prt_ctl_trc_info(charout)
+!                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
+         CALL tra_adv_tvd_crs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+         WRITE(charout, FMT="('adv2')")  ; CALL prt_ctl_trc_info(charout)
+                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
+!         CALL tra_adv_muscl ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb,      tra, jptra )          
+!         WRITE(charout, FMT="('adv3')")  ; CALL prt_ctl_trc_info(charout)
+!                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
+!         CALL tra_adv_muscl2( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+!         WRITE(charout, FMT="('adv4')")  ; CALL prt_ctl_trc_info(charout)
+!                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
+!         CALL tra_adv_ubs   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+!         WRITE(charout, FMT="('adv5')")  ; CALL prt_ctl_trc_info(charout)
+!                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
+!         CALL tra_adv_qck   ( kt, nittrc000, 'TRC', r2dt, zun, zvn, zwn, trb, trn, tra, jptra )          
+!         WRITE(charout, FMT="('adv6')")  ; CALL prt_ctl_trc_info(charout)
+                                           CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
+         !
+      END SELECT
+
+      !                                              ! print mean trends (used for debugging)
+      IF( ln_ctl )   THEN
+         WRITE(charout, FMT="('adv ')")  ;  CALL prt_ctl_trc_info(charout)
+                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      END IF
+      !
+      CALL wrk_dealloc( jpi, jpj, jpk, zun, zvn, zwn )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('trc_adv')
+      !
+   END SUBROUTINE trc_adv_crs
+
+
+   SUBROUTINE trc_adv_ctl_crs
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_adv_ctl  ***
+      !!                
+      !! ** Purpose : Control the consistency between namelist options for 
+      !!              passive tracer advection schemes and set nadv
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio
+      !!----------------------------------------------------------------------
+
+      ioptio = 0                      ! Parameter control
+      IF( ln_trcadv_cen2   )   ioptio = ioptio + 1
+      IF( ln_trcadv_tvd    )   ioptio = ioptio + 1
+      IF( ln_trcadv_muscl  )   ioptio = ioptio + 1
+      IF( ln_trcadv_muscl2 )   ioptio = ioptio + 1
+      IF( ln_trcadv_ubs    )   ioptio = ioptio + 1
+      IF( ln_trcadv_qck    )   ioptio = ioptio + 1
+      IF( lk_esopa         )   ioptio =          1
+
+      IF( ioptio /= 1 )   CALL ctl_stop( 'Choose ONE advection scheme in namelist namtrc_adv' )
+
+      !                              ! Set nadv
+      IF( ln_trcadv_cen2   )   nadv =  1
+      IF( ln_trcadv_tvd    )   nadv =  2
+      IF( ln_trcadv_muscl  )   nadv =  3
+      IF( ln_trcadv_muscl2 )   nadv =  4
+      IF( ln_trcadv_ubs    )   nadv =  5
+      IF( ln_trcadv_qck    )   nadv =  6
+      IF( lk_esopa         )   nadv = -1
+
+      IF(lwp) THEN                   ! Print the choice
+         WRITE(numout,*)
+         IF( nadv ==  1 )   WRITE(numout,*) '         2nd order scheme is used'
+         IF( nadv ==  2 )   WRITE(numout,*) '         TVD       scheme is used'
+         IF( nadv ==  3 )   WRITE(numout,*) '         MUSCL     scheme is used'
+         IF( nadv ==  4 )   WRITE(numout,*) '         MUSCL2    scheme is used'
+         IF( nadv ==  5 )   WRITE(numout,*) '         UBS       scheme is used'
+         IF( nadv ==  6 )   WRITE(numout,*) '         QUICKEST  scheme is used'
+         IF( nadv == -1 )   WRITE(numout,*) '         esopa test: use all advection scheme'
+      ENDIF
+      !
+   END SUBROUTINE trc_adv_ctl_crs
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                         Empty module
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_adv_crs( kt )
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trc_adv: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_adv_crs
+#endif
+
+  !!======================================================================
+END MODULE trcadv_crs
Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl_crs.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl_crs.F90	(revision 5105)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcbbl_crs.F90	(revision 5105)
@@ -0,0 +1,117 @@
+MODULE trcbbl_crs
+  !!======================================================================
+   !!                       ***  MODULE  trcbbl  ***
+   !! Ocean passive tracers physics :  advective and/or diffusive bottom boundary 
+   !!                                  layer scheme
+   !!======================================================================
+   !!==============================================================================
+   !! History :  OPA  !  1996-06  (L. Mortier)  Original code
+   !!            8.0  !  1997-11  (G. Madec)    Optimization
+   !!   NEMO     1.0  !  2002-08  (G. Madec)  free form + modules
+   !!             -   !  2004-01  (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl
+   !!            3.3  !  2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization 
+   !!             -   !  2010-04  (G. Madec)  Campin & Goosse advective bbl 
+   !!             -   !  2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
+   !!----------------------------------------------------------------------
+#if  defined key_top &&  defined key_trabbl 
+   !!----------------------------------------------------------------------
+   !!   'key_trabbl                      diffusive or/and adevective bottom boundary layer
+   !!----------------------------------------------------------------------
+   !!    trc_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
+   !!----------------------------------------------------------------------
+   USE oce_trc             ! ocean dynamics and active tracers variables
+   USE trc                 ! ocean passive tracers variables
+   USE trcnam_trp      ! passive tracers transport namelist variables
+   USE trabbl_crs              ! 
+   USE prtctl_trc          ! Print control for debbuging
+   USE trdmod_oce
+   USE trdtra
+
+   PUBLIC   trc_bbl_crs       !  routine called by step.F90
+
+
+   !! * Substitutions
+#  include "top_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcbbl.F90 3294 2012-01-28 16:44:18Z rblod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+
+   SUBROUTINE trc_bbl_crs( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE bbl  ***
+      !!                   
+      !! ** Purpose :   Compute the before tracer (t & s) trend associated 
+      !!     with the bottom boundary layer and add it to the general trend
+      !!     of tracer equations.
+      !!
+      !!----------------------------------------------------------------------  
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step 
+      CHARACTER (len=22) :: charout
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd
+      !!----------------------------------------------------------------------
+      !
+!      WRITE(numout,*) '          nn_bbl_ldf_test     = ' , nn_bbl_ldf
+      
+      IF( nn_timing == 1 )  CALL timing_start('trc_bbl')
+      !
+      IF( .NOT. lk_offline .AND. nn_dttrc == 1 ) THEN
+         CALL bbl_crs( kt, nittrc000, 'TRC' )      ! Online coupling with dynamics  : Computation of bbl coef and bbl transport
+         l_bbl_crs = .FALSE.                       ! Offline coupling with dynamics : Read bbl coef and bbl transport from input files
+      ENDIF
+
+      IF( l_trdtrc )  THEN
+         CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd ) ! temporary save of trends
+         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+      ENDIF
+
+      !* Diffusive bbl :
+      IF( nn_bbl_ldf == 1 ) THEN
+  
+         CALL tra_bbl_dif_crs( trb, tra, jptra )  
+         IF( ln_ctl )   THEN
+            WRITE(charout, FMT="(' bbl_dif')")  ;  CALL prt_ctl_trc_info(charout)
+            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         ENDIF
+         !
+      END IF
+
+      !* Advective bbl : bbl upstream advective trends added to the tracer trends
+      IF( nn_bbl_adv /= 0 ) THEN
+         !
+         CALL tra_bbl_adv_crs( trb, tra, jptra )  
+         IF( ln_ctl )   THEN
+            WRITE(charout, FMT="(' bbl_adv')")  ;  CALL prt_ctl_trc_info(charout)
+            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         ENDIF
+         !
+      END IF
+
+      IF( l_trdtrc )   THEN                      ! save the horizontal diffusive trends for further diagnostics
+        DO jn = 1, jptra
+           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
+           CALL trd_tra( kt, 'TRC', jn, jptra_trd_ldf, ztrtrd(:,:,:,jn) )
+        END DO
+        CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd ) ! temporary save of trends
+      ENDIF
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('trc_bbl')
+      !
+   END SUBROUTINE trc_bbl_crs
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                      No bottom boundary layer scheme
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_bbl_crs( kt )              ! Empty routine
+      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_bbl_crs
+#endif
+
+   !!======================================================================
+END MODULE trcbbl_crs
Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf_crs.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf_crs.F90	(revision 5105)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf_crs.F90	(revision 5105)
@@ -0,0 +1,277 @@
+MODULE trcldf_crs
+   !!======================================================================
+   !!                       ***  MODULE  trcldf  ***
+   !! Ocean Passive tracers : lateral diffusive trends
+   !!=====================================================================
+   !! History :  9.0  ! 2005-11 (G. Madec)  Original code
+   !!       NEMO 3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   trc_ldf     : update the tracer trend with the lateral diffusion
+   !!       ldf_ctl : initialization, namelist read, and parameters control
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! ocean dynamics and active tracers
+   USE trc             ! ocean passive tracers variables
+   USE trcnam_trp      ! passive tracers transport namelist variables
+   USE ldftra_oce      ! lateral diffusion coefficient on tracers
+   USE ldfslp          ! ???
+   USE traldf_bilapg   ! lateral mixing            (tra_ldf_bilapg routine)
+   USE traldf_bilap    ! lateral mixing            (tra_ldf_bilap routine)
+   USE traldf_iso_crs      ! lateral mixing            (tra_ldf_iso routine)
+   USE traldf_iso_grif ! lateral mixing          (tra_ldf_iso_grif routine)
+   USE traldf_lap_crs      ! lateral mixing            (tra_ldf_lap routine)
+   USE trdmod_oce
+   USE trdtra
+   USE prtctl_trc      ! Print control
+   !USE crs
+   USE trc_oce, ONLY : lk_offline  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_ldf_crs    ! called by step.F90
+   !                                                 !!: ** lateral mixing namelist (nam_trcldf) **
+   REAL(wp) ::  rldf_rat    ! ratio between active and passive tracers diffusive coefficient
+   INTEGER  ::  nldf = 0   ! type of lateral diffusion used defined from ln_trcldf_... namlist logicals)
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcldf.F90 3294 2012-01-28 16:44:18Z rblod $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE trc_ldf_crs( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_ldf  ***
+      !!
+      !! ** Purpose :   compute the lateral ocean tracer physics.
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
+      !!
+      INTEGER            :: jn
+      CHARACTER (len=22) :: charout
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('trc_ldf')
+      !
+      IF( kt == nittrc000 )   CALL ldf_ctl          ! initialisation & control of options
+
+      rldf = rldf_rat
+
+      IF( l_trdtrc )  THEN
+         CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd )
+         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+      ENDIF
+
+!      WRITE(numout,*) 'nldf', nldf
+      SELECT CASE ( nldf )                       ! compute lateral mixing trend and add it to the general trend
+      CASE ( 0 )   ;   CALL tra_ldf_lap_crs   ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level laplacian
+      CASE ( 1 )                                                                                            ! rotated laplacian
+                       IF( ln_traldf_grif ) THEN
+                          CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+                       ELSE
+                          CALL tra_ldf_iso_crs     ( kt, nittrc000, 'TRC', gtru ,gtrv , trb, tra, jptra, rn_ahtb_0 )
+                       ENDIF
+      CASE ( 2 )   ;   CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )  ! iso-level bilaplacian
+      CASE ( 3 )   ;   CALL tra_ldf_bilapg( kt, nittrc000, 'TRC',             trb, tra, jptra            )  ! s-coord. horizontal bilaplacian
+         !
+      CASE ( -1 )                                     ! esopa: test all possibility with control print
+         CALL tra_ldf_lap_crs   ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )
+         WRITE(charout, FMT="('ldf0 ')") ;  CALL prt_ctl_trc_info(charout)
+                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         IF( ln_traldf_grif ) THEN
+            CALL tra_ldf_iso_grif( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+         ELSE
+            CALL tra_ldf_iso_crs     ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra, rn_ahtb_0 )
+         ENDIF
+         WRITE(charout, FMT="('ldf1 ')") ;  CALL prt_ctl_trc_info(charout)
+                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL tra_ldf_bilap ( kt, nittrc000, 'TRC', gtru, gtrv, trb, tra, jptra            )
+         WRITE(charout, FMT="('ldf2 ')") ;  CALL prt_ctl_trc_info(charout)
+                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL tra_ldf_bilapg( kt, nittrc000, 'TRC',             trb, tra, jptra            )
+         WRITE(charout, FMT="('ldf3 ')") ;  CALL prt_ctl_trc_info(charout)
+                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      END SELECT
+      !
+      IF( l_trdtrc )   THEN                      ! save the horizontal diffusive trends for further diagnostics
+        DO jn = 1, jptra
+           ztrtrd(:,:,:,jn) = tra(:,:,:,jn) - ztrtrd(:,:,:,jn)
+           CALL trd_tra( kt, 'TRC', jn, jptra_trd_ldf, ztrtrd(:,:,:,jn) )
+        END DO
+        CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
+      ENDIF
+      !                                          ! print mean trends (used for debugging)
+      IF( ln_ctl )   THEN
+         WRITE(charout, FMT="('ldf ')") ;  CALL prt_ctl_trc_info(charout)
+                                           CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      ENDIF
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('trc_ldf')
+      !
+   END SUBROUTINE trc_ldf_crs
+
+
+   SUBROUTINE ldf_ctl
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_ctl  ***
+      !!
+      !! ** Purpose :   Choice of the operator for the lateral tracer diffusion
+      !!
+      !! ** Method  :   set nldf from the namtra_ldf logicals
+      !!      nldf == -2   No lateral diffusion
+      !!      nldf == -1   ESOPA test: ALL operators are used
+      !!      nldf ==  0   laplacian operator
+      !!      nldf ==  1   Rotated laplacian operator
+      !!      nldf ==  2   bilaplacian operator
+      !!      nldf ==  3   Rotated bilaplacian
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio, ierr         ! temporary integers
+      !!----------------------------------------------------------------------
+
+      IF (ABS(rn_aht_0) < 2._wp*TINY(1.e0)) THEN
+       !  IF (ABS(rn_ahtrc_0) < 2._wp*TINY(1.e0)) THEN  !! cc
+       !     rldf_rat = 1.0_wp
+       !  ELSE
+       !     CALL ctl_stop( 'STOP', 'ldf_ctl : cannot define rldf_rat, rn_aht_0==0, rn_ahtrc_0 /=0' )
+       !  END IF
+          rldf_rat = 1.0_wp !!cc
+      ELSE
+         rldf_rat = rn_ahtrc_0 / rn_aht_0
+      END IF
+      !  Define the lateral mixing oparator for tracers
+      ! ===============================================
+
+      !                               ! control the input
+      ioptio = 0
+      IF( ln_trcldf_lap   )   ioptio = ioptio + 1
+      IF( ln_trcldf_bilap )   ioptio = ioptio + 1
+      IF( ioptio >  1 )   CALL ctl_stop( '          use ONE or NONE of the 2 lap/bilap operator type on tracer' )
+      IF( ioptio == 0 )   nldf = -2   ! No lateral diffusion
+      ioptio = 0
+      IF( ln_trcldf_level )   ioptio = ioptio + 1
+      IF( ln_trcldf_hor   )   ioptio = ioptio + 1
+      IF( ln_trcldf_iso   )   ioptio = ioptio + 1
+      IF( ioptio /= 1 )   CALL ctl_stop( '          use only ONE direction (level/hor/iso)' )
+
+      ! defined the type of lateral diffusion from ln_trcldf_... logicals
+      ! CAUTION : nldf = 1 is used in trazdf_imp, change it carefully
+      ierr = 0
+      IF( ln_trcldf_lap ) THEN       ! laplacian operator
+         IF ( ln_zco ) THEN                ! z-coordinate
+            IF ( ln_trcldf_level )   nldf = 0      ! iso-level  (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = 0      ! horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN             ! z-coordinate
+            IF ( ln_trcldf_level )   ierr = 1      ! iso-level not allowed
+            IF ( ln_trcldf_hor   )   nldf = 0      ! horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_sco ) THEN             ! z-coordinate
+            IF ( ln_trcldf_level )   nldf = 0      ! iso-level  (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = 1      ! horizontal (   rotation)
+            IF ( ln_trcldf_iso   )   nldf = 1      ! isoneutral (   rotation)
+         ENDIF
+      ENDIF
+
+      IF( ln_trcldf_bilap ) THEN      ! bilaplacian operator
+         IF ( ln_zco ) THEN                ! z-coordinate
+            IF ( ln_trcldf_level )   nldf = 2      ! iso-level  (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = 2      ! horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_zps ) THEN             ! z-coordinate
+            IF ( ln_trcldf_level )   ierr = 1      ! iso-level not allowed
+            IF ( ln_trcldf_hor   )   nldf = 2      ! horizontal (no rotation)
+            IF ( ln_trcldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+         IF ( ln_sco ) THEN             ! z-coordinate
+            IF ( ln_trcldf_level )   nldf = 2      ! iso-level  (no rotation)
+            IF ( ln_trcldf_hor   )   nldf = 3      ! horizontal (   rotation)
+            IF ( ln_trcldf_iso   )   ierr = 2      ! isoneutral (   rotation)
+         ENDIF
+      ENDIF
+
+      IF( ierr == 1 )   CALL ctl_stop( ' iso-level in z-coordinate - partial step, not allowed' )
+      IF( ierr == 2 )   CALL ctl_stop( ' isoneutral bilaplacian operator does not exist' )
+      IF( lk_traldf_eiv .AND. .NOT.ln_trcldf_iso )   &
+           CALL ctl_stop( '          eddy induced velocity on tracers',   &
+           &              ' the eddy induced velocity on tracers requires isopycnal laplacian diffusion' )
+      IF( nldf == 1 .OR. nldf == 3 ) THEN      ! rotation
+         IF( .NOT.lk_ldfslp )   CALL ctl_stop( '          the rotation of the diffusive tensor require key_ldfslp' )
+#if defined key_offline
+         l_traldf_rot = .TRUE.                 ! needed for trazdf_imp
+#endif
+      ENDIF
+
+      IF( lk_esopa ) THEN
+         IF(lwp) WRITE(numout,*) '          esopa control: use all lateral physics options'
+         nldf = -1
+      ENDIF
+
+      IF( .NOT. ln_trcldf_diff ) THEN
+         IF(lwp) WRITE(numout,*) '          No lateral diffusion on passive tracers'
+         nldf = -2
+      ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         IF( nldf == -2 )   WRITE(numout,*) '          NO lateral diffusion'
+         IF( nldf == -1 )   WRITE(numout,*) '          ESOPA test All scheme used'
+         IF( nldf ==  0 )   WRITE(numout,*) '          laplacian operator'
+         IF( nldf ==  1 )   WRITE(numout,*) '          Rotated laplacian operator'
+         IF( nldf ==  2 )   WRITE(numout,*) '          bilaplacian operator'
+         IF( nldf ==  3 )   WRITE(numout,*) '          Rotated bilaplacian'
+      ENDIF
+
+      IF( ln_trcldf_bilap ) THEN
+         IF(lwp) WRITE(numout,*) '          biharmonic tracer diffusion'
+         IF( rn_ahtrc_0 > 0 .AND. .NOT. lk_esopa )   CALL ctl_stop( 'The horizontal diffusivity coef. rn_ahtrc_0 must be negative' )
+      ELSE
+         IF(lwp) WRITE(numout,*) '          harmonic tracer diffusion (default)'
+         IF( rn_ahtrc_0 < 0 .AND. .NOT. lk_esopa )   CALL ctl_stop('The horizontal diffusivity coef. rn_ahtrc_0 must be positive' )
+      ENDIF
+
+      ! ratio between active and passive tracers diffusive coef.
+      IF (ABS(rn_aht_0) < 2._wp*TINY(1.e0)) THEN
+       !  IF (ABS(rn_ahtrc_0) < 2._wp*TINY(1.e0)) THEN   !!cc
+       !     rldf_rat = 1.0_wp
+       !  ELSE
+       !     CALL ctl_stop( 'STOP', 'ldf_ctl : cannot define rldf_rat, rn_aht_0==0, rn_ahtrc_0 /=0' )
+       !  END IF
+          rldf_rat = 1.0_wp   !!cc
+      ELSE
+         rldf_rat = rn_ahtrc_0 / rn_aht_0
+      END IF
+      IF( rldf_rat < 0 ) THEN
+         IF( .NOT.lk_offline ) THEN 
+            CALL ctl_stop( 'Choose the same type of diffusive scheme both for active & passive tracers' )
+         ELSE
+            CALL ctl_stop( 'Change the sign of rn_aht_0 in namelist to -/+1' )
+         ENDIF 
+      ENDIF
+      !
+   END SUBROUTINE ldf_ctl
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                         Empty module
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_ldf_crs( kt )
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trc_ldf: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_ldf_crs
+#endif
+   !!======================================================================
+END MODULE trcldf_crs
Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc_crs.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc_crs.F90	(revision 5105)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc_crs.F90	(revision 5105)
@@ -0,0 +1,131 @@
+MODULE trcsbc_crs
+   !!==============================================================================
+   !!                       ***  MODULE  trcsbc  ***
+   !! Ocean passive tracers:  surface boundary condition
+   !!======================================================================
+   !! History :  8.2  !  1998-10  (G. Madec, G. Roullet, M. Imbard)  Original code
+   !!            8.2  !  2001-02  (D. Ludicone)  sea ice and free surface
+   !!            8.5  !  2002-06  (G. Madec)  F90: Free form and module
+   !!            9.0  !  2004-03  (C. Ethe)  adapted for passive tracers
+   !!                 !  2006-08  (C. Deltel) Diagnose ML trends for passive tracers
+   !!==============================================================================
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!   trc_sbc      : update the tracer trend at ocean surface
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! ocean dynamics and active tracers variables
+   USE trc             ! ocean  passive tracers variables
+   USE prtctl_trc      ! Print control for debbuging
+   USE trdmod_oce
+   USE trdtra
+!cbr   USE crs
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_sbc_crs   ! routine called by step.F90
+
+   !! * Substitutions
+#  include "top_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trcsbc.F90 3294 2012-01-28 16:44:18Z rblod $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE trc_sbc_crs ( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_sbc  ***
+      !!                   
+      !! ** Purpose :   Compute the tracer surface boundary condition trend of
+      !!      (concentration/dilution effect) and add it to the general 
+      !!       trend of tracer equations.
+      !!
+      !! ** Method :
+      !!      * concentration/dilution effect:
+      !!            The surface freshwater flux modify the ocean volume
+      !!         and thus the concentration of a tracer as :
+      !!            tra = tra + emp * trn / e3t   for k=1
+      !!         where emp, the surface freshwater budget (evaporation minus
+      !!         precipitation minus runoff) given in kg/m2/s is divided
+      !!         by 1035 kg/m3 (density of ocean water) to obtain m/s.
+      !!
+      !! ** Action  : - Update the 1st level of tra with the trend associated
+      !!                with the tracer surface boundary condition 
+      !!
+      !!----------------------------------------------------------------------
+      !
+      INTEGER, INTENT( in ) ::   kt          ! ocean time-step index
+      !
+      INTEGER  ::   ji, jj, jn           ! dummy loop indices
+      REAL(wp) ::   zsrau, zse3t   ! temporary scalars
+      CHARACTER (len=22) :: charout
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zsfx
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrtrd
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('trc_sbc')
+      !
+      ! Allocate temporary workspace
+                      CALL wrk_alloc( jpi, jpj,      zsfx  )
+      IF( l_trdtrc )  CALL wrk_alloc( jpi, jpj, jpk, ztrtrd )
+
+      IF( kt == nittrc000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'trc_sbc : Passive tracers surface boundary condition'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+      ENDIF
+
+     IF( .NOT. lk_offline .AND. lk_vvl ) THEN  ! online coupling with vvl
+         zsfx(:,:) = 0._wp
+      ELSE                                      ! online coupling free surface or offline with free surface
+         zsfx(:,:) = emp(:,:)
+      ENDIF
+
+      ! 0. initialization
+      zsrau = 1. / rau0
+      DO jn = 1, jptra
+         !
+         IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)  ! save trends
+         !                                             ! add the trend to the general tracer trend
+         DO jj = 2, jpj
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zse3t = 1. / fse3t(ji,jj,1)
+               tra(ji,jj,1,jn) = tra(ji,jj,1,jn) + zsfx(ji,jj) *  zsrau * trn(ji,jj,1,jn) * zse3t
+            END DO
+         END DO
+         
+         IF( l_trdtrc ) THEN
+            ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
+            CALL trd_tra( kt, 'TRC', jn, jptra_trd_nsr, ztrtrd )
+         END IF
+         !                                                       ! ===========
+      END DO                                                     ! tracer loop
+      !                                                          ! ===========
+      IF( ln_ctl )   THEN
+         WRITE(charout, FMT="('sbc ')") ;  CALL prt_ctl_trc_info(charout)
+                                           CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      ENDIF
+                      CALL wrk_dealloc( jpi, jpj,      zsfx  )
+      IF( l_trdtrc )  CALL wrk_dealloc( jpi, jpj, jpk, ztrtrd )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('trc_sbc')
+      !
+   END SUBROUTINE trc_sbc_crs
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                      NO passive tracer
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_sbc_crs (kt)              ! Empty routine
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trc_sbc: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_sbc_crs
+#endif
+   
+   !!======================================================================
+END MODULE trcsbc_crs
Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 5006)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 5105)
@@ -17,14 +17,23 @@
    USE trcnam_trp      ! passive tracers transport namelist variables
    USE trabbl          ! bottom boundary layer               (trc_bbl routine)
+   USE trabbl_crs      ! bottom boundary layer               (trc_bbl routine)
    USE trcbbl          ! bottom boundary layer               (trc_bbl routine)
+   USE trcbbl_crs      ! bottom boundary layer               (trc_bbl routine)
    USE zdfkpp          ! KPP non-local tracer fluxes         (trc_kpp routine)
    USE trcdmp          ! internal damping                    (trc_dmp routine)
    USE trcldf          ! lateral mixing                      (trc_ldf routine)
+   USE trcldf_crs      ! lateral mixing                      (trc_ldf routine)
    USE trcadv          ! advection                           (trc_adv routine)
+   USE trcadv_crs      ! advection                           (trc_adv routine)
    USE trczdf          ! vertical diffusion                  (trc_zdf routine)
+   USE trczdf_crs      ! vertical diffusion                  (trc_zdf routine
    USE trcnxt          ! time-stepping                       (trc_nxt routine)
    USE trcrad          ! positivity                          (trc_rad routine)
    USE trcsbc          ! surface boundary condition          (trc_sbc routine)
+   USE trcsbc_crs      ! surface boundary condition          (trc_sbc routine)
    USE zpshde          ! partial step: hor. derivative       (zps_hde routine)
+   USE zpshde_crs      ! partial step: hor. derivative       (zps_hde routine)
+   USE dom_oce , ONLY : ln_crs
+   USe crs, ONLY : jpi_crs,jpj_crs,wn_crs !cbr
 
 #if defined key_agrif
@@ -58,4 +67,6 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::  kstp  ! ocean time-step index
+      REAL(wp) :: zmin,zmax
+      INTEGER :: ji,jj,jk
       !! ---------------------------------------------------------------------
       !
@@ -64,10 +75,25 @@
       IF( .NOT. lk_c1d ) THEN
          !
-                                CALL trc_sbc( kstp )            ! surface boundary condition
-         IF( lk_trabbl )        CALL trc_bbl( kstp )            ! advective (and/or diffusive) bottom boundary layer scheme
+!         CALL test(kstp,1)
+!         IF( ln_crs ) THEN ;    CALL trc_sbc_crs( kstp )
+!         ELSE              ;    CALL trc_sbc( kstp )
+!         ENDIF
+!         CALL test(kstp,2)
+         IF( ln_crs ) THEN ;    CALL trc_bbl_crs( kstp )
+         ELSE              ;    CALL trc_bbl( kstp )
+         ENDIF
          IF( ln_trcdmp )        CALL trc_dmp( kstp )            ! internal damping trends
+!         CALL test(kstp,3)
+
+         IF( ln_crs ) THEN ;    CALL trc_adv_crs( kstp )
+         ELSE              ;    CALL trc_adv( kstp )
+         ENDIF
+!         CALL test(kstp,4)
+
          IF( ln_trcdmp_clo )    CALL trc_dmp_clo( kstp )        ! internal damping trends on closed seas only
-                                CALL trc_adv( kstp )            ! horizontal & vertical advection 
-                                CALL trc_ldf( kstp )            ! lateral mixing
+         IF( ln_crs ) THEN ;    CALL trc_ldf_crs( kstp )
+         ELSE              ;    CALL trc_ldf( kstp )
+         ENDIF
+!         CALL test(kstp,5)
          IF( .NOT. lk_offline .AND. lk_zdfkpp )    &
             &                   CALL trc_kpp( kstp )            ! KPP non-local tracer fluxes
@@ -75,6 +101,10 @@
          IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc           ! tracers sponge
 #endif
-                                CALL trc_zdf( kstp )            ! vertical mixing and after tracer fields
+         IF( ln_crs ) THEN ;    CALL trc_zdf_crs( kstp )
+         ELSE              ;    CALL trc_zdf( kstp )
+         ENDIF
+!         CALL test(kstp,6)
                                 CALL trc_nxt( kstp )            ! tracer fields at next time step     
+!         CALL test(kstp,7)
          IF( ln_trcrad )        CALL trc_rad( kstp )            ! Correct artificial negative concentrations
 
@@ -82,5 +112,10 @@
       IF( .NOT. Agrif_Root())   CALL Agrif_Update_Trc( kstp )   ! Update tracer at AGRIF zoom boundaries : children only
 #endif
-         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, sgtu=gtrui, sgtv=gtrvi )  ! Partial steps: now horizontal gradient of passive
+         IF( ln_zps    )        CALL zps_hde( kstp, jptra, trn, gtru, gtrv )  ! Partial steps: now horizontal gradient of passive
+         IF( ln_zps    )THEN
+         IF( ln_crs ) THEN ;    CALL zps_hde_crs( kstp, jptra, trn, gtru, gtrv )
+         ELSE              ;    CALL zps_hde( kstp, jptra, trn, gtru, gtrv )
+         ENDIF
+         ENDIF
                                                                 ! tracers at the bottom ocean level
          !
@@ -98,5 +133,35 @@
       !
    END SUBROUTINE trc_trp
+   SUBROUTINE test(kt,i)
+   INTEGER,INTENT(IN) :: kt,i
+   REAL(wp)::zmin,zmax
+   INTEGER :: ji,jj,jk
+   zmin=MINVAL( trb(2:jpi-1,2:jpj-1,:,1),mask=(tmask(2:jpi-1,2:jpj-1,:)==1)) ; CALL mpp_min(zmin)
+   zmax=MAXVAL( trb(2:jpi-1,2:jpj-1,:,1),mask=(tmask(2:jpi-1,2:jpj-1,:)==1)) ; CALL mpp_max(zmax)
+   IF(lwp)WRITE(numout,*)"trctrp b ",kt,i,zmin,zmax   
+   zmin=MINVAL( trn(2:jpi-1,2:jpj-1,:,1),mask=(tmask(2:jpi-1,2:jpj-1,:)==1)) ; CALL mpp_min(zmin)
+   zmax=MAXVAL( trn(2:jpi-1,2:jpj-1,:,1),mask=(tmask(2:jpi-1,2:jpj-1,:)==1)) ; CALL mpp_max(zmax)
+   IF(lwp)WRITE(numout,*)"trctrp n ",kt,i,zmin,zmax   
+   zmin=MINVAL( tra(2:jpi-1,2:jpj-1,:,1),mask=(tmask(2:jpi-1,2:jpj-1,:)==1)) ; CALL mpp_min(zmin)
+   zmax=MAXVAL( tra(2:jpi-1,2:jpj-1,:,1),mask=(tmask(2:jpi-1,2:jpj-1,:)==1)) ; CALL mpp_max(zmax)
+   IF(lwp)WRITE(numout,*)"trctrp a ",kt,i,zmin,zmax   
+   zmin=MINVAL( trn(2:jpi-1,2:jpj-1,1:jpk-1,1),mask=(tmask(2:jpi-1,2:jpj-1,1:jpk-1)==1)) ; CALL mpp_min(zmin)
+   zmax=MAXVAL( trn(2:jpi-1,2:jpj-1,1:jpk-1,1),mask=(tmask(2:jpi-1,2:jpj-1,1:jpk-1)==1)) ; CALL mpp_max(zmax)
+   IF(lwp)WRITE(numout,*)"trctrp n ",kt,i,zmin,zmax   
+   zmin=MINVAL( tra(2:jpi-1,2:jpj-1,1:jpk-1,1),mask=(tmask(2:jpi-1,2:jpj-1,1:jpk-1)==1)) ; CALL mpp_min(zmin)
+   zmax=MAXVAL( tra(2:jpi-1,2:jpj-1,1:jpk-1,1),mask=(tmask(2:jpi-1,2:jpj-1,1:jpk-1)==1)) ; CALL mpp_max(zmax)
+   IF(lwp)WRITE(numout,*)"trctrp a ",kt,i,zmin,zmax   
 
+   IF(narea==267)WRITE(narea+5000,*)"tra(17,5,74,1) = ",kt,i,tra(17,5,74,1)
+
+   DO ji=1,jpi
+   DO jj=1,jpj
+   DO jk=1,jpk
+      IF( tra(ji,jj,jk,1) .NE.  tra(ji,jj,jk,1) )WRITE(narea+200,*)"BUG7 ",ji,jj,jk, tra(ji,jj,jk,1); CALL FLUSH(narea+200)
+   ENDDO
+   ENDDO
+   ENDDO
+   
+   END SUBROUTINE test
 #else
    !!----------------------------------------------------------------------
Index: branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_crs.F90
===================================================================
--- branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_crs.F90	(revision 5105)
+++ branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trczdf_crs.F90	(revision 5105)
@@ -0,0 +1,182 @@
+MODULE trczdf_crs
+   !!==============================================================================
+   !!                 ***  MODULE  trczdf  ***
+   !! Ocean Passive tracers : vertical diffusive trends 
+   !!=====================================================================
+   !! History :  9.0  ! 2005-11 (G. Madec)  Original code
+   !!       NEMO 3.0  ! 2008-01  (C. Ethe, G. Madec)  merge TRC-TRA 
+   !!----------------------------------------------------------------------
+#if defined key_top
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                TOP models
+   !!----------------------------------------------------------------------
+   !!   trc_ldf     : update the tracer trend with the lateral diffusion
+   !!       ldf_ctl : initialization, namelist read, and parameters control
+   !!----------------------------------------------------------------------
+   USE oce_trc         ! ocean dynamics and active tracers
+   USE trc             ! ocean passive tracers variables
+   USE trcnam_trp      ! passive tracers transport namelist variables
+   USE trazdf_exp      ! vertical diffusion: explicit (tra_zdf_exp     routine)
+   USE trazdf_imp_crs      ! vertical diffusion: implicit (tra_zdf_imp     routine)
+   USE trdmod_oce
+   USE trdtra
+   USE prtctl_trc      ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_zdf_crs          ! called by step.F90 
+   PUBLIC   trc_zdf_alloc_crs    ! called by nemogcm.F90 
+
+   INTEGER ::   nzdf = 0               ! type vertical diffusion algorithm used
+      !                                ! defined from ln_zdf...  namlist logicals)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) ::  r2dt   ! vertical profile time-step, = 2 rdttra
+      !                                                 ! except at nittrc000 (=rdttra) if neuler=0
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "zdfddm_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 3.3 , NEMO Consortium (2010)
+   !! $Id: trczdf.F90 3294 2012-01-28 16:44:18Z rblod $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+   
+   INTEGER FUNCTION trc_zdf_alloc_crs()
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_zdf_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( r2dt(jpk) , STAT=trc_zdf_alloc_crs )
+      !
+      IF( trc_zdf_alloc_crs /= 0 )   CALL ctl_warn('trc_zdf_alloc_crs : failed to allocate array.')
+      !
+   END FUNCTION trc_zdf_alloc_crs
+
+
+   SUBROUTINE trc_zdf_crs( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_zdf  ***
+      !!
+      !! ** Purpose :   compute the vertical ocean tracer physics.
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::  kt      ! ocean time-step index
+      !
+      INTEGER               ::  jk, jn
+      CHARACTER (len=22)    :: charout
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztrtrd   ! 4D workspace
+      !!---------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('trc_zdf')
+      !
+      IF( kt == nittrc000 )   CALL zdf_ctl          ! initialisation & control of options
+
+#if ! defined key_pisces
+      IF( neuler == 0 .AND. kt == nittrc000 ) THEN     ! at nittrc000
+         r2dt(:) =  rdttrc(:)           ! = rdttrc (restarting with Euler time stepping)
+      ELSEIF( kt <= nittrc000 + 1 ) THEN          ! at nittrc000 or nittrc000+1
+         r2dt(:) = 2. * rdttrc(:)       ! = 2 rdttrc (leapfrog)
+      ENDIF
+#else
+      r2dt(:) =  rdttrc(:)              ! = rdttrc (for PISCES use Euler time stepping)
+#endif
+
+      IF( l_trdtrc )  THEN
+         CALL wrk_alloc( jpi, jpj, jpk, jptra, ztrtrd )
+         ztrtrd(:,:,:,:)  = tra(:,:,:,:)
+      ENDIF
+!      WRITE(numout,*) 'nzdf', nzdf
+      SELECT CASE ( nzdf )                       ! compute lateral mixing trend and add it to the general trend
+      CASE ( -1 )                                       ! esopa: test all possibility with control print
+         CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra ) 
+         WRITE(charout, FMT="('zdf1 ')") ;  CALL prt_ctl_trc_info(charout)
+                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+         CALL tra_zdf_imp_crs( kt, nittrc000, 'TRC', r2dt,                trb, tra, jptra ) 
+         WRITE(charout, FMT="('zdf2 ')") ;  CALL prt_ctl_trc_info(charout)
+                                            CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      CASE ( 0 ) ;  CALL tra_zdf_exp( kt, nittrc000, 'TRC', r2dt, nn_trczdf_exp, trb, tra, jptra )    !   explicit scheme 
+      CASE ( 1 ) ;  CALL tra_zdf_imp_crs( kt, nittrc000, 'TRC', r2dt,                trb, tra, jptra )    !   implicit scheme          
+
+      END SELECT
+
+      IF( l_trdtra )   THEN                      ! save the vertical diffusive trends for further diagnostics
+         DO jn = 1, jptra
+            DO jk = 1, jpkm1
+               ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dt(jk) ) - ztrtrd(:,:,jk,jn)
+            END DO
+            CALL trd_tra( kt, 'TRC', jn, jptra_trd_zdf, ztrtrd(:,:,:,jn) )
+         END DO
+         CALL wrk_dealloc( jpi, jpj, jpk, jptra, ztrtrd )
+      ENDIF
+      !                                          ! print mean trends (used for debugging)
+      IF( ln_ctl )   THEN
+         WRITE(charout, FMT="('zdf ')") ;  CALL prt_ctl_trc_info(charout)
+                                           CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=ctrcnm, clinfo2='trd' )
+      END IF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('trc_zdf')
+      !
+   END SUBROUTINE trc_zdf_crs
+
+
+   SUBROUTINE zdf_ctl
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE zdf_ctl  ***
+      !!
+      !! ** Purpose :   Choose the vertical mixing scheme
+      !!
+      !! ** Method  :   Set nzdf from ln_zdfexp
+      !!      nzdf = 0   explicit (time-splitting) scheme (ln_trczdf_exp=T)
+      !!           = 1   implicit (euler backward) scheme (ln_trczdf_exp=F)
+      !!      NB: The implicit scheme is required when using : 
+      !!             - rotated lateral mixing operator
+      !!             - TKE, GLS or KPP vertical mixing scheme
+      !!----------------------------------------------------------------------
+
+      !  Define the vertical tracer physics scheme
+      ! ==========================================
+
+      ! Choice from ln_zdfexp already read in namelist in zdfini module
+      IF( ln_trczdf_exp ) THEN           ! use explicit scheme
+         nzdf = 0
+      ELSE                               ! use implicit scheme
+         nzdf = 1
+      ENDIF
+
+      ! Force implicit schemes
+      IF( ln_trcldf_iso                               )   nzdf = 1      ! iso-neutral lateral physics
+      IF( ln_trcldf_hor .AND. ln_sco                  )   nzdf = 1      ! horizontal lateral physics in s-coordinate
+#if defined key_zdftke || defined key_zdfgls || defined key_zdfkpp
+                                                          nzdf = 1      ! TKE, GLS or KPP physics       
+#endif
+      IF( ln_trczdf_exp .AND. nzdf == 1 )   THEN
+         CALL ctl_stop( 'trc_zdf : If using the rotated lateral mixing operator or TKE, GLS or KPP vertical scheme ', &
+            &           '          the implicit scheme is required, set ln_trczdf_exp = .false.' )
+      ENDIF
+
+      ! Test: esopa
+      IF( lk_esopa )    nzdf = -1                      ! All schemes used
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'trc:zdf_ctl : vertical passive tracer physics scheme'
+         WRITE(numout,*) '~~~~~~~~~~~'
+         IF( nzdf == -1 )   WRITE(numout,*) '              ESOPA test All scheme used'
+         IF( nzdf ==  0 )   WRITE(numout,*) '              Explicit time-splitting scheme'
+         IF( nzdf ==  1 )   WRITE(numout,*) '              Implicit (euler backward) scheme'
+      ENDIF
+
+   END SUBROUTINE zdf_ctl
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                         Empty module
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE trc_zdf_crs( kt )
+      INTEGER, INTENT(in) :: kt  
+      WRITE(*,*) 'trc_zdf: You should not have seen this print! error?', kt
+   END SUBROUTINE trc_zdf_crs
+#endif
+   !!==============================================================================
+END MODULE trczdf_crs
