Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2374)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2375)
@@ -544,32 +544,25 @@
    rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
    rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
-   rn_ebb      =   3.75    !  coef. of the surface input of tke
-   rn_emin     =   1.e-5   !  minimum value of tke [m2/s2]
+   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
+   rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
    rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
-   rn_bshear   =   1.e-20  !  background shear (>0)
    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                            !                 = 1 bounded by the local vertical scale factor
                            !                 = 2 first vertical derivative of mixing length bounded by 1
-                           !                 = 3 same criteria as case 2 but applied in a different way
+                           !                 = 3 as =2 with distinct disspipative an mixing length scale
    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
-   ln_mxl0     = .false.   !  mixing length scale surface value as function of wind stress (T) or not (F)
-   rn_lmin     =   0.4     !  interior buoyancy lenght scale minimum value
-   rn_lmin0    =   0.4     !  surface  buoyancy lenght scale minimum value
-   nn_etau     =   0       !  exponentially deceasing penetration of tke due to internal & intertial waves
-                           !        = 0 no penetration ( O(2 km) resolution)
-                           !        = 1 additional tke source (rn_efr * en)
-                           !        = 2 additional tke source (rn_efr * en) applied only at the base of the mixed layer
-                           !        = 3 additional tke source (HF contribution: mean of stress module - module of mean stress)
-   nn_htau     =   1       !  type of exponential decrease of tke penetration
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
+   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
+   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+   nn_etau     =   0       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+                           !        = 0 no penetration
+                           !        = 1 add a tke source below the ML
+                           !        = 2 add a tke source just at the base of the ML
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                            !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m at high latitudes
-                           !        = 2  30 meters constant depth penetration
-                           !  otion used only if nn_etau = 1 or 2: 
-   rn_efr      =   0.05    !     fraction of surface tke value which penetrates inside the ocean
-                           !  otion used only if nn_etau = 3:
-   rn_addhft   =  -1.e-3   !     add offset   applied to the "mean of stress module - module of mean stress" (always kept > 0)
-   rn_sclhft   =   1.      !     scale factor applied to the "mean of stress module - module of mean stress"
-   ln_lc       = .false.   !  Langmuir cell parameterisation
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
 !------------------------------------------------------------------------
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 2374)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 2375)
@@ -544,32 +544,25 @@
    rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
    rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
-   rn_ebb      =   3.75    !  coef. of the surface input of tke
-   rn_emin     =   1.e-5   !  minimum value of tke [m2/s2]
+   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
+   rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
    rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
-   rn_bshear   =   1.e-20  !  background shear (>0)
    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                            !                 = 1 bounded by the local vertical scale factor
                            !                 = 2 first vertical derivative of mixing length bounded by 1
-                           !                 = 3 same criteria as case 2 but applied in a different way
+                           !                 = 3 as =2 with distinct disspipative an mixing length scale
    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
-   ln_mxl0     = .false.   !  mixing length scale surface value as function of wind stress (T) or not (F)
-   rn_lmin     =   0.4     !  interior buoyancy lenght scale minimum value
-   rn_lmin0    =   0.4     !  surface  buoyancy lenght scale minimum value
-   nn_etau     =   0       !  exponentially deceasing penetration of tke due to internal & intertial waves
-                           !        = 0 no penetration ( O(2 km) resolution)
-                           !        = 1 additional tke source (rn_efr * en)
-                           !        = 2 additional tke source (rn_efr * en) applied only at the base of the mixed layer
-                           !        = 3 additional tke source (HF contribution: mean of stress module - module of mean stress)
-   nn_htau     =   1       !  type of exponential decrease of tke penetration
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
+   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
+   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+   nn_etau     =   0       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+                           !        = 0 no penetration
+                           !        = 1 add a tke source below the ML
+                           !        = 2 add a tke source just at the base of the ML
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                            !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m at high latitudes
-                           !        = 2  30 meters constant depth penetration
-                           !  otion used only if nn_etau = 1 or 2: 
-   rn_efr      =   0.05    !     fraction of surface tke value which penetrates inside the ocean
-                           !  otion used only if nn_etau = 3:
-   rn_addhft   =  -1.e-3   !     add offset   applied to the "mean of stress module - module of mean stress" (always kept > 0)
-   rn_sclhft   =   1.      !     scale factor applied to the "mean of stress module - module of mean stress"
-   ln_lc       = .false.   !  Langmuir cell parameterisation
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
 !------------------------------------------------------------------------
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2374)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2375)
@@ -552,32 +552,25 @@
    rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
    rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
