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 NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/ZDF/zdftke.F90 @ 14986

Last change on this file since 14986 was 14986, checked in by sparonuz, 3 years ago

Merge trunk -r14984:HEAD
  • Property svn:keywords set to Id
File size: 50.0 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   !!            4.0  !  2017-04  (G. Madec)  remove CPP ddm key & avm at t-point only
30   !!             -   !  2017-05  (G. Madec)  add top/bottom friction as boundary condition
31   !!            4.2  !  2020-12  (G. Madec, E. Clementi) add wave coupling
32   !                  !           following Couvelard et al., 2019
33   !!----------------------------------------------------------------------
34
35   !!----------------------------------------------------------------------
36   !!   zdf_tke       : update momentum and tracer Kz from a tke scheme
37   !!   tke_tke       : tke time stepping: update tke at now time step (en)
38   !!   tke_avn       : compute mixing length scale and deduce avm and avt
39   !!   zdf_tke_init  : initialization, namelist read, and parameters control
40   !!   tke_rst       : read/write tke restart in ocean restart file
41   !!----------------------------------------------------------------------
42   USE oce            ! ocean: dynamics and active tracers variables
43   USE phycst         ! physical constants
44   USE dom_oce        ! domain: ocean
45   USE domvvl         ! domain: variable volume layer
46   USE sbc_oce        ! surface boundary condition: ocean
47   USE zdfdrg         ! vertical physics: top/bottom drag coef.
48   USE zdfmxl         ! vertical physics: mixed layer
49#if defined key_si3
50   USE ice, ONLY: hm_i, h_i
51#endif
52#if defined key_cice
53   USE sbc_ice, ONLY: h_i
54#endif
55   !
56   USE in_out_manager ! I/O manager
57   USE iom            ! I/O manager library
58   USE lib_mpp        ! MPP library
59   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
60   USE prtctl         ! Print control
61   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
62   USE sbcwave        ! Surface boundary waves
63
64   IMPLICIT NONE
65   PRIVATE
66
67   PUBLIC   zdf_tke        ! routine called in step module
68   PUBLIC   zdf_tke_init   ! routine called in opa module
69   PUBLIC   tke_rst        ! routine called in step module
70
71   !                      !!** Namelist  namzdf_tke  **
72   LOGICAL  ::   ln_mxl0   ! mixing length scale surface value as function of wind stress or not
73   LOGICAL  ::   ln_mxhsw  ! mixing length scale surface value as a fonction of wave height
74   INTEGER  ::   nn_mxlice ! type of scaling under sea-ice (=0/1/2/3)
75   REAL(wp) ::   rn_mxlice ! ice thickness value when scaling under sea-ice
76   INTEGER  ::   nn_mxl    ! type of mixing length (=0/1/2/3)
77   REAL(wp) ::   rn_mxl0   ! surface  min value of mixing length (kappa*z_o=0.4*0.1 m)  [m]
78   INTEGER  ::   nn_pdl    ! Prandtl number or not (ratio avt/avm) (=0/1)
79   REAL(wp) ::   rn_ediff  ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e)
80   REAL(wp) ::   rn_ediss  ! coefficient of the Kolmogoroff dissipation
81   REAL(wp) ::   rn_ebb    ! coefficient of the surface input of tke
82   REAL(wp) ::   rn_emin   ! minimum value of tke           [m2/s2]
83   REAL(wp) ::   rn_emin0  ! surface minimum value of tke   [m2/s2]
84   REAL(wp) ::   rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
85   INTEGER  ::   nn_etau   ! type of depth penetration of surface tke (=0/1/2/3)
86   INTEGER  ::      nn_htau   ! type of tke profile of penetration (=0/1)
87   INTEGER  ::   nn_bc_surf! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
88   INTEGER  ::   nn_bc_bot ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling
89   REAL(wp) ::      rn_efr    ! fraction of TKE surface value which penetrates in the ocean
90   LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
91   REAL(wp) ::      rn_lc     ! coef to compute vertical velocity of Langmuir cells
92   INTEGER  ::   nn_eice   ! attenutaion of langmuir & surface wave breaking under ice (=0/1/2/3)
93
94   REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
95   REAL(wp) ::   rmxl_min  ! minimum mixing length value (deduced from rn_ediff and rn_emin values)  [m]
96   REAL(wp) ::   rhftau_add = 1.e-3_wp     ! add offset   applied to HF part of taum  (nn_etau=3)
97   REAL(wp) ::   rhftau_scl = 1.0_wp       ! scale factor applied to HF part of taum  (nn_etau=3)
98
99   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau    ! depth of tke penetration (nn_htau)
100   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl   ! now mixing lenght of dissipation
101   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   apdlr   ! now mixing lenght of dissipation
102
103   !! * Substitutions
104#  include "do_loop_substitute.h90"
105#  include "domzgr_substitute.h90"
106   !!----------------------------------------------------------------------
107   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
108   !! $Id$
109   !! Software governed by the CeCILL license (see ./LICENSE)
110   !!----------------------------------------------------------------------
111CONTAINS
112
113   INTEGER FUNCTION zdf_tke_alloc()
114      !!----------------------------------------------------------------------
115      !!                ***  FUNCTION zdf_tke_alloc  ***
116      !!----------------------------------------------------------------------
117      ALLOCATE( htau(jpi,jpj) , dissl(jpi,jpj,jpk) , apdlr(jpi,jpj,jpk) ,   STAT= zdf_tke_alloc )
118      !
119      CALL mpp_sum ( 'zdftke', zdf_tke_alloc )
120      IF( zdf_tke_alloc /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_alloc: failed to allocate arrays' )
121      !
122   END FUNCTION zdf_tke_alloc
123
124
125   SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt )
126      !!----------------------------------------------------------------------
127      !!                   ***  ROUTINE zdf_tke  ***
128      !!
129      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
130      !!              coefficients using a turbulent closure scheme (TKE).
131      !!
132      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
133      !!              is computed from a prognostic equation :
134      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
135      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
136      !!                  + avt N^2                      ! stratif. destruc.
137      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
138      !!      with the boundary conditions:
139      !!         surface: en = max( rn_emin0, rn_ebb * taum )
140      !!         bottom : en = rn_emin
141      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
142      !!
143      !!        The now Turbulent kinetic energy is computed using the following
144      !!      time stepping: implicit for vertical diffusion term, linearized semi
145      !!      implicit for kolmogoroff dissipation term, and explicit forward for
146      !!      both buoyancy and shear production terms. Therefore a tridiagonal
147      !!      linear system is solved. Note that buoyancy and shear terms are
148      !!      discretized in a energy conserving form (Bruchard 2002).
149      !!
150      !!        The dissipative and mixing length scale are computed from en and
151      !!      the stratification (see tke_avn)
152      !!
153      !!        The now vertical eddy vicosity and diffusivity coefficients are
154      !!      given by:
155      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
156      !!              avt = max( avmb, pdl * avm                 )
157      !!              eav = max( avmb, avm )
158      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
159      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
160      !!
161      !! ** Action  :   compute en (now turbulent kinetic energy)
162      !!                update avt, avm (before vertical eddy coef.)
163      !!
164      !! References : Gaspar et al., JGR, 1990,
165      !!              Blanke and Delecluse, JPO, 1991
166      !!              Mellor and Blumberg, JPO 2004
167      !!              Axell, JGR, 2002
168      !!              Bruchard OM 2002
169      !!----------------------------------------------------------------------
170      INTEGER                             , INTENT(in   ) ::   kt             ! ocean time step
171      INTEGER                             , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
172      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in   ) ::   p_sh2          ! shear production term
173      REAL(wp), DIMENSION(:,:,:)          , INTENT(inout) ::   p_avm, p_avt   !  momentum and tracer Kz (w-points)
174      !!----------------------------------------------------------------------
175      !
176      CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )   ! now tke (en)
177      !
178      CALL tke_avn( Kbb, Kmm,        p_avm, p_avt )   ! now avt, avm, dissl
179      !
180  END SUBROUTINE zdf_tke
181
182
183   SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt )
184      !!----------------------------------------------------------------------
185      !!                   ***  ROUTINE tke_tke  ***
186      !!
187      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
188      !!
189      !! ** Method  : - TKE surface boundary condition
190      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
191      !!              - source term due to shear (= Kz dz[Ub] * dz[Un] )
192      !!              - Now TKE : resolution of the TKE equation by inverting
193      !!                a tridiagonal linear system by a "methode de chasse"
194      !!              - increase TKE due to surface and internal wave breaking
195      !!             NB: when sea-ice is present, both LC parameterization
196      !!                 and TKE penetration are turned off when the ice fraction
197      !!                 is smaller than 0.25
198      !!
199      !! ** Action  : - en : now turbulent kinetic energy)
200      !! ---------------------------------------------------------------------
201      USE zdf_oce , ONLY : en   ! ocean vertical physics
202      !!
203      INTEGER                              , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
204      REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(in   ) ::   p_sh2          ! shear production term
205      REAL(wp), DIMENSION(:,:,:)           , INTENT(in   ) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
206      !
207      INTEGER ::   ji, jj, jk                  ! dummy loop arguments
208      REAL(wp) ::   zetop, zebot, zmsku, zmskv ! local scalars
209      REAL(wp) ::   zrhoa  = 1.22              ! Air density kg/m3
210      REAL(wp) ::   zcdrag = 1.5e-3            ! drag coefficient
211      REAL(wp) ::   zbbrau, zbbirau, zri       ! local scalars
212      REAL(wp) ::   zfact1, zfact2, zfact3     !   -      -
213      REAL(wp) ::   ztx2  , zty2  , zcof       !   -      -
214      REAL(wp) ::   ztau  , zdif               !   -      -
215      REAL(wp) ::   zus   , zwlc  , zind       !   -      -
216      REAL(wp) ::   zzd_up, zzd_lw             !   -      -
217      REAL(wp) ::   ztaui, ztauj, z1_norm
218      INTEGER , DIMENSION(A2D(nn_hls))     ::   imlc
219      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zice_fra, zhlc, zus3, zWlc2
220      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zpelc, zdiag, zd_up, zd_lw
221      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztmp ! for diags
222      !!--------------------------------------------------------------------
223      !
224      zbbrau  = rn_ebb / rho0       ! Local constant initialisation
225      zbbirau = 3.75_wp / rho0
226      zfact1  = -.5_wp * rn_Dt
227      zfact2  = 1.5_wp * rn_Dt * rn_ediss
228      zfact3  = 0.5_wp         * rn_ediss
229      !
230      zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used
231      !
232      ! ice fraction considered for attenuation of langmuir & wave breaking
233      SELECT CASE ( nn_eice )
234      CASE( 0 )   ;   zice_fra(:,:) = 0._wp
235      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(A2D(nn_hls)) * 10._wp )
236      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(A2D(nn_hls))
237      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(A2D(nn_hls)) , 1._wp )
238      END SELECT
239      !
240      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
241      !                     !  Surface/top/bottom boundary condition on tke
242      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
243      !
244      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
245         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) )
246         zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
247         zd_lw(ji,jj,1) = 1._wp
248         zd_up(ji,jj,1) = 0._wp
249      END_2D
250      !
251      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
252      !                     !  Bottom boundary condition on tke
253      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
254      !
255      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
256      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
257      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
258      !
259      IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE
260         !
261         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )        ! bottom friction
262            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
263            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
264            !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
265            zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2  &
266               &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
267            en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
268         END_2D
269         IF( ln_isfcav ) THEN
270            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )     ! top friction
271               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
272               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
273               !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
274               zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT(  ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2  &
275                  &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
276               ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present
277               en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &
278                  &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
279            END_2D
280         ENDIF
281         !
282      ENDIF
283      !
284      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
285      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke (Axell JGR 2002)
286         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
287         !
288         !                       !* Langmuir velocity scale
289         !
290         IF ( cpl_sdrftx )  THEN       ! Surface Stokes Drift available
291            !                                ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2    with u* = (taum/rho0)^1/2
292            !                                ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s
293            !                                ! more precisely, it is the dot product that must be used :
294            !                                !     1/2  (W_lc)^2 = MAX( u* u_s + v* v_s , 0 )   only the positive part
295!!gm  ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
296!!gm  ! so we will overestimate the LC velocity....   !!gm I will do the work if !LC have an effect !
297            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
298!!XC                  zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )  )
299                  zWlc2(ji,jj) = 0.5_wp *  ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )
300            END_2D
301!
302!  Projection of Stokes drift in the wind stress direction
303!
304            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
305                  ztaui   = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
306                  ztauj   = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
307                  z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1)
308                  zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
309            END_2D
310         ELSE                          ! Surface Stokes drift deduced from surface stress
311            !                                ! Wlc = u_s   with u_s = 0.016*U_10m, the surface stokes drift  (Axell 2002, Eq.44)
312            !                                ! using |tau| = rho_air Cd |U_10m|^2 , it comes:
313            !                                ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2   and thus:
314            !                                ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag )
315            zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )      ! to convert stress in 10m wind using a constant drag
316            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
317               zWlc2(ji,jj) = zcof * taum(ji,jj)
318            END_2D
319            !
320         ENDIF
321         !
322         !                       !* Depth of the LC circulation  (Axell 2002, Eq.47)
323         !                             !- LHS of Eq.47
324         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
325            zpelc(ji,jj,1) =  MAX( rn2b(ji,jj,1), 0._wp ) * gdepw(ji,jj,1,Kmm) * e3w(ji,jj,1,Kmm)
326         END_2D
327         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpk )
328            zpelc(ji,jj,jk)  = zpelc(ji,jj,jk-1) +   &
329               &          MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
330         END_3D
331         !
332         !                             !- compare LHS to RHS of Eq.47
333         imlc(:,:) = mbkt(A2D(nn_hls)) + 1       ! Initialization to the number of w ocean point (=2 over land)
334         DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )
335            IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) )   imlc(ji,jj) = jk
336         END_3D
337         !                               ! finite LC depth
338         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
339            zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
340         END_2D
341         !
342         zcof = 0.016 / SQRT( zrhoa * zcdrag )
343         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
344            zus = SQRT( 2. * zWlc2(ji,jj) )             ! Stokes drift
345            zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
346         END_2D
347         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en
348            IF ( zus3(ji,jj) /= 0._wp ) THEN
349               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
350                  !                                           ! vertical velocity due to LC
351                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )
352                  !                                           ! TKE Langmuir circulation source term
353                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
354               ENDIF
355            ENDIF
356         END_3D
357         !
358      ENDIF
359      !
360      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
361      !                     !  Now Turbulent kinetic energy (output in en)
362      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
363      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
364      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
365      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
366      !
367      IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri )
368         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
369            !                             ! local Richardson number
370            IF (rn2b(ji,jj,jk) <= 0.0_wp) then
371                zri = 0.0_wp
372            ELSE
373                zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
374            ENDIF
375            !                             ! inverse of Prandtl number
376            apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
377         END_3D
378      ENDIF
379      !
380      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )   !* Matrix and right hand side in en
381         zcof   = zfact1 * tmask(ji,jj,jk)
382         !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
383         !                                   ! eddy coefficient (ensure numerical stability)
384         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
385            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
386         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
387            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
388         !
389         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
390         zd_lw(ji,jj,jk) = zzd_lw
391         zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
392         !
393         !                                   ! right hand side in en
394         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * (  p_sh2(ji,jj,jk)                        &   ! shear
395            &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification
396            &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
397            &                                ) * wmask(ji,jj,jk)
398      END_3D
399      !
400      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
401      !                     !  Surface boundary condition on tke if
402      !                     !  coupling with waves
403      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
404      !
405      IF ( cpl_phioc .and. ln_phioc )  THEN
406         SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
407
408         CASE ( 0 ) ! Dirichlet BC
409            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
410               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
411               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
412               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
413            END_2D
414
415         CASE ( 1 ) ! Neumann BC
416            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
417               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
418               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
419               en(ji,jj,1)    = en(ji,jj,2) + (2 * e3t(ji,jj,1,Kmm) * phioc(ji,jj)/rho0) / ( p_avm(ji,jj,1) + p_avm(ji,jj,2) )
420               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
421               zdiag(ji,jj,1) = 1._wp
422               zd_lw(ji,jj,2) = 0._wp
423            END_2D
424
425         END SELECT
426
427      ENDIF
428      !
429      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
430      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
431         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
432      END_3D
433!XC : commented to allow for neumann boundary condition
434!      DO_2D( 0, 0, 0, 0 )
435!         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
436!      END_2D
437      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
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_3D
440      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
441         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
442      END_2D
443      DO_3DS_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpk-2, 2, -1 )
444         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
445      END_3D
446      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )                ! set the minimum value of tke
447         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
448      END_3D
449      !
450      ! Kolmogorov energy of dissipation (W/kg)
451      !    ediss = Ce*sqrt(en)/L*en
452      !    dissl = sqrt(en)/L
453      IF( iom_use('ediss_k') ) THEN
454         ALLOCATE( ztmp(A2D(nn_hls),jpk) )
455         ztmp(:,:,:) = zfact3 * dissl * en * wmask
456         CALL lbc_lnk( 'zdftke', ztmp, 'W', 1._wp )
457         CALL iom_put( 'ediss_k', ztmp )
458         DEALLOCATE( ztmp )
459      ENDIF
460      !
461      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
462      !                            !  TKE due to surface and internal wave breaking
463      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
464!!gm BUG : in the exp  remove the depth of ssh !!!
465!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
466      !
467      ! penetration is partly switched off below sea-ice if nn_eice/=0
468      !
469      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
470         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
471            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
472               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
473         END_3D
474      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
475         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
476            jk = nmln(ji,jj)
477            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
478               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
479         END_2D
480      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
481         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
482            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
483            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
484            ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
485            zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
486            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
487            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
488               &                        * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
489         END_3D
490      ENDIF
491      !
492   END SUBROUTINE tke_tke
493
494
495   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt )
496      !!----------------------------------------------------------------------
497      !!                   ***  ROUTINE tke_avn  ***
498      !!
499      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
500      !!
501      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
502      !!              the tke_tke routine). First, the now mixing lenth is
503      !!      computed from en and the strafification (N^2), then the mixings
504      !!      coefficients are computed.
505      !!              - Mixing length : a first evaluation of the mixing lengh
506      !!      scales is:
507      !!                      mxl = sqrt(2*en) / N
508      !!      where N is the brunt-vaisala frequency, with a minimum value set
509      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
510      !!        The mixing and dissipative length scale are bound as follow :
511      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
512      !!                        zmxld = zmxlm = mxl
513      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
514      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
515      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
516      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
517      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
518      !!                    the surface to obtain ldown. the resulting length
519      !!                    scales are:
520      !!                        zmxld = sqrt( lup * ldown )
521      !!                        zmxlm = min ( lup , ldown )
522      !!              - Vertical eddy viscosity and diffusivity:
523      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
524      !!                      avt = max( avmb, pdlr * avm )
525      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
526      !!
527      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
528      !!----------------------------------------------------------------------
529      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics
530      !!
531      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
532      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
533      !
534      INTEGER  ::   ji, jj, jk   ! dummy loop indices
535      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
536      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
537      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      -
538      REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zmxlm, zmxld   ! 3D workspace
539      !!--------------------------------------------------------------------
540      !
541      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
542      !                     !  Mixing length
543      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
544      !
545      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
546      !
547      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
548      zmxlm(:,:,:)  = rmxl_min
549      zmxld(:,:,:)  = rmxl_min
550      !
551      IF(ln_sdw .AND. ln_mxhsw) THEN
552         zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
553         ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6
554         zcoef       = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff
555         zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
556      ELSE
557      !
558         IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
559         !
560            zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
561#if ! defined key_si3 && ! defined key_cice
562            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )                  ! No sea-ice
563               zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
564            END_2D
565#else
566            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
567            !
568            CASE( 0 )                      ! No scaling under sea-ice
569               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
570                  zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
571               END_2D
572               !
573            CASE( 1 )                      ! scaling with constant sea-ice thickness
574               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
575                  zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
576                     &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1)
577               END_2D
578               !
579            CASE( 2 )                      ! scaling with mean sea-ice thickness
580               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
581#if defined key_si3
582                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
583                     &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1)
584#elif defined key_cice
585                  zmaxice = MAXVAL( h_i(ji,jj,:) )
586                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
587                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
588#endif
589               END_2D
590               !
591            CASE( 3 )                      ! scaling with max sea-ice thickness
592               DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
593                  zmaxice = MAXVAL( h_i(ji,jj,:) )
594                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
595                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
596               END_2D
597               !
598            END SELECT
599#endif
600            !
601            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
602               zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
603            END_2D
604            !
605         ELSE
606            zmxlm(:,:,1) = rn_mxl0
607         ENDIF
608      ENDIF
609      !
610      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
611         zrn2 = MAX( rn2(ji,jj,jk), rsmall )
612         zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
613      END_3D
614      !
615      !                     !* Physical limits for the mixing length
616      !
617      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
618      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
619      !
620      SELECT CASE ( nn_mxl )
621      !
622 !!gm Not sure of that coding for ISF....
623      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
624      CASE ( 0 )           ! bounded by the distance to surface and bottom
625         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
626            zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   &
627            &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
628            ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
629            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
630               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
631            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
632               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
633         END_3D
634         !
635      CASE ( 1 )           ! bounded by the vertical scale factor
636         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
637            zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
638            zmxlm(ji,jj,jk) = zemxl
639            zmxld(ji,jj,jk) = zemxl
640         END_3D
641         !
642      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
643         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )       ! from the surface to the bottom :
644            zmxlm(ji,jj,jk) =   &
645               &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
646         END_3D
647         DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )   ! from the bottom to the surface :
648            zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
649            zmxlm(ji,jj,jk) = zemxl
650            zmxld(ji,jj,jk) = zemxl
651         END_3D
652         !
653      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
654         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )        ! from the surface to the bottom : lup
655            zmxld(ji,jj,jk) =    &
656               &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
657         END_3D
658         DO_3DS( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
659            zmxlm(ji,jj,jk) =   &
660               &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
661         END_3D
662         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
663            zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
664            zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
665            zmxlm(ji,jj,jk) = zemlm
666            zmxld(ji,jj,jk) = zemlp
667         END_3D
668         !
669      END SELECT
670      !
671      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
672      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
673      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
674      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
675         zsqen = SQRT( en(ji,jj,jk) )
676         zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
677         p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
678         p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
679         dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
680      END_3D
681      !
682      !
683      IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt
684         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
685            p_avt(ji,jj,jk)   = MAX( apdlr(ji,jj,jk) * p_avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
686         END_3D
687      ENDIF
688      !
689      IF(sn_cfctl%l_prtctl) THEN
690         CALL prt_ctl( tab3d_1=en   , clinfo1=' tke  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
691         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke  - m: ', kdim=jpk )
692      ENDIF
693      !
694   END SUBROUTINE tke_avn
695
696
697   SUBROUTINE zdf_tke_init( Kmm )
698      !!----------------------------------------------------------------------
699      !!                  ***  ROUTINE zdf_tke_init  ***
700      !!
701      !! ** Purpose :   Initialization of the vertical eddy diffivity and
702      !!              viscosity when using a tke turbulent closure scheme
703      !!
704      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
705      !!              called at the first timestep (nit000)
706      !!
707      !! ** input   :   Namlist namzdf_tke
708      !!
709      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
710      !!----------------------------------------------------------------------
711      USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag
712      !!
713      INTEGER, INTENT(in) ::   Kmm          ! time level index
714      INTEGER             ::   ji, jj, jk   ! dummy loop indices
715      INTEGER             ::   ios
716      !!
717      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  &
718         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  &
719         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
720         &                 nn_pdl  , ln_lc    , rn_lc    ,             &
721         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  ,  &
722         &                 nn_bc_surf, nn_bc_bot, ln_mxhsw
723      !!----------------------------------------------------------------------
724      !
725      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
726901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' )
727
728      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
729902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' )
730      IF(lwm) WRITE ( numond, namzdf_tke )
731      !
732      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
733      !
734      IF(lwp) THEN                    !* Control print
735         WRITE(numout,*)
736         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
737         WRITE(numout,*) '~~~~~~~~~~~~'
738         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
739         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
740         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
741         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
742         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
743         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
744         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
745         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
746         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
747         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
748         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
749         IF( ln_mxl0 ) THEN
750            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice
751            IF( nn_mxlice == 1 ) &
752            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice
753            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
754            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice'
755            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness'
756            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness'
757            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness'
758            CASE DEFAULT
759               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4')
760            END SELECT
761         ENDIF
762         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
763         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
764         IF ( cpl_phioc .and. ln_phioc )  THEN
765            SELECT CASE( nn_bc_surf)             ! Type of scaling under sea-ice
766            CASE( 0 )   ;   WRITE(numout,*) '  nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves'
767            CASE( 1 )   ;   WRITE(numout,*) '  nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves'
768            END SELECT
769         ENDIF
770         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
771         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
772         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr
773         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice
774         SELECT CASE( nn_eice )
775         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking'
776         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )'
777         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)'
778         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )'
779         CASE DEFAULT
780            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3')
781         END SELECT
782         WRITE(numout,*)
783         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
784         WRITE(numout,*)
785      ENDIF
786      !
787      IF( ln_zdfiwm ) THEN          ! Internal wave-driven mixing
788         rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
789         rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
790         IF(lwp) WRITE(numout,*) '   ==>>>   Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3'
791      ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
792         rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
793         IF(lwp) WRITE(numout,*) '   ==>>>   minimum mixing length with your parameters rmxl_min = ', rmxl_min
794      ENDIF
795      !
796      !                              ! allocate tke arrays
797      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
798      !
799      !                               !* Check of some namelist values
800      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1, 2 or 3' )
801      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1' )
802      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0 or 1' )
803      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
804      !
805      IF( ln_mxl0 ) THEN
806         IF(lwp) WRITE(numout,*)
807         IF(lwp) WRITE(numout,*) '   ==>>>   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
808         rn_mxl0 = rmxl_min
809      ENDIF
810      !                               !* depth of penetration of surface tke
811      IF( nn_etau /= 0 ) THEN
812         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
813         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
814            htau(:,:) = 10._wp
815         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
816            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )
817         END SELECT
818      ENDIF
819      !                                !* read or initialize all required files
820      CALL tke_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, dissl)
821      !
822   END SUBROUTINE zdf_tke_init
823
824
825   SUBROUTINE tke_rst( kt, cdrw )
826      !!---------------------------------------------------------------------
827      !!                   ***  ROUTINE tke_rst  ***
828      !!
829      !! ** Purpose :   Read or write TKE file (en) in restart file
830      !!
831      !! ** Method  :   use of IOM library
832      !!                if the restart does not contain TKE, en is either
833      !!                set to rn_emin or recomputed
834      !!----------------------------------------------------------------------
835      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
836      !!
837      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
838      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
839      !
840      INTEGER ::   jit, jk              ! dummy loop indices
841      INTEGER ::   id1, id2, id3, id4   ! local integers
842      !!----------------------------------------------------------------------
843      !
844      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
845         !                                   ! ---------------
846         IF( ln_rstart ) THEN                   !* Read the restart file
847            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
848            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
849            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
850            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
851            !
852            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
853               CALL iom_get( numror, jpdom_auto, 'en'   , en    )
854               CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k )
855               CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k )
856               CALL iom_get( numror, jpdom_auto, 'dissl', dissl )
857            ELSE                                          ! start TKE from rest
858               IF(lwp) WRITE(numout,*)
859               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without TKE scheme, set en to background values'
860               en   (:,:,:) = rn_emin * wmask(:,:,:)
861               dissl(:,:,:) = 1.e-12_wp
862               ! avt_k, avm_k already set to the background value in zdf_phy_init
863            ENDIF
864         ELSE                                   !* Start from rest
865            IF(lwp) WRITE(numout,*)
866            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set en to the background value'
867            en   (:,:,:) = rn_emin * wmask(:,:,:)
868            dissl(:,:,:) = 1.e-12_wp
869            ! avt_k, avm_k already set to the background value in zdf_phy_init
870         ENDIF
871         !
872      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
873         !                                   ! -------------------
874         IF(lwp) WRITE(numout,*) '---- tke_rst ----'
875         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
876         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
877         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
878         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
879         !
880      ENDIF
881      !
882   END SUBROUTINE tke_rst
883
884   !!======================================================================
885END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.