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

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 15603

Last change on this file since 15603 was 15603, checked in by mattmartin, 3 years ago

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