-   rn_ebb      =  60.      !  coef. of the surface input of tke
+   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
    rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
    rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
-   rn_bshear   =   1.e-20  !  background shear (>0)
    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                            !                 = 1 bounded by the local vertical scale factor
                            !                 = 2 first vertical derivative of mixing length bounded by 1
-                           !                 = 3 same criteria as case 2 but applied in a different way
+                           !                 = 3 as =2 with distinct disspipative an mixing length scale
    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
-   ln_mxl0     = .false.   !  mixing length scale surface value as function of wind stress (T) or not (F)
-   rn_lmin     =   0.001   !  interior buoyancy lenght scale minimum value
-   rn_lmin0    =   0.01    !  surface  buoyancy lenght scale minimum value
-   nn_etau     =   0       !  exponentially deceasing penetration of tke due to internal & intertial waves
-                           !        = 0 no penetration ( O(2 km) resolution)
-                           !        = 1 additional tke source (rn_efr * en)
-                           !        = 2 additional tke source (rn_efr * en) applied only at the base of the mixed layer
-                           !        = 3 additional tke source (HF contribution: mean of stress module - module of mean stress)
-   nn_htau     =   1       !  type of exponential decrease of tke penetration
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
+   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
+   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+   nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+                           !        = 0 no penetration
+                           !        = 1 add a tke source below the ML
+                           !        = 2 add a tke source just at the base of the ML
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                            !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m at high latitudes
-                           !        = 2  30 meters constant depth penetration
-                           !  otion used only if nn_etau = 1 or 2: 
-   rn_efr      =   0.05    !     fraction of surface tke value which penetrates inside the ocean
-                           !  otion used only if nn_etau = 3:
-   rn_addhft   =  -1.e-3   !     add offset   applied to the "mean of stress module - module of mean stress" (always kept > 0)
-   rn_sclhft   =   1.      !     scale factor applied to the "mean of stress module - module of mean stress"
-   ln_lc       = .false.   !  Langmuir cell parameterisation
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
 !------------------------------------------------------------------------
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2374)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2375)
@@ -589,32 +589,25 @@
    rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
    rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
-   rn_ebb      =  60.      !  coef. of the surface input of tke
+   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
    rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
    rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
-   rn_bshear   =   1.e-20  !  background shear (>0)
    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                            !                 = 1 bounded by the local vertical scale factor
                            !                 = 2 first vertical derivative of mixing length bounded by 1
-                           !                 = 3 same criteria as case 2 but applied in a different way
+                           !                 = 3 as =2 with distinct disspipative an mixing length scale
    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
-   ln_mxl0     = .false.   !  mixing length scale surface value as function of wind stress (T) or not (F)
-   rn_lmin     =   0.001   !  interior buoyancy lenght scale minimum value
-   rn_lmin0    =   0.01    !  surface  buoyancy lenght scale minimum value
-   nn_etau     =   0       !  exponentially deceasing penetration of tke due to internal & intertial waves
-                           !        = 0 no penetration ( O(2 km) resolution)
-                           !        = 1 additional tke source (rn_efr * en)
-                           !        = 2 additional tke source (rn_efr * en) applied only at the base of the mixed layer
-                           !        = 3 additional tke source (HF contribution: mean of stress module - module of mean stress)
-   nn_htau     =   1       !  type of exponential decrease of tke penetration
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
+   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
+   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+   nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+                           !        = 0 no penetration
+                           !        = 1 add a tke source below the ML
+                           !        = 2 add a tke source just at the base of the ML
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                            !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m at high latitudes
-                           !        = 2  30 meters constant depth penetration
-                           !  otion used only if nn_etau = 1 or 2: 
-   rn_efr      =   0.05    !     fraction of surface tke value which penetrates inside the ocean
-                           !  otion used only if nn_etau = 3:
-   rn_addhft   =  -1.e-3   !     add offset   applied to the "mean of stress module - module of mean stress" (always kept > 0)
-   rn_sclhft   =   1.      !     scale factor applied to the "mean of stress module - module of mean stress"
-   ln_lc       = .false.   !  Langmuir cell parameterisation
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
 !------------------------------------------------------------------------
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 2374)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 2375)
@@ -589,32 +589,25 @@
    rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
    rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
