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

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

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

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