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

Last change on this file since 8056 was 8056, checked in by gm, 8 years ago

wrk_OMP: 2nd step: small error correction in zdfric & zdftke
  • Property svn:keywords set to Id
File size: 41.6 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      zbbrau = rn_ebb / rau0       ! Local constant initialisation
224      zfact1 = -.5_wp * rdt 
225      zfact2 = 1.5_wp * rdt * rn_ediss
226      zfact3 = 0.5_wp       * rn_ediss
227      !
228      !
229      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
230      !                     !  Surface boundary condition on tke
231      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
232      IF ( ln_isfcav ) THEN
233         DO jj = k_Jstr, k_Jend            ! en(mikt(ji,jj))   = rn_emin
234            DO ji = k_Istr, k_Iend
235               en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
236            END DO
237         END DO
238      END IF
239      DO jj = k_Jstr, k_Jend            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
240         DO ji = k_Istr, k_Iend
241            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
242         END DO
243      END DO
244     
245!!bfr   - start commented area
246      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
247      !                     !  Bottom boundary condition on tke
248      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
249      !
250      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
251      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
252      ! The condition is coded here for completion but commented out until there is proof that the
253      ! computational cost is justified
254      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
255      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
256!!    DO jj = k_Jstr, k_Jend
257!!       DO ji = k_Jstr, k_Iend
258!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
259!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
260!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
261!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
262!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
263!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
264!!       END DO
265!!    END DO
266!!bfr   - end commented area
267      !
268      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
269      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke   !   (Axell JGR 2002)
270         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
271         !
272         !                        !* total energy produce by LC : cumulative sum over jk
273         zpelc(WRK_2D,1) =  MAX( rn2b(WRK_2D,1), 0._wp ) * pdepw(WRK_2D,1) * p_e3w(WRK_2D,1)
274         DO jk = 2, jpk
275            zpelc(WRK_2D,jk) = zpelc(WRK_2D,jk-1)   &
276               &             + MAX( rn2b(WRK_2D,jk), 0._wp )   &
277               &             * pdepw(WRK_2D,jk) * p_e3w(WRK_2D,jk)
278         END DO
279         !                        !* finite Langmuir Circulation depth
280         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
281         imlc(WRK_2D) = mbkt(WRK_2D) + 1       ! Initialization to the number of w ocean point (=2 over land)
282         DO jk = jpkm1, 2, -1
283            DO jj = k_Jstr, k_Jend               ! Last w-level at which zpelc>=0.5*us*us
284               DO ji = k_Istr, k_Iend            !      with us=0.016*wind(starting from jpk-1)
285                  zus  = zcof * taum(ji,jj)
286                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
287               END DO
288            END DO
289         END DO
290         !                               ! finite LC depth
291         DO jj = k_Jstr, k_Jend 
292            DO ji = k_Istr, k_Iend
293               zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj))
294            END DO
295         END DO
296         zcof = 0.016 / SQRT( zrhoa * zcdrag )
297         DO jk = 2, jpkm1         !* TKE Langmuir circulation source term added to en
298            DO jj = k_Jstr, k_Jend
299               DO ji = k_Istr, k_Iend
300                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
301                  !                                           ! vertical velocity due to LC
302                  zind = 0.5 - SIGN( 0.5, pdepw(ji,jj,jk) - zhlc(ji,jj) )
303                  zwlc = zind * rn_lc * zus * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) )
304                  !                                           ! TKE Langmuir circulation source term
305                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
306                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
307               END DO
308            END DO
309         END DO
310         !
311      ENDIF
312      !
313      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
314      !                     !  Now Turbulent kinetic energy (output in en)
315      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
316      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
317      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
318      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
319      !
320      IF( nn_pdl == 1 ) THEN      !* Prandtl number = F( Ri )
321         DO jk = 2, jpkm1
322            DO jj = k_Jstr, k_Jend
323               DO ji = k_Istr, k_Iend
324                  !                             ! local Richardson number
325                  zri = MAX( rn2b(ji,jj,jk), 0._wp ) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
326                  !                             ! inverse of Prandtl number
327                  apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
328               END DO
329            END DO
330         END DO
331      ENDIF
332      !         
333      DO jk = 2, jpkm1           !* Matrix and right hand side in en
334         DO jj = k_Jstr, k_Jend
335            DO ji = k_Istr, k_Iend
336               zcof   = zfact1 * tmask(ji,jj,jk)
337               !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
338               !                                   ! eddy coefficient (ensure numerical stability)
339               zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
340                  &          /    (  p_e3t(ji,jj,jk  ) * p_e3w(ji,jj,jk  )  )
341               zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
342                  &          /    (  p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk  )  )
343               !
344               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
345               zd_lw(ji,jj,jk) = zzd_lw
346               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
347               !
348               !                                   ! right hand side in en
349               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  p_sh2(ji,jj,jk)                          &   ! shear
350                  &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification
351                  &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
352                  &                                ) * wmask(ji,jj,jk)
353            END DO
354         END DO
355      END DO
356      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
357      DO jk = 3, jpkm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
358         DO jj = k_Jstr, k_Jend
359            DO ji = k_Istr, k_Iend
360               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
361            END DO
362         END DO
363      END DO
364      DO jj = k_Jstr, k_Jend                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
365         DO ji = k_Istr, k_Iend
366            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
367         END DO
368      END DO
369      DO jk = 3, jpkm1
370         DO jj = k_Jstr, k_Jend
371            DO ji = k_Istr, k_Iend
372               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)
373            END DO
374         END DO
375      END DO
376      DO jj = k_Jstr, k_Jend                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
377         DO ji = k_Istr, k_Iend
378            en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
379         END DO
380      END DO
381      DO jk = jpk-2, 2, -1
382         DO jj = k_Jstr, k_Jend
383            DO ji = k_Istr, k_Iend
384               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
385            END DO
386         END DO
387      END DO
388      DO jk = 2, jpkm1                             ! set the minimum value of tke
389         DO jj = k_Jstr, k_Jend
390            DO ji = k_Istr, k_Iend
391               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
392            END DO
393         END DO
394      END DO
395
396      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
397      !                            !  TKE due to surface and internal wave breaking
398      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
399!!gm BUG : in the exp  remove the depth of ssh !!!
400     
401     
402      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
403         DO jk = 2, jpkm1
404            DO jj = k_Jstr, k_Jend
405               DO ji = k_Istr, k_Iend
406                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
407                     &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
408               END DO
409            END DO
410         END DO
411      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
412         DO jj = k_Jstr, k_Jend
413            DO ji = k_Istr, k_Iend
414               jk = nmln(ji,jj)
415               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
416                  &                                 * MAX(0.,1._wp - 2.*fr_i(ji,jj) )  * wmask(ji,jj,jk) * tmask(ji,jj,1)
417            END DO
418         END DO
419      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
420         DO jk = 2, jpkm1
421            DO jj = k_Jstr, k_Jend
422               DO ji = k_Istr, k_Iend
423                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
424                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
425                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
426                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
427                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
428                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -pdepw(ji,jj,jk) / htau(ji,jj) )   &
429                     &                        * MAX(0.,1._wp - 2.*fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
430               END DO
431            END DO
432         END DO
433      ENDIF
434      !
435      IF( nn_timing == 1 )  CALL timing_stop('tke_tke')
436      !
437   END SUBROUTINE tke_tke
438
439
440   SUBROUTINE tke_avn( ARG_2D, pdepw, p_e3t, p_e3w, p_avm, p_avt )
441      !!----------------------------------------------------------------------
442      !!                   ***  ROUTINE tke_avn  ***
443      !!
444      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
445      !!
446      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
447      !!              the tke_tke routine). First, the now mixing lenth is
448      !!      computed from en and the strafification (N^2), then the mixings
449      !!      coefficients are computed.
450      !!              - Mixing length : a first evaluation of the mixing lengh
451      !!      scales is:
452      !!                      mxl = sqrt(2*en) / N 
453      !!      where N is the brunt-vaisala frequency, with a minimum value set
454      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
455      !!        The mixing and dissipative length scale are bound as follow :
456      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
457      !!                        zmxld = zmxlm = mxl
458      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
459      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
460      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
461      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
462      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
463      !!                    the surface to obtain ldown. the resulting length
464      !!                    scales are:
465      !!                        zmxld = sqrt( lup * ldown )
466      !!                        zmxlm = min ( lup , ldown )
467      !!              - Vertical eddy viscosity and diffusivity:
468      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
469      !!                      avt = max( avmb, pdlr * avm ) 
470      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
471      !!
472      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
473      !!----------------------------------------------------------------------
474      INTEGER                   , INTENT(in   ) ::   ARG_2D         ! inner domain start-end i-indices
475      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdepw          ! vertical eddy viscosity (w-points)
476      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   p_e3t, p_e3w     ! vertical eddy viscosity (w-points)
477      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   !  vertical eddy viscosity (w-points)
478      !
479      INTEGER  ::   ji, jj, jk   ! dummy loop indices
480      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
481      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
482      REAL(wp) ::   zemxl, zemlm, zemlp       !   -      -
483      REAL(wp), DIMENSION(WRK_3D) ::   zmxlm, zmxld
484      !!--------------------------------------------------------------------
485      !
486      IF( nn_timing == 1 )  CALL timing_start('tke_avn')
487
488      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
489      !                     !  Mixing length
490      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
491      !
492      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
493      !
494      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
495      zmxlm(WRK_3D)  = rmxl_min   
496      zmxld(WRK_3D)  = rmxl_min
497      !
498      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
499         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
500         DO jj = k_Jstr, k_Jend
501            DO ji = k_Istr, k_Iend
502               zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
503            END DO
504         END DO
505      ELSE
506         zmxlm(WRK_2D,1) = rn_mxl0
507      ENDIF
508      !
509      DO jk = 2, jpkm1              ! interior value : l=sqrt(2*e/n^2)
510         DO jj = k_Jstr, k_Jend
511            DO ji = k_Istr, k_Iend
512               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
513               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
514            END DO
515         END DO
516      END DO
517      !
518      !                     !* Physical limits for the mixing length
519      !
520      zmxld(WRK_2D, 1 ) = zmxlm(WRK_2D,1)   ! surface set to the minimum value
521      zmxld(WRK_2D,jpk) = rmxl_min          ! last level  set to the minimum value
522      !
523      SELECT CASE ( nn_mxl )
524      !
525 !!gm Not sure of that coding for ISF....
526      ! where wmask = 0 set zmxlm == p_e3w
527      CASE ( 0 )           ! bounded by the distance to surface and bottom
528         DO jk = 2, jpkm1
529            DO jj = k_Jstr, k_Jend
530               DO ji = k_Istr, k_Iend
531                  zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk),   &
532                  &            pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) )
533                  ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
534                  zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
535                  zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))
536               END DO
537            END DO
538         END DO
539         !
540      CASE ( 1 )           ! bounded by the vertical scale factor
541         DO jk = 2, jpkm1
542            DO jj = k_Jstr, k_Jend
543               DO ji = k_Istr, k_Iend
544                  zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) )
545                  zmxlm(ji,jj,jk) = zemxl
546                  zmxld(ji,jj,jk) = zemxl
547               END DO
548            END DO
549         END DO
550         !
551      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
552         DO jk = 2, jpkm1         ! from the surface to the bottom :
553            DO jj = k_Jstr, k_Jend
554               DO ji = k_Istr, k_Iend
555                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
556               END DO
557            END DO
558         END DO
559         DO jk = jpkm1, 2, -1     ! from the bottom to the surface :
560            DO jj = k_Jstr, k_Jend
561               DO ji = k_Istr, k_Iend
562                  zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
563                  zmxlm(ji,jj,jk) = zemxl
564                  zmxld(ji,jj,jk) = zemxl
565               END DO
566            END DO
567         END DO
568         !
569      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
570         DO jk = 2, jpkm1         ! from the surface to the bottom : lup
571            DO jj = k_Jstr, k_Jend
572               DO ji = k_Istr, k_Iend
573                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
574               END DO
575            END DO
576         END DO
577         DO jk = jpkm1, 2, -1     ! from the bottom to the surface : ldown
578            DO jj = k_Jstr, k_Jend
579               DO ji = k_Istr, k_Iend
580                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
581               END DO
582            END DO
583         END DO
584         DO jk = 2, jpkm1
585            DO jj = k_Jstr, k_Jend
586               DO ji = k_Istr, k_Iend
587                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
588                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
589                  zmxlm(ji,jj,jk) = zemlm
590                  zmxld(ji,jj,jk) = zemlp
591               END DO
592            END DO
593         END DO
594         !
595      END SELECT
596      !
597
598      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
599      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
600      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
601      DO jk = 1, jpkm1            !* vertical eddy viscosity & diffivity at w-points
602         DO jj = k_Jstr, k_Jend
603            DO ji = k_Istr, k_Iend
604               zsqen = SQRT( en(ji,jj,jk) )
605               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
606               p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
607               p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
608               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
609            END DO
610         END DO
611      END DO
612      !
613      !
614      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
615         DO jk = 2, jpkm1
616            DO jj = k_Jstr, k_Jend
617               DO ji = k_Istr, k_Iend
618                  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)
619              END DO
620            END DO
621         END DO
622      ENDIF
623
624      IF(ln_ctl) THEN
625         CALL prt_ctl( tab3d_1=en , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
626         CALL prt_ctl( tab3d_1=avm, clinfo1=' tke  - m: ', ovlap=1, kdim=jpk )
627      ENDIF
628      !
629      IF( nn_timing == 1 )  CALL timing_stop('tke_avn')
630      !
631   END SUBROUTINE tke_avn
632
633
634   SUBROUTINE zdf_tke_init
635      !!----------------------------------------------------------------------
636      !!                  ***  ROUTINE zdf_tke_init  ***
637      !!                     
638      !! ** Purpose :   Initialization of the vertical eddy diffivity and
639      !!              viscosity when using a tke turbulent closure scheme
640      !!
641      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
642      !!              called at the first timestep (nit000)
643      !!
644      !! ** input   :   Namlist namzdf_tke
645      !!
646      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
647      !!----------------------------------------------------------------------
648      INTEGER ::   ji, jj, jk   ! dummy loop indices
649      INTEGER ::   ios
650      !!
651      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
652         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
653         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
654         &                 nn_etau , nn_htau  , rn_efr   
655      !!----------------------------------------------------------------------
656      !
657      REWIND( numnam_ref )              ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
658      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
659901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
660
661      REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
662      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
663902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
664      IF(lwm) WRITE ( numond, namzdf_tke )
665      !
666      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
667      !
668      IF(lwp) THEN                    !* Control print
669         WRITE(numout,*)
670         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
671         WRITE(numout,*) '~~~~~~~~~~~~'
672         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
673         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
674         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
675         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
676         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
677         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
678         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
679         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
680         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
681         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
682         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
683         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
684         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
685         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
686         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
687         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
688         WRITE(numout,*)
689         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
690         WRITE(numout,*)
691      ENDIF
692      !
693      IF( ln_zdfiwm ) THEN          ! Internal wave-driven mixing
694         rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
695         rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
696         IF(lwp) WRITE(numout,*) '      Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
697      ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
698         rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
699         IF(lwp) WRITE(numout,*) '      minimum mixing length with your parameters rmxl_min = ', rmxl_min
700      ENDIF
701      !
702      !                              ! allocate tke arrays
703      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
704      !
705      !                               !* Check of some namelist values
706      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
707      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
708      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
709      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
710
711      IF( ln_mxl0 ) THEN
712         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
713         rn_mxl0 = rmxl_min
714      ENDIF
715     
716      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
717
718      !                               !* depth of penetration of surface tke
719      IF( nn_etau /= 0 ) THEN     
720         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
721         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
722            htau(:,:) = 10._wp
723         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
724            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
725         END SELECT
726      ENDIF
727      !                               !* set vertical eddy coef. to the background value
728      DO jk = 1, jpk
729         avt(:,:,jk) = avtb(jk) * wmask(:,:,jk)
730         avm(:,:,jk) = avmb(jk) * wmask(:,:,jk)
731      END DO
732      dissl(:,:,:) = 1.e-12_wp
733      !                             
734      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
735      !
736   END SUBROUTINE zdf_tke_init
737
738
739   SUBROUTINE tke_rst( kt, cdrw )
740      !!---------------------------------------------------------------------
741      !!                   ***  ROUTINE tke_rst  ***
742      !!                     
743      !! ** Purpose :   Read or write TKE file (en) in restart file
744      !!
745      !! ** Method  :   use of IOM library
746      !!                if the restart does not contain TKE, en is either
747      !!                set to rn_emin or recomputed
748      !!----------------------------------------------------------------------
749      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
750      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
751      !
752      INTEGER ::   jit, jk              ! dummy loop indices
753      INTEGER ::   id1, id2, id3, id4   ! local integers
754      !!----------------------------------------------------------------------
755      !
756      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
757         !                                   ! ---------------
758         IF( ln_rstart ) THEN                   !* Read the restart file
759            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
760            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
761            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
762            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
763            !
764            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
765               CALL iom_get( numror, jpdom_autoglo, 'en', en )
766               CALL iom_get( numror, jpdom_autoglo, 'avt_k', avt_k )
767               CALL iom_get( numror, jpdom_autoglo, 'avm_k', avm_k )
768               CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
769            ELSE                                          ! start TKE from rest
770               IF(lwp) WRITE(numout,*) '   ==>>   previous run without TKE scheme, set en to background values'
771               en(:,:,:) = rn_emin * wmask(:,:,:)
772               ! avt_k, avm_k already set to the background value in zdf_phy_init
773            ENDIF
774         ELSE                                   !* Start from rest
775            IF(lwp) WRITE(numout,*) '   ==>>   start from rest: set en to the background value'
776            en(:,:,:) = rn_emin * wmask(:,:,:)
777            ! avt_k, avm_k already set to the background value in zdf_phy_init
778         ENDIF
779         !
780      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
781         !                                   ! -------------------
782         IF(lwp) WRITE(numout,*) '---- tke-rst ----'
783         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
784         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
785         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
786         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
787         !
788      ENDIF
789      !
790   END SUBROUTINE tke_rst
791
792   !!======================================================================
793END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.