Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crs.F90	(revision 6101)
@@ -152,5 +152,5 @@
       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: hdivb_crs , hdivn_crs    
       REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: sshb_crs, sshn_crs , ssha_crs
-      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: rhop_crs,rhd_crs,rb2_crs
+      REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: rhop_crs,rhd_crs,rn2_crs,rb2_crs
       REAL(wp), DIMENSION(:,:)    , ALLOCATABLE      :: gru_crs, grv_crs
       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE      :: gtsu_crs, gtsv_crs
@@ -160,5 +160,5 @@
       REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE :: emp_crs, emp_b_crs, sfx_crs
       REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE :: fmmflx_crs
-      REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE :: utau_crs, vtau_crs
+      REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE :: utau_crs, vtau_crs, taum_crs
       REAL(wp), PUBLIC, DIMENSION(:,:)  , ALLOCATABLE :: rnf_crs
 
@@ -179,4 +179,6 @@
       ! Vertical diffusion
       REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avt_crs           !: vert. diffusivity coef. [m2/s] at w-point for temp  
+      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  en_crs            !: vert. diffusivity coef. [m2/s] at w-point for temp  
+      REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)    ::  avtb_2d_crs       !: vert. diffusivity coef. [m2/s] at w-point for temp  
 # if defined key_zdfddm
       REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:)  ::  avs_crs           !: salinity vertical diffusivity coeff. [m2/s] at w-point
@@ -261,5 +263,7 @@
          &      hdivb_crs(jpi_crs,jpj_crs,jpk) , hdivn_crs(jpi_crs,jpj_crs,jpk) , &
          &      rke_crs(jpi_crs,jpj_crs,jpk), rhop_crs(jpi_crs,jpj_crs,jpk)  , &
-         &      rb2_crs(jpi_crs,jpj_crs,jpk) ,rhd_crs(jpi_crs,jpj_crs,jpk)   , rab_crs_n(jpi_crs,jpj_crs,jpk,jpts) , &
+         &      rb2_crs(jpi_crs,jpj_crs,jpk) ,rn2_crs(jpi_crs,jpj_crs,jpk) , &
+         &      rhd_crs(jpi_crs,jpj_crs,jpk)   , rab_crs_n(jpi_crs,jpj_crs,jpk,jpts) , &
+         &      avtb_2d_crs(jpi_crs,jpj_crs), &
          &      gtsu_crs(jpi_crs,jpj_crs,jpts) ,gtsv_crs(jpi_crs,jpj_crs,jpts) , &
                 gru_crs(jpi_crs,jpj_crs) ,grv_crs(jpi_crs,jpj_crs) , STAT=ierr(11))
@@ -268,16 +272,19 @@
          &     emp_crs (jpi_crs,jpj_crs), emp_b_crs(jpi_crs,jpj_crs), &
          &     qsr_crs(jpi_crs ,jpj_crs), wndm_crs(jpi_crs,jpj_crs), utau_crs(jpi_crs,jpj_crs) , &
-         &     vtau_crs(jpi_crs,jpj_crs), rnf_crs(jpi_crs ,jpj_crs), &
+         &     vtau_crs(jpi_crs,jpj_crs), taum_crs(jpi_crs,jpj_crs), rnf_crs(jpi_crs ,jpj_crs), &
          &     fr_i_crs(jpi_crs,jpj_crs), sfx_crs(jpi_crs ,jpj_crs), fmmflx_crs(jpi_crs ,jpj_crs),  STAT=ierr(12)  )
 
 #if defined key_traldf_c3d
-      ALLOCATE( ahtt_crs(jpi_crs,jpj_crs,jpk) , ahtu_crs(jpi_crs,jpj_crs,jpk) , ahtv_crs(jpi_crs,jpj_crs,jpk) , ahtw_crs(jpi_crs,jpj_crs,jpk) , STAT=ierr(13) )
+      ALLOCATE( ahtt_crs(jpi_crs,jpj_crs,jpk) , ahtu_crs(jpi_crs,jpj_crs,jpk) , &
+              & ahtv_crs(jpi_crs,jpj_crs,jpk) , ahtw_crs(jpi_crs,jpj_crs,jpk) , STAT=ierr(13) )
 #elif defined key_traldf_c2d
-      ALLOCATE( ahtt_crs(jpi_crs,jpj_crs    ) , ahtu_crs(jpi_crs,jpj_crs    ) , ahtv_crs(jpi_crs,jpj_crs    ) , ahtw_crs(jpi_crs,jpj_crs    ) , STAT=ierr(13) )
+      ALLOCATE( ahtt_crs(jpi_crs,jpj_crs    ) , ahtu_crs(jpi_crs,jpj_crs    ) , &
+              & ahtv_crs(jpi_crs,jpj_crs    ) , ahtw_crs(jpi_crs,jpj_crs    ) , STAT=ierr(13) )
 #elif defined key_traldf_c1d
       ALLOCATE( ahtt_crs(        jpk) , ahtu_crs(        jpk) , ahtv_crs(        jpk) , ahtw_crs(        jpk) , STAT=ierr(13) )
 #endif
 
-     ALLOCATE( tsn_crs(jpi_crs,jpj_crs,jpk,jpts), tsb_crs(jpi_crs,jpj_crs,jpk,jpts), avt_crs(jpi_crs,jpj_crs,jpk),    &
+     ALLOCATE( tsn_crs(jpi_crs,jpj_crs,jpk,jpts), tsb_crs(jpi_crs,jpj_crs,jpk,jpts),  &
+               en_crs(jpi_crs,jpj_crs,jpk),   avt_crs(jpi_crs,jpj_crs,jpk),    &
 # if defined key_zdfddm
          &      avs_crs(jpi_crs,jpj_crs,jpk),    &
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsdom.F90	(revision 6101)
@@ -38,5 +38,5 @@
    USE crslbclnk
    USE lib_mpp
-   
+!cbr   USE ieee_arithmetic   
 
    IMPLICIT NONE
@@ -57,5 +57,4 @@
 #  include "domzgr_substitute.h90"
    
-   !! $Id$
 CONTAINS
 
@@ -67,4 +66,5 @@
       INTEGER  ::  iji, ijj
       REAL(wp) ::  zmask
+      INTEGER  :: ir,jr
       
       ! Initialize
@@ -130,4 +130,9 @@
                !write(narea+5000,*)"mask ",ijie,ijis,ijjs,ijje
                !ENDIF
+
+               ir=303-nimpp_crs+1 ; jr=302-njmpp_crs+1
+               IF( ji==ir .AND. jj==jr )THEN
+                   WRITE(narea+2000,*)"mask",ir,jr,ijis+nimpp-1,ijjs+njmpp-1
+               ENDIF
 
                !IF( ijje .GT. jpj )WRITE(narea+200,*)"BUG",jj,ijjs,ijje,SHAPE(tmask) ; call flush(narea+200)
@@ -215,4 +220,5 @@
       INTEGER :: ji, jj, jk                   ! dummy loop indices
       INTEGER :: ijis, ijjs
+      INTEGER  :: ir,jr
 
   
@@ -225,4 +231,9 @@
                   p_gphi_crs(ji,jj) = p_gphi(ijis,ijjs)
                   p_glam_crs(ji,jj) = p_glam(ijis,ijjs)
+                  ir=303-nimpp_crs+1 ; jr=302-njmpp_crs+1
+                  WRITE(narea+2000,*)"coordT1",ir,jr
+                  IF( ji==ir .AND. jj==jr )THEN
+                     WRITE(narea+2000,*)"coordT",ir,jr,ijis+nimpp-1,ijjs+njmpp-1
+                  ENDIF
                ENDDO
             ENDDO
@@ -530,6 +541,6 @@
       !! 
       !!  Arguments
-      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_fld   ! T, U, V or W on parent grid
-      CHARACTER(len=3),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN
+      REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)        :: p_fld   ! T, U, V or W on parent grid
+      CHARACTER(len=*),                         INTENT(in)           :: cd_op    ! Operation SUM, MAX or MIN
       CHARACTER(len=1),                         INTENT(in)           :: cd_type    ! grid type U,V 
       REAL(wp), DIMENSION(jpi,jpj,jpk),         INTENT(in)           :: p_mask    ! Parent grid T,U,V mask
@@ -547,6 +558,12 @@
       INTEGER  :: ii, ij, ijie, ijje, je_2
       REAL(wp) :: zflcrs, zsfcrs   
-      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf, zsurfmsk, zmask  
+      REAL(wp), DIMENSION(:,:,:), POINTER :: zsurf, zsurfmsk, zmask,ztabtmp
       INTEGER  :: iji, ijj
+      INTEGER  :: ir,jr
+      REAL(wp), DIMENSION(nn_factx,nn_facty):: ztmp
+      REAL(wp), DIMENSION(nn_factx*nn_facty):: ztmp1
+      REAL(wp), DIMENSION(:), ALLOCATABLE   :: ztmp2
+      INTEGER , DIMENSION(1)  :: zdim1
+      REAL(wp) :: zmin,zmax
       !!----------------------------------------------------------------  
    
@@ -554,5 +571,5 @@
 
       SELECT CASE ( cd_op )
-      
+ 
          CASE ( 'VOL' )
       
@@ -633,9 +650,10 @@
                              &     + p_fld(ji+1,jj+2,jk) * zsurfmsk(ji+1,jj+2,jk) &
                              &     + p_fld(ji+2,jj+2,jk) * zsurfmsk(ji+2,jj+2,jk) 
-
                            zsfcrs =  zsurf(ji,jj  ,jk) + zsurf(ji+1,jj  ,jk) + zsurf(ji+2,jj  ,jk) &
                              &     + zsurf(ji,jj+1,jk) + zsurf(ji+1,jj+1,jk) + zsurf(ji+2,jj+1,jk) &
                              &     + zsurf(ji,jj+2,jk) + zsurf(ji+1,jj+2,jk) + zsurf(ji+2,jj+2,jk) 
                             !
+!cbr                            IF( ieee_is_nan(p_fld_crs(ii,ij,jk))) THEN
+
                            p_fld_crs(ii,ij,jk) = zflcrs
                            IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,jk) = zflcrs / zsfcrs
@@ -648,5 +666,227 @@
 
               CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk )
