Index: /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbccpl.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbccpl.F90	(revision 2188)
+++ /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/sbccpl.F90	(revision 2189)
@@ -169,5 +169,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id: sbccpl.F90 3914 2013-06-12 20:17:48Z smasson $
+   !! $Id: sbccpl.F90 4162 2013-11-07 10:19:49Z cetlod $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -221,4 +221,5 @@
       !!
       INTEGER ::   jn   ! dummy loop index
+      INTEGER ::   ios  ! Local integer output status for namelist read
       REAL(wp), POINTER, DIMENSION(:,:) ::   zacs, zaos
       !!
@@ -236,28 +237,12 @@
       ! ================================ !
 
-      ! default definitions
-      !                    !     description       !  multiple  !    vector   !      vector          ! vector !
-      !                    !                       ! categories !  reference  !    orientation       ! grids  !
-      ! send
-      sn_snd_temp   = FLD_C( 'weighted oce and ice',    'no'    ,     ''      ,         ''           ,   ''   ) 
-      sn_snd_alb    = FLD_C( 'weighted ice'        ,    'no'    ,     ''      ,         ''           ,   ''   ) 
-      sn_snd_thick  = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''           ,   ''   ) 
-      sn_snd_crt    = FLD_C( 'none'                ,    'no'    , 'spherical' , 'eastward-northward' ,  'T'   )     
-      sn_snd_co2    = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''           ,   ''   )     
-      ! receive
-      sn_rcv_w10m   = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_taumod = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_tau    = FLD_C( 'oce only'            ,    'no'    , 'cartesian' , 'eastward-northward',  'U,V'  )  
-      sn_rcv_dqnsdt = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_qsr    = FLD_C( 'oce and ice'         ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_qns    = FLD_C( 'oce and ice'         ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_emp    = FLD_C( 'conservative'        ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_rnf    = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_cal    = FLD_C( 'coupled'             ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_iceflx = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
-      sn_rcv_co2    = FLD_C( 'none'                ,    'no'    ,     ''      ,         ''          ,   ''    )
-
-      REWIND( numnam )                    ! ... read namlist namsbc_cpl
-      READ  ( numnam, namsbc_cpl )
+      REWIND( numnam_ref )              ! Namelist namsbc_cpl in reference namelist : Variables for OASIS coupling
+      READ  ( numnam_ref, namsbc_cpl, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namsbc_cpl in configuration namelist : Variables for OASIS coupling
+      READ  ( numnam_cfg, namsbc_cpl, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_cpl in configuration namelist', lwp )
+      WRITE ( numond, namsbc_cpl )
 
       IF(lwp) THEN                        ! control print
@@ -388,5 +373,5 @@
       !
       IF( TRIM( sn_rcv_tau%cldes ) /= 'oce and ice' ) THEN        ! 'oce and ice' case ocean stress on ocean mesh used
-         srcv(jpr_itz1:jpr_itz2)%laction = .FALSE.    ! ice components not received (itx1 and ity1 used later)
+         srcv(jpr_itx1:jpr_itz2)%laction = .FALSE.    ! ice components not received
          srcv(jpr_itx1)%clgrid = 'U'                  ! ocean stress used after its transformation
          srcv(jpr_ity1)%clgrid = 'V'                  ! i.e. it is always at U- & V-points for i- & j-comp. resp.
@@ -407,5 +392,8 @@
       SELECT CASE( TRIM( sn_rcv_emp%cldes ) )
       CASE( 'oce only'      )   ;   srcv(                                 jpr_oemp   )%laction = .TRUE. 
-      CASE( 'conservative'  )   ;   srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
+      CASE( 'conservative'  )
+         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE.
+!$$ O.M.         IF ( k_ice <= 1 )  srcv(jpr_ivep)%laction = .FALSE.
+         IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE.
       CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE.
       CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' )
@@ -465,13 +453,4 @@
          CALL ctl_stop( 'sbc_cpl_init: namsbc_cpl namelist mismatch between sn_rcv_qns%cldes and sn_rcv_dqnsdt%cldes' )
       !                                                      ! ------------------------- !
-      !                                                      !    Ice Qsr penetration    !   
-      !                                                      ! ------------------------- !
-      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer 
-      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
-      ! Coupled case: since cloud cover is not received from atmosphere 
-      !               ===> defined as constant value -> definition done in sbc_cpl_init
-      IF ( ALLOCATED (fr1_i0)) fr1_i0 (:,:) = 0.18
-      IF ( ALLOCATED (fr2_i0)) fr2_i0 (:,:) = 0.82
-      !                                                      ! ------------------------- !
       !                                                      !      10m wind module      !   
       !                                                      ! ------------------------- !
@@ -508,4 +487,9 @@
       ! Allocate taum part of frcv which is used even when not received as coupling field
       IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jn)%nct) )
+      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE.
+      IF( k_ice /= 0 ) THEN
+         IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jn)%nct) )
+         IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jn)%nct) )
+      END IF
 
       ! ================================ !
@@ -939,8 +923,8 @@
       IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
 
-         !                                                      ! ======================= !
+         !                                                                                              ! ======================= !
 !AC Pour eviter un stress nul sur la glace dans le cas mixes oce-ice