-   rn_ebb      =  60.      !  coef. of the surface input of tke
+   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
    rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
    rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
-   rn_bshear   =   1.e-20  !  background shear (>0)
    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                            !                 = 1 bounded by the local vertical scale factor
                            !                 = 2 first vertical derivative of mixing length bounded by 1
-                           !                 = 3 same criteria as case 2 but applied in a different way
+                           !                 = 3 as =2 with distinct disspipative an mixing length scale
    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
-   ln_mxl0     = .false.   !  mixing length scale surface value as function of wind stress (T) or not (F)
-   rn_lmin     =   0.001   !  interior buoyancy lenght scale minimum value
-   rn_lmin0    =   0.01    !  surface  buoyancy lenght scale minimum value
-   nn_etau     =   0       !  exponentially deceasing penetration of tke due to internal & intertial waves
-                           !        = 0 no penetration ( O(2 km) resolution)
-                           !        = 1 additional tke source (rn_efr * en)
-                           !        = 2 additional tke source (rn_efr * en) applied only at the base of the mixed layer
-                           !        = 3 additional tke source (HF contribution: mean of stress module - module of mean stress)
-   nn_htau     =   1       !  type of exponential decrease of tke penetration
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
+   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
+   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+   nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+                           !        = 0 no penetration
+                           !        = 1 add a tke source below the ML
+                           !        = 2 add a tke source just at the base of the ML
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                            !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m at high latitudes
-                           !        = 2  30 meters constant depth penetration
-                           !  otion used only if nn_etau = 1 or 2: 
-   rn_efr      =   0.05    !     fraction of surface tke value which penetrates inside the ocean
-                           !  otion used only if nn_etau = 3:
-   rn_addhft   =  -1.e-3   !     add offset   applied to the "mean of stress module - module of mean stress" (always kept > 0)
-   rn_sclhft   =   1.      !     scale factor applied to the "mean of stress module - module of mean stress"
-   ln_lc       = .false.   !  Langmuir cell parameterisation
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
 !------------------------------------------------------------------------
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2374)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2375)
@@ -589,32 +589,25 @@
    rn_ediff    =   0.1     !  coef. for vertical eddy coef. (avt=rn_ediff*mxl*sqrt(e) )
    rn_ediss    =   0.7     !  coef. of the Kolmogoroff dissipation
-   rn_ebb      =  60.      !  coef. of the surface input of tke
+   rn_ebb      =  67.83    !  coef. of the surface input of tke (=67.83 suggested when ln_mxl0=T)
    rn_emin     =   1.e-6   !  minimum value of tke [m2/s2]
    rn_emin0    =   1.e-4   !  surface minimum value of tke [m2/s2]
-   rn_bshear   =   1.e-20  !  background shear (>0)
    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom
                            !                 = 1 bounded by the local vertical scale factor
                            !                 = 2 first vertical derivative of mixing length bounded by 1
-                           !                 = 3 same criteria as case 2 but applied in a different way
+                           !                 = 3 as =2 with distinct disspipative an mixing length scale
    nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm)
