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

source: branches/UKMO/dev_r5021_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 5242

Last change on this file since 5242 was 5242, checked in by davestorkey, 9 years ago

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