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

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_STOPACK/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 11400

Last change on this file since 11400 was 11400, checked in by mattmartin, 5 years ago

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