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

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 9104

Last change on this file since 9104 was 9104, checked in by gm, 6 years ago

dev_merge_2017: ZDF: timing + lnk_multi + namelist cfg ctl

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