-   ln_mxl0     = .false.   !  mixing length scale surface value as function of wind stress (T) or not (F)
-   rn_lmin     =   0.001   !  interior buoyancy lenght scale minimum value
-   rn_lmin0    =   0.01    !  surface  buoyancy lenght scale minimum value
-   nn_etau     =   0       !  exponentially deceasing penetration of tke due to internal & intertial waves
-                           !        = 0 no penetration ( O(2 km) resolution)
-                           !        = 1 additional tke source (rn_efr * en)
-                           !        = 2 additional tke source (rn_efr * en) applied only at the base of the mixed layer
-                           !        = 3 additional tke source (HF contribution: mean of stress module - module of mean stress)
-   nn_htau     =   1       !  type of exponential decrease of tke penetration
+   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F)
+   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value
+   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002)
+   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+   nn_etau     =   1       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+                           !        = 0 no penetration
+                           !        = 1 add a tke source below the ML
+                           !        = 2 add a tke source just at the base of the ML
+                           !        = 3 as = 1 applied on HF part of the stress    ("key_coupled")
+   rn_efr      =   0.05    !  fraction of surface tke value which penetrates below the ML (nn_etau=1 or 2)
+   nn_htau     =   1       !  type of exponential decrease of tke penetration below the ML
                            !        = 0  constant 10 m length scale
-                           !        = 1  0.5m at the equator to 30m at high latitudes
-                           !        = 2  30 meters constant depth penetration
-                           !  otion used only if nn_etau = 1 or 2: 
-   rn_efr      =   0.05    !     fraction of surface tke value which penetrates inside the ocean
-                           !  otion used only if nn_etau = 3:
-   rn_addhft   =  -1.e-3   !     add offset   applied to the "mean of stress module - module of mean stress" (always kept > 0)
-   rn_sclhft   =   1.      !     scale factor applied to the "mean of stress module - module of mean stress"
-   ln_lc       = .false.   !  Langmuir cell parameterisation
-   rn_lc       =   0.15    !  coef. associated to Langmuir cells
+                           !        = 1  0.5m at the equator to 30m poleward of 40 degrees
 /
 !------------------------------------------------------------------------
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 2374)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 2375)
@@ -61,34 +61,35 @@
 
 #if defined key_c1d
-   !                                                                !!* 1D cfg only *
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers
+   !                                                           !!** 1D cfg only  **   ('key_c1d')
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
 #endif
 
-   !                                       !!! ** Namelist  namzdf_tke  **
-   LOGICAL  ::   ln_mxl0   = .FALSE.         ! mixing length scale surface value as function of wind stress or not
-   INTEGER  ::   nn_mxl    =  2              ! type of mixing length (=0/1/2/3)
-   REAL(wp) ::   rn_lmin0  = 0.4_wp          ! surface  min value of mixing length   [m]
-   REAL(wp) ::   rn_lmin   = 0.1_wp          ! interior min value of mixing length   [m]
-   INTEGER  ::   nn_pdl    =  1              ! Prandtl number or not (ratio avt/avm) (=0/1)
-   REAL(wp) ::   rn_ediff  = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
-   REAL(wp) ::   rn_ediss  = 0.7_wp          ! coefficient of the Kolmogoroff dissipation 
-   REAL(wp) ::   rn_ebb    = 3.75_wp         ! coefficient of the surface input of tke
-   REAL(wp) ::   rn_emin   = 0.7071e-6_wp    ! minimum value of tke           [m2/s2]
-   REAL(wp) ::   rn_emin0  = 1.e-4_wp        ! surface minimum value of tke   [m2/s2]
-   REAL(wp) ::   rn_bshear = 1.e-20          ! background shear (>0)
-   INTEGER  ::   nn_etau   = 0               ! type of depth penetration of surface tke (=0/1/2)
-   INTEGER  ::   nn_htau   = 0               ! type of tke profile of penetration (=0/1)
-   REAL(wp) ::   rn_efr    = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean
-   REAL(wp) ::   rn_addhft = 0.0_wp          ! add offset   applied to HF tau
-   REAL(wp) ::   rn_sclhft = 1.0_wp          ! scale factor applied to HF tau
-   LOGICAL  ::   ln_lc     = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
-   REAL(wp) ::   rn_lc     = 0.15_wp         ! 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), DIMENSION(jpi,jpj)     ::   htau      ! depth of tke penetration (nn_htau)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en        ! now turbulent kinetic energy   [m2/s2]
-   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dissl     ! now mixing lenght of dissipation
+   !                                      !!** Namelist  namzdf_tke  **
+   LOGICAL  ::   ln_mxl0   = .FALSE.       ! mixing length scale surface value as function of wind stress or not
+   INTEGER  ::   nn_mxl    =  2            ! type of mixing length (=0/1/2/3)
+   REAL(wp) ::   rn_mxl0   = 0.04_wp       ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
+   INTEGER  ::   nn_pdl    =  1            ! Prandtl number or not (ratio avt/avm) (=0/1)
+   REAL(wp) ::   rn_ediff  = 0.1_wp        ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
+   REAL(wp) ::   rn_ediss  = 0.7_wp        ! coefficient of the Kolmogoroff dissipation 
+   REAL(wp) ::   rn_ebb    = 3.75_wp       ! coefficient of the surface input of tke
+   REAL(wp) ::   rn_emin   = 0.7071e-6_wp  ! minimum value of tke           [m2/s2]
+   REAL(wp) ::   rn_emin0  = 1.e-4_wp      ! surface minimum value of tke   [m2/s2]
+   REAL(wp) ::   rn_bshear = 1.e-20_wp     ! background shear (>0) currently a numerical threshold (do not change it)
+   INTEGER  ::   nn_etau   = 0             ! type of depth penetration of surface tke (=0/1/2/3)
+   INTEGER  ::   nn_htau   = 0             ! type of tke profile of penetration (=0/1)
+   REAL(wp) ::   rn_efr    = 1.0_wp        ! fraction of TKE surface value which penetrates in the ocean
+   LOGICAL  ::   ln_lc     = .FALSE.       ! Langmuir cells (LC) as a source term of TKE or not
+   REAL(wp) ::   rn_lc     = 0.15_wp       ! 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), DIMENSION(jpi,jpj,jpk), PUBLIC ::   en   ! now turbulent kinetic energy   [m2/s2]
+   
+   REAL(wp), DIMENSION(jpi,jpj)     ::   htau    ! depth of tke penetration (nn_htau)
+   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dissl   ! now mixing lenght of dissipation
 
    !! * Substitutions
