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 @ 7525

Last change on this file since 7525 was 7525, checked in by mocavero, 7 years ago

changes after review
  • Property svn:keywords set to Id
File size: 51.9 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      INTEGER             ::   jk, jj, ji 
174      !!----------------------------------------------------------------------
175      !
176#if defined key_agrif 
177      ! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
178      IF( .NOT.Agrif_Root() )   CALL Agrif_Tke
179#endif
180      !
181      IF( kt /= nit000 ) THEN   ! restore before value to compute tke
182!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
183         DO jk = 1, jpk
184            DO jj = 1, jpj
185               DO ji = 1, jpi
186                  avt (ji,jj,jk) = avt_k (ji,jj,jk) 
187                  avm (ji,jj,jk) = avm_k (ji,jj,jk) 
188                  avmu(ji,jj,jk) = avmu_k(ji,jj,jk) 
189                  avmv(ji,jj,jk) = avmv_k(ji,jj,jk) 
190               END DO
191            END DO
192         END DO
193      ENDIF 
194      !
195      CALL tke_tke      ! now tke (en)
196      !
197      CALL tke_avn      ! now avt, avm, avmu, avmv
198      !
199!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
200         DO jk = 1, jpk
201            DO jj = 1, jpj
202               DO ji = 1, jpi
203                  avt_k (ji,jj,jk) = avt (ji,jj,jk) 
204                  avm_k (ji,jj,jk) = avm (ji,jj,jk) 
205                  avmu_k(ji,jj,jk) = avmu(ji,jj,jk) 
206                  avmv_k(ji,jj,jk) = avmv(ji,jj,jk) 
207               END DO
208            END DO
209         END DO
210      !
211#if defined key_agrif
212      ! Update child grid f => parent grid
213      IF( .NOT.Agrif_Root() )   CALL Agrif_Update_Tke( kt )      ! children only
214#endif     
215     !
216  END SUBROUTINE zdf_tke
217
218
219   SUBROUTINE tke_tke
220      !!----------------------------------------------------------------------
221      !!                   ***  ROUTINE tke_tke  ***
222      !!
223      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
224      !!
225      !! ** Method  : - TKE surface boundary condition
226      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
227      !!              - source term due to shear (saved in avmu, avmv arrays)
228      !!              - Now TKE : resolution of the TKE equation by inverting
229      !!                a tridiagonal linear system by a "methode de chasse"
230      !!              - increase TKE due to surface and internal wave breaking
231      !!
232      !! ** Action  : - en : now turbulent kinetic energy)
233      !!              - avmu, avmv : production of TKE by shear at u and v-points
234      !!                (= Kz dz[Ub] * dz[Un] )
235      !! ---------------------------------------------------------------------
236      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
237!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
238!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
239      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
240      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
241      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
242      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
243      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
244      REAL(wp) ::   ztau  , zdif                    !    -         -
245      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
246      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
247!!bfr      REAL(wp) ::   zebot                           !    -         -
248      INTEGER , POINTER, DIMENSION(:,:  ) ::   imlc
249      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zhlc
250      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv
251      REAL(wp)                            ::   zri  !   local Richardson number
252      !!--------------------------------------------------------------------
253      !
254      IF( nn_timing == 1 )  CALL timing_start('tke_tke')
255      !
256      CALL wrk_alloc( jpi,jpj,       imlc )    ! integer
257      CALL wrk_alloc( jpi,jpj,       zhlc ) 
258      CALL wrk_alloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
259      !
260      zbbrau = rn_ebb / rau0       ! Local constant initialisation
261      zfact1 = -.5_wp * rdt 
262      zfact2 = 1.5_wp * rdt * rn_ediss
263      zfact3 = 0.5_wp       * rn_ediss
264      !
265      !
266      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
267      !                     !  Surface boundary condition on tke
268      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
269      IF ( ln_isfcav ) THEN
270!$OMP PARALLEL DO schedule(static) private(jj, ji)
271         DO jj = 2, jpjm1            ! en(mikt(ji,jj))   = rn_emin
272            DO ji = fs_2, fs_jpim1   ! vector opt.
273               en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
274            END DO
275         END DO
276      END IF
277!$OMP PARALLEL DO schedule(static) private(jj, ji)
278      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
279         DO ji = fs_2, fs_jpim1   ! vector opt.
280            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
281         END DO
282      END DO
283     
284!!bfr   - start commented area
285      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
286      !                     !  Bottom boundary condition on tke
287      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
288      !
289      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
290      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
291      ! The condition is coded here for completion but commented out until there is proof that the
292      ! computational cost is justified
293      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
294      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
295!!    DO jj = 2, jpjm1
296!!       DO ji = fs_2, fs_jpim1   ! vector opt.
297!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
298!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
299!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
300!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
301!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
302!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
303!!       END DO
304!!    END DO
305!!bfr   - end commented area
306      !
307      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
308      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
309         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
310         !
311         !                        !* total energy produce by LC : cumulative sum over jk
312!$OMP PARALLEL
313!$OMP DO schedule(static) private(jj, ji)
314         DO jj =1, jpj
315            DO ji=1, jpi
316               zpelc(ji,jj,1) =  MAX( rn2b(ji,jj,1), 0._wp ) * gdepw_n(ji,jj,1) * e3w_n(ji,jj,1)
317            END DO
318         END DO
319         DO jk = 2, jpk
320!$OMP DO schedule(static) private(jj, ji)
321            DO jj =1, jpj
322               DO ji=1, jpi
323                  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)
324               END DO
325            END DO
326         END DO
327         !                        !* finite Langmuir Circulation depth
328         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
329!$OMP DO schedule(static) private(jj,ji)
330            DO jj = 1, jpj
331               DO ji = 1, jpi
332                  imlc(ji,jj) = mbkt(ji,jj) + 1       ! Initialization to the number of w ocean point (=2 over land)
333               END DO
334            END DO
335         DO jk = jpkm1, 2, -1
336!$OMP DO schedule(static) private(jj, ji, zus)
337            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
338               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
339                  zus  = zcof * taum(ji,jj)
340                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
341               END DO
342            END DO
343         END DO
344         !                               ! finite LC depth
345!$OMP DO schedule(static) private(jj, ji)
346         DO jj = 1, jpj 
347            DO ji = 1, jpi
348               zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))
349            END DO
350         END DO
351         zcof = 0.016 / SQRT( zrhoa * zcdrag )
352!$OMP DO schedule(static) private(jk, jj, ji, zus, zind, zwlc)
353         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
354            DO jj = 2, jpjm1
355               DO ji = fs_2, fs_jpim1   ! vector opt.
356                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
357                  !                                           ! vertical velocity due to LC
358                  zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) )
359                  zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )
360                  !                                           ! TKE Langmuir circulation source term
361                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
362                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
363               END DO
364            END DO
365         END DO
366!$OMP END DO NOWAIT
367!$OMP END PARALLEL
368         !
369      ENDIF
370      !
371      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
372      !                     !  Now Turbulent kinetic energy (output in en)
373      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
374      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
375      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
376      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
377      !
378!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
379      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
380         DO jj = 1, jpjm1
381            DO ji = 1, fs_jpim1   ! vector opt.
382               z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji+1,jj,jk) )   &
383                  &                 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
384                  &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
385                  &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
386               z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk  ) + avm(ji,jj+1,jk) )   &
387                  &                 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
388                  &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
389                  &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
390            END DO
391         END DO
392      END DO
393      !
394      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: compute apdlr
395         ! Note that zesh2 is also computed in the next loop.
396         ! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops
397!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zesh2, zri)
398         DO jk = 2, jpkm1
399            DO jj = 2, jpjm1
400               DO ji = fs_2, fs_jpim1   ! vector opt.
401                  !                                          ! shear prod. at w-point weightened by mask
402                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
403                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
404                  !                                          ! local Richardson number
405                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear )
406                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
407                 
408               END DO
409            END DO
410         END DO
411         !
412      ENDIF
413      !         
414!$OMP PARALLEL
415!$OMP DO schedule(static) private(jk, jj, ji, zcof, zzd_up, zzd_lw, zesh2)
416      DO jk = 2, jpkm1           !* Matrix and right hand side in en
417         DO jj = 2, jpjm1
418            DO ji = fs_2, fs_jpim1   ! vector opt.
419               zcof   = zfact1 * tmask(ji,jj,jk)
420# if defined key_zdftmx_new
421               ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
422               zzd_up = zcof * MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp )   &  ! upper diagonal
423                  &          / (  e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  )  )
424               zzd_lw = zcof * MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp )   &  ! lower diagonal
425                  &          / (  e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  )  )
426# else
427               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
428                  &          / ( e3t_n(ji,jj,jk  ) * e3w_n(ji,jj,jk  ) )
429               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
430                  &          / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk  ) )
431# endif
432               !                                   ! shear prod. at w-point weightened by mask
433               zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
434                  &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
435               !
436               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
437               zd_lw(ji,jj,jk) = zzd_lw
438               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
439               !
440               !                                   ! right hand side in en
441               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
442                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) &
443                  &                                 * wmask(ji,jj,jk)
444            END DO
445         END DO
446      END DO
447      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
448      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
449!$OMP DO schedule(static) private(jj, ji)
450         DO jj = 2, jpjm1
451            DO ji = fs_2, fs_jpim1    ! vector opt.
452               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(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                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
458         DO ji = fs_2, fs_jpim1   ! vector opt.
459            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
460         END DO
461      END DO
462      DO jk = 3, jpkm1
463!$OMP DO schedule(static) private(jj, ji)
464         DO jj = 2, jpjm1
465            DO ji = fs_2, fs_jpim1    ! vector opt.
466               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)
467            END DO
468         END DO
469      END DO
470!$OMP DO schedule(static) private(jj, ji)
471      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
472         DO ji = fs_2, fs_jpim1   ! vector opt.
473            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
474         END DO
475      END DO
476      DO jk = jpk-2, 2, -1
477!$OMP DO schedule(static) private(jj, ji)
478         DO jj = 2, jpjm1
479            DO ji = fs_2, fs_jpim1    ! vector opt.
480               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
481            END DO
482         END DO
483      END DO
484!$OMP DO schedule(static) private(jk,jj, ji)
485      DO jk = 2, jpkm1                             ! set the minimum value of tke
486         DO jj = 2, jpjm1
487            DO ji = fs_2, fs_jpim1   ! vector opt.
488               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
489            END DO
490         END DO
491      END DO
492!$OMP END DO NOWAIT
493!$OMP END PARALLEL
494
495      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
496      !                            !  TKE due to surface and internal wave breaking
497      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
498!!gm BUG : in the exp  remove the depth of ssh !!!
499     
500     
501      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
502!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
503         DO jk = 2, jpkm1
504            DO jj = 2, jpjm1
505               DO ji = fs_2, fs_jpim1   ! vector opt.
506                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
507                     &                                 * ( 1._wp - fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
508               END DO
509            END DO
510         END DO
511      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
512!$OMP PARALLEL DO schedule(static) private(jj, ji, jk)
513         DO jj = 2, jpjm1
514            DO ji = fs_2, fs_jpim1   ! vector opt.
515               jk = nmln(ji,jj)
516               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * 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      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
521!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ztx2, zty2, ztau, zdif)
522         DO jk = 2, jpkm1
523            DO jj = 2, jpjm1
524               DO ji = fs_2, fs_jpim1   ! vector opt.
525                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
526                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
527                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
528                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
529                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
530                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) )   &
531                     &                        * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
532               END DO
533            END DO
534         END DO
535      ENDIF
536      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
537      !
538      CALL wrk_dealloc( jpi,jpj,       imlc )    ! integer
539      CALL wrk_dealloc( jpi,jpj,       zhlc ) 
540      CALL wrk_dealloc( jpi,jpj,jpk,   zpelc, zdiag, zd_up, zd_lw, z3du, z3dv ) 
541      !
542      IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
543      !
544   END SUBROUTINE tke_tke
545
546
547   SUBROUTINE tke_avn
548      !!----------------------------------------------------------------------
549      !!                   ***  ROUTINE tke_avn  ***
550      !!
551      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
552      !!
553      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
554      !!              the tke_tke routine). First, the now mixing lenth is
555      !!      computed from en and the strafification (N^2), then the mixings
556      !!      coefficients are computed.
557      !!              - Mixing length : a first evaluation of the mixing lengh
558      !!      scales is:
559      !!                      mxl = sqrt(2*en) / N 
560      !!      where N is the brunt-vaisala frequency, with a minimum value set
561      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
562      !!        The mixing and dissipative length scale are bound as follow :
563      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
564      !!                        zmxld = zmxlm = mxl
565      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
566      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
567      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
568      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
569      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
570      !!                    the surface to obtain ldown. the resulting length
571      !!                    scales are:
572      !!                        zmxld = sqrt( lup * ldown )
573      !!                        zmxlm = min ( lup , ldown )
574      !!              - Vertical eddy viscosity and diffusivity:
575      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
576      !!                      avt = max( avmb, pdlr * avm ) 
577      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
578      !!
579      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
580      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
581      !!----------------------------------------------------------------------
582      INTEGER  ::   ji, jj, jk   ! dummy loop indices
583      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
584      REAL(wp) ::   zdku, zri, zsqen            !   -      -
585      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
586      REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
587      !!--------------------------------------------------------------------
588      !
589      IF( nn_timing == 1 )  CALL timing_start('tke_avn')
590
591      CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
592
593      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
594      !                     !  Mixing length
595      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
596      !
597      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
598      !
599      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
600!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
601      DO jk = 1, jpk
602         DO jj = 1, jpj
603            DO ji = 1, jpi
604               zmxlm(ji,jj,jk)  = rmxl_min   
605               zmxld(ji,jj,jk)  = rmxl_min
606            END DO
607         END DO
608      END DO
609      !
610      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
611!$OMP PARALLEL DO schedule(static) private(jj, ji, zraug)
612         DO jj = 2, jpjm1
613            DO ji = fs_2, fs_jpim1
614               zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
615               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
616            END DO
617         END DO
618      ELSE 
619!$OMP PARALLEL DO schedule(static) private(jj,ji)
620         DO jj = 1, jpj
621            DO ji = 1, jpi
622               zmxlm(ji,jj,1) = rn_mxl0
623            END DO
624         END DO
625      ENDIF
626      !
627!$OMP PARALLEL
628!$OMP DO schedule(static) private(jk, jj, ji, zrn2)
629      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
630         DO jj = 2, jpjm1
631            DO ji = fs_2, fs_jpim1   ! vector opt.
632               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
633               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
634            END DO
635         END DO
636      END DO
637      !
638      !                     !* Physical limits for the mixing length
639      !
640!$OMP DO schedule(static) private(jj,ji)
641      DO jj = 1, jpj
642         DO ji = 1, jpi
643            zmxld(ji,jj, 1 ) = zmxlm(ji,jj,1)   ! surface set to the minimum value
644            zmxld(ji,jj,jpk) = rmxl_min       ! last level  set to the minimum value
645         END DO
646      END DO
647!$OMP END PARALLEL
648      !
649      SELECT CASE ( nn_mxl )
650      !
651 !!gm Not sure of that coding for ISF....
652      ! where wmask = 0 set zmxlm == e3w_n
653      CASE ( 0 )           ! bounded by the distance to surface and bottom
654!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
655         DO jk = 2, jpkm1
656            DO jj = 2, jpjm1
657               DO ji = fs_2, fs_jpim1   ! vector opt.
658                  zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
659                  &            gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) )
660                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
661                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
662                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
663               END DO
664            END DO
665         END DO
666         !
667      CASE ( 1 )           ! bounded by the vertical scale factor
668!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
669         DO jk = 2, jpkm1
670            DO jj = 2, jpjm1
671               DO ji = fs_2, fs_jpim1   ! vector opt.
672                  zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )
673                  zmxlm(ji,jj,jk) = zemxl
674                  zmxld(ji,jj,jk) = zemxl
675               END DO
676            END DO
677         END DO
678         !
679      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
680!$OMP PARALLEL
681         DO jk = 2, jpkm1         ! from the surface to the bottom :
682!$OMP DO schedule(static) private(jj, ji)
683            DO jj = 2, jpjm1
684               DO ji = fs_2, fs_jpim1   ! vector opt.
685                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
686               END DO
687            END DO
688         END DO
689         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
690!$OMP DO schedule(static) private(jj, ji, zemxl)
691            DO jj = 2, jpjm1
692               DO ji = fs_2, fs_jpim1   ! vector opt.
693                  zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
694                  zmxlm(ji,jj,jk) = zemxl
695                  zmxld(ji,jj,jk) = zemxl
696               END DO
697            END DO
698         END DO
699!$OMP END PARALLEL
700         !
701      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
702!$OMP PARALLEL
703         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
704!$OMP DO schedule(static) private(jj, ji)
705            DO jj = 2, jpjm1
706               DO ji = fs_2, fs_jpim1   ! vector opt.
707                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
708               END DO
709            END DO
710         END DO
711         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
712!$OMP DO schedule(static) private(jj, ji)
713            DO jj = 2, jpjm1
714               DO ji = fs_2, fs_jpim1   ! vector opt.
715                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
716               END DO
717            END DO
718         END DO
719!$OMP DO schedule(static) private(jk, jj, ji, zemlm, zemlp)
720         DO jk = 2, jpkm1
721            DO jj = 2, jpjm1
722               DO ji = fs_2, fs_jpim1   ! vector opt.
723                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
724                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
725                  zmxlm(ji,jj,jk) = zemlm
726                  zmxld(ji,jj,jk) = zemlp
727               END DO
728            END DO
729         END DO
730!$OMP END PARALLEL
731         !
732      END SELECT
733      !
734# if defined key_c1d
735!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
736      DO jk = 1, jpk
737         DO jj = 1, jpj
738            DO ji = 1, jpi
739               e_dis(ji,jj,jk) = zmxld(ji,jj,jk)      ! c1d configuration : save mixing and dissipation turbulent length scales
740               e_mix(ji,jj,jk) = zmxlm(ji,jj,jk)
741            END DO
742         END DO
743      END DO
744# endif
745
746      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
747      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
748      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
749!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zsqen, zav)
750      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
751         DO jj = 2, jpjm1
752            DO ji = fs_2, fs_jpim1   ! vector opt.
753               zsqen = SQRT( en(ji,jj,jk) )
754               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
755               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
756               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
757               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
758            END DO
759         END DO
760      END DO
761      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
762      !
763!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
764      DO jk = 2, jpkm1            !* vertical eddy viscosity at wu- and wv-points
765         DO jj = 2, jpjm1
766            DO ji = fs_2, fs_jpim1   ! vector opt.
767               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
768               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
769            END DO
770         END DO
771      END DO
772      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
773      !
774      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
775!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
776         DO jk = 2, jpkm1
777            DO jj = 2, jpjm1
778               DO ji = fs_2, fs_jpim1   ! vector opt.
779                  avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
780# if defined key_c1d
781                  e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
782!!gm bug NO zri here....
783!!gm remove the specific diag for c1d !
784                  e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk)                            ! c1d config. : save Ri
785# endif
786              END DO
787            END DO
788         END DO
789      ENDIF
790      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
791
792      IF(ln_ctl) THEN
793         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
794         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
795            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
796      ENDIF
797      !
798      CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 
799      !
800      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
801      !
802   END SUBROUTINE tke_avn
803
804
805   SUBROUTINE zdf_tke_init
806      !!----------------------------------------------------------------------
807      !!                  ***  ROUTINE zdf_tke_init  ***
808      !!                     
809      !! ** Purpose :   Initialization of the vertical eddy diffivity and
810      !!              viscosity when using a tke turbulent closure scheme
811      !!
812      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
813      !!              called at the first timestep (nit000)
814      !!
815      !! ** input   :   Namlist namzdf_tke
816      !!
817      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
818      !!----------------------------------------------------------------------
819      INTEGER ::   ji, jj, jk   ! dummy loop indices
820      INTEGER ::   ios
821      !!
822      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
823         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
824         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
825         &                 nn_etau , nn_htau  , rn_efr   
826      !!----------------------------------------------------------------------
827      !
828      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
829      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
830901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
831
832      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
833      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
834902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
835      IF(lwm) WRITE ( numond, namzdf_tke )
836      !
837      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
838# if defined key_zdftmx_new
839      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
840      rn_emin  = 1.e-10_wp
841      rmxl_min = 1.e-03_wp
842      IF(lwp) THEN                  ! Control print
843         WRITE(numout,*)
844         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
845         WRITE(numout,*) '~~~~~~~~~~~~'
846      ENDIF
847# else
848      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
849# endif
850      !
851      IF(lwp) THEN                    !* Control print
852         WRITE(numout,*)
853         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
854         WRITE(numout,*) '~~~~~~~~~~~~'
855         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
856         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
857         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
858         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
859         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
860         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
861         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
862         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
863         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
864         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
865         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
866         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
867         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
868         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
869         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
870         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
871         WRITE(numout,*)
872         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
873      ENDIF
874      !
875      !                              ! allocate tke arrays
876      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
877      !
878      !                               !* Check of some namelist values
879      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
880      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
881      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
882      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
883
884      IF( ln_mxl0 ) THEN
885         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
886         rn_mxl0 = rmxl_min
887      ENDIF
888     
889      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
890
891      !                               !* depth of penetration of surface tke
892      IF( nn_etau /= 0 ) THEN     
893         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
894         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
895            htau(:,:) = 10._wp
896         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
897            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
898         END SELECT
899      ENDIF
900      !                               !* set vertical eddy coef. to the background value
901!$OMP PARALLEL
902!$OMP DO schedule(static) private(jk,jj,ji)
903      DO jk = 1, jpk
904         DO jj = 1, jpj
905            DO ji = 1, jpi
906               avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk)
907               avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk)
908               avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk)
909               avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk)
910            END DO
911         END DO
912      END DO
913!$OMP END DO NOWAIT
914!$OMP DO schedule(static) private(jk,jj,ji)
915      DO jk = 1, jpk
916         DO jj = 1, jpj
917            DO ji = 1, jpi
918               dissl(ji,jj,jk) = 1.e-12_wp
919            END DO
920         END DO
921      END DO
922!$OMP END PARALLEL
923      !                             
924      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
925      !
926   END SUBROUTINE zdf_tke_init
927
928
929   SUBROUTINE tke_rst( kt, cdrw )
930     !!---------------------------------------------------------------------
931     !!                   ***  ROUTINE tke_rst  ***
932     !!                     
933     !! ** Purpose :   Read or write TKE file (en) in restart file
934     !!
935     !! ** Method  :   use of IOM library
936     !!                if the restart does not contain TKE, en is either
937     !!                set to rn_emin or recomputed
938     !!----------------------------------------------------------------------
939     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
940     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
941     !
942     INTEGER ::   jit, jk, jj, ji   ! dummy loop indices
943     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
944     !!----------------------------------------------------------------------
945     !
946     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
947        !                                   ! ---------------
948        IF( ln_rstart ) THEN                   !* Read the restart file
949           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
950           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
951           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
952           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
953           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
954           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
955           !
956           IF( id1 > 0 ) THEN                       ! 'en' exists
957              CALL iom_get( numror, jpdom_autoglo, 'en', en )
958              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
959                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
960                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
961                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
962                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
963                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
964              ELSE                                                 ! one at least array is missing
965                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
966              ENDIF
967           ELSE                                     ! No TKE array found: initialisation
968              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
969!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
970              DO jk = 1, jpk
971                 DO jj = 1, jpj
972                    DO ji = 1, jpi
973                       en (ji,jj,jk) = rn_emin * tmask(ji,jj,jk)
974                    END DO
975                 END DO
976              END DO
977              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
978              !
979!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
980              DO jk = 1, jpk
981                 DO jj = 1, jpj
982                    DO ji = 1, jpi
983                       avt_k (ji,jj,jk) = avt (ji,jj,jk)
984                       avm_k (ji,jj,jk) = avm (ji,jj,jk)
985                       avmu_k(ji,jj,jk) = avmu(ji,jj,jk)
986                       avmv_k(ji,jj,jk) = avmv(ji,jj,jk)
987                    END DO
988                 END DO
989              END DO
990              !
991              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
992           ENDIF
993        ELSE                                   !* Start from rest
994!$OMP PARALLEL
995!$OMP DO schedule(static) private(jk,jj,ji)
996           DO jk = 1, jpk
997              DO jj = 1, jpj
998                 DO ji = 1, jpi
999                    en(ji,jj,jk) = rn_emin * tmask(ji,jj,jk)
1000                 END DO
1001              END DO
1002           END DO
1003!$OMP END DO NOWAIT
1004!$OMP DO schedule(static) private(jk)
1005           DO jk = 1, jpk                           ! set the Kz to the background value
1006              avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
1007              avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
1008              avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
1009              avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
1010           END DO
1011!$OMP END DO NOWAIT
1012!$OMP END PARALLEL
1013        ENDIF
1014        !
1015     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1016        !                                   ! -------------------
1017        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
1018        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
1019        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
1020        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
1021        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
1022        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
1023        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
1024        !
1025     ENDIF
1026     !
1027   END SUBROUTINE tke_rst
1028
1029#else
1030   !!----------------------------------------------------------------------
1031   !!   Dummy module :                                        NO TKE scheme
1032   !!----------------------------------------------------------------------
1033   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
1034CONTAINS
1035   SUBROUTINE zdf_tke_init           ! Dummy routine
1036   END SUBROUTINE zdf_tke_init
1037   SUBROUTINE zdf_tke( kt )          ! Dummy routine
1038      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
1039   END SUBROUTINE zdf_tke
1040   SUBROUTINE tke_rst( kt, cdrw )
1041     CHARACTER(len=*) ::   cdrw
1042     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
1043   END SUBROUTINE tke_rst
1044#endif
1045
1046   !!======================================================================
1047END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.