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

source: branches/UKMO/dev_r5518_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 7066

Last change on this file since 7066 was 7066, checked in by till, 8 years ago

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