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.
zdftke2.F90 in trunk/NEMO/OPA_SRC/ZDF – NEMO

source: trunk/NEMO/OPA_SRC/ZDF/zdftke2.F90 @ 1518

Last change on this file since 1518 was 1492, checked in by ctlod, 15 years ago

style and cosmetic changes
File size: 43.2 KB
Line 
1MODULE zdftke2
2   !!======================================================================
3   !!                       ***  MODULE  zdftke2  ***
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   !!             -   !  2002-08  (G. Madec)  rn_cri and Free form, F90
17   !!             -   !  2004-10  (C. Ethe )  1D configuration
18   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
19   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
20   !!                 !           - tke penetration (wind steering)
21   !!                 !           - suface condition for tke & mixing length
22   !!                 !           - Langmuir cells
23   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
24   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
25   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
26   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
27   !!                 !                                + cleaning of the parameters + bugs correction
28   !!----------------------------------------------------------------------
29#if defined key_zdftke2   ||   defined key_esopa
30   !!----------------------------------------------------------------------
31   !!   'key_zdftke2'                                  TKE vertical physics
32   !!----------------------------------------------------------------------
33   !!   zdf_tke2      : update momentum and tracer Kz from a tke scheme
34   !!   tke_tke       : tke time stepping: update tke at now time step (en)
35   !!   tke_avn       : compute mixing length scale and deduce avm and avt
36   !!   tke_init      : initialization, namelist read, and parameters control
37   !!   tke2_rst      : read/write tke restart in ocean restart file
38   !!----------------------------------------------------------------------
39   USE oce            ! ocean dynamics and active tracers
40   USE dom_oce        ! ocean space and time domain
41   USE domvvl         ! ocean space and time domain : variable volume layer
42   USE zdf_oce        ! ocean vertical physics
43   USE sbc_oce        ! surface boundary condition: ocean
44   USE phycst         ! physical constants
45   USE zdfmxl         ! mixed layer
46   USE restart        ! only for lrst_oce
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
52   IMPLICIT NONE
53   PRIVATE
54
55   PUBLIC   zdf_tke2   ! routine called in step module
56   PUBLIC   tke2_rst   ! routine called in step module
57
58   LOGICAL , PUBLIC, PARAMETER              ::   lk_zdftke2 = .TRUE.  !: TKE vertical mixing flag
59
60#if defined key_c1d
61   !                                                                !!* 1D cfg only *
62   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_dis, e_mix        !: dissipation and mixing turbulent lengh scales
63   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   e_pdl, e_ric        !: prandl and local Richardson numbers
64#endif
65
66!!gm  - variables to be suppressed from the namelist:
67   LOGICAL  ::   ln_rstke = .FALSE.         ! =T restart with tke from a run without tke
68   INTEGER  ::   nn_itke  = 50              ! number of restart iterative loops
69   INTEGER  ::   nn_ave   =  1              ! horizontal average or not on avt, avmu, avmv (=0/1)
70   REAL(wp) ::   rn_cri   = 2._wp / 9._wp   ! critic Richardson number
71   REAL(wp) ::   rn_efave = 1._wp           ! coefficient for ave : ave=rn_efave*avm
72!!gm end
73
74!!gm  - variables to be added in the namelist (with a larger default value)
75   REAL(wp) ::   rn_bshear = 1.e-20   ! backgrounf shear (>0)
76!!gm
77
78   !                                       !!! ** Namelist  namtke  **
79   LOGICAL  ::   ln_mxl0  = .FALSE.         ! mixing length scale surface value as function of wind stress or not
80   INTEGER  ::   nn_mxl   =  2              ! type of mixing length (=0/1/2/3)
81!!gm to be change: the ref value of lmin0 and lmin
82   REAL(wp) ::   rn_lmin0 = 0.4_wp          ! surface  min value of mixing length   [m]
83   REAL(wp) ::   rn_lmin  = 0.1_wp          ! interior min value of mixing length   [m]
84   INTEGER  ::   nn_pdl   =  1              ! Prandtl number or not (ratio avt/avm) (=0/1)
85   INTEGER  ::   nn_avb   =  0              ! constant or profile background on avt (=0/1)
86   REAL(wp) ::   rn_ediff = 0.1_wp          ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
87   REAL(wp) ::   rn_ediss = 0.7_wp          ! coefficient of the Kolmogoroff dissipation
88   REAL(wp) ::   rn_ebb   = 3.75_wp         ! coefficient of the surface input of tke
89   REAL(wp) ::   rn_emin  = 0.7071e-6_wp    ! minimum value of tke           [m2/s2]
90   REAL(wp) ::   rn_emin0 = 1.e-4_wp        ! surface minimum value of tke   [m2/s2]
91   INTEGER  ::   nn_havtb = 1               ! horizontal shape or not for avtb (=0/1)
92   INTEGER  ::   nn_etau  = 0               ! type of depth penetration of surface tke (=0/1/2)
93   INTEGER  ::   nn_htau  = 0               ! type of tke profile of penetration (=0/1/2)
94   REAL(wp) ::   rn_efr   = 1.0_wp          ! fraction of TKE surface value which penetrates in the ocean
95   LOGICAL  ::   ln_lc    = .FALSE.         ! Langmuir cells (LC) as a source term of TKE or not
96   REAL(wp) ::   rn_lc    = 0.15_wp         ! coef to compute vertical velocity of Langmuir cells
97
98   REAL(wp) ::   ri_cri   ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
99
100   REAL(wp), DIMENSION(jpi,jpj)     ::   htau      ! depth of tke penetration (nn_htau)
101   REAL(wp), DIMENSION(jpi,jpj)     ::   avtb_2d   ! set in tke2_init, for other modif than ice
102   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   en        ! now turbulent kinetic energy   [m2/s2]
103   REAL(wp), DIMENSION(jpi,jpj,jpk) ::   dissl     ! now mixing lenght of dissipation
104
105   !! * Substitutions
106#  include "domzgr_substitute.h90"
107#  include "vectopt_loop_substitute.h90"
108   !!----------------------------------------------------------------------
109   !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
110   !! $Id: zdftke2.F90 1201 2008-09-24 13:24:21Z rblod $
111   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
112   !!----------------------------------------------------------------------
113
114CONTAINS
115
116   SUBROUTINE zdf_tke2( kt )
117      !!----------------------------------------------------------------------
118      !!                   ***  ROUTINE zdf_tke2  ***
119      !!
120      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
121      !!              coefficients using a turbulent closure scheme (TKE).
122      !!
123      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
124      !!              is computed from a prognostic equation :
125      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
126      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
127      !!                  + avt N^2                      ! stratif. destruc.
128      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
129      !!      with the boundary conditions:
130      !!         surface: en = max( rn_emin0, rn_ebb sqrt(utau^2 + vtau^2) )
131      !!         bottom : en = rn_emin
132      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
133      !!
134      !!        The now Turbulent kinetic energy is computed using the following
135      !!      time stepping: implicit for vertical diffusion term, linearized semi
136      !!      implicit for kolmogoroff dissipation term, and explicit forward for
137      !!      both buoyancy and shear production terms. Therefore a tridiagonal
138      !!      linear system is solved. Note that buoyancy and shear terms are
139      !!      discretized in a energy conserving form (Bruchard 2002).
140      !!
141      !!        The dissipative and mixing length scale are computed from en and
142      !!      the stratification (see tke_avn)
143      !!
144      !!        The now vertical eddy vicosity and diffusivity coefficients are
145      !!      given by:
146      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
147      !!              avt = max( avmb, pdl * avm                 ) 
148      !!              eav = max( avmb, avm )
149      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
150      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
151      !!
152      !! ** Action  :   compute en (now turbulent kinetic energy)
153      !!                update avt, avmu, avmv (before vertical eddy coef.)
154      !!
155      !! References : Gaspar et al., JGR, 1990,
156      !!              Blanke and Delecluse, JPO, 1991
157      !!              Mellor and Blumberg, JPO 2004
158      !!              Axell, JGR, 2002
159      !!              Bruchard OM 2002
160      !!----------------------------------------------------------------------
161      INTEGER, INTENT(in) ::   kt   ! ocean time step
162      !!----------------------------------------------------------------------
163      !
164      IF( kt == nit000 )   CALL tke_init     ! initialisation
165                           !
166                           CALL tke_tke      ! now tke (en)
167                           !
168                           CALL tke_avn      ! now avt, avm, avmu, avmv
169      !
170   END SUBROUTINE zdf_tke2
171
172
173   SUBROUTINE tke_tke
174      !!----------------------------------------------------------------------
175      !!                   ***  ROUTINE tke_tke  ***
176      !!
177      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
178      !!
179      !! ** Method  : - TKE surface boundary condition
180      !!              - source term due to Langmuir cells (ln_lc=T)
181      !!              - source term due to shear (saved in avmu, avmv arrays)
182      !!              - Now TKE : resolution of the TKE equation by inverting
183      !!                a tridiagonal linear system by a "methode de chasse"
184      !!              - increase TKE due to surface and internal wave breaking
185      !!
186      !! ** Action  : - en : now turbulent kinetic energy)
187      !!              - avmu, avmv : production of TKE by shear at u and v-points
188      !!                (= Kz dz[Ub] * dz[Un] )
189      !! ---------------------------------------------------------------------
190      USE oce,   zdiag  =>   ua   ! use ua as workspace
191      USE oce,   zd_up  =>   va   ! use va as workspace
192      USE oce,   zd_lw  =>   ta   ! use ta as workspace
193      !!
194      INTEGER  ::   ji, jj, jk                ! dummy loop arguments
195      REAL(wp) ::   zbbrau, zesurf, zesh2     ! temporary scalars
196      REAL(wp) ::   zfact1, zfact2, zfact3    !    -         -
197      REAL(wp) ::   ztx2  , zty2  , zcof      !    -         -
198      REAL(wp) ::   zus   , zwlc  , zind      !    -         -
199      REAL(wp) ::   zzd_up, zzd_lw            !    -         -
200      INTEGER , DIMENSION(jpi,jpj)     ::   imlc    ! 2D workspace
201      REAL(wp), DIMENSION(jpi,jpj)     ::   zhlc    !  -      -
202      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc   ! 3D workspace
203      !!--------------------------------------------------------------------
204      !
205      zbbrau =  .5 * rn_ebb / rau0       ! Local constant initialisation
206      zfact1 = -.5 * rdt 
207      zfact2 = 1.5 * rdt * rn_ediss
208      zfact3 = 0.5       * rn_ediss
209      !
210      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
211      !                     !  Surface boundary condition on tke
212      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
213!CDIR NOVERRCHK
214      DO jj = 2, jpjm1         ! en(1)   = rn_ebb sqrt(utau^2+vtau^2) / rau0  (min value rn_emin0)
215!CDIR NOVERRCHK
216         DO ji = fs_2, fs_jpim1   ! vector opt.
217            ztx2   = utau(ji-1,jj  ) + utau(ji,jj)
218            zty2   = vtau(ji  ,jj-1) + vtau(ji,jj)
219            zesurf = zbbrau * SQRT( ztx2 * ztx2 + zty2 * zty2 )
220            en(ji,jj,1) = MAX( zesurf, rn_emin0 ) * tmask(ji,jj,1)
221         END DO
222      END DO
223      !
224      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
225      !                     !  Bottom boundary condition on tke
226      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
227!!gm to be added soon
228      !
229      !
230      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
231      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke
232         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
233         !
234         !                        !* total energy produce by LC : cumulative sum over jk
235         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0. ) * fsdepw(:,:,1) * fse3w(:,:,1)
236         DO jk = 2, jpk
237            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0. ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
238         END DO
239         !                        !* finite Langmuir Circulation depth
240         imlc(:,:) = mbathy(:,:)         ! Initialization to the number of w ocean point mbathy
241         DO jk = jpkm1, 2, -1
242            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
243               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
244                  zus  = 0.000128 * wndm(ji,jj) * wndm(ji,jj)
245                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
246               END DO
247            END DO
248         END DO
249         !                               ! finite LC depth
250# if defined key_vectopt_loop
251         DO jj = 1, 1
252            DO ji = 1, jpij   ! vector opt. (forced unrolling)
253# else
254         DO jj = 1, jpj 
255            DO ji = 1, jpi
256# endif
257               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
258            END DO
259         END DO
260# if defined key_c1d
261         hlc(:,:) = zhlc(:,:) * tmask(:,:,1)      ! c1d configuration: save finite Langmuir Circulation depth
262# endif
263!CDIR NOVERRCHK
264         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
265!CDIR NOVERRCHK
266            DO jj = 2, jpjm1
267!CDIR NOVERRCHK
268               DO ji = fs_2, fs_jpim1   ! vector opt.
269!!gm replace here wndn by a formulation with the stress module
270                  zus  = 0.016 * wndm(ji,jj)                  ! Stokes drift
271                  !                                           ! vertical velocity due to LC
272                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
273                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
274                  !                                           ! TKE Langmuir circulation source term
275                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
276               END DO
277            END DO
278         END DO
279         !
280      ENDIF
281      !
282      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
283      !                     !  Now Turbulent kinetic energy (output in en)
284      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
285      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
286      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
287      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
288      !
289      DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
290         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
291            DO ji = 1, jpi
292               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
293                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
294                  &           / (  fse3uw_n(ji,jj,jk)         &
295                  &              * fse3uw_b(ji,jj,jk) )
296               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
297                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
298                  &                            / (  fse3vw_n(ji,jj,jk)               &
299                  &                              *  fse3vw_b(ji,jj,jk)  )
300            END DO
301         END DO
302      END DO
303      !
304      DO jk = 2, jpkm1           !* Matrix and right hand side in en
305         DO jj = 2, jpjm1
306            DO ji = fs_2, fs_jpim1   ! vector opt.
307               zcof   = zfact1 * tmask(ji,jj,jk)
308               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
309                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
310               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
311                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
312                  !                                                           ! shear prod. at w-point weightened by mask
313               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1.e0 , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
314                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1.e0 , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
315                  !
316               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
317               zd_lw(ji,jj,jk) = zzd_lw
318               zdiag(ji,jj,jk) = 1.e0 - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
319               !
320               !                                   ! right hand side in en
321               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
322                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk)
323            END DO
324         END DO
325      END DO
326      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
327      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
328         DO jj = 2, jpjm1
329            DO ji = fs_2, fs_jpim1    ! vector opt.
330               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
331            END DO
332         END DO
333      END DO
334      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
335         DO ji = fs_2, fs_jpim1    ! vector opt.
336            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
337         END DO
338      END DO
339      DO jk = 3, jpkm1
340         DO jj = 2, jpjm1
341            DO ji = fs_2, fs_jpim1    ! vector opt.
342               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)
343            END DO
344         END DO
345      END DO
346      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
347         DO ji = fs_2, fs_jpim1    ! vector opt.
348            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
349         END DO
350      END DO
351      DO jk = jpk-2, 2, -1
352         DO jj = 2, jpjm1
353            DO ji = fs_2, fs_jpim1    ! vector opt.
354               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
355            END DO
356         END DO
357      END DO
358      DO jk = 2, jpkm1                             ! set the minimum value of tke
359         DO jj = 2, jpjm1
360            DO ji = fs_2, fs_jpim1   ! vector opt.
361               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
362            END DO
363         END DO
364      END DO
365
366      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
367      !                            !  TKE due to surface and internal wave breaking
368      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
369      IF( nn_etau == 1 ) THEN           !*  penetration throughout the water column
370         DO jk = 2, jpkm1
371            DO jj = 2, jpjm1
372               DO ji = fs_2, fs_jpim1   ! vector opt.
373                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
374                     &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
375               END DO
376            END DO
377         END DO
378      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)
379         DO jj = 2, jpjm1
380            DO ji = fs_2, fs_jpim1   ! vector opt.
381               jk = nmln(ji,jj)
382               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
383                  &                                               * ( 1.e0 - fr_i(ji,jj) )  * tmask(ji,jj,jk)
384            END DO
385         END DO
386      ENDIF
387      !
388      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
389      !
390   END SUBROUTINE tke_tke
391
392
393   SUBROUTINE tke_avn
394      !!----------------------------------------------------------------------
395      !!                   ***  ROUTINE tke_avn  ***
396      !!
397      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
398      !!
399      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
400      !!              the tke_tke routine). First, the now mixing lenth is
401      !!      computed from en and the strafification (N^2), then the mixings
402      !!      coefficients are computed.
403      !!              - Mixing length : a first evaluation of the mixing lengh
404      !!      scales is:
405      !!                      mxl = sqrt(2*en) / N 
406      !!      where N is the brunt-vaisala frequency, with a minimum value set
407      !!      to rn_lmin (rn_lmin0) in the interior (surface) ocean.
408      !!        The mixing and dissipative length scale are bound as follow :
409      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
410      !!                        zmxld = zmxlm = mxl
411      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
412      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
413      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
414      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
415      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
416      !!                    the surface to obtain ldown. the resulting length
417      !!                    scales are:
418      !!                        zmxld = sqrt( lup * ldown )
419      !!                        zmxlm = min ( lup , ldown )
420      !!              - Vertical eddy viscosity and diffusivity:
421      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
422      !!                      avt = max( avmb, pdlr * avm ) 
423      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
424      !!
425      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
426      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
427      !!----------------------------------------------------------------------
428      USE oce,     zmpdl  =>   ua   ! use ua as workspace
429      USE oce,     zmxlm  =>   va   ! use va as workspace
430      USE oce,     zmxld  =>   ta   ! use ta as workspace
431      !!
432      INTEGER  ::   ji, jj, jk            ! dummy loop arguments
433      REAL(wp) ::   zrn2, zraug           ! temporary scalars
434      REAL(wp) ::   ztx2, zdku            !    -         -
435      REAL(wp) ::   zty2, zdkv            !    -         -
436      REAL(wp) ::   zcoef, zav            !    -         -
437      REAL(wp) ::   zpdlr, zri, zsqen     !    -         -
438      REAL(wp) ::   zemxl, zemlm, zemlp   !    -         -
439      !!--------------------------------------------------------------------
440
441      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
442      !                     !  Mixing length
443      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
444      !
445      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
446      !
447      IF( ln_mxl0 ) THEN         ! surface mixing length = F(stress) : l=vkarmn*2.e5*sqrt(utau^2 + vtau^2)/(rau0*g)
448!!gm  this should be useless
449         zmxlm(:,:,1) = 0.e0
450!!gm end
451         zraug = 0.5 * vkarmn * 2.e5 / ( rau0 * grav )
452         DO jj = 2, jpjm1
453!CDIR NOVERRCHK
454            DO ji = fs_2, fs_jpim1   ! vector opt.
455               ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
456               zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
457               zmxlm(ji,jj,1) = MAX(  rn_lmin0,  zraug * SQRT( ztx2 * ztx2 + zty2 * zty2 )  )
458            END DO
459         END DO
460      ELSE                       ! surface set to the minimum value
461         zmxlm(:,:,1) = rn_lmin0
462      ENDIF
463      zmxlm(:,:,jpk) = rn_lmin   ! bottom set to the interior minium value
464      !
465!CDIR NOVERRCHK
466      DO jk = 2, jpkm1           ! interior value : l=sqrt(2*e/n**2)
467!CDIR NOVERRCHK
468         DO jj = 2, jpjm1
469!CDIR NOVERRCHK
470            DO ji = fs_2, fs_jpim1   ! vector opt.
471               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
472               zmxlm(ji,jj,jk) = MAX(  rn_lmin,  SQRT( 2. * en(ji,jj,jk) / zrn2 )  )
473            END DO
474         END DO
475      END DO
476      !
477!!gm  CAUTION: to be added here the bottom boundary condition on zmxlm
478      !
479      !                     !* Physical limits for the mixing length
480      !
481      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
482      zmxld(:,:,jpk) = rn_lmin        ! bottom  set to the minimum value
483      !
484      SELECT CASE ( nn_mxl )
485      !
486      CASE ( 0 )           ! bounded by the distance to surface and bottom
487         DO jk = 2, jpkm1
488            DO jj = 2, jpjm1
489               DO ji = fs_2, fs_jpim1   ! vector opt.
490                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
491                  &            fsdepw(ji,jj,mbathy(ji,jj)) - fsdepw(ji,jj,jk) )
492                  zmxlm(ji,jj,jk) = zemxl
493                  zmxld(ji,jj,jk) = zemxl
494               END DO
495            END DO
496         END DO
497         !
498      CASE ( 1 )           ! bounded by the vertical scale factor
499         DO jk = 2, jpkm1
500            DO jj = 2, jpjm1
501               DO ji = fs_2, fs_jpim1   ! vector opt.
502                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
503                  zmxlm(ji,jj,jk) = zemxl
504                  zmxld(ji,jj,jk) = zemxl
505               END DO
506            END DO
507         END DO
508         !
509      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
510         DO jk = 2, jpkm1         ! from the surface to the bottom :
511            DO jj = 2, jpjm1
512               DO ji = fs_2, fs_jpim1   ! vector opt.
513                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
514               END DO
515            END DO
516         END DO
517         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
518            DO jj = 2, jpjm1
519               DO ji = fs_2, fs_jpim1   ! vector opt.
520                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
521                  zmxlm(ji,jj,jk) = zemxl
522                  zmxld(ji,jj,jk) = zemxl
523               END DO
524            END DO
525         END DO
526         !
527      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
528         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
529            DO jj = 2, jpjm1
530               DO ji = fs_2, fs_jpim1   ! vector opt.
531                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
532               END DO
533            END DO
534         END DO
535         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
536            DO jj = 2, jpjm1
537               DO ji = fs_2, fs_jpim1   ! vector opt.
538                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
539               END DO
540            END DO
541         END DO
542!CDIR NOVERRCHK
543         DO jk = 2, jpkm1
544!CDIR NOVERRCHK
545            DO jj = 2, jpjm1
546!CDIR NOVERRCHK
547               DO ji = fs_2, fs_jpim1   ! vector opt.
548                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
549                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
550                  zmxlm(ji,jj,jk) = zemlm
551                  zmxld(ji,jj,jk) = zemlp
552               END DO
553            END DO
554         END DO
555         !
556      END SELECT
557      !
558# if defined key_c1d
559      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
560      e_mix(:,:,:) = zmxlm(:,:,:)
561# endif
562
563      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
564      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
565      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
566!CDIR NOVERRCHK
567      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
568!CDIR NOVERRCHK
569         DO jj = 2, jpjm1
570!CDIR NOVERRCHK
571            DO ji = fs_2, fs_jpim1   ! vector opt.
572               zsqen = SQRT( en(ji,jj,jk) )
573               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
574               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
575               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
576               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
577            END DO
578         END DO
579      END DO
580      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
581      !
582      DO jk = 2, jpkm1            !* vertical eddy viscosity at u- and v-points
583         DO jj = 2, jpjm1
584            DO ji = fs_2, fs_jpim1   ! vector opt.
585               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
586               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
587            END DO
588         END DO
589      END DO
590      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
591      !
592      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
593         DO jk = 2, jpkm1
594            DO jj = 2, jpjm1
595               DO ji = fs_2, fs_jpim1   ! vector opt.
596                  zcoef = 0.5 / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
597                  !                                          ! shear
598                  zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) )   &
599                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
600                  zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) )   &
601                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
602                  !                                          ! local Richardson number
603                  zri   = MAX( rn2b(ji,jj,jk), 0. ) * avm(ji,jj,jk) / (zcoef * (zdku + zdkv + rn_bshear ) )
604                  zpdlr = MAX(  0.1,  0.2 / MAX( 0.2 , zri )  )
605!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
606!!gm              zpdlr = MAX(  0.1,  ri_crit / MAX( ri_crit , zri )  )
607                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
608# if defined key_c1d
609                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
610                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
611# endif
612              END DO
613            END DO
614         END DO
615      ENDIF
616      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
617
618      IF(ln_ctl) THEN
619         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
620         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
621            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
622      ENDIF
623      !
624   END SUBROUTINE tke_avn
625
626
627   SUBROUTINE tke_init
628      !!----------------------------------------------------------------------
629      !!                  ***  ROUTINE zdf_tke2_init  ***
630      !!                     
631      !! ** Purpose :   Initialization of the vertical eddy diffivity and
632      !!              viscosity when using a tke turbulent closure scheme
633      !!
634      !! ** Method  :   Read the namtke namelist and check the parameters
635      !!              called at the first timestep (nit000)
636      !!
637      !! ** input   :   Namlist namtke
638      !!
639      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
640      !!----------------------------------------------------------------------
641      INTEGER ::   ji, jj, jk   ! dummy loop indices
642      !!
643      NAMELIST/namtke/ ln_rstke, rn_ediff, rn_ediss, rn_ebb  , rn_efave, rn_emin,   &
644         &             rn_emin0, rn_cri  , nn_itke , nn_mxl  , nn_pdl  , nn_ave ,   &
645         &             nn_avb  , ln_mxl0 , rn_lmin , rn_lmin0, nn_havtb, nn_etau,   &
646         &             nn_htau , rn_efr  , ln_lc   , rn_lc 
647      !!----------------------------------------------------------------------
648
649      REWIND ( numnam )               !* Read Namelist namtke : Turbulente Kinetic Energy
650      READ   ( numnam, namtke )
651     
652      ri_cri = 2. / ( 2. + rn_ediss / rn_ediff )      ! resulting critical Richardson number
653
654      IF(lwp) THEN                    !* Control print
655         WRITE(numout,*)
656         WRITE(numout,*) 'zdf_tke2 : tke turbulent closure scheme - initialisation'
657         WRITE(numout,*) '~~~~~~~~'
658         WRITE(numout,*) '          Namelist namtke : set tke mixing parameters'
659         WRITE(numout,*) '             restart with tke from no tke              ln_rstke = ', ln_rstke
660         WRITE(numout,*) '             coef. to compute avt                      rn_ediff = ', rn_ediff
661         WRITE(numout,*) '             Kolmogoroff dissipation coef.             rn_ediss = ', rn_ediss
662         WRITE(numout,*) '             tke surface input coef.                   rn_ebb   = ', rn_ebb
663         WRITE(numout,*) '             minimum value of tke                      rn_emin  = ', rn_emin
664         WRITE(numout,*) '             surface minimum value of tke              rn_emin0 = ', rn_emin0
665         WRITE(numout,*) '             mixing length type                        nn_mxl   = ', nn_mxl
666         WRITE(numout,*) '             prandl number flag                        nn_pdl   = ', nn_pdl
667         WRITE(numout,*) '             constant background or profile            nn_avb   = ', nn_avb
668         WRITE(numout,*) '             surface mixing length = F(stress) or not  ln_mxl0  = ', ln_mxl0
669         WRITE(numout,*) '             surface  mixing length minimum value      rn_lmin0 = ', rn_lmin0
670         WRITE(numout,*) '             interior mixing length minimum value      rn_lmin0 = ', rn_lmin0
671         WRITE(numout,*) '             horizontal variation for avtb             nn_havtb = ', nn_havtb
672         WRITE(numout,*) '             test param. to add tke induced by wind    nn_etau  = ', nn_etau
673         WRITE(numout,*) '             flag for computation of exp. tke profile  nn_htau  = ', nn_htau
674         WRITE(numout,*) '             % of rn_emin0 which pene. the thermocline rn_efr   = ', rn_efr
675         WRITE(numout,*) '             flag to take into acc.  Langmuir circ.    ln_lc    = ', ln_lc
676         WRITE(numout,*) '             coef to compute verticla velocity of LC   rn_lc    = ', rn_lc
677         WRITE(numout,*)
678         WRITE(numout,*) '             critical Richardson nb with your choice of coefs.  = ', ri_cri
679      ENDIF
680
681      !                               !* Check of some namelist values
682      IF( nn_mxl  < 0    .OR. nn_mxl  > 3   )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
683      IF( nn_pdl  < 0    .OR. nn_pdl  > 1   )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
684      IF( nn_htau < 0    .OR. nn_htau > 2   )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
685      IF( rn_lc   < 0.15 .OR. rn_lc   > 0.2 )   CALL ctl_stop( 'bad value: rn_lc must be between 0.15 and 0.2 ' )
686
687      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
688
689      !                               !* Background eddy viscosity and diffusivity profil
690      IF( nn_avb == 0 ) THEN                ! Define avmb, avtb from namelist parameter
691         avmb(:) = avm0
692         avtb(:) = avt0
693      ELSE                                  ! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
694         avmb(:) = avm0
695         avtb(:) = avt0 + ( 3.0e-4 - 2 * avt0 ) * 1.0e-4 * gdepw_0(:)   ! m2/s
696         IF(ln_sco .AND. lwp)   CALL ctl_warn( '          avtb profile not valid in sco' )
697      ENDIF
698      !
699      !                                     ! 2D shape of the avtb
700      avtb_2d(:,:) = 1.e0                        ! uniform
701      !
702      IF( nn_havtb == 1 ) THEN                   ! decrease avtb in the equatorial band
703           !  -15S -5S : linear decrease from avt0 to avt0/10.
704           !  -5S  +5N : cst value avt0/10.
705           !   5N  15N : linear increase from avt0/10, to avt0
706           WHERE(-15. <= gphit .AND. gphit < -5 )   avtb_2d = (1.  - 0.09 * (gphit + 15.))
707           WHERE( -5. <= gphit .AND. gphit <  5 )   avtb_2d =  0.1
708           WHERE(  5. <= gphit .AND. gphit < 15 )   avtb_2d = (0.1 + 0.09 * (gphit -  5.))
709      ENDIF
710
711      !                               !* depth of penetration of surface tke
712      IF( nn_etau /= 0 ) THEN     
713         SELECT CASE( nn_htau )           ! Choice of the depth of penetration
714         CASE( 0 )                                    ! constant depth penetration (here 10 meters)
715            htau(:,:) = 10.e0
716         CASE( 1 )                                    ! F(latitude) : 5m to 40m at high lat.
717            DO jj = 1, jpj
718               DO ji = 1, jpi 
719                  htau(ji,jj) = MAX( 5., MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
720               END DO
721            END DO
722         CASE( 2 )                                    ! F(latitude) : 5m to 60m at high lat.
723            DO jj = 1, jpj
724               DO ji = 1, jpi
725                  htau(ji,jj) =  MAX( 5.,6./4.* MIN( 40., 60.*ABS( SIN( rpi/180. * gphit(ji,jj) ) ) )   )
726               END DO
727            END DO
728         END SELECT
729      ENDIF
730
731      !                               !* set vertical eddy coef. to the background value
732      DO jk = 1, jpk
733         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
734         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
735         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
736         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
737      END DO
738      dissl(:,:,:) = 1.e-12
739      !                               !* read or initialize all required files
740      CALL tke2_rst( nit000, 'READ' )
741      !
742   END SUBROUTINE tke_init
743
744
745   SUBROUTINE tke2_rst( kt, cdrw )
746     !!---------------------------------------------------------------------
747     !!                   ***  ROUTINE ts_rst  ***
748     !!                     
749     !! ** Purpose :   Read or write TKE file (en) in restart file
750     !!
751     !! ** Method  :   use of IOM library
752     !!                if the restart does not contain TKE, en is either
753     !!                set to rn_emin or recomputed (nn_itke/=0)
754     !!----------------------------------------------------------------------
755     INTEGER         , INTENT(in) ::   kt         ! ocean time-step
756     CHARACTER(len=*), INTENT(in) ::   cdrw       ! "READ"/"WRITE" flag
757     !
758     INTEGER ::   jit, jk   ! dummy loop indices
759     INTEGER ::   id1, id2, id3, id4, id5, id6
760     !!----------------------------------------------------------------------
761     !
762     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
763        !                                   ! ---------------
764        IF( ln_rstart ) THEN                   !* Read the restart file
765           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
766           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
767           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
768           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
769           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
770           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
771           !
772           IF( id1 > 0 ) THEN                       ! 'en' exists
773              CALL iom_get( numror, jpdom_autoglo, 'en', en )
774              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
775                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
776                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
777                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
778                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
779                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
780              ELSE                                                 ! one at least array is missing
781                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
782              ENDIF
783           ELSE                                     ! No TKE array found: initialisation
784              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
785              en (:,:,:) = rn_emin * tmask(:,:,:)
786              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
787              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke2( jit )   ;   END DO
788           ENDIF
789        ELSE                                   !* Start from rest
790           en(:,:,:) = rn_emin * tmask(:,:,:)
791           DO jk = 1, jpk                           ! set the Kz to the background value
792              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
793              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
794              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
795              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
796           END DO
797        ENDIF
798        !
799     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
800        !                                   ! -------------------
801        IF(lwp) WRITE(numout,*) '---- tke2-rst ----'
802        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
803        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt    )
804        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm    )
805        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu   )
806        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv   )
807        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
808        !
809     ENDIF
810     !
811   END SUBROUTINE tke2_rst
812
813#else
814   !!----------------------------------------------------------------------
815   !!   Dummy module :                                        NO TKE scheme
816   !!----------------------------------------------------------------------
817   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke2 = .FALSE.   !: TKE flag
818CONTAINS
819   SUBROUTINE zdf_tke2( kt )          ! Empty routine
820      WRITE(*,*) 'zdf_tke2: You should not have seen this print! error?', kt
821   END SUBROUTINE zdf_tke2
822   SUBROUTINE tke2_rst( kt, cdrw )
823     CHARACTER(len=*) ::   cdrw
824     WRITE(*,*) 'tke2_rst: You should not have seen this print! error?', kt, cdwr
825   END SUBROUTINE tke2_rst
826#endif
827
828   !!======================================================================
829END MODULE zdftke2
Note: See TracBrowser for help on using the repository browser.