-
+         CASE ( 'LOGVOL' )
+
+            CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk, ztabtmp )
+
+            zmin=MINVAL(p_fld) ; zmax=MAXVAL(p_fld);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"p_fld",zmin,zmax; CALL flush(numout)
+
+            ztabtmp(:,:,:)=0._wp
+            WHERE(p_fld* p_mask .NE. 0._wp ) ztabtmp =  LOG10(p_fld * p_mask)*p_mask
+            zmin=MINVAL(ztabtmp) ; zmax=MAXVAL(ztabtmp);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"LOG()",zmin,zmax; CALL flush(numout)
+            ztabtmp = ztabtmp * p_mask
+            zmin=MINVAL(ztabtmp) ; zmax=MAXVAL(ztabtmp);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"LOG()*tmask",zmin,zmax; CALL flush(numout)
+
+            SELECT CASE ( cd_type )
+
+               CASE( 'T', 'W' )
+                     DO jk = 1, jpk
+                        zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk)
+                        zsurfmsk(:,:,jk) =  zsurf(:,:,jk)
+                    ENDDO
+
+                  IF( nldj_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
+                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
+                        je_2 = mje_crs(2)
+                        DO jk = 1, jpk
+                           DO ji = nistr, niend, nn_factx
+                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2
+                              zflcrs =  ztabtmp(ji  ,je_2,jk) * zsurfmsk(ji  ,je_2,jk)   &
+                                &     + ztabtmp(ji+1,je_2,jk) * zsurfmsk(ji+1,je_2,jk)   &
+                                &     + ztabtmp(ji+2,je_2,jk) * zsurfmsk(ji+2,je_2,jk)
+
+                              zsfcrs =  zsurf(ji,je_2,jk) + zsurf(ji+1,je_2,jk) + zsurf(ji+2,je_2,jk)
+                              !
+                              p_fld_crs(ii,2,jk) = 0._wp
+                              IF( zsfcrs /= 0.0 )  p_fld_crs(ii,2,jk) = zflcrs / zsfcrs
+                              p_fld_crs(ii,2,jk) = 10 ** ( p_fld_crs(ii,2,jk) * p_mask_crs(ii,2,jk) ) * p_mask_crs(ii,2,jk)
+                           ENDDO
+                        ENDDO
+                     ENDIF
+                  ELSE
+                     je_2 = mjs_crs(2)
+                     DO jk = 1, jpk
+                        DO ji = nistr, niend, nn_factx
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2
+                           zflcrs =  ztabtmp(ji  ,je_2  ,jk) * zsurfmsk(ji  ,je_2  ,jk) &
+                             &     + ztabtmp(ji+1,je_2  ,jk) * zsurfmsk(ji+1,je_2  ,jk) &
+                             &     + ztabtmp(ji+2,je_2  ,jk) * zsurfmsk(ji+2,je_2  ,jk) &
+                             &     + ztabtmp(ji  ,je_2+1,jk) * zsurfmsk(ji  ,je_2+1,jk) &
+                             &     + ztabtmp(ji+1,je_2+1,jk) * zsurfmsk(ji+1,je_2+1,jk) &
+                             &     + ztabtmp(ji+2,je_2+1,jk) * zsurfmsk(ji+2,je_2+1,jk) &
+                             &     + ztabtmp(ji  ,je_2+2,jk) * zsurfmsk(ji  ,je_2+2,jk) &
+                             &     + ztabtmp(ji+1,je_2+2,jk) * zsurfmsk(ji+1,je_2+2,jk) &
+                             &     + ztabtmp(ji+2,je_2+2,jk) * zsurfmsk(ji+2,je_2+2,jk)
+
+                           zsfcrs =  zsurf(ji,je_2  ,jk) + zsurf(ji+1,je_2  ,jk) + zsurf(ji+2,je_2  ,jk) &
+                             &     + zsurf(ji,je_2+1,jk) + zsurf(ji+1,je_2+1,jk) + zsurf(ji+2,je_2+1,jk) &
+                             &     + zsurf(ji,je_2+2,jk) + zsurf(ji+1,je_2+2,jk) + zsurf(ji+2,je_2+2,jk)
+                            !
+                            p_fld_crs(ii,2,jk) = 0._wp
+                            IF( zsfcrs /= 0.0 )  p_fld_crs(ii,2,jk) = zflcrs / zsfcrs
+                            p_fld_crs(ii,2,jk) = 10 ** ( p_fld_crs(ii,2,jk) * p_mask_crs(ii,2,jk) ) * p_mask_crs(ii,2,jk)
+                        ENDDO
+                     ENDDO
+                  ENDIF
+                  !
+                  DO jk = 1, jpk
+                     DO jj  = njstr, njend, nn_facty
+                        DO ji = nistr, niend, nn_factx
+                           ii = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                           ij = ( jj - njstr ) * rfacty_r + 3
+                           zflcrs =  ztabtmp(ji  ,jj  ,jk) * zsurfmsk(ji  ,jj  ,jk) &
+                             &     + ztabtmp(ji+1,jj  ,jk) * zsurfmsk(ji+1,jj  ,jk) &
+                             &     + ztabtmp(ji+2,jj  ,jk) * zsurfmsk(ji+2,jj  ,jk) &
+                             &     + ztabtmp(ji  ,jj+1,jk) * zsurfmsk(ji  ,jj+1,jk) &
+                             &     + ztabtmp(ji+1,jj+1,jk) * zsurfmsk(ji+1,jj+1,jk) &
+                             &     + ztabtmp(ji+2,jj+1,jk) * zsurfmsk(ji+2,jj+1,jk) &
+                             &     + ztabtmp(ji  ,jj+2,jk) * zsurfmsk(ji  ,jj+2,jk) &
+                             &     + ztabtmp(ji+1,jj+2,jk) * zsurfmsk(ji+1,jj+2,jk) &
+                             &     + ztabtmp(ji+2,jj+2,jk) * zsurfmsk(ji+2,jj+2,jk)
+                           zsfcrs =  zsurf(ji,jj  ,jk) + zsurf(ji+1,jj  ,jk) + zsurf(ji+2,jj  ,jk) &
+                             &     + zsurf(ji,jj+1,jk) + zsurf(ji+1,jj+1,jk) + zsurf(ji+2,jj+1,jk) &
+                             &     + zsurf(ji,jj+2,jk) + zsurf(ji+1,jj+2,jk) + zsurf(ji+2,jj+2,jk)
+                            !
+                           p_fld_crs(ii,ij,jk) = 0._wp
+                           IF( zsfcrs /= 0.0 )  p_fld_crs(ii,ij,jk) = zflcrs / zsfcrs
+                           p_fld_crs(ii,ij,jk) = 10 ** ( p_fld_crs(ii,ij,jk) *  p_mask_crs(ii,ij,jk) ) * p_mask_crs(ii,ij,jk)
+                        ENDDO
+                     ENDDO
+                  ENDDO
+               CASE DEFAULT
+                    STOP
+               END SELECT
+
+
+              !WHERE( p_fld .NE. 0._wp ) p_fld=10**(p_fld)
+              !zmin=MINVAL(p_fld) ; zmax=MAXVAL(p_fld);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld)",zmin,zmax ; CALL flush(numout)
+              !p_fld = p_fld * p_mask
+              !zmin=MINVAL(p_fld) ; zmax=MAXVAL(p_fld);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld)*tmask",zmin,zmax ; CALL flush(numout)
+
+              zmin=MINVAL(p_fld_crs) ; zmax=MAXVAL(p_fld_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"p_fld_crs",zmin,zmax; CALL flush(numout)
+              !p_fld_crs=10**(p_fld_crs*p_mask_crs)
+              !zmin=MINVAL(p_fld_crs) ; zmax=MAXVAL(p_fld_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld_crs)",zmin,zmax; CALL flush(numout)
+              !p_fld_crs=p_fld_crs*p_mask_crs
+              !zmin=MINVAL(p_fld_crs) ; zmax=MAXVAL(p_fld_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"10**(p_fld_crs)*tmask",zmin,zmax; CALL flush(numout)
+
+              CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk ,ztabtmp )
+         CASE ( 'MED' )
+
+            CALL wrk_alloc( jpi, jpj, jpk, zsurf, zsurfmsk )
+
+            SELECT CASE ( cd_type )
+
+               CASE( 'T', 'W' )
+                     DO jk = 1, jpk
+                        zsurf   (:,:,jk) =  p_e12(:,:) * p_e3(:,:,jk) *  p_mask(:,:,jk)
+                        zsurfmsk(:,:,jk) =  zsurf(:,:,jk)
+                    ENDDO
+
+                  IF( nldj_crs == 1 .AND. ( ( mje_crs(2) - mjs_crs(2) ) < 2 ) ) THEN     !!cc bande du sud style ORCA2
+                     IF( mje_crs(2) - mjs_crs(2) == 1 ) THEN
+                        je_2 = mje_crs(2)
+                        DO jk = 1, jpk
+                           DO ji = nistr, niend, nn_factx
+                              ii   = ( ji - mis_crs(2) ) * rfactx_r + 2
+
+                              ztmp1(:) = 0._wp
+                              ztmp1(1:3) =  p_fld(ji:ji+2,je_2,jk)
+                              CALL PIKSRT(nn_factx*nn_facty,ztmp1)
+                              ir=0
+                              jr=1
+                              DO WHILE( jr .LE. nn_factx*nn_facty )
+                                 IF( ztmp1(jr) == 0. )THEN
+                                    ir=jr
+                                    jr=jr+1
+                                 ELSE
+                                    EXIT
+                                 ENDIF
+                              ENDDO
+                              IF( ir .LE. nn_factx*nn_facty-1 )THEN
+                                 ALLOCATE( ztmp2(nn_factx*nn_facty-ir) )
+                                 ztmp2(1:nn_factx*nn_facty-ir) = ztmp1(1+ir:nn_factx*nn_facty)
+                                 jr=INT( 0.5 * REAL(nn_factx*nn_facty-ir,wp) )+1
+                                 p_fld_crs(ii,2,jk) = ztmp2(jr)
+                                 DEALLOCATE( ztmp2 )
+                              ELSE
+                                 p_fld_crs(ii,ij,jk) = 0._wp
+                              ENDIF
+
+                           ENDDO
+                        ENDDO
+                     ENDIF
+                  ELSE
+                     je_2 = mjs_crs(2)
+                     DO jk = 1, jpk
+                        DO ji = nistr, niend, nn_factx
+                           ii   = ( ji - mis_crs(2) ) * rfactx_r + 2
+                           
+                           ztmp(:,:)= p_fld(ji:ji+2,je_2:je_2+2,jk)
+                           zdim1(1)=nn_factx*nn_facty
+                           ztmp1(:) = RESHAPE( ztmp(:,:) , zdim1 )
+                           CALL PIKSRT(nn_factx*nn_facty,ztmp1)
+                           ir=0
+                           jr=1
+                           DO WHILE( jr .LE. nn_factx*nn_facty )
+                              IF( ztmp1(jr) == 0. ) THEN
+                                 ir=jr
+                                 jr=jr+1
+                              ELSE
+                                 EXIT
+                              ENDIF
+                           ENDDO
+                           IF( ir .LE. nn_factx*nn_facty-1 )THEN
+                              ALLOCATE( ztmp2(nn_factx*nn_facty-ir) )
+                              ztmp2(1:nn_factx*nn_facty-ir) = ztmp1(1+ir:nn_factx*nn_facty)
+                              jr=INT( 0.5 * REAL(nn_factx*nn_facty-ir,wp) )+1
+                              p_fld_crs(ii,2,jk) = ztmp2(jr)
+                              DEALLOCATE( ztmp2 )
+                           ELSE
+                           p_fld_crs(ii,ij,jk) = 0._wp
+                           ENDIF
+
+                        ENDDO
+                     ENDDO
+                  ENDIF
+                  !
+                  DO jk = 1, jpk
+                     DO jj  = njstr, njend, nn_facty
+                        DO ji = nistr, niend, nn_factx
+                           ii = ( ji - mis_crs(2) ) * rfactx_r + 2                 ! cordinate in parent grid
+                           ij = ( jj - njstr ) * rfacty_r + 3
+
+                           ztmp(:,:)= p_fld(ji:ji+2,jj:jj+2,jk) 
+                           zdim1(1)=nn_factx*nn_facty
+                           ztmp1(:) = RESHAPE( ztmp(:,:) , zdim1 )
+                           CALL PIKSRT(nn_factx*nn_facty,ztmp1)
+                           ir=0
+                           jr=1
+                           DO WHILE( jr .LE. nn_factx*nn_facty )
+                              IF( ztmp1(jr) == 0. ) THEN
+                                 ir=jr
+                                 jr=jr+1
+                              ELSE
+                                 EXIT
+                              ENDIF
+                           ENDDO
+                           IF( ir .LE. nn_factx*nn_facty-1 )THEN
+                              ALLOCATE( ztmp2(nn_factx*nn_facty-ir) )
+                              ztmp2(1:nn_factx*nn_facty-ir) = ztmp1(1+ir:nn_factx*nn_facty)
+                              jr=INT( 0.5 * REAL(nn_factx*nn_facty-ir,wp) )+1
+                              p_fld_crs(ii,ij,jk) = ztmp2(jr)
+                              DEALLOCATE( ztmp2 )
+                           ELSE
+                              p_fld_crs(ii,ij,jk) = 0._wp
+                           ENDIF
+
+                        ENDDO
+                     ENDDO
+                  ENDDO
+               CASE DEFAULT
+                    STOP
+               END SELECT
+
+              CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk )
+ 
          CASE ( 'SUM' )
          
