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

source: branches/UKMO/dev_r5518_GC3p0_package/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 6457

Last change on this file since 6457 was 6457, checked in by dancopsey, 9 years ago

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