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

source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 13191

Last change on this file since 13191 was 13191, checked in by jwhile, 4 years ago

Updates for 1d runnig
File size: 53.4 KB
Line 
1MODULE zdftke
2   !!======================================================================
3   !!                       ***  MODULE  zdftke  ***
4   !! Ocean physics:  vertical mixing coefficient computed from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !! History :  OPA  !  1991-03  (b. blanke)  Original code
8   !!            7.0  !  1991-11  (G. Madec)   bug fix
9   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
10   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
11   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
12   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
13   !!            8.0  !  1997-07  (G. Madec)   lbc
14   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
15   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
16   !!             -   !  2004-10  (C. Ethe )  1D configuration
17   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
18   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
19   !!                 !           - tke penetration (wind steering)
20   !!                 !           - suface condition for tke & mixing length
21   !!                 !           - Langmuir cells
22   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
23   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
24   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
25   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
26   !!                 !                                + cleaning of the parameters + bugs correction
27   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
28   !!            3.6  !  2014-11  (P. Mathiot) add ice shelf capability
29   !!----------------------------------------------------------------------
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) 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      ENDIF
850# else
851      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
852# endif
853      !
854      IF(lwp) THEN                    !* Control print
855         WRITE(numout,*)
856         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
857         WRITE(numout,*) '~~~~~~~~~~~~'
858         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
859         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
860         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
861         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
862         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
863         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
864         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
865         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
866         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
867         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
868         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
869         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
870         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
871         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
872         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
873         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
874         WRITE(numout,*) '      fraction of TKE added within the mixed layer by nn_etau rn_c    = ', rn_c
875         WRITE(numout,*)
876         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
877      ENDIF
878      !
879      !                              ! allocate tke arrays
880      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
881      !
882      !                               !* Check of some namelist values
883      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
884      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
885      IF( nn_htau < 0  .OR.  nn_htau > 5 )   CALL ctl_stop( 'bad flag: nn_htau is 0 to 5    ' )
886      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
887
888      IF( ln_mxl0 ) THEN
889         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
890         rn_mxl0 = rmxl_min
891      ENDIF
892
893      ALLOCATE( rn_lc0   (jpi,jpj) ) ; rn_lc0    = rn_lc
894      ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff
895      ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss
896      ALLOCATE( rn_ebb0  (jpi,jpj) ) ; rn_ebb0   = rn_ebb
897      ALLOCATE( rn_efr0  (jpi,jpj) ) ; rn_efr0   = rn_efr
898
899      IF( nn_etau == 2  ) THEN
900          ierr = zdf_mxl_alloc()
901          nmln(:,:) = nlb10           ! Initialization of nmln
902      ENDIF
903
904      IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN
905          ierr = zdf_mxl_alloc()
906          nmln(:,:) = nlb10           ! Initialization of nmln
907      ENDIF
908
909      !                               !* depth of penetration of surface tke
910      IF( nn_etau /= 0 ) THEN
911         htau(:,:) = 0._wp
912         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
913         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
914            htau(:,:) = 10._wp
915         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
916            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
917         CASE( 2 )                                 ! fraction of depth-integrated TKE within mixed-layer
918            rhtau = -1._wp / LOG( 1._wp - rn_c )
919         CASE( 3 )                                 ! F(latitude) : 0.5m to 15m poleward of 20 degrees
920            htau(:,:) = MAX(  0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
921         CASE( 4 )                                 ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south
922            DO jj = 2, jpjm1
923               DO ji = fs_2, fs_jpim1   ! vector opt.
924                  IF( gphit(ji,jj) <= 0._wp ) THEN
925                     htau(ji,jj) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
926                  ELSE
927                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
928                  ENDIF
929               END DO
930            END DO
931         CASE ( 5 )                                ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south,
932            DO jj = 2, jpjm1                       !               10m to 30m between 30/45 degrees south
933               DO ji = fs_2, fs_jpim1   ! vector opt.
934                  IF( gphit(ji,jj) <= -30._wp ) THEN
935                     htau(ji,jj) = MAX(  10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) )   )
936                  ELSE
937                     htau(ji,jj) = MAX(  0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) )   )
938                  ENDIF
939               END DO
940            END DO
941         END SELECT
942         !
943         IF( nn_etau == 4 .AND. nn_htau /= 2 ) THEN            ! efr dependant on constant htau
944            DO jj = 2, jpjm1
945               DO ji = fs_2, fs_jpim1   ! vector opt.
946                  efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) )
947               END DO
948            END DO
949         ENDIF
950      ENDIF
951      !                               !* set vertical eddy coef. to the background value
952      DO jk = 1, jpk
953         avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)
954         avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)
955         avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)
956         avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)
957      END DO
958      dissl(:,:,:) = 1.e-12_wp
959      !
960      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
961      !
962   END SUBROUTINE zdf_tke_init
963
964
965   SUBROUTINE tke_rst( kt, cdrw )
966     !!---------------------------------------------------------------------
967     !!                   ***  ROUTINE tke_rst  ***
968     !!
969     !! ** Purpose :   Read or write TKE file (en) in restart file
970     !!
971     !! ** Method  :   use of IOM library
972     !!                if the restart does not contain TKE, en is either
973     !!                set to rn_emin or recomputed
974     !!----------------------------------------------------------------------
975     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
976     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
977     !
978     INTEGER ::   jit, jk   ! dummy loop indices
979     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
980     !!----------------------------------------------------------------------
981     !
982     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
983        !                                   ! ---------------
984        IF( ln_rstart ) THEN                   !* Read the restart file
985           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
986           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
987           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
988           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
989           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
990           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
991           !
992           IF( id1 > 0 ) THEN                       ! 'en' exists
993              CALL iom_get( numror, jpdom_autoglo, 'en', en )
994              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
995                 IF(nn_timing == 2)  CALL timing_start('iom_rstget')
996                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
997                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
998                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
999                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
1000                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
1001                 IF(nn_timing == 2)  CALL timing_stop('iom_rstget')
1002              ELSE                                                 ! one at least array is missing
1003                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
1004              ENDIF
1005           ELSE                                     ! No TKE array found: initialisation
1006              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
1007              en (:,:,:) = rn_emin * tmask(:,:,:)
1008              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
1009              !
1010              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
1011           ENDIF
1012        ELSE                                   !* Start from rest
1013           en(:,:,:) = rn_emin * tmask(:,:,:)
1014        ENDIF
1015        !
1016        avt_k (:,:,:) = avt (:,:,:)
1017        avm_k (:,:,:) = avm (:,:,:)
1018        avmu_k(:,:,:) = avmu(:,:,:)
1019        avmv_k(:,:,:) = avmv(:,:,:)
1020
1021     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1022        !                                   ! -------------------
1023        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
1024        IF(nn_timing == 2)  CALL timing_start('iom_rstput')
1025        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en     )
1026        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt_k  )
1027        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm_k  )
1028        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
1029        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
1030        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl  )
1031        IF(nn_timing == 2)  CALL timing_stop('iom_rstput')
1032        !
1033     ENDIF
1034     !
1035   END SUBROUTINE tke_rst
1036
1037#else
1038   !!----------------------------------------------------------------------
1039   !!   Dummy module :                                        NO TKE scheme
1040   !!----------------------------------------------------------------------
1041   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
1042CONTAINS
1043   SUBROUTINE zdf_tke_init           ! Dummy routine
1044   END SUBROUTINE zdf_tke_init
1045   SUBROUTINE zdf_tke( kt )          ! Dummy routine
1046      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
1047   END SUBROUTINE zdf_tke
1048   SUBROUTINE tke_rst( kt, cdrw )
1049     CHARACTER(len=*) ::   cdrw
1050     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
1051   END SUBROUTINE tke_rst
1052#endif
1053
1054   !!======================================================================
1055END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.