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

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

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