@@ -98,7 +99,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -165,5 +165,5 @@
       !!
       !! ** Method  : - TKE surface boundary condition
-      !!              - source term due to Langmuir cells (ln_lc=T)
+      !!              - 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
@@ -197,7 +197,7 @@
       !
       zbbrau = rn_ebb / rau0       ! Local constant initialisation
-      zfact1 = -.5 * rdt 
-      zfact2 = 1.5 * rdt * rn_ediss
-      zfact3 = 0.5       * rn_ediss
+      zfact1 = -.5_wp * rdt 
+      zfact2 = 1.5_wp * rdt * rn_ediss
+      zfact3 = 0.5_wp       * rn_ediss
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -245,11 +245,11 @@
       !
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke
+      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. ) * fsdepw(:,:,1) * fse3w(:,:,1)
+         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
          DO jk = 2, jpk
-            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
+            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
          END DO
          !                        !* finite Langmuir Circulation depth
@@ -328,10 +328,10 @@
                   &          / ( 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.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
-                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )    
+               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.e0 - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
+               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
                !
                !                                   ! right hand side in en
@@ -384,22 +384,22 @@
       !                            !  TKE due to surface and internal wave breaking
       !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( nn_etau == 1 ) THEN           !* penetration throughout the water column
+      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.e0 - 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)
+                     &                                               * ( 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.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
-            END DO
-         END DO
-      ELSEIF( nn_etau == 3 ) THEN       !* penetration throughout the water column
+                  &                                               * ( 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
@@ -410,9 +410,9 @@
                   ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
                   zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
-                  ztau = 0.5 * SQRT( ztx2 * ztx2 + zty2 * zty2 )      ! module of the mean stress 
-                  zdif = taum(ji,jj) - ztau                           ! mean of the module - module of the mean 
-                  zdif = rn_sclhft * MAX( 0.e0, zdif + rn_addhft )    ! apply some modifications...
+                  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.e0 - fr_i(ji,jj) ) * tmask(ji,jj,jk)
+                     &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
                END DO
             END DO
@@ -439,5 +439,5 @@
       !!                      mxl = sqrt(2*en) / N  
       !!      where N is the brunt-vaisala frequency, with a minimum value set
-      !!      to rn_lmin (rn_lmin0) in the interior (surface) ocean.
+      !!      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.
@@ -479,21 +479,14 @@
       !                     !* 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)