@@ -2390,5 +2630,5 @@
       CALL crs_lbc_lnk( p_surf_crs_msk, cd_type, 1.0, pval=1.0 )
 
-      CALL wrk_dealloc( jpi, jpj, jpk, zsurfmsk, zsurf )
+      CALL wrk_dealloc( jpi, jpj, jpk, zsurf, zsurfmsk )
 
    END SUBROUTINE crs_dom_sfc
@@ -2893,4 +3133,6 @@
       ENDDO
      
+      CALL wrk_alloc( jpi_crs, jpj_crs, zmbk )
+
       zmbk(:,:) = 0.0
       zmbk(:,:) = REAL( mbathy_crs(:,:), wp ) ;   CALL crs_lbc_lnk(zmbk,'T',1.0)   ;   mbathy_crs(:,:) = INT( zmbk(:,:) )
@@ -2921,4 +3163,25 @@
    END SUBROUTINE crs_dom_bat
 
+   SUBROUTINE PIKSRT(N,ARR)
+     INTEGER                  ,INTENT(IN) :: N
+     REAL(kind=8),DIMENSION(N),INTENT(INOUT) :: ARR
+
+     INTEGER      :: i,j
+     REAL(kind=8) :: a
+     !!----------------------------------------------------------------
+
+     DO j=2, N
+       a=ARR(j)
+       DO i=j-1,1,-1
+          IF(ARR(i)<=a) goto 10
+          ARR(i+1)=ARR(i)
+       ENDDO
+       i=0
+10     ARR(i+1)=a
+     ENDDO
+     RETURN
+
+   END SUBROUTINE PIKSRT
+
 
 END MODULE crsdom
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsfld.F90	(revision 6101)
@@ -27,4 +27,10 @@
    USE iom
    USE zdfmxl_crs
+   USE eosbn2
+   USE zdfevd_crs
+   USE zdftke
+   USE zdftke_crs
+
+!   USE ieee_arithmetic
 
    IMPLICIT NONE
@@ -40,5 +46,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
+   !! $Id $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -65,12 +71,18 @@
       !!
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zfse3t, zfse3u, zfse3v, zfse3w ! 3D workspace for e3
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zs 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt, zs , ztmp
       REAL(wp), POINTER, DIMENSION(:,:)   :: z2d,z2d_crs
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_crs, zs_crs !
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zt_crs, zs_crs, zerr_crs,zmax_crs
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ztmp_crs
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) :: avte_crs
       REAL(wp)       :: z2dcrsu, z2dcrsv
-      REAL(wp)       :: zmin,zmax
+      REAL(wp)       :: zmin,zmax,icnt1,icnt2
       INTEGER :: i,j,ijis,ijie,ijjs,ijje
       REAL(wp)       :: zw,zwp1,zum1,zu,zvm1,zv,zsnm,zsm,z
+      REAL(wp)       :: zerr, zerr0, zerr1, zmean
+      INTEGER,DIMENSION(4,3) :: ind
+      REAL(wp),DIMENSION(4) :: zwgt
       INTEGER ::  iji,ijj
+      INTEGER :: jl,jm,jn
       !!
       !!----------------------------------------------------------------------
@@ -81,8 +93,10 @@
       CALL wrk_alloc( jpi, jpj, jpk, zfse3t, zfse3w )
       CALL wrk_alloc( jpi, jpj, jpk, zfse3u, zfse3v )
-      CALL wrk_alloc( jpi, jpj, jpk, zt, zs         )
+      CALL wrk_alloc( jpi, jpj, jpk, zt, zs , ztmp        )
       CALL wrk_alloc( jpi, jpj,      z2d            )
       !
-      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs )
+      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs, zerr_crs,zmax_crs )
+      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, ztmp_crs )
+      CALL wrk_alloc( jpi_crs, jpj_crs, jpk, 4, avte_crs )
       CALL wrk_alloc( jpi_crs, jpj_crs, z2d_crs     )
 
@@ -129,5 +143,10 @@
       CALL iom_put( "sst" , tsn_crs(:,:,1,jp_tem) )    ! sst
 
-      
+      !n2 before
+      zt(:,:,:) = rn2b(:,:,:)  ;      zt_crs(:,:,:) = 0._wp
+      CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
+      rb2_crs(:,:,:) = zt_crs(:,:,:)
+      CALL iom_put("rb2_crs",rb2_crs)
+
       !  Salinity
       zs(:,:,:) = tsb(:,:,:,jp_sal)  ;      zs_crs(:,:,:) = 0._wp
