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

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 3558

Last change on this file since 3558 was 3558, checked in by rblod, 11 years ago

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