Index: anches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl_crs.F90
===================================================================
--- /branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl_crs.F90	(revision 7222)
+++ 	(revision )
@@ -1,135 +1,0 @@
-MODULE domvvl_crs
-   !!======================================================================
-   !!                       ***  MODULE domvvl   ***
-   !! Ocean : 
-   !!======================================================================
-   !! History :  2.0  !  2006-06  (B. Levier, L. Marie)  original code
-   !!            3.1  !  2009-02  (G. Madec, M. Leclair, R. Benshila)  pure z* coordinate
-   !!            3.3  !  2011-10  (M. Leclair) totally rewrote domvvl:
-   !!                                          vvl option includes z_star and z_tilde coordinates
-   !!----------------------------------------------------------------------
-   !!   'key_vvl'                              variable volume
-   !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   dom_vvl_init     : define initial vertical scale factors, depths and column thickness
-   !!   dom_vvl_sf_nxt   : Compute next vertical scale factors
-   !!   dom_vvl_sf_swp   : Swap vertical scale factors and update the vertical grid
-   !!   dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another
-   !!   dom_vvl_rst      : read/write restart file
-   !!   dom_vvl_ctl      : Check the vvl options
-   !!   dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors 
-   !!                    : to account for manual changes to e[1,2][u,v] in some Straits 
-   !!----------------------------------------------------------------------
-   !! * Modules used
-!   USE oce             ! ocean dynamics and tracers
-   USE crs             ! ocean space and time domain
-   USE sbc_oce         ! ocean surface boundary condition
-   USE in_out_manager  ! I/O manager
-   USE iom             ! I/O manager library
-   USE restart         ! ocean restart
-   USE lib_mpp         ! distributed memory computing library
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE wrk_nemo        ! Memory allocation
-   USE timing          ! Timing
-   USE crslbclnk
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Routine accessibility
-   PUBLIC  dom_vvl_interpol_crs   ! called by dynnxt.F90
-
-   !                                                                                            ! conservation: not used yet
-
-   !! * Module variables
-!#  include "domzgr_substitute.h90"
-!#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 
-   !! $Id: domvvl.F90 4998 2014-12-22 16:16:45Z mocavero $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE dom_vvl_interpol_crs( pe3_in, pe3_out, pout )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE dom_vvl__interpol  ***
-      !!
-      !! ** Purpose :   interpolate scale factors from one grid point to another
-      !!
-      !! ** Method  :   e3_out = e3_0 + interpolation(e3_in - e3_0)
-      !!                - horizontal interpolation: grid cell surface averaging
-      !!                - vertical interpolation: simple averaging
-      !!----------------------------------------------------------------------
-      !! * Arguments
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT( in    ) ::  pe3_in     ! input e3 to be interpolated
-      REAL(wp), DIMENSION(jpi_crs,jpj_crs,jpk), INTENT( inout ) ::  pe3_out    ! output interpolated e3
-      CHARACTER(LEN=*), INTENT( in )                    ::  pout       ! grid point of out scale factors
-      !                                                                !   =  'U', 'V', 'W, 'F', 'UW' or 'VW'
-      !! * Local declarations
-      INTEGER ::   ji, jj, jk                                          ! dummy loop indices
-      REAL(wp)::   zz 
-      !!----------------------------------------------------------------------
-      IF( nn_timing == 1 )  CALL timing_start('dom_vvl_interpol')
-         !
-      SELECT CASE ( pout )
-         !               ! ------------------------------------- !
-      CASE( 'U' )        ! interpolation from T-point to U-point !
-         !               ! ------------------------------------- !
-         ! horizontal surface weighted interpolation
-         DO jk = 1, jpk
-            DO jj = 1, jpj_crs-1
-               DO ji = 1, jpi_crs-1   ! vector opt.
-                  zz=1._wp/(e1u_crs(ji,jj)*e2u_crs(ji,jj))
-                  pe3_out(ji,jj,jk) = 0.5_wp * umask_crs(ji,jj,jk) * zz                                           &
-                     &                       * (   e1t_crs(ji  ,jj) * e2t_crs(ji  ,jj) * ( pe3_in(ji  ,jj,jk) - e3t_0_crs(ji  ,jj,jk) )     &
-                     &                           + e1t_crs(ji+1,jj) * e2t_crs(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0_crs(ji+1,jj,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! boundary conditions
-         CALL crs_lbc_lnk( pe3_out(:,:,:), 'U', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0_crs(:,:,:)
-         !               ! ------------------------------------- !
-      CASE( 'V' )        ! interpolation from T-point to V-point !
-         !               ! ------------------------------------- !
-         ! horizontal surface weighted interpolation
-         DO jk = 1, jpk
-            DO jj = 1, jpj_crs
-               DO ji = 1, jpi_crs   ! vector opt.
-                  zz=1._wp/(e1v_crs(ji,jj)*e2v_crs(ji,jj))
-                  pe3_out(ji,jj,jk) = 0.5_wp * vmask_crs(ji,jj,jk) * zz                                               &
-                     &                       * (   e1t_crs(ji,jj  ) * e2t_crs(ji,jj  ) * ( pe3_in(ji,jj  ,jk) - e3t_0(ji,jj  ,jk) )     &
-                     &                           + e1t_crs(ji,jj+1) * e2t_crs(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) )
-               END DO
-            END DO
-         END DO
-         !
-         ! boundary conditions
-         CALL crs_lbc_lnk( pe3_out(:,:,:), 'V', 1._wp )
-         pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0_crs(:,:,:)
-         !               ! ------------------------------------- !
-      CASE( 'W' )        ! interpolation from T-point to W-point !
-         !               ! ------------------------------------- !
-         ! vertical simple interpolation
-         pe3_out(:,:,1) = e3w_0_crs(:,:,1) + pe3_in(:,:,1) - e3t_0_crs(:,:,1)
-         ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing
-         DO jk = 2, jpk
-            pe3_out(:,:,jk) = e3w_0_crs(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask_crs(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0_crs(:,:,jk-1) )   &
-               &                                +            0.5_wp * tmask_crs(:,:,jk)   * ( pe3_in(:,:,jk  ) - e3t_0_crs(:,:,jk  ) )
-         END DO
-      END SELECT
-      !
-
-      IF( nn_timing == 1 )  CALL timing_stop('dom_vvl_interpol')
-
-   END SUBROUTINE dom_vvl_interpol_crs
-
-
-   !!======================================================================
-END MODULE domvvl_crs
-
-
-
Index: anches/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 7222)
+++ 	(revision )
@@ -1,126 +1,0 @@
-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: anches/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 7222)
+++ 	(revision )
@@ -1,380 +1,0 @@
-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) + fse3t_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) + fse3t_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 * fse3w_crs(ji,jj,jk) * fse3w_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
