Changeset 1060


Ignore:
Timestamp:
2008-06-05T11:47:45+02:00 (13 years ago)
Author:
ctlod
Message:

trunk: follow DOCTOR rules for TKE namelist parameters names and zdftke.F90 style code review, see ticket: #188

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/NEMO/OPA_SRC/ZDF/zdftke.F90

    r1050 r1060  
    22   !!====================================================================== 
    33   !!                       ***  MODULE  zdftke  *** 
    4    !! Ocean physics:  vertical mixing coefficient compute from the tke  
     4   !! Ocean physics:  vertical mixing coefficient computed from the tke  
    55   !!                 turbulent closure parameterization 
    66   !!===================================================================== 
    7    !! History :   6.0  !  91-03  (b. blanke)  Original code 
    8    !!             7.0  !  91-11  (G. Madec)   bug fix 
    9    !!             7.1  !  92-10  (G. Madec)   new mixing length and eav 
    10    !!             7.2  !  93-03  (M. Guyon)   symetrical conditions 
    11    !!             7.3  !  94-08  (G. Madec, M. Imbard)   npdl flag 
    12    !!             7.5  !  96-01  (G. Madec)   s-coordinates 
    13    !!             8.0  !  97-07  (G. Madec)   lbc 
    14    !!             8.1  !  99-01  (E. Stretta) new option for the mixing length 
    15    !!             8.5  !  02-06  (G. Madec) add zdf_tke_init routine 
    16    !!             8.5  !  02-08  (G. Madec)  ri_c and Free form, F90 
    17    !!             9.0  !  04-10  (C. Ethe )  1D configuration 
    18    !!             9.0  !  02-08  (G. Madec)  autotasking optimization 
    19    !!             9.0  !  06-07  (S. Masson)  distributed restart using iom 
    20    !!              -   !  2005-07  (C. Ethe,  G.Madec) : update TKE physics: 
     7   !! History :   OPA  !  1991-03  (b. blanke)  Original code 
     8   !!             7.0  !  1991-11  (G. Madec)   bug fix 
     9   !!             7.1  !  1992-10  (G. Madec)   new mixing length and eav 
     10   !!             7.2  !  1993-03  (M. Guyon)   symetrical conditions 
     11   !!             7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag 
     12   !!             7.5  !  1996-01  (G. Madec)   s-coordinates 
     13   !!             8.0  !  1997-07  (G. Madec)   lbc 
     14   !!             8.1  !  1999-01  (E. Stretta) new option for the mixing length 
     15   !!  NEMO       1.0  !  2002-06  (G. Madec) add zdf_tke_init routine 
     16   !!              -   !  2002-08  (G. Madec)  rn_cri and Free form, F90 
     17   !!              -   !  2004-10  (C. Ethe )  1D configuration 
     18   !!             2.0  !  2006-07  (S. Masson)  distributed restart using iom 
     19   !!             3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics: 
    2120   !!                              - tke penetration (wind steering) 
    2221   !!                              - suface condition for tke & mixing length 
    2322   !!                              - Langmuir cells 
    24    !!             3.0  !  2007-11  (G. Madec)  avtb_2d, nn_avtb_2d  
     23   !!              -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb 
     24   !!              -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters 
    2525   !!---------------------------------------------------------------------- 
    2626#if defined key_zdftke   ||   defined key_esopa 
     
    3939   USE phycst          ! physical constants 
    4040   USE zdfmxl          ! mixed layer 
     41   USE restart         ! only for lrst_oce 
    4142   USE lbclnk          ! ocean lateral boundary conditions (or mpp link) 
    4243   USE prtctl          ! Print control 
    4344   USE in_out_manager  ! I/O manager 
    44    USE iom 
    45    USE restart         ! only for lrst_oce 
     45   USE iom             ! I/O manager library 
    4646 
    4747   IMPLICIT NONE 
    4848   PRIVATE 
    4949 
    50    PUBLIC   zdf_tke        ! routine called in step module 
     50   PUBLIC   zdf_tke    ! routine called in step module 
    5151 
    5252   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag 
     
    5454   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   en                  !: now turbulent kinetic energy 
    5555# if defined key_vectopt_memory 
     56   !                                                                !!! key_vectopt_memory 
    5657   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   etmean              !: coefficient used for horizontal smoothing 
    5758   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   eumean, evmean      !: at t-, u- and v-points 
    5859# endif 
    59  
    60    !! * Namelist (namtke) 
    61    LOGICAL , PUBLIC ::   ln_rstke = .FALSE.         !: =T restart with tke from a run without tke with  
    62      !                                              !  a none zero initial value for en 
    63    INTEGER , PUBLIC ::   nitke = 50 ,            &  !: number of restart iterative loops 
    64       &                  nmxl  =  2 ,            &  !: = 0/1/2/3 flag for the type of mixing length used 
    65       &                  npdl  =  1 ,            &  !: = 0/1/2 flag on prandtl number on vert. eddy coeff. 
    66       &                  nave  =  1 ,            &  !: = 0/1 flag for horizontal average on avt, avmu, avmv 
    67       &                  navb  =  0                 !: = 0/1 flag for constant or profile background avt 
    68    REAL(wp), PUBLIC ::   ediff = 0.1_wp       ,  &  !: coeff. for vertical eddy coef.; avt=ediff*mxl*sqrt(e) 
    69       &                  ediss = 0.7_wp       ,  &  !: coef. of the Kolmogoroff dissipation  
    70       &                  ebb   = 3.75_wp      ,  &  !: coef. of the surface input of tke 
    71       &                  efave = 1._wp        ,  &  !: coef. for the tke vert. diff. coeff.; avtke=efave*avm 
    72       &                  emin  = 0.7071e-6_wp ,  &  !: minimum value of tke (m2/s2) 
    73       &                  emin0 = 1.e-4_wp     ,  &  !: surface minimum value of tke (m2/s2) 
    74       &                  ri_c  = 2._wp / 9._wp      !: critic Richardson number 
    75    !                                      !!! ** Namelist (namtke) ** 
    76    INTEGER  ::   nn_avtb_2d = 1            !: = 0/1 horizontal shape for avtb 
    77    INTEGER  ::   nn_etau    = 0            !: = 0/1/2 tke depth penetration 
    78    INTEGER  ::   nn_htau    = 0            !: = 0/1/2 type of tke profile of penetration 
    79    REAL(wp) ::   rn_lmin    = 0.4_wp       !: surface min value of mixing turbulent length scale 
    80    REAL(wp) ::   rn_efr     = 1.0_wp       !: fraction of TKE surface value which penetrates in the ocean 
    81    LOGICAL  ::   ln_lsfc    = .FALSE.      !: mixing length scale surface value as function of wind stress or not 
    82    LOGICAL  ::   ln_lc      = .FALSE.      !: Langmuir cells (LC) as a source term of TKE or not 
    83    REAL(wp) ::   rn_lc      = 0.15_wp      !: coef to compute vertical velocity of LC 
    84  
    85    REAL(wp), DIMENSION (jpi,jpj)  :: avtb_2d ! set in tke_init, for other modif than ice 
    86  
    87 # if defined key_c1d 
    88    !                                                                   ! 1D cfg only 
    89    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix,      &  ! dissipation and mixing turbulent lengh scales 
    90       &                                          e_pdl, e_ric          ! prandl and local Richardson numbers 
     60#if defined key_c1d 
     61   !                                                                !!! 1D cfg only 
     62   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales 
     63   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers 
    9164#endif 
     65 
     66   !                                       !!! ** Namelist  namtke  ** 
     67   LOGICAL  ::   ln_rstke = .FALSE.         ! =T restart with tke from a run without tke 
     68   LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not 
     69   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not 
     70   INTEGER  ::   nn_itke  = 50              ! number of restart iterative loops 
     71   INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3) 
     72   INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1) 
     73   INTEGER  ::   nn_ave   =  1              ! horizontal average or not on avt, avmu, avmv (=0/1) 
     74   INTEGER  ::   nn_avb   =  0              ! constant or profile background on avt (=0/1) 
     75   REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 
     76   REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation  
     77   REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke 
     78   REAL(wp) ::   rn_efave = 1._wp           ! coefficient for ave : ave=rn_efave*avm 
     79   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke (m2/s2) 
     80   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke (m2/s2) 
     81   REAL(wp) ::   rn_cri   = 2._wp / 9._wp   ! critic Richardson number 
     82   INTEGER  ::   nn_havtb = 1               ! horizontal shape or not for avtb (=0/1) 
     83   INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2) 
     84   INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1/2) 
     85   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length 
     86   REAL(wp) ::   rn_lmin  = 0.4_wp          ! interior min value of mixing length 
     87   REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean 
     88   REAL(wp) ::   rn_lc    = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells 
     89 
     90   REAL(wp), DIMENSION (jpi,jpj) ::   avtb_2d   ! set in tke_init, for other modif than ice 
    9291 
    9392   !! * Substitutions 
     
    9594#  include "vectopt_loop_substitute.h90" 
    9695   !!---------------------------------------------------------------------- 
    97    !!   OPA 9.0 , LOCEAN-IPSL (2006)  
     96   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)  
    9897   !! $Id$ 
    9998   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     
    111110      !! ** Method  :   The time evolution of the turbulent kinetic energy  
    112111      !!      (tke) is computed from a prognostic equation : 
    113       !!         d(en)/dt = eboost eav (d(u)/dz)**2       ! shear production 
    114       !!                  + d( efave eav d(en)/dz )/dz    ! diffusion of tke 
    115       !!                  + grav/rau0 pdl eav d(rau)/dz   ! stratif. destruc. 
    116       !!                  - ediss / emxl en**(2/3)        ! dissipation 
     112      !!         d(en)/dt = eboost eav (d(u)/dz)**2         ! shear production 
     113      !!                  + d( rn_efave eav d(en)/dz )/dz   ! diffusion of tke 
     114      !!                  + grav/rau0 pdl eav d(rau)/dz     ! stratif. destruc. 
     115      !!                  - rn_ediss / emxl en**(2/3)       ! dissipation 
    117116      !!      with the boundary conditions: 
    118       !!         surface: en = max( emin0,ebb sqrt(utau^2 + vtau^2) ) 
    119       !!         bottom : en = emin 
     117      !!         surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) ) 
     118      !!         bottom : en = rn_emin 
    120119      !!      -1- The dissipation and mixing turbulent lengh scales are computed 
    121       !!      from the usual diagnostic buoyancy length scale:   
    122       !!         mxl= 1/(sqrt(en)/N)  WHERE N is the brunt-vaisala frequency 
     120      !!         from the usual diagnostic buoyancy length scale:   
     121      !!         mxl= sqrt(2*en)/N  where N is the brunt-vaisala frequency 
     122      !!         with mxl = rn_lmin at the bottom minimum value of 0.4 
    123123      !!      Four cases :  
    124       !!         nmxl=0 : mxl bounded by the distance to surface and bottom. 
     124      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom. 
    125125      !!                  zmxld = zmxlm = mxl 
    126       !!         nmxl=1 : mxl bounded by the vertical scale factor. 
     126      !!         nn_mxl=1 : mxl bounded by the vertical scale factor. 
    127127      !!                  zmxld = zmxlm = mxl 
    128       !!         nmxl=2 : mxl bounded such that the vertical derivative of mxl  
     128      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl  
    129129      !!                  is less than 1 (|d/dz(xml)|<1).  
    130130      !!                  zmxld = zmxlm = mxl 
    131       !!         nmxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface  
     131      !!         nn_mxl=3 : lup = mxl bounded using |d/dz(xml)|<1 from the surface  
    132132      !!                        to the bottom 
    133133      !!                  ldown = mxl bounded using |d/dz(xml)|<1 from the bottom  
     
    140140      !!      solved. 
    141141      !!         Note that - the shear production is multiplied by eboost in order 
    142       !!      to set the critic richardson number to ri_c (namelist parameter) 
     142      !!      to set the critic richardson number to rn_cri (namelist parameter) 
    143143      !!                   - the destruction by stratification term is multiplied 
    144144      !!      by the Prandtl number (defined by an empirical funtion of the local  
    145       !!      Richardson number) if npdl=1 (namelist parameter) 
     145      !!      Richardson number) if nn_pdl=1 (namelist parameter) 
    146146      !!      coefficient (zesh2): 
    147147      !!      -3- Compute the now vertical eddy vicosity and diffusivity 
    148148      !!      coefficients from en (before the time stepping) and zmxlm: 
    149       !!              avm = max( avtb, ediff*zmxlm*en^1/2 ) 
    150       !!              avt = max( avmb, pdl*avm )  (pdl=1 if npdl=0) 
     149      !!              avm = max( avtb, rn_ediff*zmxlm*en^1/2 ) 
     150      !!              avt = max( avmb, pdl*avm )  (pdl=1 if nn_pdl=0) 
    151151      !!              eav = max( avmb, avm ) 
    152152      !!      avt and avm are horizontally averaged to avoid numerical insta- 
     
    159159      !!                update avt, avmu, avmv (before vertical eddy coef.) 
    160160      !! 
    161       !! References : Gaspar et al., jgr, 95, 1990, 
    162       !!              Blanke and Delecluse, jpo, 1991 
     161      !! References : Gaspar et al., JGR, 1990, 
     162      !!              Blanke and Delecluse, JPO, 1991 
     163      !!              Mellor and Blumberg, JPO 2004 
     164      !!              Axell, JGR, 2002 
    163165      !!---------------------------------------------------------------------- 
    164       USE oce     , zwd   => ua,  &  ! use ua as workspace 
    165          &          zmxlm => ta,  &  ! use ta as workspace 
    166          &          zmxld => sa      ! use sa as workspace 
    167       USE oce     , ztkelc => va     ! use va as workspace 
    168       ! 
    169       INTEGER, INTENT(in) ::   kt ! ocean time step 
    170       ! 
    171       INTEGER  ::   ji, jj, jk                  ! dummy loop arguments 
    172       REAL(wp) ::   zmlmin, zbbrau,          &  ! temporary scalars 
    173          &          zfact1, zfact2, zfact3,  &  ! 
    174          &          zrn2, zesurf,            &  ! 
    175          &          ztx2, zty2, zav,         &  ! 
    176          &          zcoef, zcof, zsh2,       &  ! 
    177          &          zdku, zdkv, zpdl, zri,   &  ! 
    178          &          zsqen, zesh2,            &  ! 
    179          &          zemxl, zemlm, zemlp 
     166      USE oce,     zwd    =>   ua   ! use ua as workspace 
     167      USE oce,     zmxlm  =>   va   ! use va as workspace 
     168      USE oce,     zmxld  =>   ta   ! use ta as workspace 
     169      USE oce,     ztkelc =>   sa   ! use sa as workspace 
     170      ! 
     171      INTEGER, INTENT(in) ::   kt   ! ocean time step 
     172      ! 
     173      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments 
     174      REAL(wp) ::   zbbrau, zrn2, zesurf            ! temporary scalars 
     175      REAL(wp) ::   zfact1, ztx2, zdku              !    -         - 
     176      REAL(wp) ::   zfact2, zty2, zdkv              !    -         - 
     177      REAL(wp) ::   zfact3, zcoef, zcof, zav        !    -         - 
     178      REAL(wp) ::   zsh2, zpdl, zri, zsqen, zesh2   !    -         - 
     179      REAL(wp) ::   zemxl, zemlm, zemlp             !    -         - 
     180      REAL(wp) ::   zraug, zus, zwlc, zind          !    -         - 
    180181      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace 
    181       REAL(wp) ::   zraug, zus, zwlc, zind     ! temporary scalar 
    182       REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   ! 2D workspace 
    183       REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    ! 2D workspace 
     182      REAL(wp), DIMENSION(jpi,jpj)     ::   zhtau   !  -      - 
     183      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      - 
    184184      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace 
    185185      !!-------------------------------------------------------------------- 
     
    188188 
    189189      !                                            ! Local constant initialization 
    190       zmlmin = 0.4 
    191       zbbrau =  .5 * ebb / rau0 
    192       zfact1 = -.5 * rdt * efave 
    193       zfact2 = 1.5 * rdt * ediss 
    194       zfact3 = 0.5 * rdt * ediss 
     190      zbbrau =  .5 * rn_ebb / rau0 
     191      zfact1 = -.5 * rdt * rn_efave 
     192      zfact2 = 1.5 * rdt * rn_ediss 
     193      zfact3 = 0.5 * rdt * rn_ediss 
    195194 
    196195      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     
    200199      ! Buoyancy length scale: l=sqrt(2*e/n**2) 
    201200      ! --------------------- 
    202       IF( ln_lsfc ) THEN      ! lsurf is function of wind stress : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g) 
     201      IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=2.e5*sqrt(utau^2 + vtau^2)/(rau0*g) 
     202!!gm  this should be useless 
    203203         zmxlm(:,:,1) = 0.e0 
     204!!gm end 
    204205         zraug = 0.5 * 2.e5 / ( rau0 * grav ) 
    205206         DO jj = 2, jpjm1 
     
    208209               ztx2 = utau(ji-1,jj  ) + utau(ji,jj) 
    209210               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    210                zmxlm(ji,jj,1  ) = zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 
    211                ! set the surface minimum value to rn_lmin 
    212                zmxlm(ji,jj,1  ) = MAX( zmxlm(ji,jj,1) , rn_lmin ) 
    213             END DO 
    214          END DO 
    215       ELSE                    ! surface set to the minimum value 
    216          zmxlm(:,:,1) = rn_lmin   
     211               zmxlm(ji,jj,1) = MAX(  rn_lmin0,  zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )  ) 
     212            END DO 
     213         END DO 
     214      ELSE                       ! surface set to the minimum value 
     215         zmxlm(:,:,1) = rn_lmin0 
    217216      ENDIF 
    218       zmxlm(:,:,jpk) = rn_lmin   ! bottom  set to the minimum value 
    219 !CDIR NOVERRCHK 
    220       DO jk = 2, jpkm1 
     217      zmxlm(:,:,jpk) = rn_lmin   ! bottom set to the interior minium value 
     218      ! 
     219!CDIR NOVERRCHK 
     220      DO jk = 2, jpkm1           ! interior value : l=sqrt(2*e/n**2) 
    221221!CDIR NOVERRCHK 
    222222         DO jj = 2, jpjm1 
     
    224224            DO ji = fs_2, fs_jpim1   ! vector opt. 
    225225               zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 
    226                zmxlm(ji,jj,jk) = MAX( SQRT( 2. * en(ji,jj,jk) / zrn2 ), zmlmin  ) 
     226               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  ) 
    227227            END DO 
    228228         END DO 
     
    232232      ! ------------------------------------- 
    233233      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value 
    234       zmxld(:,:,jpk) = rn_lmin           ! bottom  set to the minimum value 
    235  
    236       SELECT CASE ( nmxl ) 
    237  
     234      zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value 
     235 
     236      SELECT CASE ( nn_mxl ) 
     237      ! 
    238238      CASE ( 0 )           ! bounded by the distance to surface and bottom 
    239239         DO jk = 2, jpkm1 
     
    247247            END DO 
    248248         END DO 
    249  
     249         ! 
    250250      CASE ( 1 )           ! bounded by the vertical scale factor 
    251251         DO jk = 2, jpkm1 
     
    258258            END DO 
    259259         END DO 
    260  
     260         ! 
    261261      CASE ( 2 )           ! |dk[xml]| bounded by e3t : 
    262262         DO jk = 2, jpkm1         ! from the surface to the bottom : 
     
    276276            END DO 
    277277         END DO 
    278  
     278         ! 
    279279      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t : 
    280280         DO jk = 2, jpkm1         ! from the surface to the bottom : lup 
     
    305305            END DO 
    306306         END DO 
    307  
     307         ! 
    308308      END SELECT 
    309309 
    310310# if defined key_c1d 
    311       ! save mixing and dissipation turbulent length scales 
     311      ! c1d configuration : save mixing and dissipation turbulent length scales 
    312312      e_dis(:,:,:) = zmxld(:,:,:) 
    313313      e_mix(:,:,:) = zmxlm(:,:,:) 
     
    346346         ! 
    347347# if defined key_c1d 
    348          hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! save finite Langmuir Circulation depth 
     348         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth 
    349349# endif 
    350350         ! 
     
    382382            DO ji = fs_2, fs_jpim1   ! vector opt. 
    383383               zsqen = SQRT( en(ji,jj,jk) ) 
    384                zav   = ediff * zmxlm(ji,jj,jk) * zsqen 
     384               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen 
    385385               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 
    386386               zmxlm(ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk) 
     
    392392      ! 2. Surface boundary condition on tke and its eddy viscosity (zmxlm) 
    393393      ! ------------------------------------------------- 
    394       ! en(1)   = ebb sqrt(utau^2+vtau^2) / rau0  (min value emin0) 
     394      ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0) 
    395395      ! zmxlm(1) = avmb(1) and zmxlm(jpk) = 0. 
    396396!CDIR NOVERRCHK 
     
    401401            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj) 
    402402            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 ) 
    403             en (ji,jj,1) = MAX( zesurf, emin0 ) * tmask(ji,jj,1) 
    404             zav  =  ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) )* tmask(ji,jj,1) 
    405             zmxlm(ji,jj,1) = MAX( zav, avmb(1) ) * tmask(ji,jj,1) 
    406             avt  (ji,jj,1) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1) 
     403            en (ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1) 
     404            zav  =  rn_ediff * zmxlm(ji,jj,1) * SQRT( en(ji,jj,1) ) 
     405            zmxlm(ji,jj,1  ) = MAX( zav, avmb(1)                  ) * tmask(ji,jj,1) 
     406            avt  (ji,jj,1  ) = MAX( zav, avtb(1) * avtb_2d(ji,jj) ) * tmask(ji,jj,1) 
    407407            zmxlm(ji,jj,jpk) = 0.e0 
    408408         END DO 
     
    412412      ! ------------------------------- 
    413413      ! Resolution of a tridiagonal linear system by a "methode de chasse" 
    414       ! computation from level 2 to jpkm1  (e(1) already computed and 
    415       ! e(jpk)=0 ).  
    416  
    417       SELECT CASE ( npdl ) 
    418  
     414      ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).  
     415 
     416      SELECT CASE ( nn_pdl ) 
     417      ! 
    419418      CASE ( 0 )           ! No Prandtl number 
    420419         DO jk = 2, jpkm1 
    421420            DO jj = 2, jpjm1 
    422421               DO ji = fs_2, fs_jpim1   ! vector opt. 
    423                   ! zesh2 = eboost * (du/dz)^2 - N^2 
     422                  !                                                             ! shear prod. - stratif. destruction 
    424423                  zcoef = 0.5 / fse3w(ji,jj,jk) 
    425                   ! shear 
    426                   zdku = zcoef * (   ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   & 
    427                   &                - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  ) 
    428                   zdkv = zcoef * (   vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   & 
    429                   &                - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  ) 
    430                   ! coefficient (zesh2) 
    431                   zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk) 
    432  
    433                   ! Matrix 
     424                  zdku = zcoef * (  ub(ji-1, jj ,jk-1) + ub(ji,jj,jk-1)   &          ! shear 
     425                     &            - ub(ji-1, jj ,jk  ) - ub(ji,jj,jk  )  ) 
     426                  zdkv = zcoef * (  vb( ji ,jj-1,jk-1) + vb(ji,jj,jk-1)   & 
     427                     &            - vb( ji ,jj-1,jk  ) - vb(ji,jj,jk  )  ) 
     428                  zesh2 =  eboost * ( zdku*zdku + zdkv*zdkv ) - rn2(ji,jj,jk)        ! coefficient (zesh2) 
     429                  ! 
     430                  !                                                             ! Matrix 
    434431                  zcof = zfact1 * tmask(ji,jj,jk) 
    435                   ! lower diagonal 
     432                  !                                                                  ! lower diagonal 
    436433                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   & 
    437                                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) ) 
    438                   ! upper diagonal 
     434                     &                  / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) ) 
     435                  !                                                                  ! upper diagonal 
    439436                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   & 
    440                                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) ) 
    441                   ! diagonal 
     437                     &                  / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) ) 
     438                  !                                                                  ! diagonal 
    442439                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 
    443                   ! right hand side in en  
    444                   IF( .NOT. ln_lc ) THEN 
     440                  ! 
     441                  !                                                             ! right hand side in en  
     442                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells 
    445443                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   & 
    446444                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2 
    447                   ELSE 
     445                  ELSE                                                               ! Langmuir cell source term 
    448446                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en(ji,jj,jk)   & 
    449447                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2          & 
     
    453451            END DO 
    454452         END DO 
    455  
     453         ! 
    456454      CASE ( 1 )           ! Prandtl number 
    457455         DO jk = 2, jpkm1 
    458456            DO jj = 2, jpjm1 
    459457               DO ji = fs_2, fs_jpim1   ! vector opt. 
    460                   ! zesh2 =  eboost * (du/dz)^2 - pdl * N^2 
    461                   zcoef = 0.5 / fse3w(ji,jj,jk) 
    462                   ! shear 
    463                   zdku = zcoef * (   ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1) & 
    464                   &                - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )   ) 
    465                   zdkv = zcoef * (   vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1) & 
    466                   &                - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )   ) 
    467                   ! square of vertical shear 
    468                   zsh2 = zdku * zdku + zdkv * zdkv 
    469                   ! local Richardson number 
    470                   zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 ) 
     458                  !                                                             ! shear prod. - stratif. destruction 
     459                  zcoef = 0.5 / fse3w(ji,jj,jk)                                      
     460                  zdku = zcoef * (  ub(ji-1,jj  ,jk-1) + ub(ji,jj,jk-1)   &          ! shear 
     461                  &               - ub(ji-1,jj  ,jk  ) - ub(ji,jj,jk  )  ) 
     462                  zdkv = zcoef * (  vb(ji  ,jj-1,jk-1) + vb(ji,jj,jk-1)   & 
     463                  &               - vb(ji  ,jj-1,jk  ) - vb(ji,jj,jk  )  ) 
     464                  zsh2 = zdku * zdku + zdkv * zdkv                                   ! square of shear 
     465                  zri  = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + 1.e-20 )                ! local Richardson number 
    471466# if defined key_c1d 
    472                   ! save masked local Richardson number in zmxlm array 
    473                   e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) 
     467                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri  
    474468# endif 
    475                   ! Prandtl number 
    476                   zpdl = 1.0 
    477                   IF( zri >= 0.2 ) zpdl = 0.2 / zri 
     469                  zpdl = 1.0                                                         ! Prandtl number 
     470                  IF( zri >= 0.2 )   zpdl = 0.2 / zri 
    478471                  zpdl = MAX( 0.1, zpdl ) 
    479                   ! coefficient (esh2) 
    480                   zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk) 
    481  
    482                   ! Matrix 
     472                  zesh2 = eboost * zsh2 - zpdl * rn2(ji,jj,jk)                       ! coefficient (esh2) 
     473                  ! 
     474                  !                                                             ! Matrix 
    483475                  zcof = zfact1 * tmask(ji,jj,jk) 
    484                   ! lower diagonal 
     476                  !                                                                 ! lower diagonal 
    485477                  avmv(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk  ) + zmxlm(ji,jj,jk-1) )   & 
    486478                  &                     / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) ) 
    487                   ! upper diagonal 
     479                  !                                                                 ! upper diagonal 
    488480                  avmu(ji,jj,jk) = zcof * ( zmxlm(ji,jj,jk+1) + zmxlm(ji,jj,jk  ) )   & 
    489481                  &                     / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk) ) 
    490                   ! diagonal 
     482                  !                                                                 ! diagonal 
    491483                  zwd(ji,jj,jk) = 1. - avmv(ji,jj,jk) - avmu(ji,jj,jk) + zfact2 * zmxld(ji,jj,jk) 
    492                   ! right hand side in en  
    493                   IF( .NOT. ln_lc ) THEN 
     484                  ! 
     485                  !                                                             ! right hand side in en  
     486                  IF( .NOT. ln_lc ) THEN                                             ! No Langmuir cells 
    494487                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   & 
    495488                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2 
    496                   ELSE 
     489                  ELSE                                                               ! Langmuir cell source term 
    497490                     en(ji,jj,jk) = en(ji,jj,jk) + zfact3 * zmxld (ji,jj,jk) * en   (ji,jj,jk)   & 
    498491                        &                        +   rdt  * zmxlm (ji,jj,jk) * zesh2             & 
    499492                        &                        +   rdt  * ztkelc(ji,jj,jk) 
    500493                  ENDIF 
    501                   ! save masked Prandlt number in zmxld array 
    502                   zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk) 
    503                END DO 
    504             END DO 
    505          END DO 
    506  
     494                  zmxld(ji,jj,jk) = zpdl * tmask(ji,jj,jk)                      ! store masked Prandlt number in zmxld array 
     495               END DO 
     496            END DO 
     497         END DO 
     498         ! 
    507499      END SELECT 
    508500 
    509501# if defined key_c1d 
    510       !  save masked Prandlt number 
    511       e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1) 
     502      e_pdl(:,:,2:jpkm1) = zmxld(:,:,2:jpkm1)      ! c1d configuration : save masked Prandlt number 
    512503      e_pdl(:,:,      1) = e_pdl(:,:,      2) 
    513504      e_pdl(:,:,    jpk) = e_pdl(:,:,  jpkm1)       
     
    516507      ! 4. Matrix inversion from level 2 (tke prescribed at level 1) 
    517508      !!-------------------------------- 
    518       ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    519       DO jk = 3, jpkm1 
     509      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 
    520510         DO jj = 2, jpjm1 
    521511            DO ji = fs_2, fs_jpim1    ! vector opt. 
     
    524514         END DO 
    525515      END DO 
    526  
    527       ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    528       DO jj = 2, jpjm1 
     516      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 
    529517         DO ji = fs_2, fs_jpim1    ! vector opt. 
    530518            avmv(ji,jj,2) = en(ji,jj,2) - avmv(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke 
     
    538526         END DO 
    539527      END DO 
    540  
    541       ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    542       DO jj = 2, jpjm1 
     528      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 
    543529         DO ji = fs_2, fs_jpim1    ! vector opt. 
    544530            en(ji,jj,jpkm1) = avmv(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 
     
    552538         END DO 
    553539      END DO 
    554  
    555       ! Save the result in en and set minimum value of tke : emin 
    556       DO jk = 2, jpkm1 
     540      DO jk = 2, jpkm1                             ! set the minimum value of tke 
    557541         DO jj = 2, jpjm1 
    558542            DO ji = fs_2, fs_jpim1   ! vector opt. 
    559                en(ji,jj,jk) = MAX( en(ji,jj,jk), emin ) * tmask(ji,jj,jk) 
     543               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 
    560544            END DO 
    561545         END DO 
    562546      END DO 
    563547      
    564       ! Modify the value of TKE by an exponential assumption 
    565       ! en(ji,jj,jk)=en(ji,jj,jk)+en(ji,jj,1)*EXP(-fsdepw(ji,jj,jk)/ zhtau(ji,jj) ) 
    566  
    567       SELECT CASE( nn_htau )      ! Choice of H value 
    568       ! 
    569       CASE( 0 ) 
    570          DO jj = 2, jpjm1 
    571             DO ji = fs_2, fs_jpim1   ! vector opt. 
    572                zhtau(ji,jj) = 10. 
    573             END DO 
    574          END DO 
    575          ! 
    576       CASE( 1 ) 
    577          DO jj = 2, jpjm1 
    578             DO ji = fs_2, fs_jpim1   ! vector opt. 
    579                zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   ) 
    580             END DO 
    581          END DO 
    582          ! 
    583       CASE( 2 ) 
    584          DO jj = 2, jpjm1 
    585             DO ji = fs_2, fs_jpim1   ! vector opt. 
    586                zhtau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   ) 
    587             END DO 
    588          END DO 
    589          ! 
    590       END SELECT 
    591  
    592       IF( nn_etau == 1 ) THEN 
    593          DO jk = 2, jpkm1 
    594             DO jj = 2, jpjm1 
    595                DO ji = fs_2, fs_jpim1   ! vector opt. 
    596                   en(ji,jj,jk) = en(ji,jj,jk)   & 
     548      ! 5. Add extra TKE due to surface and internal wave breaking (nn_etau /= 0) 
     549      !!---------------------------------------------------------- 
     550      IF( nn_etau /= 0 ) THEN        ! extra tke : en = en + rn_efr * en(1) * exp( -z/zhtau ) 
     551         ! 
     552         SELECT CASE( nn_htau )           ! Choice of the depth of penetration 
     553         CASE( 0 )                                    ! constant depth penetration (here 10 meters) 
     554            DO jj = 2, jpjm1 
     555               DO ji = fs_2, fs_jpim1   ! vector opt. 
     556                  zhtau(ji,jj) = 10. 
     557               END DO 
     558            END DO 
     559         CASE( 1 )                                    ! meridional profile  1 
     560            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 40 m at high lat.) 
     561               DO ji = fs_2, fs_jpim1   ! vector opt. 
     562                  zhtau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   ) 
     563               END DO 
     564            END DO 
     565         CASE( 2 )                                    ! meridional profile 2 
     566            DO jj = 2, jpjm1                          ! ( 5m in the tropics to a maximum of 60 m at high lat.) 
     567               DO ji = fs_2, fs_jpim1   ! vector opt. 
     568                  zhtau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   ) 
     569               END DO 
     570            END DO 
     571         END SELECT 
     572         !  
     573         IF( nn_etau == 1 ) THEN          ! extra term throughout the water column 
     574            DO jk = 2, jpkm1 
     575               DO jj = 2, jpjm1 
     576                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     577                     en(ji,jj,jk) = en(ji,jj,jk)   & 
     578                        &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 
     579                        &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     580                  END DO 
     581               END DO 
     582            END DO 
     583         ELSEIF( nn_etau == 2 ) THEN      ! extra term only at the base of the mixed layer 
     584            DO jj = 2, jpjm1 
     585               DO ji = fs_2, fs_jpim1   ! vector opt. 
     586                  jk = nmln(ji,jj) 
     587                  en(ji,jj,jk) = en(ji,jj,jk)    & 
    597588                     &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 
    598                      &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
    599                END DO 
    600             END DO 
    601          END DO 
    602       ELSEIF( nn_etau == 2 ) THEN     !  only at the base of the mixed layer 
    603          DO jj = 2, jpjm1 
    604             DO ji = fs_2, fs_jpim1   ! vector opt. 
    605                jk = nmln(ji,jj) 
    606                en(ji,jj,jk) = en(ji,jj,jk)    & 
    607                   &         + rn_efr * en(ji,jj,1)*EXP( -fsdepw(ji,jj,jk) / zhtau(ji,jj) ) & 
    608                   &                  * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
    609             END DO 
    610          END DO 
     589                     &                   * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk) 
     590               END DO 
     591            END DO 
     592         ENDIF 
     593         ! 
    611594      ENDIF 
    612595 
    613596 
    614       ! Lateral boundary conditions on ( avt, en )   (sign unchanged) 
    615       CALL lbc_lnk( en , 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. ) 
     597      ! Lateral boundary conditions (sign unchanged) 
     598      CALL lbc_lnk( en, 'W', 1. )   ;   CALL lbc_lnk( avt, 'W', 1. ) 
     599 
    616600 
    617601      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    618602      ! IV.  Before vertical eddy vicosity and diffusivity coefficients 
    619603      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 
    620  
    621       SELECT CASE ( nave ) 
    622           
    623       CASE ( 0 )                ! no horizontal average 
    624  
    625          ! Vertical eddy viscosity 
    626  
    627          DO jk = 2, jpkm1                                 ! Horizontal slab 
     604      ! 
     605      SELECT CASE ( nn_ave ) 
     606      CASE ( 0 )                                     ! no horizontal average  
     607         DO jk = 2, jpkm1                                 ! only vertical eddy viscosity  
    628608            DO jj = 2, jpjm1 
    629609               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    635615            END DO 
    636616         END DO 
    637  
    638          ! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged) 
    639          CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. ) 
    640           
    641       CASE ( 1 )                ! horizontal average 
    642  
    643          !                                                ( 1/2  1/2 ) 
    644          ! Eddy viscosity: horizontal average: avmu = 1/4 ( 1    1   ) 
    645          !                      ( 1/2  1 1/2 )            ( 1/2  1/2 ) 
    646          !           avmv = 1/4 ( 1/2  1 1/2 )       
    647           
    648 !! caution vectopt_memory change the solution (last digit of the solver stat) 
     617         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions 
     618         ! 
     619      CASE ( 1 )                                     ! horizontal average                         ( 1/2  1/2 ) 
     620         !                                                ! Vertical eddy viscosity    avmu = 1/4 ( 1    1   ) 
     621         !                                                !                                       ( 1/2  1/2 ) 
     622         !                                                ! 
     623         !                                                !                                       ( 1/2  1   1/2 )             
     624         !                                                !                            avmv = 1/4 ( 1/2  1   1/2 )       
     625         DO jk = 2, jpkm1                                 
     626            DO jj = 2, jpjm1 
     627               DO ji = fs_2, fs_jpim1   ! vector opt. 
     628# if defined key_vectopt_memory 
     629                  !                                       ! caution : vectopt_memory change the solution 
     630                  !                                       !           (last digit of the solver stat) 
     631                  avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   & 
     632                     &              +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   & 
     633                     &                   +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk) 
     634 
     635                  avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   & 
     636                     &              +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   & 
     637                     &                   +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk) 
     638# else 
     639                  avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   & 
     640                     &           +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   & 
     641                     &                +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)   & 
     642                     &    / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   & 
     643                     &           +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   & 
     644                     &                +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  ) 
     645 
     646                  avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   & 
     647                     &           +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   & 
     648                     &                +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)   & 
     649                     &   /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   & 
     650                     &           +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   & 
     651                     &                +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  ) 
     652# endif 
     653               END DO 
     654            END DO 
     655         END DO 
     656         CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions  
     657         ! 
     658         !                                                ! Vertical eddy diffusivity             (1 2 1) 
     659         !                                                !                            avt = 1/16 (2 4 2) 
     660         !                                                !                                       (1 2 1) 
     661         DO jk = 2, jpkm1          
     662            DO jj = 2, jpjm1 
     663               DO ji = fs_2, fs_jpim1   ! vector opt. 
    649664#  if defined key_vectopt_memory 
    650          DO jk = 2, jpkm1                                 ! Horizontal slab 
    651             DO jj = 2, jpjm1 
    652                DO ji = fs_2, fs_jpim1   ! vector opt. 
    653                   avmu(ji,jj,jk) = (      avt(ji,jj  ,jk) + avt(ji+1,jj  ,jk)   & 
    654                   &                 +.5*( avt(ji,jj-1,jk) + avt(ji+1,jj-1,jk)   & 
    655                   &                      +avt(ji,jj+1,jk) + avt(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk) 
    656  
    657                   avmv(ji,jj,jk) = (      avt(ji  ,jj,jk) + avt(ji  ,jj+1,jk)   & 
    658                   &                 +.5*( avt(ji-1,jj,jk) + avt(ji-1,jj+1,jk)   & 
    659                   &                      +avt(ji+1,jj,jk) + avt(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk) 
    660                END DO 
    661             END DO 
    662          END DO 
     665                  avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)                          & 
     666                     &            + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk) 
    663667#  else 
    664          DO jk = 2, jpkm1                                 ! Horizontal slab 
    665             DO jj = 2, jpjm1 
    666                DO ji = fs_2, fs_jpim1   ! vector opt. 
    667                   avmu(ji,jj,jk) = (   avt  (ji,jj  ,jk) + avt  (ji+1,jj  ,jk)   & 
    668                   &              +.5*( avt  (ji,jj-1,jk) + avt  (ji+1,jj-1,jk)   & 
    669                   &                   +avt  (ji,jj+1,jk) + avt  (ji+1,jj+1,jk) ) ) * umask(ji,jj,jk)  & 
    670                   &       / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   & 
    671                   &              +.5*( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   & 
    672                   &                   +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  ) 
    673  
    674                   avmv(ji,jj,jk) = (   avt  (ji  ,jj,jk) + avt  (ji  ,jj+1,jk)   & 
    675                   &              +.5*( avt  (ji-1,jj,jk) + avt  (ji-1,jj+1,jk)   & 
    676                   &                   +avt  (ji+1,jj,jk) + avt  (ji+1,jj+1,jk) ) ) * vmask(ji,jj,jk)  & 
    677                   &      /  MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   & 
    678                   &              +.5*( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   & 
    679                   &                   +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  ) 
    680                END DO 
    681             END DO 
    682          END DO 
     668                  avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)                        & 
     669                     &            + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   & 
     670                     &  / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)                        & 
     671                     &            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  ) 
    683672#  endif 
    684  
    685          ! Lateral boundary conditions (avmu,avmv) (sign unchanged) 
    686          CALL lbc_lnk( avmu, 'U', 1. )   ;    CALL lbc_lnk( avmv, 'V', 1. ) 
    687  
    688          ! Vertical eddy diffusivity 
    689          ! ------------------------------ 
    690          !                                (1 2 1) 
    691          ! horizontal average  avt = 1/16 (2 4 2) 
    692          !                                (1 2 1) 
    693          DO jk = 2, jpkm1                                 ! Horizontal slab 
    694 #  if defined key_vectopt_memory 
    695             DO jj = 2, jpjm1 
    696                DO ji = fs_2, fs_jpim1   ! vector opt. 
    697                   avt(ji,jj,jk) = ( avmu(ji,jj,jk) + avmu(ji-1,jj  ,jk)    & 
    698                   &               + avmv(ji,jj,jk) + avmv(ji  ,jj-1,jk)  ) * etmean(ji,jj,jk) 
    699                END DO 
    700             END DO 
    701 #  else 
    702             DO jj = 2, jpjm1 
    703                DO ji = fs_2, fs_jpim1   ! vector opt. 
    704                   avt(ji,jj,jk) = ( avmu (ji,jj,jk) + avmu (ji-1,jj  ,jk)   & 
    705                   &               + avmv (ji,jj,jk) + avmv (ji  ,jj-1,jk)  ) * tmask(ji,jj,jk)   & 
    706                   &     / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   & 
    707                   &               + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk)  ) 
    708                END DO 
    709             END DO 
    710 #  endif 
    711          END DO 
    712  
     673               END DO 
     674            END DO 
     675         END DO 
     676         ! 
    713677      END SELECT 
    714  
    715       ! multiplied by the Prandtl number (npdl>1) 
    716       ! ---------------------------------------- 
    717       IF( npdl == 1 ) THEN 
    718          DO jk = 2, jpkm1                                 ! Horizontal slab 
     678      ! 
     679      IF( nn_pdl == 1 ) THEN                            ! Ponderation by the Prandtl number (nn_pdl=1) 
     680         DO jk = 2, jpkm1  
    719681            DO jj = 2, jpjm1 
    720682               DO ji = fs_2, fs_jpim1   ! vector opt. 
     
    725687         END DO 
    726688      ENDIF 
    727  
    728       ! Minimum value on the eddy viscosity 
    729       ! ---------------------------------------- 
    730       DO jk = 2, jpkm1                                 ! Horizontal slab 
    731          DO jj = 1, jpj 
    732             DO ji = 1, jpi 
    733                avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk) 
    734                avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk) 
    735             END DO 
    736          END DO 
    737       END DO 
    738  
    739       ! Lateral boundary conditions on avt  (sign unchanged) 
    740       ! ------------------------------===== 
    741       CALL lbc_lnk( avt, 'W', 1. ) 
    742  
    743       ! write en in restart file 
    744       ! ------------------------ 
    745       IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' ) 
     689      ! 
     690      DO jk = 2, jpkm1                                  ! Minimum value on the eddy viscosity 
     691         avmu(:,:,jk) = MAX( avmu(:,:,jk), avmb(jk) ) * umask(:,:,jk) 
     692         avmv(:,:,jk) = MAX( avmv(:,:,jk), avmb(jk) ) * vmask(:,:,jk) 
     693      END DO 
     694 
     695      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged) 
     696 
     697      IF( lrst_oce )   CALL tke_rst( kt, 'WRITE' )      ! write en in restart file 
    746698 
    747699      IF(ln_ctl) THEN 
    748          CALL prt_ctl(tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt , clinfo2=' t: ', ovlap=1, kdim=jpk) 
    749          CALL prt_ctl(tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask, & 
    750             &         tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk) 
     700         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk) 
     701         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                  & 
     702            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk ) 
    751703      ENDIF 
    752  
     704      ! 
    753705   END SUBROUTINE zdf_tke 
    754706 
     
    773725      ! 
    774726# if defined key_vectopt_memory 
    775       ! caution vectopt_memory change the solution (last digit of the solver stat) 
    776727      INTEGER ::   ji, jj, jk   ! dummy loop indices 
    777728# else 
     
    779730# endif 
    780731      !! 
    781       NAMELIST/namtke/ ln_rstke, ediff, ediss, ebb, efave, emin, emin0,   & 
    782                        ri_c, nitke, nmxl, npdl, nave, navb,    & 
    783                        ln_lsfc, rn_lmin, nn_avtb_2d, nn_etau, nn_htau, rn_efr, & 
    784                        ln_lc , rn_lc  
     732      NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb  , rn_efave, rn_emin,   & 
     733         &             rn_emin0, rn_cri  , nn_itke , nn_mxl  , nn_pdl  , nn_ave ,   & 
     734         &             nn_avb  , ln_mxl0 , rn_lmin , rn_lmin0, nn_havtb, nn_etau,  & 
     735         &             nn_htau , rn_efr  , ln_lc  , rn_lc  
    785736      !!---------------------------------------------------------------------- 
    786737 
     
    791742 
    792743      ! Compute boost associated with the Richardson critic 
    793       !     (control values: ri_c = 0.3   ==> eboost=1.25 for npdl=1 or 2) 
    794       !     (                ri_c = 0.222 ==> eboost=1.                  ) 
    795       eboost = ri_c * ( 2. + ediss / ediff ) / 2. 
     744      !     (control values: rn_cri = 0.3   ==> eboost=1.25 for nn_pdl=1) 
     745      !     (                rn_cri = 0.222 ==> eboost=1.               ) 
     746      eboost = rn_cri * ( 2. + rn_ediss / rn_ediff ) / 2. 
     747 
    796748 
    797749 
     
    804756         WRITE(numout,*) '~~~~~~~~~~~~' 
    805757         WRITE(numout,*) '          Namelist namtke : set tke mixing parameters' 
    806          WRITE(numout,*) '             restart with tke from no tke ln_rstke = ', ln_rstke 
    807          WRITE(numout,*) '             coef. to compute avt           ediff  = ', ediff 
    808          WRITE(numout,*) '             Kolmogoroff dissipation coef.  ediss  = ', ediss 
    809          WRITE(numout,*) '             tke surface input coef.        ebb    = ', ebb 
    810          WRITE(numout,*) '             tke diffusion coef.            efave  = ', efave 
    811          WRITE(numout,*) '             minimum value of tke           emin   = ', emin 
    812          WRITE(numout,*) '             surface minimum value of tke   emin0  = ', emin0 
    813          WRITE(numout,*) '             number of restart iter loops   nitke  = ', nitke 
    814          WRITE(numout,*) '             mixing length type             nmxl   = ', nmxl 
    815          WRITE(numout,*) '             prandl number flag             npdl   = ', npdl 
    816          WRITE(numout,*) '             horizontal average flag        nave   = ', nave 
    817          WRITE(numout,*) '             critic Richardson nb           ri_c   = ', ri_c 
    818          WRITE(numout,*) '                and its associated coeff.   eboost = ', eboost 
    819          WRITE(numout,*) '             constant background or profile navb   = ', navb 
    820          WRITE(numout,*) '             flag for compu.of bls as fun. of wind     ln_lsfc    = ', ln_lsfc 
    821          WRITE(numout,*) '             buoyancy lenght scale surface min value   rn_lmin    = ', rn_lmin 
    822          WRITE(numout,*) '             horizontal variation for avtb             nn_avtb_2d = ', nn_avtb_2d 
    823          WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau    = ', nn_etau 
    824          WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau    = ', nn_htau 
    825          WRITE(numout,*) '             % of emin which pene. the thermocline     rn_efr     = ', rn_efr 
    826          WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc      = ', ln_lc 
    827          WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc      = ', rn_lc 
     758         WRITE(numout,*) '             restart with tke from no tke              ln_rstke = ', ln_rstke 
     759         WRITE(numout,*) '             coef. to compute avt                      rn_ediff = ', rn_ediff 
     760         WRITE(numout,*) '             Kolmogoroff dissipation coef.             rn_ediss = ', rn_ediss 
     761         WRITE(numout,*) '             tke surface input coef.                   rn_ebb   = ', rn_ebb 
     762         WRITE(numout,*) '             tke diffusion coef.                       rn_efave = ', rn_efave 
     763         WRITE(numout,*) '             minimum value of tke                      rn_emin  = ', rn_emin 
     764         WRITE(numout,*) '             surface minimum value of tke              rn_emin0 = ', rn_emin0 
     765         WRITE(numout,*) '             number of restart iter loops              nn_itke  = ', nn_itke 
     766         WRITE(numout,*) '             mixing length type                        nn_mxl   = ', nn_mxl 
     767         WRITE(numout,*) '             prandl number flag                        nn_pdl   = ', nn_pdl 
     768         WRITE(numout,*) '             horizontal average flag                   nn_ave   = ', nn_ave 
     769         WRITE(numout,*) '             critic Richardson nb                      rn_cri   = ', rn_cri 
     770         WRITE(numout,*) '                and its associated coeff.              eboost   = ', eboost 
     771         WRITE(numout,*) '             constant background or profile            nn_avb   = ', nn_avb 
     772         WRITE(numout,*) '             surface mixing length = F(stress) or not  ln_mxl0  = ', ln_mxl0 
     773         WRITE(numout,*) '             surface  mixing length minimum value      rn_lmin0 = ', rn_lmin0 
     774         WRITE(numout,*) '             interior mixing length minimum value      rn_lmin0 = ', rn_lmin0 
     775         WRITE(numout,*) '             horizontal variation for avtb             nn_havtb = ', nn_havtb 
     776         WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau  = ', nn_etau 
     777         WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau  = ', nn_htau 
     778         WRITE(numout,*) '             % of rn_emin0 which pene. the thermocline rn_efr   = ', rn_efr 
     779         WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc    = ', ln_lc 
     780         WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc    = ', rn_lc 
    828781         WRITE(numout,*) 
    829782      ENDIF 
    830783 
    831       ! Check nmxl and npdl values 
    832       IF( nmxl < 0 .OR. nmxl > 3 ) CALL ctl_stop( '          bad flag: nmxl is < 0 or > 3 ' ) 
    833       IF( npdl < 0 .OR. npdl > 1 ) CALL ctl_stop( '          bad flag: npdl is < 0 or > 1 ' ) 
    834       IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is < 0 or > 3 ' ) 
    835       IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be > 0.15 and < 0.2 ' ) 
     784      ! Check of some namelist values 
     785      IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' ) 
     786      IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' ) 
     787      IF( nn_ave  < 0    .OR. nn_ave  > 1   )   CALL ctl_stop( 'bad flag: nn_ave is  0 or 1    ' ) 
     788      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 
     789      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 ' ) 
    836790 
    837791      IF( nn_etau == 2 )   CALL zdf_mxl( nit000 )      ! Initialization of nmln  
     
    839793      ! Horizontal average : initialization of weighting arrays  
    840794      ! ------------------- 
    841        
    842       SELECT CASE ( nave ) 
    843  
     795      ! 
     796      SELECT CASE ( nn_ave ) 
     797      ! Notice: mean arrays have not to by defined at local domain boundaries due to the vertical nature of TKE 
     798      ! 
    844799      CASE ( 0 )                ! no horizontal average 
    845800         IF(lwp) WRITE(numout,*) '          no horizontal average on avt, avmu, avmv' 
     
    854809         eumean(:,:,:) = 0.e0 
    855810         evmean(:,:,:) = 0.e0 
    856           
     811         ! 
    857812         DO jk = 1, jpkm1 
    858813            DO jj = 2, jpjm1 
    859814               DO ji = fs_2, fs_jpim1   ! vector opt. 
    860                   etmean(ji,jj,jk) = tmask(ji,jj,jk)                     & 
    861                   &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   & 
    862                   &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  ) 
    863                    
    864                   eumean(ji,jj,jk) = umask(ji,jj,jk)                     & 
    865                   &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk)  ) 
    866  
    867                   evmean(ji,jj,jk) = vmask(ji,jj,jk)                     & 
    868                   &  / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk)  ) 
     815                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,  umask(ji,jj,jk) + umask(ji-1,jj  ,jk)   & 
     816                  &                                            + vmask(ji,jj,jk) + vmask(ji  ,jj-1,jk) ) 
     817                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji+1,jj  ,jk) ) 
     818                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,  tmask(ji,jj,jk) + tmask(ji  ,jj+1,jk) ) 
    869819               END DO 
    870820            END DO 
    871821         END DO 
    872822# endif 
    873  
     823         ! 
    874824      CASE ( 1 )                ! horizontal average  
    875825         IF(lwp) WRITE(numout,*) '          horizontal average on avt, avmu, avmv' 
     
    883833         eumean(:,:,:) = 0.e0 
    884834         evmean(:,:,:) = 0.e0 
    885           
     835         ! 
    886836         DO jk = 1, jpkm1 
    887837            DO jj = 2, jpjm1 
    888838               DO ji = fs_2, fs_jpim1   ! vector opt. 
    889                   etmean(ji,jj,jk) = tmask(ji,jj,jk)                     & 
    890                   &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   & 
    891                   &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  ) 
    892                    
    893                   eumean(ji,jj,jk) = umask(ji,jj,jk)                        & 
    894                   &  / MAX( 1.,   tmask(ji,jj  ,jk) + tmask(ji+1,jj  ,jk)   & 
    895                   &       +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk)   & 
    896                   &              +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) )  ) 
    897  
    898                   evmean(ji,jj,jk) = vmask(ji,jj,jk)                        & 
    899                   &  / MAX( 1.,   tmask(ji  ,jj,jk) + tmask(ji  ,jj+1,jk)   & 
    900                   &       +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk)   & 
    901                   &              +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) )  ) 
     839                  etmean(ji,jj,jk) = tmask(ji,jj,jk) / MAX( 1.,   umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk)   & 
     840                  &                                       +       vmask(ji  ,jj-1,jk) + vmask(ji  ,jj  ,jk) ) 
     841                  eumean(ji,jj,jk) = umask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji+1,jj  ,jk)   & 
     842                  &                                       +.5 * ( tmask(ji  ,jj-1,jk) + tmask(ji+1,jj-1,jk)   & 
     843                  &                                             + tmask(ji  ,jj+1,jk) + tmask(ji+1,jj+1,jk) )  ) 
     844                  evmean(ji,jj,jk) = vmask(ji,jj,jk) / MAX( 1.,   tmask(ji  ,jj  ,jk) + tmask(ji  ,jj+1,jk)   & 
     845                  &                                       +.5 * ( tmask(ji-1,jj  ,jk) + tmask(ji-1,jj+1,jk)   & 
     846                  &                                             + tmask(ji+1,jj  ,jk) + tmask(ji+1,jj+1,jk) )  ) 
    902847               END DO 
    903848            END DO 
    904849         END DO 
    905850# endif 
    906  
    907       CASE DEFAULT 
    908          WRITE(ctmp1,*) '          bad flag value for nave = ', nave 
    909          CALL ctl_stop( ctmp1 ) 
    910  
     851         ! 
    911852      END SELECT 
    912853 
     
    914855      ! Background eddy viscosity and diffusivity profil 
    915856      ! ------------------------------------------------ 
    916       IF( navb == 0 ) THEN      ! Define avmb, avtb from namelist parameter 
     857      IF( nn_avb == 0 ) THEN      ! Define avmb, avtb from namelist parameter 
    917858         avmb(:) = avm0 
    918859         avtb(:) = avt0 
    919860      ELSE                      ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)  
    920861         avmb(:) = avm0 
    921 !!bug  this is not valide neither in scoord 
    922          IF(ln_sco .AND. lwp) WRITE(numout,cform_war) 
    923          IF(ln_sco .AND. lwp) WRITE(numout,*) '          avtb profile not valid in sco' 
    924  
    925862         avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s 
     863         IF(ln_sco .AND. lwp)   CALL ctl_warn( '          avtb profile not valid in sco' ) 
    926864      ENDIF 
    927  
     865      ! 
    928866      !                         ! 2D shape of the avtb 
    929867      avtb_2d(:,:) = 1.e0             ! uniform  
    930868      ! 
    931       IF( nn_avtb_2d == 1 ) THEN      ! decrease avtb in the equatorial band 
     869      IF( nn_havtb == 1 ) THEN      ! decrease avtb in the equatorial band 
    932870           !  -15S -5S : linear decrease from avt0 to avt0/10. 
    933871           !  -5S  +5N : cst value avt0/10. 
     
    938876      ENDIF 
    939877 
     878!!gm  the increase is only required when using cen2 advection scheme  
     879!!gm  for the equatorial upwelling.       ===>  use of TVD in ORCA  so this have to be suppressed 
    940880      !   Increase the background in the surface layers 
    941 !!gm  the increase is only required when using cen2 advection scheme  
    942 !!gm  for the equatorial upwelling. 
    943881      avmb(1) = 10.  * avmb(1)      ;      avtb(1) = 10.  * avtb(1) 
    944882      avmb(2) = 10.  * avmb(2)      ;      avtb(2) = 10.  * avtb(2) 
    945883      avmb(3) =  5.  * avmb(3)      ;      avtb(3) =  5.  * avtb(3) 
    946884      avmb(4) =  2.5 * avmb(4)      ;      avtb(4) =  2.5 * avtb(4) 
     885!!gm end 
    947886 
    948887 
     
    967906     !!                   ***  ROUTINE ts_rst  *** 
    968907     !!                      
    969      !! ** Purpose : Read or write filtered free surface arrays in restart file 
     908     !! ** Purpose :   Read or write TKE file (en) in restart file 
    970909     !! 
    971      !! ** Method  :  
    972      !! 
     910     !! ** Method  :   use of IOM library 
     911     !!                if the restart does not contain TKE, en is either  
     912     !!                set to rn_emin or recomputed (nn_itke/=0) 
    973913     !!---------------------------------------------------------------------- 
    974914     INTEGER         , INTENT(in) ::   kt         ! ocean time-step 
     
    986926                 &                       WRITE(numout,*) ' ===>>>> : previous run without tke scheme' 
    987927              IF( lwp .AND. ln_rstke )   WRITE(numout,*) ' ===>>>> : We do not use en from the restart file' 
    988               IF( lwp                )   WRITE(numout,*) ' ===>>>> : en set by iterative loop' 
    989               IF( lwp                )   WRITE(numout,*) ' =======             =========' 
    990               en (:,:,:) = emin * tmask(:,:,:) 
    991               DO jit = 2, nitke+1 
     928              IF( lwp                )   WRITE(numout,*) ' ===>>>> : en is set by iterative loop' 
     929              en (:,:,:) = rn_emin * tmask(:,:,:) 
     930              DO jit = 2, nn_itke + 1 
    992931                 CALL zdf_tke( jit ) 
    993932              END DO 
    994933           ENDIF 
    995934        ELSE 
    996            en(:,:,:) = emin * tmask(:,:,:)      ! no restart: en set to emin 
     935           en(:,:,:) = rn_emin * tmask(:,:,:)      ! no restart: en set to its minimum value 
    997936        ENDIF 
    998937     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 
Note: See TracChangeset for help on using the changeset viewer.