@@ -252,4 +271,148 @@
          CASE ( 2 )
             CALL crs_dom_ope( avt, 'MIN', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
+         CASE ( 3 )
+            CALL crs_dom_ope( avt, 'LOGVOL', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, p_mask_crs=tmask_crs, psgn=1.0 )
+         CASE ( 4 )
+            CALL crs_dom_ope( avt, 'MED', 'W', tmask, avt_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
+         CASE ( 5 )
+            CALL crs_dom_ope( en , 'VOL', 'W', tmask, en_crs, p_e12=e1e2t, p_e3=zfse3w, psgn=1.0 )
+            CALL crs_dom_ope( taum , 'SUM', 'T', tmask, taum_crs , p_e12=e1e2t, p_surf_crs=e1e2t_crs, psgn=1.0 )
+            CALL crs_dom_ope( rn2(:,:,:), 'VOL', 'W', tmask, rn2_crs, p_e12=e1e2t, p_e3=zfse3t, psgn=1.0 )
+            IF( kt==nit000 )CALL tke_avn_ini_crs
+            CALL tke_avn_crs
+            CALL zdf_evd_crs(kt)
+         CASE ( 6 )
+
+            avte_crs(:,:,:,:) = 0._wp
+            ztmp(:,:,:)=1.
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk 
+               WHERE( fse3w(:,:,jk) .NE. 0._wp) zs(:,:,jk) = ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) / fse3w(:,:,jk)
+               zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,1) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,1));zmax=MAXVAL(avte_crs(:,:,:,1));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax 
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk 
+               WHERE( fse3w(:,:,jk) .NE. 0._wp) zs(:,:,jk) = ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) / fse3w(:,:,jk)
+               zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte1_crs",zmin,zmax
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,2) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,2));zmax=MAXVAL(avte_crs(:,:,:,2));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte2_crs",zmin,zmax
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk
+                WHERE( fse3w(:,:,jk) .NE. 0._wp ) zs(:,:,jk)= ( rn_a0 * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) +  &
+                                                                  &   rn_b0 * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) )/ fse3w(:,:,jk) 
+                zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)                                                           
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,3) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,3));zmax=MAXVAL(avte_crs(:,:,:,3));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte3_crs",zmin,zmax
+
+            zt(:,:,:) = 0._wp
+            zs(:,:,:) = 0._wp
+            DO jk=2,jpk
+                WHERE( fse3w(:,:,jk) .NE. 0._wp ) zs(:,:,jk)= ( rn_a0 * ( tsn(:,:,jk-1,jp_tem) - tsn(:,:,jk,jp_tem) ) -  &
+                                                                  &   rn_b0 * ( tsn(:,:,jk-1,jp_sal) - tsn(:,:,jk,jp_sal) ) )/ fse3w(:,:,jk) 
+                zt(:,:,jk)=  avt(:,:,jk) *  zs(:,:,jk)
+            ENDDO
+            CALL crs_dom_ope( zt, 'VOL', 'W', tmask, zt_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            CALL crs_dom_ope( zs, 'VOL', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=ztmp, psgn=1.0 )
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax 
+            zmin=MINVAL(zs_crs);zmax=MAXVAL(zs_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax 
+            zt_crs=tmask_crs*zt_crs
+            zmin=MINVAL(zt_crs);zmax=MAXVAL(zt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax
+            WHERE( zs_crs .NE. 0._wp ) avte_crs(:,:,:,4) = zt_crs / zs_crs
+            zmin=MINVAL(avte_crs(:,:,:,4));zmax=MAXVAL(avte_crs(:,:,:,4));CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avte4_crs",zmin,zmax
+
+            CALL iom_put( "avte_crs1", avte_crs(:,:,:,1) )   !  Kz
+            CALL iom_put( "avte_crs2", avte_crs(:,:,:,2) )   !  Kz
+            CALL iom_put( "avte_crs3", avte_crs(:,:,:,3) )   !  Kz
+            CALL iom_put( "avte_crs4", avte_crs(:,:,:,4) )   !  Kz
+!--------------------- 
+            CALL crs_dom_ope( avt, 'MED', 'W', tmask, zs_crs, p_e12=e1e2t, p_e3=zfse3w, p_mask_crs=tmask_crs, psgn=1.0 )
+!?            zmin=MINVAL(zs_crs*tmask_crs);zmax=MAXVAL(zs_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"logvol zs_crs*tmask ",zmin,zmax ; call flush(numout)
+            CALL iom_put( "zs_crs", zs_crs )   !  Kzlogvol
+!--------------------- ok
+
+            CALL crs_dom_ope( avt, 'VOL', 'W', tmask, zmax_crs, p_e12=e1e2t, p_e3=zfse3w,  psgn=1.0 )
+            WRITE(narea+200,*)"zmax_crs ",SHAPE(zmax_crs) ; call flush(narea+200)
+            CALL iom_put( "zmax_crs", zmax_crs )   !  Kzlogvol
+            zmin=MINVAL(zmax_crs*tmask_crs);zmax=MAXVAL(zmax_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"vol zmax_crs*tmask ",zmin,zmax ; call flush(numout)
+!-------------------------nok
+            avt_crs=zs_crs
+
+
+            zerr0=0.01
+
+            icnt1=0
+            icnt2=0
+
+            zt_crs(:,:,:)=0._wp
+            zerr_crs(:,:,:)=0._wp
+            DO ji=1,jpi_crs 
+            DO jj=1,jpj_crs 
+            DO jk=1,jpk
+
+ 
+!--------------
+               zwgt(1:4)=0._wp
+               DO jm=1,4 ; IF( avte_crs(ji,jj,jk,jm)  .GE. 0._wp .AND.  avte_crs(ji,jj,jk,jm)  .LE. zmax_crs(ji,jj,jk)  ) zwgt(jm) = 1._wp ; ENDDO
+!--------------
+               IF( SUM(zwgt(1:4)) .NE. 0._wp ) THEN  
+                  zmean = SUM( zwgt(1:4)*avte_crs(ji,jj,jk,1:4)) / SUM(zwgt(1:4) )
+                  zerr  = SQRT(SUM( zwgt(1:4)*(avte_crs(ji,jj,jk,1:4)-zmean)*(avte_crs(ji,jj,jk,1:4)-zmean) ) / SUM(zwgt(1:4) ) )
+               ELSE
+                  zmean=0._wp
+                 zerr=1.e10
+               ENDIF
+!--------------
+
+               zerr_crs(ji,jj,jk)=zerr
+
+               IF( zerr .LE. zerr0 .AND. zmean .GT. 0._wp )zt_crs(ji,jj,jk)=zmean
+               IF( zerr .LE. zerr0 .AND. zmean .GT. 0._wp )avt_crs(ji,jj,jk)=zmean
+
+               IF( tmask_crs(ji,jj,jk) == 1 ) icnt1=icnt1+1
+               IF( tmask_crs(ji,jj,jk) == 1 .AND.  zerr .LE. zerr0 .AND. zmean .GT. 0._wp ) icnt2=icnt2+1
+
+               !IF( ieee_is_nan(  zt_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANMEANEFF ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
+               !IF( ieee_is_nan(  zs_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANLOG ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
+               !IF( ieee_is_nan( avt_crs(ji,jj,jk))   )WRITE(narea+200,*)"NANAVT ",ji,jj,jk,tmask_crs(ji,jj,jk)  ; call flush(narea+200)
+            ENDDO
+            ENDDO
+            ENDDO
+            zmin=MINVAL(avt_crs);zmax=MAXVAL(avt_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avt_crs ",zmin,zmax  ; call flush(numout)
+            zmin=MINVAL(avt_crs*tmask_crs);zmax=MAXVAL(avt_crs*tmask_crs);CALL mpp_min(zmin);CALL mpp_max(zmax);IF(lwp)WRITE(numout,*)"avt_crs*tmask ",zmin,zmax  ; call flush(numout)
+
+            CALL mpp_sum(icnt1)
+            CALL mpp_sum(icnt2)
+            IF(lwp)WRITE(numout,*)"TOTO",kt,icnt1,icnt2
+            CALL iom_put( "zt_crs", zt_crs )   !  Kz
+            CALL iom_put( "zerr_crs", zerr_crs )   !  Kz
+
       END SELECT
       !
@@ -293,7 +456,9 @@
       CALL wrk_dealloc( jpi, jpj, jpk, zfse3t, zfse3w )
       CALL wrk_dealloc( jpi, jpj, jpk, zfse3u, zfse3v )
-      CALL wrk_dealloc( jpi, jpj, jpk, zt, zs         )
+      CALL wrk_dealloc( jpi, jpj, jpk, zt, zs, ztmp   )
       CALL wrk_dealloc( jpi, jpj, z2d                 )
-      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs )
+      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, zt_crs, zs_crs, zerr_crs,zmax_crs )
+      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, ztmp_crs )
+      CALL wrk_dealloc( jpi_crs, jpj_crs, jpk, 4, avte_crs )
       CALL wrk_dealloc( jpi_crs, jpj_crs, z2d_crs     )
       !
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/CRS/crsini.F90	(revision 6101)
@@ -256,5 +256,4 @@
      CALL crs_dom_ope( gdepw_0, 'MAX', 'W', tmask, gdepw_crs, p_e3=zfse3w, psgn=1.0 )
 
-
      !---------------------------------------------------------
      ! 4.  Coarse grid ocean volume and averaging weights
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 6101)
@@ -22,5 +22,6 @@
                        gphiv, gphif, &
                        agrif_root, agrif_cfixed, lk_agrif, &
-                       rdt,rdttra, gdept_0, ln_crs, gdepw_0,  adatrj, fjulday
+                       rdt,rdttra, gdept_0, ln_crs, gdepw_0,  adatrj, fjulday, &
+                       mikt
 
    USE c1d             ! 1D vertical configuration
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp_crs.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp_crs.F90	(revision 6101)
@@ -493,6 +493,9 @@
          !
       ELSE                             ! Madec operator : slopes at u-, v-, and w-points
-         ALLOCATE( uslp_crs(jpi_crs,jpj_crs,jpk) , vslp_crs(jpi_crs,jpj_crs,jpk) , wslpi_crs(jpi_crs,jpj_crs,jpk) , wslpj_crs(jpi_crs,jpj_crs,jpk) ,  &
-            &   omlmask(jpi_crs,jpj_crs,jpk) , uslpml(jpi_crs,jpj_crs)   , vslpml(jpi_crs,jpj_crs)    , wslpiml(jpi_crs,jpj_crs)   , wslpjml(jpi_crs,jpj_crs) , STAT=ierr )
+         ALLOCATE( uslp_crs(jpi_crs,jpj_crs,jpk) , vslp_crs(jpi_crs,jpj_crs,jpk) , & 
+                 & wslpi_crs(jpi_crs,jpj_crs,jpk) , wslpj_crs(jpi_crs,jpj_crs,jpk) ,  &
+                 & omlmask(jpi_crs,jpj_crs,jpk) ,  &
+                 & uslpml(jpi_crs,jpj_crs)   ,  vslpml(jpi_crs,jpj_crs)  , & 
+                 & wslpiml(jpi_crs,jpj_crs)   , wslpjml(jpi_crs,jpj_crs) , STAT=ierr )
          IF( ierr > 0 )   CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' )
 
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c2d_crs.h90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c2d_crs.h90	(revision 6101)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_c2d_crs.h90	(revision 6101)
@@ -0,0 +1,166 @@
+   !!----------------------------------------------------------------------
+   !!                      ***  ldftra_c2d.h90  ***
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: ldftra_c2d.h90 4147 2013-11-04 11:51:55Z cetlod $ 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE ldf_tra_c2d_crs( ld_print )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldftra_c2d  ***
+      !!              
+      !! ** Purpose :   initializations of horizontally non uniform eddy 
+      !!      diffusivity coefficients
+      !!
+      !! ** Method :
+      !!       biharmonic operator    : ahtt = defined at T-level
+      !!                                ahtu,ahtv,ahtw never used
+      !!       harmonic operator (ahtt never used)
+      !!           iso-model level   : ahtu, ahtv defined at u-, v-points
+      !!         isopycnal         : ahtu, ahtv, ahtw defined at u-, v-, w-pts
+      !!         or geopotential   
+      !!       eddy induced velocity
+      !!           always harmonic   : aeiu, aeiv, aeiw defined at u-, v-, w-pts
+      !!----------------------------------------------------------------------
+      LOGICAL, INTENT (in) ::   ld_print   ! If true, print arrays in numout
+      !
+      INTEGER ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   za00, zd_max, zeumax, zevmax, zetmax
+      !!----------------------------------------------------------------------
+
+      IF( lk_traldf_eiv ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' ldf_tra_c2d : 2D eddy diffusivity and eddy'
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --  induced velocity coefficients'
+      ELSE
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) ' ldf_tra2d : 2D eddy diffusivity coefficient'
+         IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~   --'
+      ENDIF
+
+      zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+      IF( lk_mpp ) CALL mpp_max( zd_max )   ! max over the global domain
+
+      ! harmonic operator : (U-, V-, W-points)
+      ! ==================
+      IF( ln_traldf_lap ) THEN
+         !
+         za00 = aht0 / zd_max
+         !
+         DO jj = 1, jpj 
+            DO ji = 1, jpi 
+               zeumax = MAX( e1u(ji,jj), e2u(ji,jj) ) 
+               zevmax = MAX( e1v(ji,jj), e2v(ji,jj) ) 
+               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
+               ahtu(ji,jj) = za00 * zeumax ! set ahtu = ahtv at u- and v-points, 
+               ahtv(ji,jj) = za00 * zevmax ! and ahtw at w-point (idem T-point) 
+               ahtw(ji,jj) = za00 * zetmax ! 
+            END DO
+         END DO
+
+         CALL crs_lbc_lnk( ahtu, 'U', 1. )   ! Lateral boundary conditions
+         CALL crs_lbc_lnk( ahtv, 'V', 1. )   ! (no change of sign)
+         CALL crs_lbc_lnk( ahtw, 'W', 1. )
+         ! Special case for ORCA R2 and R4 configurations (overwrite the value of ahtu ahtv ahtw)
+         ! ==============================================
+         IF( cp_cfg == "orca" .AND. ( jp_cfg == 2 .OR. jp_cfg == 4 ) )   THEN
+            ahtu(:,:) = aht0              ! set ahtu = ahtv at u- and v-points,
+            ahtv(:,:) = aht0              ! and ahtw at w-point
+            ahtw(:,:) = aht0              ! (here : no space variation)
+            IF(lwp) WRITE(numout,*) '               ORCA R2 or R4 case'
+            IF(lwp) WRITE(numout,*) '               Constant values used for eddy diffusivity coefficients'
+            IF(lwp) WRITE(numout,*) '               Variation lat/lon only for eddy induced velocity coefficients'
+         ENDIF
+
+         ! Control print
+         IF( lwp .AND. ld_print ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'inildf: ahtu array'
+            CALL prihre( ahtu, jpi, jpj, 1, jpi, 1,   &
+               &                         1, jpj, 1, 1.e-3, numout )
+            WRITE(numout,*)
+            WRITE(numout,*) 'inildf: ahtv array'
+            CALL prihre( ahtv, jpi, jpj, 1, jpi, 1,   &
+               &                         1, jpj, 1, 1.e-3, numout )
+            WRITE(numout,*)
+            WRITE(numout,*) 'inildf: ahtw array'
+            CALL prihre( ahtw, jpi, jpj, 1, jpi, 1,   &
+               &                         1, jpj, 1, 1.e-3, numout )
+         ENDIF
+      ENDIF
+      
+      ! biharmonic operator : (T-point)
+      ! ====================
+      IF( ln_traldf_bilap ) THEN
+         ! (USER: modify ahtt following your desiderata)
+         ! Here: ahm is proportional to the cube of the maximum of the gridspacing
+         !       in the to horizontal direction
+
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
+
+         za00 = aht0 / ( zd_max * zd_max * zd_max )
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zetmax = MAX( e1t(ji,jj), e2t(ji,jj) )
+               ahtt(ji,jj) = za00 * zetmax * zetmax * zetmax      ! set ahtt at T-point
+            END DO
+         END DO
+
+         CALL crs_lbc_lnk( ahtt, 'T', 1. )   ! Lateral boundary conditions on ( ahtt )
+         ! Control print
+         IF( lwp .AND. ld_print ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'inildf: 2D ahtt array'
+            CALL prihre( ahtt, jpi, jpj, 1, jpi, 1,   &
+               &                         1, jpj, 1, 1.e-3, numout )
+         ENDIF
+      ENDIF
+
+
+# if defined key_traldf_eiv
+      ! set aeiu = aeiv at u- and v-points, and aeiw at w-point (idem T-point)
+      ! (here no space variation)
+      aeiu(:,:) = aeiv0
+      aeiv(:,:) = aeiv0
+      aeiw(:,:) = aeiv0
+      
+      IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN
+         !                                 ! Cancel eiv in Gibraltar strait
+         aeiu( mi0(68):mi1(71) , mj0(50):mj1(53) ) = 0.e0
+         aeiv( mi0(68):mi1(71) , mj0(50):mj1(53) ) = 0.e0
+         aeiw( mi0(68):mi1(71) , mj0(50):mj1(53) ) = 0.e0
+         !                                 ! Cancel eiv in Mediterrannean sea
+         aeiu( mi0(70):mi1(90) , mj0(49):mj1(56) ) = 0.e0
+         aeiv( mi0(70):mi1(90) , mj0(49):mj1(56) ) = 0.e0
+         aeiw( mi0(70):mi1(90) , mj0(49):mj1(56) ) = 0.e0
+      ENDIF
+
+      ! Lateral boundary conditions on ( aeiu, aeiv, aeiw )
+#if defined key_crs
+      CALL crs_lbc_lnk( aeiu, 'U', 1. )
+      CALL crs_lbc_lnk( aeiu, 'U', 1. )
+      CALL crs_lbc_lnk( aeiv, 'V', 1. )
+#else
+      CALL lbc_lnk( aeiw, 'W', 1. )
+      CALL lbc_lnk( aeiv, 'V', 1. )
+      CALL lbc_lnk( aeiw, 'W', 1. )
+#endif
+      ! Control print
+      IF( lwp .AND. ld_print ) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'inildf: aeiu array'
+         CALL prihre(aeiu,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'inildf: aeiv array'
+         CALL prihre(aeiv,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+         WRITE(numout,*)
+         WRITE(numout,*) 'inildf: aeiw array'
+         CALL prihre(aeiw,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout)
+      ENDIF
+# endif
+      !
+   END SUBROUTINE ldf_tra_c2d_crs
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_crs.F90	(revision 6101)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/LDF/ldftra_crs.F90	(revision 6101)
@@ -0,0 +1,90 @@
+MODULE ldftra_crs
+   !!======================================================================
+   !!                       ***  MODULE  ldftra  ***
+   !! Ocean physics:  lateral diffusivity coefficient 
+   !!=====================================================================
+   !! History :        ! 1997-07  (G. Madec)  from inimix.F split in 2 routines
+   !!   NEMO      1.0  ! 2002-09  (G. Madec)  F90: Free form and module
+   !!             2.0  ! 2005-11  (G. Madec)  
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   ldf_tra_init : initialization, namelist read, and parameters control
+   !!   ldf_tra_c3d   : 3D eddy viscosity coefficient initialization
+   !!   ldf_tra_c2d   : 2D eddy viscosity coefficient initialization
+   !!   ldf_tra_c1d   : 1D eddy viscosity coefficient initialization
+   !!----------------------------------------------------------------------
+!   USE oce             ! ocean dynamics and tracers
+!   USE dom_oce         ! ocean space and time domain
+!   USE phycst          ! physical constants
+   USE ldftra_oce,ONLY : lk_traldf_eiv,ln_traldf_lap,ln_traldf_bilap     ! ocean tracer   lateral physics
+!   USE ldfslp          ! ???
+!   USE in_out_manager  ! I/O manager
+!   USE ioipsl
+!   USE lib_mpp         ! distribued memory computing library
+!   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE oce_trc
+   USE crslbclnk , ONLY : crs_lbc_lnk
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ldf_tra_crs_init   ! called by opa.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: ldftra.F90 4624 2014-04-28 12:09:03Z acc $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ldf_tra_crs_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ldf_tra_init  ***
+      !! 
+      !! ** Purpose :   initializations of the tracer lateral mixing coeff.
+      !!
+      !! ** Method  :   the Eddy diffusivity and eddy induced velocity ceoff.
+      !!      are defined as follows:
+      !!         default option   : constant coef. aht0, aeiv0 (namelist)
+      !!        'key_traldf_c1d': depth dependent coef. defined in 
+      !!                            in ldf_tra_c1d routine
+      !!        'key_traldf_c2d': latitude and longitude dependent coef.
+      !!                            defined in ldf_tra_c2d routine
+      !!        'key_traldf_c3d': latitude, longitude, depth dependent coef.
+      !!                            defined in ldf_tra_c3d routine
+      !!
+      !!      N.B. User defined include files.  By default, 3d and 2d coef.
+      !!      are set to a constant value given in the namelist and the 1d
+      !!      coefficients are initialized to a hyperbolic tangent vertical
+      !!      profile.
+      !!----------------------------------------------------------------------
+      INTEGER ::   ioptio               ! temporary integer
+      INTEGER ::   ios                  ! temporary integer
+      LOGICAL ::   ll_print = .FALSE.   ! =T print eddy coef. in numout
+      !!----------------------------------------------------------------------
+
+!#if defined key_traldf_c3d
+!      CALL ldf_tra_c3d( ll_print )      ! aht = 3D coef. = F( longitude, latitude, depth )
+!#elif defined key_traldf_c2d
+      CALL ldf_tra_c2d_crs( ll_print )      ! aht = 2D coef. = F( longitude, latitude )
+!#elif defined key_traldf_c1d
+!      CALL ldf_tra_c1d( ll_print )      ! aht = 1D coef. = F( depth )
+!#endif
+                                        ! Constant coefficients
+      !
+   END SUBROUTINE ldf_tra_crs_init
+
+!#if defined key_traldf_c3d
+!#   include "ldftra_c3d.h90"
+!#elif defined key_traldf_c2d
+#   include "ldftra_c2d_crs.h90"
+!#elif defined key_traldf_c1d
+!#   include "ldftra_c1d.h90"
+!#endif
+
+   !!======================================================================
+END MODULE ldftra_crs
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/eosbn2.F90	(revision 6101)
@@ -80,6 +80,6 @@
    !                                   !!!  simplified eos coefficients
    ! default value: Vallis 2006
-   REAL(wp) ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff. 
-   REAL(wp) ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff. 
+   REAL(wp), PUBLIC ::   rn_a0      = 1.6550e-1_wp     ! thermal expansion coeff. 
+   REAL(wp), PUBLIC ::   rn_b0      = 7.6554e-1_wp     ! saline  expansion coeff. 
    REAL(wp) ::   rn_lambda1 = 5.9520e-2_wp     ! cabbeling coeff. in T^2        
    REAL(wp) ::   rn_lambda2 = 5.4914e-4_wp     ! cabbeling coeff. in S^2        
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd_crs.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd_crs.F90	(revision 6101)
@@ -178,7 +178,7 @@
          END IF
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) )
-           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) )
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) )
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) )
          ENDIF
 
@@ -240,7 +240,7 @@
          END IF
          !                                 ! "Poleward" heat and salt transports (contribution of upstream fluxes)
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN  
-           IF( jn == jp_tem )  htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
-           IF( jn == jp_sal )  str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN  
+           IF( jn == jp_tem )  htr_adv(:) = ptr_sj( zwy(:,:,:) ) + htr_adv(:)
+           IF( jn == jp_sal )  str_adv(:) = ptr_sj( zwy(:,:,:) ) + str_adv(:)
          ENDIF
          !
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_crs.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_iso_crs.F90	(revision 6101)
@@ -219,7 +219,7 @@
          !
          ! "Poleward" diffusive heat or salt transports (T-S case only)
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
-            IF( jn == jp_tem)   htr_ldf(:) = ptr_vj( zftv(:,:,:) )
-            IF( jn == jp_sal)   str_ldf(:) = ptr_vj( zftv(:,:,:) )
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
+            IF( jn == jp_tem)   htr_ldf(:) = ptr_sj( zftv(:,:,:) )
+            IF( jn == jp_sal)   str_ldf(:) = ptr_sj( zftv(:,:,:) )
          ENDIF
  
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_crs.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/TRA/traldf_lap_crs.F90	(revision 6101)
@@ -150,7 +150,7 @@
          !
          ! "Poleward" diffusive heat or salt transports
-         IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
-            IF( jn  == jp_tem)   htr_ldf(:) = ptr_vj( ztv(:,:,:) )
-            IF( jn  == jp_sal)   str_ldf(:) = ptr_vj( ztv(:,:,:) )
+         IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
+            IF( jn  == jp_tem)   htr_ldf(:) = ptr_sj( ztv(:,:,:) )
+            IF( jn  == jp_sal)   str_ldf(:) = ptr_sj( ztv(:,:,:) )
          ENDIF
          !                                                  ! ==================
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd_crs.F90	(revision 6101)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd_crs.F90	(revision 6101)
@@ -0,0 +1,126 @@
+MODULE zdfevd_crs
+   !!======================================================================
+   !!                       ***  MODULE  zdfevd  ***
+   !! Ocean physics: parameterization of convection through an enhancement
+   !!                of vertical eddy mixing coefficient
+   !!======================================================================
+   !! History :  OPA  !  1997-06  (G. Madec, A. Lazar)  Original code
+   !!   NEMO     1.0  !  2002-06  (G. Madec)  F90: Free form and module
+   !!             -   !  2005-06  (C. Ethe) KPP parameterization
+   !!            3.2  !  2009-03  (M. Leclair, G. Madec, R. Benshila) test on both before & after
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   zdf_evd      : increase the momentum and tracer Kz at the location of
+   !!                  statically unstable portion of the water column (ln_zdfevd=T)
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain variables
+   USE zdf_oce         ! ocean vertical physics variables
+   USE zdfkpp          ! KPP vertical mixing
+   USE in_out_manager  ! I/O manager
+   USE iom             ! for iom_put
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE timing          ! Timing
+   USE wrk_nemo       ! work arrays
+   USE crs
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_evd_crs    ! called by step.F90
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: zdfevd.F90 4990 2014-12-15 16:42:49Z timgraham $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE zdf_evd_crs( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_evd  ***
+      !!                   
+      !! ** Purpose :   Local increased the vertical eddy viscosity and diffu-
+      !!      sivity coefficients when a static instability is encountered.
+      !!
+      !! ** Method  :   avt, avm, and the 4 neighbouring avmu, avmv coefficients
+      !!      are set to avevd (namelist parameter) if the water column is 
+      !!      statically unstable (i.e. if rn2 < -1.e-12 )
+      !!
+      !! ** Action  :   avt, avm, avmu, avmv updted in static instability cases
+      !!
+      !! References :   Lazar, A., these de l'universite Paris VI, France, 1997
+      !!----------------------------------------------------------------------
+      !USE oce,   zavt_evd => ua , zavm_evd => va  ! (ua,va) used ua workspace
+      !
+      INTEGER, INTENT( in ) ::   kt   ! ocean time-step indexocean time step
+      !
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zavt_evd_crs,zavm_evd_crs
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zdf_evd_crs')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'zdf_evd_crs : Enhanced Vertical Diffusion (evd)'
+         IF(lwp) WRITE(numout,*) '~~~~~~~ '
+         IF(lwp) WRITE(numout,*)
+      ENDIF
+
+      CALL wrk_alloc( jpi_crs,jpj_crs,jpk, zavt_evd_crs, zavm_evd_crs )
+
+      zavt_evd_crs(:,:,:) = avt_crs(:,:,:)           ! set avt prior to evd application
+
+      SELECT CASE ( nn_evdm )
+      !
+      !CASE ( 1 )           ! enhance vertical eddy viscosity and diffusivity (if rn2<-1.e-12)
+      !   !
+      !   zavm_evd(:,:,:) = avm_crs(:,:,:)           ! set avm prior to evd application
+      !   !
+      !   DO jk = 1, jpkm1 
+      !      DO jj = 2, jpj_crs             ! no vector opt.
+      !         DO ji = 2, jpi_crs
+      !            IF(  MIN( rn2_crs(ji,jj,jk), rb2_crs(ji,jj,jk) ) <= -1.e-12 ) THEN
+      !               avt_crs (ji  ,jj  ,jk) = rn_avevd * tmask_crs(ji  ,jj  ,jk)
+      !               avm_crs (ji  ,jj  ,jk) = rn_avevd * tmask_crs(ji  ,jj  ,jk)
+      !               avmu_crs(ji  ,jj  ,jk) = rn_avevd * umask_crs(ji  ,jj  ,jk)
+      !               avmu_crs(ji-1,jj  ,jk) = rn_avevd * umask_crs(ji-1,jj  ,jk)
+      !               avmv_crs(ji  ,jj  ,jk) = rn_avevd * vmask_crs(ji  ,jj  ,jk)
+      !               avmv_crs(ji  ,jj-1,jk) = rn_avevd * vmask_crs(ji  ,jj-1,jk)
+      !            ENDIF
+      !         END DO
+      !      END DO
+      !   END DO 
+      !   CALL crs_lbc_lnk( avt_crs , 'W', 1. )   ;   CALL crs_lbc_lnk( avm_crs , 'W', 1. )   ! Lateral boundary conditions
+      !   CALL crs_lbc_lnk( avmu_crs, 'U', 1. )   ;   CALL crs_lbc_lnk( avmv_crs, 'V', 1. )
+      !   !
+      !   zavm_evd_crs(:,:,:) = avm_crs(:,:,:) - zavm_evd_crs(:,:,:)   ! change in avm due to evd
+      !   CALL iom_put( "avm_evd", zavm_evd_crs )              ! output this change
+      !   !
+      CASE DEFAULT         ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 
+         DO jk = 1, jpkm1
+            DO jj = 1, jpj_crs             ! loop over the whole domain (no lbc_lnk call)
+               DO ji = 1, jpi_crs
+                  IF(  MIN( rn2_crs(ji,jj,jk), rb2_crs(ji,jj,jk) ) <= -1.e-12 )   &
+                     avt_crs(ji,jj,jk) = rn_avevd * tmask_crs(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+      END SELECT 
+
+      zavt_evd_crs(:,:,:) = avt_crs(:,:,:) - zavt_evd_crs(:,:,:)   ! change in avt due to evd
+      CALL iom_put( "avt_evd_crs", zavt_evd_crs )              ! output this change
+      !
+      CALL wrk_dealloc( jpi_crs,jpj_crs,jpk, zavt_evd_crs, zavm_evd_crs )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zdf_evd_crs')
+      !
+   END SUBROUTINE zdf_evd_crs
+
+   !!======================================================================
+END MODULE zdfevd_crs
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90	(revision 6101)
@@ -134,4 +134,6 @@
          END DO
       END DO
+      CALL iom_put("hmlpt",hmlpt)
+
       IF( .NOT.lk_offline ) THEN            ! no need to output in offline mode
          CALL iom_put( "mldr10_1", hmlp )   ! mixed layer depth
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl_crs.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl_crs.F90	(revision 6101)
@@ -64,4 +64,5 @@
       REAL(wp) ::   zN2_c        ! local scalar
       INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
+      REAL(wp), POINTER, DIMENSION(:,:) ::  z2d   ! 2D workspace
       !!----------------------------------------------------------------------
       !
@@ -69,4 +70,5 @@
       !
       CALL wrk_alloc( jpi_crs,jpj_crs, imld )
+      CALL wrk_alloc( jpi_crs,jpj_crs, z2d )
 
       IF( kt == nit000 ) THEN
@@ -98,5 +100,10 @@
       END DO
       !
+      z2d=REAL(nmln_crs,wp)
+      CALL iom_put("nmln_crs",z2d)
+      CALL iom_put("hmlpt_crs",hmlpt_crs)
+      !
       CALL wrk_dealloc( jpi_crs,jpj_crs, imld )
+      CALL wrk_dealloc( jpi_crs,jpj_crs, z2d )
       !
       IF( nn_timing == 1 )  CALL timing_stop('zdf_mxl_crs')
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 6101)
@@ -64,14 +64,14 @@
 
    !                      !!** Namelist  namzdf_tke  **
-   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
-   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
-   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
-   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
-   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
+   LOGICAL  , PUBLIC  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
+   INTEGER  , PUBLIC  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
+   REAL(wp) , PUBLIC  ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
+   INTEGER  , PUBLIC  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
+   REAL(wp) , PUBLIC  ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
    REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
    REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
    REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
    REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
-   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
+   REAL(wp) , PUBLIC  ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
    INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
    INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
@@ -81,5 +81,5 @@
 
    REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
-   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
+   REAL(wp) , PUBLIC  ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
    REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
    REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke_crs.F90	(revision 6101)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke_crs.F90	(revision 6101)
@@ -0,0 +1,380 @@
+MODULE zdftke_crs
+   !!======================================================================
+   !!                       ***  MODULE  zdftke  ***
+   !! Ocean physics:  vertical mixing coefficient computed from the tke 
+   !!                 turbulent closure parameterization
+   !!=====================================================================
+   !! History :  OPA  !  1991-03  (b. blanke)  Original code
+   !!            7.0  !  1991-11  (G. Madec)   bug fix
+   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
+   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
+   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
+   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
+   !!            8.0  !  1997-07  (G. Madec)   lbc
+   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
+   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
+   !!             -   !  2004-10  (C. Ethe )  1D configuration
+   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
+   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
+   !!                 !           - tke penetration (wind steering)
+   !!                 !           - suface condition for tke & mixing length
+   !!                 !           - Langmuir cells
+   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
+   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
+   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term 
+   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl 
+   !!                 !                                + cleaning of the parameters + bugs correction
+   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!----------------------------------------------------------------------
+#if defined key_zdftke   ||   defined key_esopa
+   !!----------------------------------------------------------------------
+   !!   'key_zdftke'                                   TKE vertical physics
+   !!----------------------------------------------------------------------
+   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
+   !!   tke_tke       : tke time stepping: update tke at now time step (en)
+   !!   tke_avn       : compute mixing length scale and deduce avm and avt
+   !!   zdf_tke_init  : initialization, namelist read, and parameters control
+   !!   tke_rst       : read/write tke restart in ocean restart file
+   !!----------------------------------------------------------------------
+   USE crs
+   USE zdf_oce , ONLY : avtb, avmb, avtb_2d
+   USE zdftke
+   USE crslbclnk
+   !USE oce            ! ocean: dynamics and active tracers variables
+   USE phycst         ! physical constants
+   !USE dom_oce        ! domain: ocean
+   !USE domvvl         ! domain: variable volume layer
+   !USE sbc_oce        ! surface boundary condition: ocean
+   !USE zdf_oce        ! vertical physics: ocean variables
+   !USE zdfmxl         ! vertical physics: mixed layer
+   !USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   !USE prtctl         ! Print control
+   !USE in_out_manager ! I/O manager
+   !USE iom            ! I/O manager library
+   !USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   !USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC tke_avn_crs
+   PUBLIC tke_avn_ini_crs
+
+   !LOGICAL , PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
+
+   !                      !!** Namelist  namzdf_tke  **
+   !LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
+   !INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
+   !REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
+   !INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
+   !REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
+   !REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation 
+   !REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
+   !REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
+   !REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
+   !REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
+   !INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
+   !INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
+   !REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean
+   !LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
+   !REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
+
+   !REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
+   !REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
+   !REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
+   !REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
+
+   !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
+   !REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
+   !REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
+   !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avt_k , avm_k  ! not enhanced Kz
+   !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   avmu_k, avmv_k ! not enhanced Kz
+#if defined key_c1d
+   !                                                                        !!** 1D cfg only  **   ('key_c1d')
+   !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales
+   !REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
+#endif
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id: zdftke.F90 4990 2014-12-15 16:42:49Z timgraham $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+!   INTEGER FUNCTION zdf_tke_alloc()
+!      !!----------------------------------------------------------------------
+!      !!                ***  FUNCTION zdf_tke_alloc  ***
+!      !!----------------------------------------------------------------------
+!      ALLOCATE(                                                                    &
+!#if defined key_c1d
+!         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
+!         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
+!#endif
+!         &      en    (jpi,jpj,jpk) , htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     & 
+!         &      avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk),                          &
+!         &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      )
+!         !
+!      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
+!      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
+!      !
+!   END FUNCTION zdf_tke_alloc
+
+
+   SUBROUTINE tke_avn_crs
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tke_avn  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!
+      !! ** Method  :   At this stage, en, the now TKE, is known (computed in 
+      !!              the tke_tke routine). First, the now mixing lenth is 
+      !!      computed from en and the strafification (N^2), then the mixings
+      !!      coefficients are computed.
+      !!              - Mixing length : a first evaluation of the mixing lengh
+      !!      scales is:
+      !!                      mxl = sqrt(2*en) / N  
+      !!      where N is the brunt-vaisala frequency, with a minimum value set
+      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
+      !!        The mixing and dissipative length scale are bound as follow : 
+      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
+      !!                        zmxld = zmxlm = mxl
+      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
+      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 
+      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
+      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
+      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to 
+      !!                    the surface to obtain ldown. the resulting length 
+      !!                    scales are:
+      !!                        zmxld = sqrt( lup * ldown ) 
+      !!                        zmxlm = min ( lup , ldown )
+      !!              - Vertical eddy viscosity and diffusivity:
+      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
+      !!                      avt = max( avmb, pdlr * avm )  
+      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
+      !!
+      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
+      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
+      REAL(wp) ::   zdku, zdkv, zpdlr, zri, zsqen     !   -      -
+      REAL(wp) ::   zemlm, zemlp !   -      -
+      !REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxlm_crs,zmxld_crs
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: avm_crs, avmu_crs,avmv_crs
+      !!--------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tke_avn')
+
+      !CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      CALL wrk_alloc( jpi_crs,jpj_crs,jpk, zmxlm_crs, zmxld_crs, avm_crs, avmu_crs, avmv_crs  ) 
+
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Mixing length
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
+      !
+      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
+         DO jj = 2, jpj_crs-1
+            DO ji = 2, jpi_crs-1
+               !IF (mikt(ji,jj) .GT. 1) THEN
+               !   zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min
+               !ELSE
+                  zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+                  zmxlm_crs(ji,jj,1) = MAX( rn_mxl0, zraug * taum_crs(ji,jj) )
+               !END IF
+            END DO
+         END DO
+      !ELSE 
+      !   DO jj = 2, jpjm1
+      !      DO ji = fs_2, fs_jpim1                         ! surface set to the minimum value
+      !         zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min)
+      !      END DO
+      !   END DO
+      ENDIF
+      zmxlm_crs(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value
+      !
+      DO jj = 2, jpj_crs-1
+         DO ji = 2, jpi_crs-1   ! vector opt.
+            DO jk = 2, jpk-1              ! interior value : l=sqrt(2*e/n^2)
+               zrn2 = MAX( rn2_crs(ji,jj,jk), rsmall )
+               zmxlm_crs(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en_crs(ji,jj,jk) / zrn2 )  )
+            END DO
+            zmxld_crs(ji,jj,1) = zmxlm_crs(ji,jj,1)   ! surface set to the minimum value 
+         END DO
+      END DO
+      !
+      !                     !* Physical limits for the mixing length
+      !
+      zmxld_crs(:,:, 1 ) = zmxlm_crs(:,:,1)   ! surface set to the zmxlm   value
+      zmxld_crs(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
+      !
+      SELECT CASE ( nn_mxl )
+      !
+      !CASE ( 0 )           ! bounded by the distance to surface and bottom
+      !   DO jj = 2, jpjm1
+      !      DO ji = fs_2, fs_jpim1   ! vector opt.
+      !         DO jk = mikt(ji,jj)+1, jpkm1
+      !            zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
+      !            &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
+      !            zmxlm(ji,jj,jk) = zemxl
+      !            zmxld(ji,jj,jk) = zemxl
+      !         END DO
+      !      END DO
+      !   END DO
+      !   !
+      !CASE ( 1 )           ! bounded by the vertical scale factor
+      !   DO jj = 2, jpjm1
+      !      DO ji = fs_2, fs_jpim1   ! vector opt.
+      !         DO jk = mikt(ji,jj)+1, jpkm1
+      !            zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
+      !            zmxlm(ji,jj,jk) = zemxl
+      !            zmxld(ji,jj,jk) = zemxl
+      !         END DO
+      !      END DO
+      !   END DO
+      !   !
+      !CASE ( 2 )           ! |dk[xml]| bounded by e3t :
+      !   DO jj = 2, jpjm1
+      !      DO ji = fs_2, fs_jpim1   ! vector opt.
+      !         DO jk = mikt(ji,jj)+1, jpkm1         ! from the surface to the bottom :
+      !            zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+      !         END DO
+      !         DO jk = jpkm1, mikt(ji,jj)+1, -1     ! from the bottom to the surface :
+      !            zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+      !            zmxlm(ji,jj,jk) = zemxl
+      !            zmxld(ji,jj,jk) = zemxl
+      !         END DO
+      !      END DO
+      !   END DO
+      !   !
+      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
+         DO jj = 2, jpj_crs-1
+            DO ji = 2, jpi_crs-1   ! vector opt.
+               DO jk = 2, jpkm1         ! from the surface to the bottom : lup
+                  zmxld_crs(ji,jj,jk) = MIN( zmxld_crs(ji,jj,jk-1) + e3t_crs(ji,jj,jk-1), zmxlm_crs(ji,jj,jk) )
+               END DO
+               DO jk = jpkm1, 2 , -1     ! from the bottom to the surface : ldown
+                  zmxlm_crs(ji,jj,jk) = MIN( zmxlm_crs(ji,jj,jk+1) + e3t_crs(ji,jj,jk+1), zmxlm_crs(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         DO jk = 2, jpkm1
+            DO jj = 2, jpj_crs-1
+               DO ji = 2, jpi_crs-1   ! vector opt.
+                  zemlm = MIN ( zmxld_crs(ji,jj,jk),  zmxlm_crs(ji,jj,jk) )
+                  zemlp = SQRT( zmxld_crs(ji,jj,jk) * zmxlm_crs(ji,jj,jk) )
+                  zmxlm_crs(ji,jj,jk) = zemlm
+                  zmxld_crs(ji,jj,jk) = zemlp
+               END DO
+            END DO
+         END DO
+         !
+      END SELECT
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
+         DO jj = 2, jpj_crs-1
+            DO ji = 2, jpi_crs-1   ! vector opt.
+               zsqen = SQRT( en_crs(ji,jj,jk) )
+               zav   = rn_ediff * zmxlm_crs(ji,jj,jk) * zsqen
+               avm_crs(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask_crs(ji,jj,jk)
+               avt_crs(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask_crs(ji,jj,jk)
+               !dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL crs_lbc_lnk( avm_crs, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      DO jj = 2, jpj_crs-1
+         DO ji = 2, jpi_crs-1   ! vector opt.
+            DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
+               avmu_crs(ji,jj,jk) = 0.5 * ( avm_crs(ji,jj,jk) + avm_crs(ji+1,jj  ,jk) ) * umask_crs(ji,jj,jk)
+            END DO
+            DO jk = 2, jpkm1
+               avmv_crs(ji,jj,jk) = 0.5 * ( avm_crs(ji,jj,jk) + avm_crs(ji  ,jj+1,jk) ) * vmask_crs(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL crs_lbc_lnk( avmu_crs, 'U', 1. )   ;   CALL crs_lbc_lnk( avmv_crs, 'V', 1. )      ! Lateral boundary conditions
+      !
+      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
+         DO jj = 2, jpj_crs-1
+            DO ji = 2, jpi_crs-1   ! vector opt.
+               DO jk = 2, jpkm1
+                  zcoef = avm_crs(ji,jj,jk) * 2._wp * e3w_crs(ji,jj,jk) * e3w_crs(ji,jj,jk)
+                  !                                          ! shear
+                  zdku = avmu_crs(ji-1,jj,jk) * ( un_crs(ji-1,jj,jk-1) - un_crs(ji-1,jj,jk) ) * ( ub_crs(ji-1,jj,jk-1) - ub_crs(ji-1,jj,jk) )   &
+                    &  + avmu_crs(ji  ,jj,jk) * ( un_crs(ji  ,jj,jk-1) - un_crs(ji  ,jj,jk) ) * ( ub_crs(ji  ,jj,jk-1) - ub_crs(ji  ,jj,jk) )
+                  zdkv = avmv_crs(ji,jj-1,jk) * ( vn_crs(ji,jj-1,jk-1) - vn_crs(ji,jj-1,jk) ) * ( vb_crs(ji,jj-1,jk-1) - vb_crs(ji,jj-1,jk) )   &
+                    &  + avmv_crs(ji,jj  ,jk) * ( vn_crs(ji,jj  ,jk-1) - vn_crs(ji,jj  ,jk) ) * ( vb_crs(ji,jj  ,jk-1) - vb_crs(ji,jj  ,jk) )
+                  !                                          ! local Richardson number
+                  zri   = MAX( rb2_crs(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
+                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
+                  avt_crs(ji,jj,jk)   = MAX( zpdlr * avt_crs(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask_crs(ji,jj,jk)
+              END DO
+            END DO
+         END DO
+      ENDIF
+      CALL crs_lbc_lnk( avt_crs, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
+
+      !IF(ln_ctl) THEN
+      !   CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
+      !   CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
+      !      &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
+      !ENDIF
+      !
+      CALL wrk_dealloc( jpi_crs,jpj_crs,jpk, zmxlm_crs, zmxld_crs, avm_crs, avmu_crs, avmv_crs  ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
+      !
+   END SUBROUTINE tke_avn_crs
+
+   SUBROUTINE tke_avn_ini_crs
+      !!----------------------------------------------------------------------
+      !!
+      !!
+      !!
+      !!
+      !!----------------------------------------------------------------------
+      !IF( nn_avb == 0 ) THEN                ! Define avmb, avtb from namelist parameter
+      !   avtb(:) = rn_avt0
+      !ELSE                                  ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
+      !   avmb(:) = rn_avm0
+      !   avtb(:) = rn_avt0 + ( 3.e-4_wp - 2._wp * rn_avt0 ) * 1.e-4_wp * gdepw_1d(:)   ! m2/s
+      !   IF(ln_sco .AND. lwp)   CALL ctl_warn( 'avtb profile not valid in sco' )
+      !ENDIF
+
+
+
+      avtb_2d_crs(:,:) = 1.e0
+      !IF( nn_havtb == 1 ) THEN                   ! decrease avtb in the equatorial band
+      !     !  -15S -5S : linear decrease from avt0 to avt0/10.
+      !     !  -5S  +5N : cst value avt0/10.
+      !     !   5N  15N : linear increase from avt0/10, to avt0
+      !     WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
+      !     WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
+      !     WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
+      !ENDIF
+
+
+
+   END SUBROUTINE tke_avn_ini_crs
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                                        NO TKE scheme
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
+#endif
+
+   !!======================================================================
+END MODULE zdftke_crs
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 6101)
@@ -33,5 +33,5 @@
    USE step_oce         ! time stepping definition modules
    USE iom
-   USE crs
+   use wrk_nemo
 
    IMPLICIT NONE
@@ -76,4 +76,6 @@
       INTEGER ::   indic    ! error indicator if < 0
       INTEGER ::   kcall    ! optional integer argument (dom_vvl_sf_nxt)
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zrhop,zrhd
+
       !! ---------------------------------------------------------------------
 
@@ -117,6 +119,6 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !  THERMODYNAMICS
-                         !cbr not used CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
-                         !cbr not used CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
+                         CALL eos_rab( tsb, rab_b )       ! before local thermal/haline expension ratio at T-points
+                         CALL eos_rab( tsn, rab_n )       ! now    local thermal/haline expension ratio at T-points
                          CALL bn2    ( tsb, rab_b, rn2b ) ! before Brunt-Vaisala frequency
                          CALL bn2    ( tsn, rab_n, rn2  ) ! now    Brunt-Vaisala frequency
@@ -169,4 +171,5 @@
                          CALL iom_put("rhd",rhd)
                          CALL iom_put("rn2b",rn2b)
+                         CALL iom_put("rn2",rn2)
                          CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
          ENDIF
@@ -240,4 +243,11 @@
       IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
 
+      CALL wrk_alloc( jpi, jpj, jpk, zrhop, zrhd )
+      CALL eos( tsn, zrhd, zrhop, fsdept_n(:,:,:) )                 ! now in situ and potential density
+      zrhop(:,:,jpk) = 0._wp
+      CALL iom_put( 'rhop', zrhop )
+
+      CALL wrk_dealloc( jpi, jpj, jpk, zrhop, zrhd )
+
 #if defined key_top
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -247,7 +257,9 @@
                          CALL dom_grid_crs
 
-                         CALL eos_rab_crs( tsn_crs, rab_crs_n )       ! now    local thermal/haline expension ratio at T-points
-                         CALL bn2_crs    ( tsn_crs, rab_crs_n, rb2_crs  ) ! now    Brunt-Vaisala frequency
+                         !CALL eos_rab_crs( tsn_crs, rab_crs_n )       ! now    local thermal/haline expension ratio at T-points
+                         !CALL bn2_crs    ( tsn_crs, rab_crs_n, rb2_crs  ) ! now    Brunt-Vaisala frequency
                          CALL eos_crs ( tsn_crs, rhd_crs, rhop_crs, gdept_crs(:,:,:) ) ! now in situ density for hpg computation
+                         CALL iom_put("rhop_crs",rhop_crs)
+                         CALL iom_put("rhd_crs",rhd_crs)
 
          IF( ln_zps )    CALL zps_hde_crs( kstp, jpts, tsb_crs, gtsu_crs, gtsv_crs, rhd_crs, gru_crs, grv_crs )
@@ -265,4 +277,5 @@
 
       IF( ln_crs_top )   CALL dom_grid_glo
+
 #endif
 
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/stpctl.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/stpctl.F90	(revision 6101)
@@ -136,5 +136,5 @@
             WRITE(numout,*) '          output of last fields in numwso'
          ENDIF
-         kindic = -3
+         WHERE( tsn(:,:,:,jp_sal) .LE. 0. )  tsn(:,:,:,jp_sal) = 0.1
       ENDIF
 9500  FORMAT (' kt=',i6,' min SSS: ',1pg11.4,', i j: ',2i5)
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcdmp.F90	(revision 6101)
@@ -24,5 +24,5 @@
    USE trdtra
    USE trd_oce
-   USE iom
+   USE iom , ONLY : iom_open, iom_get, iom_close, jpdom_autoglo
 
    IMPLICIT NONE
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 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcldf_crs.F90	(revision 6101)
@@ -219,9 +219,4 @@
          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
 
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 6101)
@@ -19,5 +19,5 @@
    USE trc             ! ocean  passive tracers variables
    USE prtctl_trc      ! Print control for debbuging
-   USE iom
+   USE iom, ONLY : iom_varid, iom_get, iom_rstput,jpdom_autoglo
    USE trd_oce
    USE trdtra
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 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90	(revision 6101)
@@ -34,8 +34,7 @@
    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 dom_oce , ONLY : ln_crs, ln_isfcav
    USE crs     , ONLY : jpi_crs,jpj_crs,wn_crs,ln_crs_top !cbr
    USE ldfslp_crs
-
 #if defined key_agrif
    USE agrif_top_sponge ! tracers sponges
@@ -111,9 +110,10 @@
                CALL zps_hde_crs( kstp, jptra, trn, gtru, gtrv )
             ELSE
-               IF( ln_isfcav)        &
+               IF( ln_isfcav )THEN
                   CALL zps_hde_isf( kstp, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! Partial steps: now horizontal gradient of passive
                ELSE
                   CALL zps_hde    ( kstp, jptra, trn, gtru, gtrv )   ! Partial steps: now horizontal gradient of passive
                ENDIF
+            ENDIF
          ENDIF
                                                                 ! tracers at the bottom ocean level
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 6101)
@@ -76,4 +76,5 @@
    USE dom_oce , ONLY : ln_sco => ln_sco
    USE dom_oce , ONLY : neuler => neuler
+   USE dom_oce , ONLY : mikt      => mikt          !: f-points (m)
 
    USE crs,  ONLY : mi0 => mi0 
@@ -282,4 +283,5 @@
    USE dom_oce , ONLY :   mbku      => mbku          !: f-points (m)
    USE dom_oce , ONLY :   mbkv      => mbkv          !: f-points (m)
+   USE dom_oce , ONLY :   mikt      => mikt          !: f-points (m)
 
    !* IO manager *
@@ -401,4 +403,10 @@
 #endif
 
+
+   USE dom_oce , ONLY : ndastp
+   USE sbc_oce , ONLY : nn_ice_embd
+   USE sbc_oce , ONLY : ln_cpl
+   USE sbc_oce , ONLY : ncpl_qsr_freq
+
 #else
    !!----------------------------------------------------------------------
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcdta.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcdta.F90	(revision 6101)
@@ -19,4 +19,5 @@
    USE oce_trc       !  shared variables between ocean and passive tracers
    USE trc           !  passive tracers common variables
+!   USE iom           !  I/O manager
    USE lib_mpp       !  MPP library
    USE fldread       !  read input fields
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcini.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcini.F90	(revision 6101)
@@ -25,13 +25,13 @@
    USE trcini_my_trc   ! MY_TRC   initialisation
    USE trcdta          ! initialisation from files
-   USE zpshde,ONLY: zps_hde    ! partial step: hor. derivative   (zps_hde routine)
+   USE zpshde,ONLY: zps_hde, zps_hde_isf    ! partial step: hor. derivative   (zps_hde routine)
    USE zpshde_crs      ! partial step: hor. derivative   (zps_hde routine)
    USE prtctl_trc      ! Print control passive tracers (prt_ctl_trc_init routine)
    USE trcsub          ! variables to substep passive tracers
    USE lib_mpp         ! distribued memory computing library
-   USE sbc_oce, ONLY : ltrcdm2dc
    USE crs    , ONLY : ln_crs
-   USE dom_oce, ONLY : nn_cla
+   USE dom_oce, ONLY : nn_cla, ln_isfcav
    USE trcice          ! tracers in sea ice
+   USE sbc_oce
  
    IMPLICIT NONE
@@ -148,13 +148,11 @@
  
       tra(:,:,:,:) = 0._wp
-      IF( ln_zps .AND. .NOT. lk_c1d )THEN              ! Partial steps: before horizontal gradient of passive
-         IF( ln_crs )  THEN
-            CALL zps_hde_crs( nit000, jptra, trn, gtru, gtrv )
-         ELSE
-            IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav )   &              ! Partial steps: before horizontal gradient of passive
-            &    CALL zps_hde    ( nit000, jptra, trn, gtru, gtrv  )  ! Partial steps: before horizontal gradient
-            IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav )   &
-            &    CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )       ! tracers at the bottom ocean level
-         ENDIF
+      IF( ln_crs )  THEN
+         CALL zps_hde_crs( nit000, jptra, trn, gtru, gtrv )
+      ELSE
+         IF( ln_zps .AND. .NOT. lk_c1d .AND. .NOT. ln_isfcav )   &              ! Partial steps: before horizontal gradient of passive
+         &    CALL zps_hde    ( nit000, jptra, trn, gtru, gtrv  )  ! Partial steps: before horizontal gradient
+         IF( ln_zps .AND. .NOT. lk_c1d .AND.       ln_isfcav )   &
+         &    CALL zps_hde_isf( nit000, jptra, trn, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )       ! tracers at the bottom ocean level
       ENDIF
 
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcstp.F90	(revision 6101)
@@ -11,5 +11,4 @@
    !!----------------------------------------------------------------------
    USE oce_trc          ! ocean dynamics and active tracers variables
-   USE sbc_oce , ONLY : ltrcdm2dc,qsr_mean
    USE trc
    USE trctrp           ! passive tracers transport
@@ -24,5 +23,6 @@
    USE in_out_manager
    USE trcsub
-   USE dom_oce, ONLY : nday, nmonth, nyear
+   USE dom_oce, ONLY : nday, nmonth, nyear, nsec1jan000, nsec_year
+   !USE sbc_oce
 
    IMPLICIT NONE
Index: /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcwri.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 6100)
+++ /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/TOP_SRC/trcwri.F90	(revision 6101)
@@ -15,5 +15,5 @@
    USE oce_trc     ! shared variables between ocean and passive tracers
    USE trc         ! passive tracers common variables 
-   USE iom         ! I/O manager
+!   USE iom         ! I/O manager
    USE dianam      ! Output file name
    USE trcwri_pisces
