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/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/2017/wrk_OMP_test_for_Silvia/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 8279

Last change on this file since 8279 was 8279, checked in by mocavero, 7 years ago

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