-         IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN                      !   ice stress received   !
-            !                                                   ! ======================= !
+         IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN               !   ice stress received   !
+            !                                                                                           ! ======================= !
             !  
             IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
@@ -1331,4 +1315,13 @@
       END SELECT
 
+      !    Ice Qsr penetration used (only?)in lim2 or lim3 
+      ! fraction of net shortwave radiation which is not absorbed in the thin surface layer 
+      ! and penetrates inside the ice cover ( Maykut and Untersteiner, 1971 ; Elbert anbd Curry, 1993 )
+      ! Coupled case: since cloud cover is not received from atmosphere 
+      !               ===> defined as constant value -> definition done in sbc_cpl_init
+      fr1_i0(:,:) = 0.18
+      fr2_i0(:,:) = 0.82
+
+
       CALL wrk_dealloc( jpi,jpj, zcptn, ztmp, zicefr )
       !
Index: /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/trc_oce.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/trc_oce.F90	(revision 2189)
+++ /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/trc_oce.F90	(revision 2189)
@@ -0,0 +1,285 @@
+MODULE trc_oce
+   !!======================================================================
+   !!                      ***  MODULE  trc_oce  ***
+   !! Ocean passive tracer  :  share SMS/Ocean variables
+   !!======================================================================
+   !! History :  1.0  !  2004-03  (C. Ethe)  Original code
+   !!----------------------------------------------------------------------
+
+   !!----------------------------------------------------------------------
+   !!   trc_oce_rgb : tabulated attenuation coefficients for RGB light penetration         
+   !!----------------------------------------------------------------------
+   USE par_oce
+   USE in_out_manager  ! I/O manager
+   USE dom_oce         ! ocean space and time domain
+   USE lib_mpp         ! MPP library
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   trc_oce_rgb        ! routine called by traqsr.F90
+   PUBLIC   trc_oce_rgb_read   ! routine called by traqsr.F90
+   PUBLIC   trc_oce_ext_lev    ! function called by traqsr.F90 at least
+   PUBLIC   trc_oce_alloc      ! function called by nemogcm.F90
+
+   INTEGER , PUBLIC                                      ::   nn_dttrc      !: frequency of step on passive tracers
+   REAL(wp), PUBLIC                                      ::   r_si2         !: largest depth of extinction (blue & 0.01 mg.m-3)  (RGB)
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   etot3         !: light absortion coefficient
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) ::   facvol        !: volume for degraded regions
+
+#if defined key_top 
+   !!----------------------------------------------------------------------
+   !!   'key_top'                                                 bio-model          
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_qsr_bio = .TRUE.   !: bio-model light absorption flag
+#else
+   !!----------------------------------------------------------------------
+   !! Default option                          No bio-model light absorption      
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_qsr_bio = .FALSE.   !: bio-model light absorption flag
+#endif
+
+#if defined key_offline
+   !!----------------------------------------------------------------------
+   !!   'key_offline'                                     OFFLINE mode          
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_offline = .TRUE.   !: offline flag
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                   NO  OFFLINE mode          
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_offline = .FALSE.   !: offline flag
+#endif
+#if defined key_degrad
+   !!----------------------------------------------------------------------
+   !!   'key_degrad'                                     Degradation mode          
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_degrad = .TRUE.   !: degradation flag
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option                                   NO  Degradation mode          
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_degrad = .FALSE.   !: degradation flag
+#endif
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: trc_oce.F90 3770 2013-02-04 10:51:28Z cetlod $ 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION trc_oce_alloc()
+      !!----------------------------------------------------------------------
+      !!                  ***  trc_oce_alloc  ***
+      !!----------------------------------------------------------------------
+      INTEGER ::   ierr(2)        ! Local variables
+      !!----------------------------------------------------------------------
+      ierr(:) = 0
+                     ALLOCATE( etot3 (jpi,jpj,jpk), STAT=ierr(1) )
+      IF( lk_degrad) ALLOCATE( facvol(jpi,jpj,jpk), STAT=ierr(2) )
+      trc_oce_alloc  = MAXVAL( ierr )
+      !
+      IF( trc_oce_alloc /= 0 )   CALL ctl_warn('trc_oce_alloc: failed to allocate etot3 array')
+   END FUNCTION trc_oce_alloc
+
+
+   SUBROUTINE trc_oce_rgb( prgb )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_opt_init  ***
+      !!
+      !! ** Purpose :   Initialization of of the optical scheme
+      !!
+      !! ** Method  :   Set a look up table for the optical coefficients
+      !!                i.e. the attenuation coefficient for R-G-B light 
+      !!                tabulated in Chlorophyll class (from JM Andre)
+      !!
+      !! ** Action  :   prgb(3,61) tabulated R-G-B attenuation coef. 
+      !!
+      !! Reference  : Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(3,61), INTENT(out) ::   prgb   ! tabulated attenuation coefficient
+      !!
+      INTEGER  ::   jc     ! dummy loop indice
+      INTEGER  ::   irgb   ! temporary integer
+      REAL(wp) ::   zchl   ! temporary scalar
+      REAL(wp), DIMENSION(4,61) ::   zrgb   ! tabulated attenuation coefficient (formerly read in 'kRGB61.txt')
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'trc_oce_rgb : Initialisation of the optical look-up table'
+         WRITE(numout,*) '~~~~~~~~~~~ '
+      ENDIF
+      !
+      !  Chlorophyll        !     Blue attenuation     !     Green attenuation    !     Red attenuation      !
+      zrgb(1, 1) =  0.010   ;   zrgb(2, 1) = 0.01618   ;   zrgb(3, 1) = 0.07464   ;   zrgb(4, 1) = 0.37807
+      zrgb(1, 2) =  0.011   ;   zrgb(2, 2) = 0.01654   ;   zrgb(3, 2) = 0.07480   ;   zrgb(4, 2) = 0.37823
+      zrgb(1, 3) =  0.013   ;   zrgb(2, 3) = 0.01693   ;   zrgb(3, 3) = 0.07499   ;   zrgb(4, 3) = 0.37840
+      zrgb(1, 4) =  0.014   ;   zrgb(2, 4) = 0.01736   ;   zrgb(3, 4) = 0.07518   ;   zrgb(4, 4) = 0.37859
+      zrgb(1, 5) =  0.016   ;   zrgb(2, 5) = 0.01782   ;   zrgb(3, 5) = 0.07539   ;   zrgb(4, 5) = 0.37879
+      zrgb(1, 6) =  0.018   ;   zrgb(2, 6) = 0.01831   ;   zrgb(3, 6) = 0.07562   ;   zrgb(4, 6) = 0.37900
+      zrgb(1, 7) =  0.020   ;   zrgb(2, 7) = 0.01885   ;   zrgb(3, 7) = 0.07586   ;   zrgb(4, 7) = 0.37923
+      zrgb(1, 8) =  0.022   ;   zrgb(2, 8) = 0.01943   ;   zrgb(3, 8) = 0.07613   ;   zrgb(4, 8) = 0.37948
+      zrgb(1, 9) =  0.025   ;   zrgb(2, 9) = 0.02005   ;   zrgb(3, 9) = 0.07641   ;   zrgb(4, 9) = 0.37976
+      zrgb(1,10) =  0.028   ;   zrgb(2,10) = 0.02073   ;   zrgb(3,10) = 0.07672   ;   zrgb(4,10) = 0.38005
+      zrgb(1,11) =  0.032   ;   zrgb(2,11) = 0.02146   ;   zrgb(3,11) = 0.07705   ;   zrgb(4,11) = 0.38036
+      zrgb(1,12) =  0.035   ;   zrgb(2,12) = 0.02224   ;   zrgb(3,12) = 0.07741   ;   zrgb(4,12) = 0.38070
+      zrgb(1,13) =  0.040   ;   zrgb(2,13) = 0.02310   ;   zrgb(3,13) = 0.07780   ;   zrgb(4,13) = 0.38107
+      zrgb(1,14) =  0.045   ;   zrgb(2,14) = 0.02402   ;   zrgb(3,14) = 0.07821   ;   zrgb(4,14) = 0.38146
+      zrgb(1,15) =  0.050   ;   zrgb(2,15) = 0.02501   ;   zrgb(3,15) = 0.07866   ;   zrgb(4,15) = 0.38189
+      zrgb(1,16) =  0.056   ;   zrgb(2,16) = 0.02608   ;   zrgb(3,16) = 0.07914   ;   zrgb(4,16) = 0.38235
+      zrgb(1,17) =  0.063   ;   zrgb(2,17) = 0.02724   ;   zrgb(3,17) = 0.07967   ;   zrgb(4,17) = 0.38285
+      zrgb(1,18) =  0.071   ;   zrgb(2,18) = 0.02849   ;   zrgb(3,18) = 0.08023   ;   zrgb(4,18) = 0.38338
+      zrgb(1,19) =  0.079   ;   zrgb(2,19) = 0.02984   ;   zrgb(3,19) = 0.08083   ;   zrgb(4,19) = 0.38396
+      zrgb(1,20) =  0.089   ;   zrgb(2,20) = 0.03131   ;   zrgb(3,20) = 0.08149   ;   zrgb(4,20) = 0.38458
+      zrgb(1,21) =  0.100   ;   zrgb(2,21) = 0.03288   ;   zrgb(3,21) = 0.08219   ;   zrgb(4,21) = 0.38526
+      zrgb(1,22) =  0.112   ;   zrgb(2,22) = 0.03459   ;   zrgb(3,22) = 0.08295   ;   zrgb(4,22) = 0.38598
+      zrgb(1,23) =  0.126   ;   zrgb(2,23) = 0.03643   ;   zrgb(3,23) = 0.08377   ;   zrgb(4,23) = 0.38676
+      zrgb(1,24) =  0.141   ;   zrgb(2,24) = 0.03842   ;   zrgb(3,24) = 0.08466   ;   zrgb(4,24) = 0.38761
+      zrgb(1,25) =  0.158   ;   zrgb(2,25) = 0.04057   ;   zrgb(3,25) = 0.08561   ;   zrgb(4,25) = 0.38852
+      zrgb(1,26) =  0.178   ;   zrgb(2,26) = 0.04289   ;   zrgb(3,26) = 0.08664   ;   zrgb(4,26) = 0.38950
+      zrgb(1,27) =  0.200   ;   zrgb(2,27) = 0.04540   ;   zrgb(3,27) = 0.08775   ;   zrgb(4,27) = 0.39056
+      zrgb(1,28) =  0.224   ;   zrgb(2,28) = 0.04811   ;   zrgb(3,28) = 0.08894   ;   zrgb(4,28) = 0.39171
+      zrgb(1,29) =  0.251   ;   zrgb(2,29) = 0.05103   ;   zrgb(3,29) = 0.09023   ;   zrgb(4,29) = 0.39294
+      zrgb(1,30) =  0.282   ;   zrgb(2,30) = 0.05420   ;   zrgb(3,30) = 0.09162   ;   zrgb(4,30) = 0.39428
+      zrgb(1,31) =  0.316   ;   zrgb(2,31) = 0.05761   ;   zrgb(3,31) = 0.09312   ;   zrgb(4,31) = 0.39572
+      zrgb(1,32) =  0.355   ;   zrgb(2,32) = 0.06130   ;   zrgb(3,32) = 0.09474   ;   zrgb(4,32) = 0.39727
+      zrgb(1,33) =  0.398   ;   zrgb(2,33) = 0.06529   ;   zrgb(3,33) = 0.09649   ;   zrgb(4,33) = 0.39894
+      zrgb(1,34) =  0.447   ;   zrgb(2,34) = 0.06959   ;   zrgb(3,34) = 0.09837   ;   zrgb(4,34) = 0.40075
+      zrgb(1,35) =  0.501   ;   zrgb(2,35) = 0.07424   ;   zrgb(3,35) = 0.10040   ;   zrgb(4,35) = 0.40270
+      zrgb(1,36) =  0.562   ;   zrgb(2,36) = 0.07927   ;   zrgb(3,36) = 0.10259   ;   zrgb(4,36) = 0.40480
+      zrgb(1,37) =  0.631   ;   zrgb(2,37) = 0.08470   ;   zrgb(3,37) = 0.10495   ;   zrgb(4,37) = 0.40707
+      zrgb(1,38) =  0.708   ;   zrgb(2,38) = 0.09056   ;   zrgb(3,38) = 0.10749   ;   zrgb(4,38) = 0.40952
+      zrgb(1,39) =  0.794   ;   zrgb(2,39) = 0.09690   ;   zrgb(3,39) = 0.11024   ;   zrgb(4,39) = 0.41216
+      zrgb(1,40) =  0.891   ;   zrgb(2,40) = 0.10374   ;   zrgb(3,40) = 0.11320   ;   zrgb(4,40) = 0.41502
+      zrgb(1,41) =  1.000   ;   zrgb(2,41) = 0.11114   ;   zrgb(3,41) = 0.11639   ;   zrgb(4,41) = 0.41809
+      zrgb(1,42) =  1.122   ;   zrgb(2,42) = 0.11912   ;   zrgb(3,42) = 0.11984   ;   zrgb(4,42) = 0.42142
+      zrgb(1,43) =  1.259   ;   zrgb(2,43) = 0.12775   ;   zrgb(3,43) = 0.12356   ;   zrgb(4,43) = 0.42500
+      zrgb(1,44) =  1.413   ;   zrgb(2,44) = 0.13707   ;   zrgb(3,44) = 0.12757   ;   zrgb(4,44) = 0.42887
+      zrgb(1,45) =  1.585   ;   zrgb(2,45) = 0.14715   ;   zrgb(3,45) = 0.13189   ;   zrgb(4,45) = 0.43304
+      zrgb(1,46) =  1.778   ;   zrgb(2,46) = 0.15803   ;   zrgb(3,46) = 0.13655   ;   zrgb(4,46) = 0.43754
+      zrgb(1,47) =  1.995   ;   zrgb(2,47) = 0.16978   ;   zrgb(3,47) = 0.14158   ;   zrgb(4,47) = 0.44240
+      zrgb(1,48) =  2.239   ;   zrgb(2,48) = 0.18248   ;   zrgb(3,48) = 0.14701   ;   zrgb(4,48) = 0.44765
+      zrgb(1,49) =  2.512   ;   zrgb(2,49) = 0.19620   ;   zrgb(3,49) = 0.15286   ;   zrgb(4,49) = 0.45331
+      zrgb(1,50) =  2.818   ;   zrgb(2,50) = 0.21102   ;   zrgb(3,50) = 0.15918   ;   zrgb(4,50) = 0.45942
+      zrgb(1,51) =  3.162   ;   zrgb(2,51) = 0.22703   ;   zrgb(3,51) = 0.16599   ;   zrgb(4,51) = 0.46601
+      zrgb(1,52) =  3.548   ;   zrgb(2,52) = 0.24433   ;   zrgb(3,52) = 0.17334   ;   zrgb(4,52) = 0.47313
+      zrgb(1,53) =  3.981   ;   zrgb(2,53) = 0.26301   ;   zrgb(3,53) = 0.18126   ;   zrgb(4,53) = 0.48080
+      zrgb(1,54) =  4.467   ;   zrgb(2,54) = 0.28320   ;   zrgb(3,54) = 0.18981   ;   zrgb(4,54) = 0.48909
+      zrgb(1,55) =  5.012   ;   zrgb(2,55) = 0.30502   ;   zrgb(3,55) = 0.19903   ;   zrgb(4,55) = 0.49803
+      zrgb(1,56) =  5.623   ;   zrgb(2,56) = 0.32858   ;   zrgb(3,56) = 0.20898   ;   zrgb(4,56) = 0.50768
+      zrgb(1,57) =  6.310   ;   zrgb(2,57) = 0.35404   ;   zrgb(3,57) = 0.21971   ;   zrgb(4,57) = 0.51810
+      zrgb(1,58) =  7.079   ;   zrgb(2,58) = 0.38154   ;   zrgb(3,58) = 0.23129   ;   zrgb(4,58) = 0.52934
+      zrgb(1,59) =  7.943   ;   zrgb(2,59) = 0.41125   ;   zrgb(3,59) = 0.24378   ;   zrgb(4,59) = 0.54147
+      zrgb(1,60) =  8.912   ;   zrgb(2,60) = 0.44336   ;   zrgb(3,60) = 0.25725   ;   zrgb(4,60) = 0.55457
+      zrgb(1,61) = 10.000   ;   zrgb(2,61) = 0.47804   ;   zrgb(3,61) = 0.27178   ;   zrgb(4,61) = 0.56870
+      !
+      prgb(:,:) = zrgb(2:4,:)
+      !
+      r_si2 = 1.e0 / zrgb(2, 1)        ! blue with the smallest chlorophyll concentration)
+      !!! O.M : essai type NEMO 3.2
+      IF(lwp) WRITE(numout,*) '     RGB longest depth of extinction    r_si2 set back to NEMO 3.2 value'
+      r_si2 = 62.0
+      IF(lwp) WRITE(numout,*) '      RGB longest depth of extinction    r_si2 = ', r_si2
+      !
+      DO jc = 1, 61                         ! check
+         zchl = zrgb(1,jc)
+         irgb = NINT( 41 + 20.* LOG10( zchl ) + 1.e-15 )
+         IF(lwp) WRITE(numout,*) '    jc =', jc, '  Chl = ', zchl, '  irgb = ', irgb
+         IF( irgb /= jc ) THEN
+            IF(lwp) WRITE(numout,*) '    jc =', jc, '  Chl = ', zchl, '  Chl class = ', irgb
+            CALL ctl_stop( 'trc_oce_rgb : inconsistency in Chl tabulated attenuation coeff.' )
+         ENDIF
+      END DO
+      !
+   END SUBROUTINE trc_oce_rgb
+
+
+   SUBROUTINE trc_oce_rgb_read( prgb )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE p4z_opt_init  ***
+      !!
+      !! ** Purpose :   Initialization of of the optical scheme
+      !!
+      !! ** Method  :   read the look up table for the optical coefficients
+      !!
+      !! ** input   :   xkrgb(61) precomputed array corresponding to the  
+      !!                          attenuation coefficient (from JM Andre)
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(3,61), INTENT(out) ::   prgb   ! tabulated attenuation coefficient
+      !!
+      INTEGER  ::   jc, jb ! dummy loop indice
+      INTEGER  ::   irgb   ! temporary integer
+      REAL(wp) ::   zchl   ! temporary scalar
+      INTEGER  ::   numlight
+      !!----------------------------------------------------------------------
+      !
+      IF(lwp) THEN                         ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) ' trc_oce_rgb_read : optical look-up table read in kRGB61.txt file'
+         WRITE(numout,*) ' ~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) 
+      ENDIF
+      !
+      CALL ctl_opn( numlight, 'kRGB61.txt', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )
+      DO jc = 1, 61
+         READ(numlight,*) zchl, ( prgb(jb,jc), jb = 1, 3 )
+         irgb = NINT( 41 + 20.* LOG10( zchl ) + 1.e-15 )   
+         IF(lwp) WRITE(numout,*) '    jc =', jc, '  Chl = ', zchl, '  irgb = ', irgb  
+         IF( irgb /= jc ) THEN  
+            IF(lwp) WRITE(numout,*) '    jc =', jc, '  Chl = ', zchl, '  Chl class = ', irgb
+            CALL ctl_stop( 'trc_oce_rgb_read : inconsistency in Chl tabulated attenuation coeff.' )
+         ENDIF
+      END DO
+      CLOSE( numlight )
+      !
+      r_si2 = 1.e0 / prgb(1, 1)      ! blue with the smallest chlorophyll concentration)
+      IF(lwp) WRITE(numout,*) '      RGB longest depth of extinction    r_si2 = ', r_si2
+      !
+   END SUBROUTINE trc_oce_rgb_read
+
+
+   FUNCTION trc_oce_ext_lev( prldex, pqsr_frc ) RESULT( pjl )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE trc_oce_ext_lev  ***
+      !!       
+      !! ** Purpose :   compute max. level for light penetration
+      !!          
+      !! ** Method  :   the function provides the level at which irradiance 
+      !!                becomes negligible (i.e. = 1.e-15 W/m2) for 3 or 2 bands light
+      !!                penetration: I(z) = pqsr_frc * EXP(hext/prldex) = 1.e-15 W/m2
+      !!                # prldex is the longest depth of extinction:
+      !!                   - prldex = 23 m (2 bands case)
+      !!                   - prldex = 62 m (3 bands case: blue waveband & 0.01 mg/m2 for the chlorophyll)
+      !!                # pqsr_frc is the fraction of solar radiation which penetrates,
+      !!                considering Qsr=240 W/m2 and rn_abs = 0.58:
+      !!                   - pqsr_frc = Qsr * (1-rn_abs)   = 1.00e2 W/m2 (2 bands case)
+      !!                   - pqsr_frc = Qsr * (1-rn_abs)/3 = 0.33e2 W/m2 (3 bands case & equi-partition)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), INTENT(in) ::   prldex    ! longest depth of extinction
+      REAL(wp), INTENT(in) ::   pqsr_frc  ! frac. solar radiation which penetrates 
+      !!
+      INTEGER  ::   jk, pjl            ! levels
+      REAL(wp) ::   zhext              ! deepest level till which light penetrates
+      REAL(wp) ::   zprec = 15._wp     ! precision to reach -LOG10(1.e-15)
+      REAL(wp) ::   zem                ! temporary scalar 
+      !!----------------------------------------------------------------------
+      !
+      ! It is not necessary to compute anything bellow the following depth
+      zhext = prldex * ( LOG(10._wp) * zprec + LOG(pqsr_frc) )
+      !
+      ! Level of light extinction
+      pjl = jpkm1
+      DO jk = jpkm1, 1, -1
+         zem = MAXVAL( fsdepw(:,:,jk+1) * tmask(:,:,jk) )
+         IF( zem >= zhext )   pjl = jk                       ! last T-level reached by Qsr
+      END DO
+      !
+   END FUNCTION trc_oce_ext_lev
+
+   !!======================================================================
+END MODULE trc_oce
Index: /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90
===================================================================
--- /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90	(revision 2189)
+++ /CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90	(revision 2189)
@@ -0,0 +1,868 @@
+MODULE zdftke
+   !!======================================================================
+   !!                       ***  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 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   zdf_tke        ! routine called in step module
+   PUBLIC   zdf_tke_init   ! routine called in opa module
+   PUBLIC   tke_rst        ! routine called in step module
+
+   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 4147 2013-11-04 11:51:55Z cetlod $
+   !! 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 zdf_tke( kt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_tke  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!              coefficients using a turbulent closure scheme (TKE).
+      !!
+      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
+      !!              is computed from a prognostic equation :
+      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
+      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
+      !!                  + avt N^2                      ! stratif. destruc.
+      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
+      !!      with the boundary conditions:
+      !!         surface: en = max( rn_emin0, rn_ebb * taum )
+      !!         bottom : en = rn_emin
+      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 
+      !!
+      !!        The now Turbulent kinetic energy is computed using the following 
+      !!      time stepping: implicit for vertical diffusion term, linearized semi
+      !!      implicit for kolmogoroff dissipation term, and explicit forward for 
+      !!      both buoyancy and shear production terms. Therefore a tridiagonal 
+      !!      linear system is solved. Note that buoyancy and shear terms are
+      !!      discretized in a energy conserving form (Bruchard 2002).
+      !!
+      !!        The dissipative and mixing length scale are computed from en and
+      !!      the stratification (see tke_avn)
+      !!
+      !!        The now vertical eddy vicosity and diffusivity coefficients are
+      !!      given by: 
+      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
+      !!              avt = max( avmb, pdl * avm                 )  
+      !!              eav = max( avmb, avm )
+      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
+      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1 
+      !!
+      !! ** Action  :   compute en (now turbulent kinetic energy)
+      !!                update avt, avmu, avmv (before vertical eddy coef.)
+      !!
+      !! References : Gaspar et al., JGR, 1990,
+      !!              Blanke and Delecluse, JPO, 1991
+      !!              Mellor and Blumberg, JPO 2004
+      !!              Axell, JGR, 2002
+      !!              Bruchard OM 2002
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time step
+      !!----------------------------------------------------------------------
+      !
+      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
+         avt (:,:,:) = avt_k (:,:,:) 
+         avm (:,:,:) = avm_k (:,:,:) 
+         avmu(:,:,:) = avmu_k(:,:,:) 
+         avmv(:,:,:) = avmv_k(:,:,:) 
+      ENDIF 
+      !
+      CALL tke_tke      ! now tke (en)
+      !
+      CALL tke_avn      ! now avt, avm, avmu, avmv
+      !
+      avt_k (:,:,:) = avt (:,:,:) 
+      avm_k (:,:,:) = avm (:,:,:) 
+      avmu_k(:,:,:) = avmu(:,:,:) 
+      avmv_k(:,:,:) = avmv(:,:,:) 
+      !
+   END SUBROUTINE zdf_tke
+
+
+   SUBROUTINE tke_tke
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE tke_tke  ***
+      !!
+      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
+      !!
+      !! ** Method  : - TKE surface boundary condition
+      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
+      !!              - source term due to shear (saved in avmu, avmv arrays)
+      !!              - Now TKE : resolution of the TKE equation by inverting
+      !!                a tridiagonal linear system by a "methode de chasse"
+      !!              - increase TKE due to surface and internal wave breaking
+      !!
+      !! ** Action  : - en : now turbulent kinetic energy)
+      !!              - avmu, avmv : production of TKE by shear at u and v-points
+      !!                (= Kz dz[Ub] * dz[Un] )
+      !! ---------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
+!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
+!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
+      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
+      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
+      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
+      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
+      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
+      REAL(wp) ::   ztau  , zdif                    !    -         -
+      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
+      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
+!!bfr      REAL(wp) ::   zebot                           !    -         -
+      INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
+      REAL(wp), POINTER, DIMENSION(:,:  ) :: zhlc
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw
+      !!--------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tke_tke')
+      !
+      CALL wrk_alloc( jpi,jpj, imlc )    ! integer
+      CALL wrk_alloc( jpi,jpj, zhlc ) 
+      CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      !
+      zbbrau = rn_ebb / rau0       ! Local constant initialisation
+      zfact1 = -.5_wp * rdt 
+      zfact2 = 1.5_wp * rdt * rn_ediss
+      zfact3 = 0.5_wp       * rn_ediss
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Surface boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+         END DO
+      END DO
+      
+!!bfr   - start commented area
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Bottom boundary condition on tke
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
+      ! The condition is coded here for completion but commented out until there is proof that the
+      ! computational cost is justified
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
+!CDIR NOVERRCHK
+!!    DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+!!       DO ji = fs_2, fs_jpim1   ! vector opt.
+!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
+!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
+!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
+!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
+!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
+!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
+!!       END DO
+!!    END DO
+!!bfr   - end commented area
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
+         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+         !
+         !                        !* total energy produce by LC : cumulative sum over jk
+         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
+         DO jk = 2, jpk
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
+         END DO
+         !                        !* finite Langmuir Circulation depth
+         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
+         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
+         DO jk = jpkm1, 2, -1
+            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us 
+               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
+                  zus  = zcof * taum(ji,jj)
+                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
+               END DO
+            END DO
+         END DO
+         !                               ! finite LC depth
+# if defined key_vectopt_loop
+         DO jj = 1, 1
+            DO ji = 1, jpij   ! vector opt. (forced unrolling)
+# else
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+# endif
+               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
+            END DO
+         END DO
+         zcof = 0.016 / SQRT( zrhoa * zcdrag )
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
+                  !                                           ! vertical velocity due to LC
+                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
+                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
+                  !                                           ! TKE Langmuir circulation source term
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Now Turbulent kinetic energy (output in en)
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
+      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
+      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
+      !
+      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
+         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
+            DO ji = 1, jpi
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
+                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
+                  &           / (  fse3uw_n(ji,jj,jk)         &
+                  &              * fse3uw_b(ji,jj,jk) )
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
+                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
+                  &                            / (  fse3vw_n(ji,jj,jk)               &
+                  &                              *  fse3vw_b(ji,jj,jk)  )
+            END DO
+         END DO
+      END DO
+      !
+      DO jk = 2, jpkm1           !* Matrix and right hand side in en
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zcof   = zfact1 * tmask(ji,jj,jk)
+               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
+                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
+               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
+                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
+                  !                                                           ! shear prod. at w-point weightened by mask
+               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+                  !
+               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
+               zd_lw(ji,jj,jk) = zzd_lw
+               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
+               !
+               !                                   ! right hand side in en
+               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
+                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
+      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
+         DO ji = fs_2, fs_jpim1    ! vector opt.
+            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
+         END DO
+      END DO
+      DO jk = 3, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
+            END DO
+         END DO
+      END DO
+      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
+         DO ji = fs_2, fs_jpim1    ! vector opt.
+            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
+         END DO
+      END DO
+      DO jk = jpk-2, 2, -1
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1    ! vector opt.
+               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      DO jk = 2, jpkm1                             ! set the minimum value of tke
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                            !  TKE due to surface and internal wave breaking
+      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               jk = nmln(ji,jj)
+               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                  &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
+            END DO
+         END DO
+      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
+                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
+                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )    ! module of the mean stress 
+                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean 
+                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
+                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
+                     &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      ENDIF
+      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
+      CALL wrk_dealloc( jpi,jpj, zhlc ) 
+      CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
+      !
+   END SUBROUTINE tke_tke
+
+
+   SUBROUTINE tke_avn
+      !!----------------------------------------------------------------------
+      !!                   ***  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, zpdlr, zri, zsqen     !   -      -
+      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
+      !!--------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tke_avn')
+
+      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  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)
+         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
+         zmxlm(:,:,1) = MAX(  rn_mxl0,  zraug * taum(:,:)  )
+      ELSE                          ! surface set to the minimum value
+         zmxlm(:,:,1) = rn_mxl0
+      ENDIF
+      zmxlm(:,:,jpk)  = rmxl_min     ! last level set to the interior minium value
+      !
+!CDIR NOVERRCHK
+      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
+               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
+            END DO
+         END DO
+      END DO
+      !
+      !                     !* Physical limits for the mixing length
+      !
+      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the zmxlm   value
+      zmxld(:,:,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 jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemxl = MIN( fsdepw(ji,jj,jk), 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 jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  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 jk = 2, jpkm1         ! from the surface to the bottom :
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  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 jk = 2, jpkm1         ! from the surface to the bottom : lup
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
+               END DO
+            END DO
+         END DO
+!CDIR NOVERRCHK
+         DO jk = 2, jpkm1
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
+                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
+                  zmxlm(ji,jj,jk) = zemlm
+                  zmxld(ji,jj,jk) = zemlp
+               END DO
+            END DO
+         END DO
+         !
+      END SELECT
+      !
+# if defined key_c1d
+      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
+      e_mix(:,:,:) = zmxlm(:,:,:)
+# endif
+
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
+      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+!CDIR NOVERRCHK
+      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               zsqen = SQRT( en(ji,jj,jk) )
+               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
+               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
+               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
+               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      !
+      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
+               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
+      !
+      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
+                  !                                          ! shear
+                  zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   &
+                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
+                  zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   &
+                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
+                  !                                          ! local Richardson number
+                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
+                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
+!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
+!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
+                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
+# if defined key_c1d
+                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
+                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
+# endif
+              END DO
+            END DO
+         END DO
+      ENDIF
+      CALL lbc_lnk( avt, '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,jpj,jpk, zmpdl, zmxlm, zmxld ) 
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
+      !
+   END SUBROUTINE tke_avn
+
+
+   SUBROUTINE zdf_tke_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_tke_init  ***
+      !!                     
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and 
+      !!              viscosity when using a tke turbulent closure scheme
+      !!
+      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
+      !!              called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namlist namzdf_tke
+      !!
+      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
+      !!----------------------------------------------------------------------
+      INTEGER ::   ji, jj, jk   ! dummy loop indices
+      INTEGER ::   ios
+      !!
+      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
+         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
+         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
+         &                 nn_etau , nn_htau  , rn_efr   
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
+      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
+      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
+      WRITE ( numond, namzdf_tke )
+      !
+      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
+      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
+      !!! O.M. 
+      rmxl_min = 1.0E-3
+      WRITE(numout,*) 'Set rmxl_min to IPSLCM5 value : ', rmxl_min
+      !
+      IF(lwp) THEN                    !* Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
+         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
+         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
+         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
+         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
+         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
+         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
+         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
+         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
+         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
+         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
+         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
+         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
+         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
+         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
+         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
+         WRITE(numout,*)
+         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
+      ENDIF
+      !
+      !                              ! allocate tke arrays
+      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
+      !
+      !                               !* Check of some namelist values
+      IF( nn_mxl  < 0  .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
+      IF( nn_pdl  < 0  .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
+      IF( nn_htau < 0  .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
+#if ! key_coupled
+      IF( nn_etau == 3 )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
+#endif
+
+      IF( ln_mxl0 ) THEN
+         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
+         rn_mxl0 = rmxl_min
+      ENDIF
+      
+      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln 
+
+      !                               !* depth of penetration of surface tke
+      IF( nn_etau /= 0 ) THEN      
+         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
+         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
+            htau(:,:) = 10._wp
+         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
+            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )            
+         END SELECT
+      ENDIF
+      !                               !* set vertical eddy coef. to the background value
+      DO jk = 1, jpk
+         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
+         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
+         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
+      END DO
+      dissl(:,:,:) = 1.e-12_wp
+      !                              
+      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
+      !
+   END SUBROUTINE zdf_tke_init
+
+
+   SUBROUTINE tke_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE tke_rst  ***
+     !!                     
+     !! ** Purpose :   Read or write TKE file (en) in restart file
+     !!
+     !! ** Method  :   use of IOM library
+     !!                if the restart does not contain TKE, en is either 
+     !!                set to rn_emin or recomputed 
+     !!----------------------------------------------------------------------
+     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
+     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+     !
+     INTEGER ::   jit, jk   ! dummy loop indices
+     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
+     !!----------------------------------------------------------------------
+     !
+     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise 
+        !                                   ! ---------------
+        IF( ln_rstart ) THEN                   !* Read the restart file
+           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
+           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
+           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
+           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
+           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
+           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
+           !
+           IF( id1 > 0 ) THEN                       ! 'en' exists
+              CALL iom_get( numror, jpdom_autoglo, 'en', en )
+              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
+                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
+                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
+                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
+              ELSE                                                 ! one at least array is missing
+                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
+              ENDIF
+           ELSE                                     ! No TKE array found: initialisation
+              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
+              en (:,:,:) = rn_emin * tmask(:,:,:)
+              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
+              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
+           ENDIF
+        ELSE                                   !* Start from rest
+           en(:,:,:) = rn_emin * tmask(:,:,:)
+           DO jk = 1, jpk                           ! set the Kz to the background value
+              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
+              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
+              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
+           END DO
+        ENDIF
+        !
+     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+        !                                   ! -------------------
+        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
+        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
+        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
+        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
+        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
+        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
+        !
+     ENDIF
+     !
+   END SUBROUTINE tke_rst
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Dummy module :                                        NO TKE scheme
+   !!----------------------------------------------------------------------
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
+CONTAINS
+   SUBROUTINE zdf_tke_init           ! Dummy routine
+   END SUBROUTINE zdf_tke_init
+   SUBROUTINE zdf_tke( kt )          ! Dummy routine
+      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
+   END SUBROUTINE zdf_tke
+   SUBROUTINE tke_rst( kt, cdrw )
+     CHARACTER(len=*) ::   cdrw
+     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
+   END SUBROUTINE tke_rst
+#endif
+
+   !!======================================================================
+END MODULE zdftke
