New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdftke.F90 in branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago

ORCA2_LIM_PISCES hybrid version update
  • Property svn:keywords set to Id
File size: 50.0 KB
Line 
1MODULE zdftke
2   !!======================================================================
3   !!                       ***  MODULE  zdftke  ***
4   !! Ocean physics:  vertical mixing coefficient computed from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
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 tke_init routine
16   !!             -   !  2004-10  (C. Ethe )  1D configuration
17   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
18   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
19   !!                 !           - tke penetration (wind steering)
20   !!                 !           - suface condition for tke & mixing length
21   !!                 !           - Langmuir cells
22   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
23   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
24   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
25   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
26   !!                 !                                + cleaning of the parameters + bugs correction
27   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
29   !!----------------------------------------------------------------------
30#if defined key_zdftke
31   !!----------------------------------------------------------------------
32   !!   'key_zdftke'                                   TKE vertical physics
33   !!----------------------------------------------------------------------
34   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
35   !!   tke_tke       : tke time stepping: update tke at now time step (en)
36   !!   tke_avn       : compute mixing length scale and deduce avm and avt
37   !!   zdf_tke_init  : initialization, namelist read, and parameters control
38   !!   tke_rst       : read/write tke restart in ocean restart file
39   !!----------------------------------------------------------------------
40   USE oce            ! ocean: dynamics and active tracers variables
41   USE phycst         ! physical constants
42   USE dom_oce        ! domain: ocean
43   USE domvvl         ! domain: variable volume layer
44   USE sbc_oce        ! surface boundary condition: ocean
45   USE zdf_oce        ! vertical physics: ocean variables
46   USE zdfmxl         ! vertical physics: mixed layer
47   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
48   USE prtctl         ! Print control
49   USE in_out_manager ! I/O manager
50   USE iom            ! I/O manager library
51   USE lib_mpp        ! MPP library
52   USE wrk_nemo       ! work arrays
53   USE timing         ! Timing
54   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
55#if defined key_agrif
56   USE agrif_opa_interp
57   USE agrif_opa_update
58#endif
59
60   IMPLICIT NONE
61   PRIVATE
62
63   PUBLIC   zdf_tke        ! routine called in step module
64   PUBLIC   zdf_tke_init   ! routine called in opa module
65   PUBLIC   tke_rst        ! routine called in step module
66
67   LOGICAL , PUBLIC, PARAMETER ::   lk_zdftke = .TRUE.  !: TKE vertical mixing flag
68
69   !                      !!** Namelist  namzdf_tke  **
70   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
71   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
72   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
73   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
74   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
75   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
76   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
77   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
78   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
79   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
80   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
81   INTEGER  ::   nn_htau   ! type of tke profile of penetration (=0/1)
82   REAL(wp) ::   rn_efr    ! fraction of TKE surface value which penetrates in the ocean
83   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
84   REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
85
86   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
87   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
88   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
89   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
90
91   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
92   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
93   REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr          ! now mixing lenght of dissipation
94#if defined key_c1d
95   !                                                                        !!** 1D cfg only  **   ('key_c1d')
96   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_dis, e_mix   !: dissipation and mixing turbulent lengh scales
97   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   e_pdl, e_ric   !: prandl and local Richardson numbers
98#endif
99
100   !! * Substitutions
101#  include "vectopt_loop_substitute.h90"
102   !!----------------------------------------------------------------------
103   !! NEMO/OPA 3.7 , NEMO Consortium (2015)
104   !! $Id$
105   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
106   !!----------------------------------------------------------------------
107CONTAINS
108
109   INTEGER FUNCTION zdf_tke_alloc()
110      !!----------------------------------------------------------------------
111      !!                ***  FUNCTION zdf_tke_alloc  ***
112      !!----------------------------------------------------------------------
113      ALLOCATE(                                                                    &
114#if defined key_c1d
115         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
116         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
117#endif
118         &      htau  (jpi,jpj)    , dissl(jpi,jpj,jpk) ,     & 
119         &      apdlr(jpi,jpj,jpk) ,                                           STAT= zdf_tke_alloc      )
120         !
121      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
122      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
123      !
124   END FUNCTION zdf_tke_alloc
125
126
127   SUBROUTINE zdf_tke( kt )
128      !!----------------------------------------------------------------------
129      !!                   ***  ROUTINE zdf_tke  ***
130      !!
131      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
132      !!              coefficients using a turbulent closure scheme (TKE).
133      !!
134      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
135      !!              is computed from a prognostic equation :
136      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
137      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
138      !!                  + avt N^2                      ! stratif. destruc.
139      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
140      !!      with the boundary conditions:
141      !!         surface: en = max( rn_emin0, rn_ebb * taum )
142      !!         bottom : en = rn_emin
143      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
144      !!
145      !!        The now Turbulent kinetic energy is computed using the following
146      !!      time stepping: implicit for vertical diffusion term, linearized semi
147      !!      implicit for kolmogoroff dissipation term, and explicit forward for
148      !!      both buoyancy and shear production terms. Therefore a tridiagonal
149      !!      linear system is solved. Note that buoyancy and shear terms are
150      !!      discretized in a energy conserving form (Bruchard 2002).
151      !!
152      !!        The dissipative and mixing length scale are computed from en and
153      !!      the stratification (see tke_avn)
154      !!
155      !!        The now vertical eddy vicosity and diffusivity coefficients are
156      !!      given by:
157      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
158      !!              avt = max( avmb, pdl * avm                 ) 
159      !!              eav = max( avmb, avm )
160      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
161      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
162      !!
163      !! ** Action  :   compute en (now turbulent kinetic energy)
164      !!                update avt, avmu, avmv (before vertical eddy coef.)
165      !!
166      !! References : Gaspar et al., JGR, 1990,
167      !!              Blanke and Delecluse, JPO, 1991
168      !!              Mellor and Blumberg, JPO 2004
169      !!              Axell, JGR, 2002
170      !!              Bruchard OM 2002
171      !!----------------------------------------------------------------------
172      INTEGER, INTENT(in) ::   kt   ! ocean time step
173      !!----------------------------------------------------------------------
174      !
175#if defined key_agrif 
176      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
177      IF( .NOT.Agrif_Root() )   CALL Agrif_Tke
178#endif
179      !
180      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
181!$OMP PARALLEL WORKSHARE
182         avt (:,:,:) = avt_k (:,:,:) 
183         avm (:,:,:) = avm_k (:,:,:) 
184         avmu(:,:,:) = avmu_k(:,:,:) 
185         avmv(:,:,:) = avmv_k(:,:,:) 
186!$OMP END PARALLEL WORKSHARE
187      ENDIF 
188      !
189      CALL tke_tke      ! now tke (en)
190      !
191      CALL tke_avn      ! now avt, avm, avmu, avmv
192      !
193!$OMP PARALLEL WORKSHARE
194      avt_k (:,:,:) = avt (:,:,:) 
195      avm_k (:,:,:) = avm (:,:,:) 
196      avmu_k(:,:,:) = avmu(:,:,:) 
197      avmv_k(:,:,:) = avmv(:,:,:) 
198!$OMP END PARALLEL WORKSHARE
199      !
200#if defined key_agrif
201      ! Update child grid f => parent grid
202      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
203#endif     
204     !
205  END SUBROUTINE zdf_tke
206
207
208   SUBROUTINE tke_tke
209      !!----------------------------------------------------------------------
210      !!                   ***  ROUTINE tke_tke  ***
211      !!
212      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
213      !!
214      !! ** Method  : - TKE surface boundary condition
215      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
216      !!              - source term due to shear (saved in avmu, avmv arrays)
217      !!              - Now TKE : resolution of the TKE equation by inverting
218      !!                a tridiagonal linear system by a "methode de chasse"
219      !!              - increase TKE due to surface and internal wave breaking
220      !!
221      !! ** Action  : - en : now turbulent kinetic energy)
222      !!              - avmu, avmv : production of TKE by shear at u and v-points
223      !!                (= Kz dz[Ub] * dz[Un] )
224      !! ---------------------------------------------------------------------
225      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
226!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
227!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
228      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
229      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
230      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
231      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
232      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
233      REAL(wp) ::   ztau  , zdif                    !    -         -
234      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
235      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
236!!bfr      REAL(wp) ::   zebot                           !    -         -
237      INTEGER , POINTER, DIMENSION(:,:  ) ::   imlc
238      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc
239      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv
240      REAL(wp)                            ::   zri  !   local Richardson number
241      !!--------------------------------------------------------------------
242      !
243      IF( nn_timing == 1 )  CALL timing_start('tke_tke')
244      !
245      CALL wrk_alloc( jpi,jpj,       imlc )    ! integer
246      CALL wrk_alloc( jpi,jpj,       zhlc ) 
247      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
248      !
249      zbbrau = rn_ebb / rau0       ! Local constant initialisation
250      zfact1 = -.5_wp * rdt 
251      zfact2 = 1.5_wp * rdt * rn_ediss
252      zfact3 = 0.5_wp       * rn_ediss
253      !
254      !
255      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
256      !                     !  Surface boundary condition on tke
257      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
258      IF ( ln_isfcav ) THEN
259!$OMP PARALLEL DO schedule(static) private(jj, ji)
260         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
261            DO ji = fs_2, fs_jpim1   ! vector opt.
262               en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
263            END DO
264         END DO
265      END IF
266!$OMP PARALLEL DO schedule(static) private(jj, ji)
267      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
268         DO ji = fs_2, fs_jpim1   ! vector opt.
269            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
270         END DO
271      END DO
272     
273!!bfr   - start commented area
274      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
275      !                     !  Bottom boundary condition on tke
276      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
277      !
278      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
279      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
280      ! The condition is coded here for completion but commented out until there is proof that the
281      ! computational cost is justified
282      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
283      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
284!!    DO jj = 2, jpjm1
285!!       DO ji = fs_2, fs_jpim1   ! vector opt.
286!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
287!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
288!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
289!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
290!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
291!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
292!!       END DO
293!!    END DO
294!!bfr   - end commented area
295      !
296      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
297      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
298         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
299         !
300         !                        !* total energy produce by LC : cumulative sum over jk
301!$OMP PARALLEL
302!$OMP DO schedule(static) private(jj, ji)
303         DO jj =1, jpj
304            DO ji=1, jpi
305               zpelc(ji,jj,1) =  MAX( rn2b(ji,jj,1), 0._wp ) * gdepw_n(ji,jj,1) * e3w_n(ji,jj,1)
306            END DO
307         END DO
308         DO jk = 2, jpk
309!$OMP DO schedule(static) private(jj, ji)
310            DO jj =1, jpj
311               DO ji=1, jpi
312                  zpelc(ji,jj,jk)  = zpelc(ji,jj,jk-1) + MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
313               END DO
314            END DO
315         END DO
316         !                        !* finite Langmuir Circulation depth
317         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
318!$OMP WORKSHARE
319         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
320!$OMP END WORKSHARE
321         DO jk = jpkm1, 2, -1
322!$OMP DO schedule(static) private(jj, ji, zus)
323            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
324               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
325                  zus  = zcof * taum(ji,jj)
326                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
327               END DO
328            END DO
329         END DO
330         !                               ! finite LC depth
331!$OMP DO schedule(static) private(jj, ji)
332         DO jj = 1, jpj 
333            DO ji = 1, jpi
334               zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))
335            END DO
336         END DO
337         zcof = 0.016 / SQRT( zrhoa * zcdrag )
338!$OMP DO schedule(static) private(jk, jj, ji, zus, zind, zwlc)
339         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
340            DO jj = 2, jpjm1
341               DO ji = fs_2, fs_jpim1   ! vector opt.
342                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
343                  !                                           ! vertical velocity due to LC
344                  zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) )
345                  zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )
346                  !                                           ! TKE Langmuir circulation source term
347                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
348                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
349               END DO
350            END DO
351         END DO
352!$OMP END DO NOWAIT
353!$OMP END PARALLEL
354         !
355      ENDIF
356      !
357      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
358      !                     !  Now Turbulent kinetic energy (output in en)
359      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
360      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
361      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
362      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
363      !
364!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
365      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
366         DO jj = 1, jpjm1
367            DO ji = 1, fs_jpim1   ! vector opt.
368               z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   &
369                  &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
370                  &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
371                  &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
372               z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   &
373                  &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
374                  &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
375                  &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
376            END DO
377         END DO
378      END DO
379      !
380      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: compute apdlr
381         ! Note that zesh2 is also computed in the next loop.
382         ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops
383!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zesh2, zri)
384         DO jk = 2, jpkm1
385            DO jj = 2, jpjm1
386               DO ji = fs_2, fs_jpim1   ! vector opt.
387                  !                                          ! shear prod. at w-point weightened by mask
388                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
389                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
390                  !                                          ! local Richardson number
391                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear )
392                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
393                 
394               END DO
395            END DO
396         END DO
397         !
398      ENDIF
399      !         
400!$OMP PARALLEL
401!$OMP DO schedule(static) private(jk, jj, ji, zcof, zzd_up, zzd_lw, zesh2)
402      DO jk = 2, jpkm1           !* Matrix and right hand side in en
403         DO jj = 2, jpjm1
404            DO ji = fs_2, fs_jpim1   ! vector opt.
405               zcof   = zfact1 * tmask(ji,jj,jk)
406# if defined key_zdftmx_new
407               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
408               zzd_up = zcof * MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp )   &  ! upper diagonal
409                  &          / (  e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  )  )
410               zzd_lw = zcof * MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp )   &  ! lower diagonal
411                  &          / (  e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  )
412# else
413               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
414                  &          / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  ) )
415               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
416                  &          / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) )
417# endif
418               !                                   ! shear prod. at w-point weightened by mask
419               zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
420                  &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
421               !
422               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
423               zd_lw(ji,jj,jk) = zzd_lw
424               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
425               !
426               !                                   ! right hand side in en
427               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
428                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
429                  &                                 * wmask(ji,jj,jk)
430            END DO
431         END DO
432      END DO
433      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
434      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
435!$OMP DO schedule(static) private(jj, ji)
436         DO jj = 2, jpjm1
437            DO ji = fs_2, fs_jpim1    ! vector opt.
438               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
439            END DO
440         END DO
441      END DO
442!$OMP DO schedule(static) private(jj, ji)
443      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
444         DO ji = fs_2, fs_jpim1   ! vector opt.
445            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
446         END DO
447      END DO
448      DO jk = 3, jpkm1
449!$OMP DO schedule(static) private(jj, ji)
450         DO jj = 2, jpjm1
451            DO ji = fs_2, fs_jpim1    ! vector opt.
452               zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
453            END DO
454         END DO
455      END DO
456!$OMP DO schedule(static) private(jj, ji)
457      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
458         DO ji = fs_2, fs_jpim1   ! vector opt.
459            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
460         END DO
461      END DO
462      DO jk = jpk-2, 2, -1
463!$OMP DO schedule(static) private(jj, ji)
464         DO jj = 2, jpjm1
465            DO ji = fs_2, fs_jpim1    ! vector opt.
466               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
467            END DO
468         END DO
469      END DO
470!$OMP DO schedule(static) private(jk,jj, ji)
471      DO jk = 2, jpkm1                             ! set the minimum value of tke
472         DO jj = 2, jpjm1
473            DO ji = fs_2, fs_jpim1   ! vector opt.
474               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
475            END DO
476         END DO
477      END DO
478!$OMP END DO NOWAIT
479!$OMP END PARALLEL
480
481      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
482      !                            !  TKE due to surface and internal wave breaking
483      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
484!!gm BUG : in the exp  remove the depth of ssh !!!
485     
486     
487      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
488!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
489         DO jk = 2, jpkm1
490            DO jj = 2, jpjm1
491               DO ji = fs_2, fs_jpim1   ! vector opt.
492                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
493                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
494               END DO
495            END DO
496         END DO
497      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
498!$OMP PARALLEL DO schedule(static) private(jj, ji, jk)
499         DO jj = 2, jpjm1
500            DO ji = fs_2, fs_jpim1   ! vector opt.
501               jk = nmln(ji,jj)
502               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
503                  &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
504            END DO
505         END DO
506      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
507!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ztx2, zty2, ztau, zdif)
508         DO jk = 2, jpkm1
509            DO jj = 2, jpjm1
510               DO ji = fs_2, fs_jpim1   ! vector opt.
511                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
512                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
513                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
514                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
515                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
516                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
517                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
518               END DO
519            END DO
520         END DO
521      ENDIF
522      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
523      !
524      CALL wrk_dealloc( jpi,jpj,       imlc )    ! integer
525      CALL wrk_dealloc( jpi,jpj,       zhlc ) 
526      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
527      !
528      IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
529      !
530   END SUBROUTINE tke_tke
531
532
533   SUBROUTINE tke_avn
534      !!----------------------------------------------------------------------
535      !!                   ***  ROUTINE tke_avn  ***
536      !!
537      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
538      !!
539      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
540      !!              the tke_tke routine). First, the now mixing lenth is
541      !!      computed from en and the strafification (N^2), then the mixings
542      !!      coefficients are computed.
543      !!              - Mixing length : a first evaluation of the mixing lengh
544      !!      scales is:
545      !!                      mxl = sqrt(2*en) / N 
546      !!      where N is the brunt-vaisala frequency, with a minimum value set
547      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
548      !!        The mixing and dissipative length scale are bound as follow :
549      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
550      !!                        zmxld = zmxlm = mxl
551      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
552      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
553      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
554      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
555      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
556      !!                    the surface to obtain ldown. the resulting length
557      !!                    scales are:
558      !!                        zmxld = sqrt( lup * ldown )
559      !!                        zmxlm = min ( lup , ldown )
560      !!              - Vertical eddy viscosity and diffusivity:
561      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
562      !!                      avt = max( avmb, pdlr * avm ) 
563      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
564      !!
565      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
566      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
567      !!----------------------------------------------------------------------
568      INTEGER  ::   ji, jj, jk   ! dummy loop indices
569      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
570      REAL(wp) ::   zdku, zri, zsqen            !   -      -
571      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
572      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
573      !!--------------------------------------------------------------------
574      !
575      IF( nn_timing == 1 )  CALL timing_start('tke_avn')
576
577      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
578
579      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
580      !                     !  Mixing length
581      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
582      !
583      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
584      !
585      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
586!$OMP PARALLEL WORKSHARE
587      zmxlm(:,:,:)  = rmxl_min   
588      zmxld(:,:,:)  = rmxl_min
589!$OMP END PARALLEL WORKSHARE
590      !
591      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
592!$OMP PARALLEL DO schedule(static) private(jj, ji, zraug)
593         DO jj = 2, jpjm1
594            DO ji = fs_2, fs_jpim1
595               zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
596               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
597            END DO
598         END DO
599      ELSE 
600!$OMP PARALLEL WORKSHARE
601         zmxlm(:,:,1) = rn_mxl0
602!$OMP END PARALLEL WORKSHARE
603      ENDIF
604      !
605!$OMP PARALLEL
606!$OMP DO schedule(static) private(jk, jj, ji, zrn2)
607      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
608         DO jj = 2, jpjm1
609            DO ji = fs_2, fs_jpim1   ! vector opt.
610               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
611               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
612            END DO
613         END DO
614      END DO
615      !
616      !                     !* Physical limits for the mixing length
617      !
618!$OMP WORKSHARE
619      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
620      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
621!$OMP END WORKSHARE NOWAIT
622!$OMP END PARALLEL
623      !
624      SELECT CASE ( nn_mxl )
625      !
626 !!gm Not sure of that coding for ISF....
627      ! where wmask = 0 set zmxlm == e3w_n
628      CASE ( 0 )           ! bounded by the distance to surface and bottom
629!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
630         DO jk = 2, jpkm1
631            DO jj = 2, jpjm1
632               DO ji = fs_2, fs_jpim1   ! vector opt.
633                  zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
634                  &            gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) )
635                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
636                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
637                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
638               END DO
639            END DO
640         END DO
641         !
642      CASE ( 1 )           ! bounded by the vertical scale factor
643!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
644         DO jk = 2, jpkm1
645            DO jj = 2, jpjm1
646               DO ji = fs_2, fs_jpim1   ! vector opt.
647                  zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )
648                  zmxlm(ji,jj,jk) = zemxl
649                  zmxld(ji,jj,jk) = zemxl
650               END DO
651            END DO
652         END DO
653         !
654      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
655!$OMP PARALLEL
656         DO jk = 2, jpkm1         ! from the surface to the bottom :
657!$OMP DO schedule(static) private(jj, ji)
658            DO jj = 2, jpjm1
659               DO ji = fs_2, fs_jpim1   ! vector opt.
660                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
661               END DO
662            END DO
663         END DO
664         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
665!$OMP DO schedule(static) private(jj, ji, zemxl)
666            DO jj = 2, jpjm1
667               DO ji = fs_2, fs_jpim1   ! vector opt.
668                  zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
669                  zmxlm(ji,jj,jk) = zemxl
670                  zmxld(ji,jj,jk) = zemxl
671               END DO
672            END DO
673         END DO
674!$OMP END PARALLEL
675         !
676      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
677!$OMP PARALLEL
678         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
679!$OMP DO schedule(static) private(jj, ji)
680            DO jj = 2, jpjm1
681               DO ji = fs_2, fs_jpim1   ! vector opt.
682                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
683               END DO
684            END DO
685         END DO
686         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
687!$OMP DO schedule(static) private(jj, ji)
688            DO jj = 2, jpjm1
689               DO ji = fs_2, fs_jpim1   ! vector opt.
690                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
691               END DO
692            END DO
693         END DO
694!$OMP DO schedule(static) private(jk, jj, ji, zemlm, zemlp)
695         DO jk = 2, jpkm1
696            DO jj = 2, jpjm1
697               DO ji = fs_2, fs_jpim1   ! vector opt.
698                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
699                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
700                  zmxlm(ji,jj,jk) = zemlm
701                  zmxld(ji,jj,jk) = zemlp
702               END DO
703            END DO
704         END DO
705!$OMP END PARALLEL
706         !
707      END SELECT
708      !
709# if defined key_c1d
710!$OMP PARALLEL WORKSHARE
711      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
712      e_mix(:,:,:) = zmxlm(:,:,:)
713!$OMP END PARALLEL WORKSHARE
714# endif
715
716      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
717      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
718      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
719!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zsqen, zav)
720      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
721         DO jj = 2, jpjm1
722            DO ji = fs_2, fs_jpim1   ! vector opt.
723               zsqen = SQRT( en(ji,jj,jk) )
724               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
725               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
726               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
727               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
728            END DO
729         END DO
730      END DO
731      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
732      !
733!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
734      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
735         DO jj = 2, jpjm1
736            DO ji = fs_2, fs_jpim1   ! vector opt.
737               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
738               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
739            END DO
740         END DO
741      END DO
742      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
743      !
744      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
745!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
746         DO jk = 2, jpkm1
747            DO jj = 2, jpjm1
748               DO ji = fs_2, fs_jpim1   ! vector opt.
749                  avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
750# if defined key_c1d
751                  e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
752!!gm bug NO zri here....
753!!gm remove the specific diag for c1d !
754                  e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk)                            ! c1d config. : save Ri
755# endif
756              END DO
757            END DO
758         END DO
759      ENDIF
760      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
761
762      IF(ln_ctl) THEN
763         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
764         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
765            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
766      ENDIF
767      !
768      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
769      !
770      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
771      !
772   END SUBROUTINE tke_avn
773
774
775   SUBROUTINE zdf_tke_init
776      !!----------------------------------------------------------------------
777      !!                  ***  ROUTINE zdf_tke_init  ***
778      !!                     
779      !! ** Purpose :   Initialization of the vertical eddy diffivity and
780      !!              viscosity when using a tke turbulent closure scheme
781      !!
782      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
783      !!              called at the first timestep (nit000)
784      !!
785      !! ** input   :   Namlist namzdf_tke
786      !!
787      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
788      !!----------------------------------------------------------------------
789      INTEGER ::   ji, jj, jk   ! dummy loop indices
790      INTEGER ::   ios
791      !!
792      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
793         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
794         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
795         &                 nn_etau , nn_htau  , rn_efr   
796      !!----------------------------------------------------------------------
797      !
798      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
799      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
800901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
801
802      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
803      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
804902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
805      IF(lwm) WRITE ( numond, namzdf_tke )
806      !
807      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
808# if defined key_zdftmx_new
809      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
810      rn_emin  = 1.e-10_wp
811      rmxl_min = 1.e-03_wp
812      IF(lwp) THEN                  ! Control print
813         WRITE(numout,*)
814         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
815         WRITE(numout,*) '~~~~~~~~~~~~'
816      ENDIF
817# else
818      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
819# endif
820      !
821      IF(lwp) THEN                    !* Control print
822         WRITE(numout,*)
823         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
824         WRITE(numout,*) '~~~~~~~~~~~~'
825         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
826         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
827         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
828         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
829         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
830         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
831         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
832         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
833         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
834         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
835         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
836         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
837         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
838         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
839         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
840         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
841         WRITE(numout,*)
842         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
843      ENDIF
844      !
845      !                              ! allocate tke arrays
846      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
847      !
848      !                               !* Check of some namelist values
849      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
850      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
851      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
852      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
853
854      IF( ln_mxl0 ) THEN
855         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
856         rn_mxl0 = rmxl_min
857      ENDIF
858     
859      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
860
861      !                               !* depth of penetration of surface tke
862      IF( nn_etau /= 0 ) THEN     
863         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
864         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
865            htau(:,:) = 10._wp
866         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
867            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
868         END SELECT
869      ENDIF
870      !                               !* set vertical eddy coef. to the background value
871!$OMP PARALLEL
872!$OMP DO schedule(static) private(jk)
873      DO jk = 1, jpk
874         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
875         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
876         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
877         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
878      END DO
879!$OMP END DO NOWAIT
880!$OMP WORKSHARE
881      dissl(:,:,:) = 1.e-12_wp
882!$OMP END WORKSHARE
883!$OMP END PARALLEL
884      !                             
885      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
886      !
887   END SUBROUTINE zdf_tke_init
888
889
890   SUBROUTINE tke_rst( kt, cdrw )
891     !!---------------------------------------------------------------------
892     !!                   ***  ROUTINE tke_rst  ***
893     !!                     
894     !! ** Purpose :   Read or write TKE file (en) in restart file
895     !!
896     !! ** Method  :   use of IOM library
897     !!                if the restart does not contain TKE, en is either
898     !!                set to rn_emin or recomputed
899     !!----------------------------------------------------------------------
900     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
901     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
902     !
903     INTEGER ::   jit, jk   ! dummy loop indices
904     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
905     !!----------------------------------------------------------------------
906     !
907     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
908        !                                   ! ---------------
909        IF( ln_rstart ) THEN                   !* Read the restart file
910           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
911           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
912           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
913           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
914           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
915           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
916           !
917           IF( id1 > 0 ) THEN                       ! 'en' exists
918              CALL iom_get( numror, jpdom_autoglo, 'en', en )
919              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
920                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
921                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
922                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
923                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
924                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
925              ELSE                                                 ! one at least array is missing
926                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
927              ENDIF
928           ELSE                                     ! No TKE array found: initialisation
929              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
930!$OMP PARALLEL WORKSHARE
931              en (:,:,:) = rn_emin * tmask(:,:,:)
932!$OMP END PARALLEL WORKSHARE
933              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
934              !
935!$OMP PARALLEL WORKSHARE
936              avt_k (:,:,:) = avt (:,:,:)
937              avm_k (:,:,:) = avm (:,:,:)
938              avmu_k(:,:,:) = avmu(:,:,:)
939              avmv_k(:,:,:) = avmv(:,:,:)
940!$OMP END PARALLEL WORKSHARE
941              !
942              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
943           ENDIF
944        ELSE                                   !* Start from rest
945!$OMP PARALLEL
946!$OMP WORKSHARE
947           en(:,:,:) = rn_emin * tmask(:,:,:)
948!$OMP END WORKSHARE NOWAIT
949!$OMP DO schedule(static) private(jk)
950           DO jk = 1, jpk                           ! set the Kz to the background value
951              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
952              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
953              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
954              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
955           END DO
956!$OMP END DO NOWAIT
957!$OMP END PARALLEL
958        ENDIF
959        !
960     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
961        !                                   ! -------------------
962        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
963        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
964        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
965        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
966        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
967        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
968        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
969        !
970     ENDIF
971     !
972   END SUBROUTINE tke_rst
973
974#else
975   !!----------------------------------------------------------------------
976   !!   Dummy module :                                        NO TKE scheme
977   !!----------------------------------------------------------------------
978   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
979CONTAINS
980   SUBROUTINE zdf_tke_init           ! Dummy routine
981   END SUBROUTINE zdf_tke_init
982   SUBROUTINE zdf_tke( kt )          ! Dummy routine
983      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
984   END SUBROUTINE zdf_tke
985   SUBROUTINE tke_rst( kt, cdrw )
986     CHARACTER(len=*) ::   cdrw
987     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
988   END SUBROUTINE tke_rst
989#endif
990
991   !!======================================================================
992END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.