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

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 8568

Last change on this file since 8568 was 8568, checked in by gm, 7 years ago

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