-!!gm  this should be useless
-         zmxlm(:,:,1) = 0.e0
-!!gm end
-         zraug = vkarmn * 2.e5 / ( rau0 * grav )
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               zmxlm(ji,jj,1) = MAX(  rn_lmin0,  zraug * taum(ji,jj)  )
-            END DO
-         END DO
-      ELSE                       ! surface set to the minimum value
-         zmxlm(:,:,1) = rn_lmin0
+      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(ji,jj)  )
+      ELSE                          ! surface set to the minimum value
+         zmxlm(:,:,1) = rn_mxl0
       ENDIF
-      zmxlm(:,:,jpk) = rn_lmin   ! bottom set to the interior minium value
-      !
-!CDIR NOVERRCHK
-      DO jk = 2, jpkm1           ! interior value : l=sqrt(2*e/n**2)
+      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
@@ -501,15 +494,14 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                zrn2 = MAX( rn2(ji,jj,jk), rsmall )
-               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
-            END DO
-         END DO
-      END DO
-      !
-!!gm  CAUTION: to be added here the bottom boundary condition on zmxlm
+               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 minimum value
-      zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value
+      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 )
@@ -625,5 +617,5 @@
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1   ! vector opt.
-                  zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
+                  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) )   &
@@ -632,8 +624,8 @@
                     &  + 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. ) * avm(ji,jj,jk) / (zcoef * (zdku + zdkv + rn_bshear ) )
-                  zpdlr = MAX(  0.1,  0.2 / MAX( 0.2 , zri )  )
+                  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,  ri_crit / MAX( ri_crit , zri )  )
+!!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
@@ -672,9 +664,8 @@
       INTEGER ::   ji, jj, jk   ! dummy loop indices
       !!
-      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,   &
-         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,   &
-         &                 rn_lmin , rn_lmin0 , nn_pdl   , nn_etau  ,   &
-         &                 nn_htau , rn_efr   , rn_addhft, rn_sclhft,   &
-         &                 ln_lc   , rn_lc 
+      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   
       !!----------------------------------------------------------------------
 
@@ -682,5 +673,6 @@
       READ   ( numnam, namzdf_tke )
       
-      ri_cri = 2. / ( 2. + rn_ediss / rn_ediff )      ! resulting critical Richardson number
+      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
 
       IF(lwp) THEN                    !* Control print
@@ -698,13 +690,10 @@
          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_lmin0  = ', rn_lmin0
-         WRITE(numout,*) '      interior mixing length minimum value        rn_lmin0  = ', rn_lmin0
+         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,*) '      add offset   applied to HF tau              rn_addhft = ', rn_addhft
-         WRITE(numout,*) '      scale factor applied to HF tau              rn_sclhft = ', rn_sclhft
-         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,*)
          WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
@@ -712,11 +701,16 @@
 
       !                               !* 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 > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
-      IF( rn_lc   < 0.15  .OR.  rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' )
-      IF( rn_sclhft == 0. .AND. nn_etau == 3  )   CALL ctl_stop( 'force null HF tau to penetrate the thermocline...' )
-      IF( .NOT. lhftau    .AND. nn_etau == 3  )   CALL ctl_stop( 'bad flag: nn_etau == 3 must be used with HF tau' )
-
+      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 
 
@@ -724,14 +718,8 @@
       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.e0
-         CASE( 1 )                                    ! F(latitude) : 0.5m to 30m at high lat.
-            DO jj = 1, jpj
-               DO ji = 1, jpi  
-                  htau(ji,jj) = MAX( 0.5, 3./4. * MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
-               END DO
-            END DO
-         CASE( 2 )                                    ! constant depth penetration (here 30 meters)
-            htau(:,:) = 30.e0
+         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
@@ -744,5 +732,5 @@
          avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
       END DO
-      dissl(:,:,:) = 1.e-12
+      dissl(:,:,:) = 1.e-12_wp
       !                               !* read or initialize all required files 
       CALL tke_rst( nit000, 'READ' )
