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

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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 2382

Last change on this file since 2382 was 2382, checked in by gm, 13 years ago

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