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

Last change on this file since 1417 was 1268, checked in by ctlod, 15 years ago

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