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/trunk/src/OCE/ZDF – NEMO

source: NEMO/trunk/src/OCE/ZDF/zdftke.F90 @ 14007

Last change on this file since 14007 was 14007, checked in by emanuelaclementi, 4 years ago

merging branch dev_r12702_ASINTER-02_emanuelaclementi_Waves

  • Property svn:keywords set to Id
File size: 48.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   !!            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(:,:,:), 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(:,:,:) , 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(jpi,jpj)     ::   imlc
219      REAL(wp), DIMENSION(jpi,jpj)     ::   zice_fra, zhlc, zus3, zWlc2
220      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zpelc, zdiag, zd_up, zd_lw
221      !!--------------------------------------------------------------------
222      !
223      zbbrau  = rn_ebb / rho0       ! Local constant initialisation
224      zbbirau = 3.75_wp / rho0
225      zfact1  = -.5_wp * rn_Dt 
226      zfact2  = 1.5_wp * rn_Dt * rn_ediss
227      zfact3  = 0.5_wp         * rn_ediss
228      !
229      zpelc(:,:,:) = 0._wp ! need to be initialised in case ln_lc is not used
230      !
231      ! ice fraction considered for attenuation of langmuir & wave breaking
232      SELECT CASE ( nn_eice )
233      CASE( 0 )   ;   zice_fra(:,:) = 0._wp
234      CASE( 1 )   ;   zice_fra(:,:) =        TANH( fr_i(:,:) * 10._wp )
235      CASE( 2 )   ;   zice_fra(:,:) =              fr_i(:,:)
236      CASE( 3 )   ;   zice_fra(:,:) = MIN( 4._wp * fr_i(:,:) , 1._wp )
237      END SELECT
238      !
239      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
240      !                     !  Surface/top/bottom boundary condition on tke
241      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
242      !
243      DO_2D( 0, 0, 0, 0 )
244         en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
245         zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)
246         zd_lw(ji,jj,1) = 1._wp 
247         zd_up(ji,jj,1) = 0._wp
248      END_2D
249      !
250      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
251      !                     !  Bottom boundary condition on tke
252      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
253      !
254      !   en(bot)   = (ebb0/rho0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
255      ! where ebb0 does not includes surface wave enhancement (i.e. ebb0=3.75)
256      ! Note that stress averaged is done using an wet-only calculation of u and v at t-point like in zdfsh2
257      !
258      IF( .NOT.ln_drg_OFF ) THEN    !== friction used as top/bottom boundary condition on TKE
259         !
260         DO_2D( 0, 0, 0, 0 )        ! bottom friction
261            zmsku = ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
262            zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
263            !                       ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0)
264            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  &
265               &                                           + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  )
266            en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj)
267         END_2D
268         IF( ln_isfcav ) THEN
269            DO_2D( 0, 0, 0, 0 )     ! top friction
270               zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) )
271               zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) )
272               !                             ! where 0.001875 = (rn_ebb0/rho0) * 0.5 = 3.75*0.5/1000.  (CAUTION CdU<0)
273               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  &
274                  &                                           + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2  )
275               ! (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj) = 1 where ice shelves are present
276               en(ji,jj,mikt(ji,jj)) = en(ji,jj,1)           * tmask(ji,jj,1) &
277                  &                  + MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) * ssmask(ji,jj)
278            END_2D
279         ENDIF
280         !
281      ENDIF
282      !
283      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
284      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke (Axell JGR 2002)
285         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
286         !
287         !                       !* Langmuir velocity scale
288         !
289         IF ( cpl_sdrftx )  THEN       ! Surface Stokes Drift available
290            !                                ! Craik-Leibovich velocity scale Wlc = ( u* u_s )^1/2    with u* = (taum/rho0)^1/2
291            !                                ! associated kinetic energy : 1/2 (Wlc)^2 = u* u_s
292            !                                ! more precisely, it is the dot product that must be used :
293            !                                !     1/2  (W_lc)^2 = MAX( u* u_s + v* v_s , 0 )   only the positive part
294!!gm  ! PS: currently we don't have neither the 2 stress components at t-point !nor the angle between u* and u_s
295!!gm  ! so we will overestimate the LC velocity....   !!gm I will do the work if !LC have an effect !
296            DO_2D( 0, 0, 0, 0 )
297!!XC                  zWlc2(ji,jj) = 0.5_wp * SQRT( taum(ji,jj) * r1_rho0 * ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )  )
298                  zWlc2(ji,jj) = 0.5_wp *  ( ut0sd(ji,jj)**2 +vt0sd(ji,jj)**2 )
299            END_2D
300!
301!  Projection of Stokes drift in the wind stress direction
302!
303            DO_2D( 0, 0, 0, 0 )
304                  ztaui   = 0.5_wp * ( utau(ji,jj) + utau(ji-1,jj) )
305                  ztauj   = 0.5_wp * ( vtau(ji,jj) + vtau(ji,jj-1) )
306                  z1_norm = 1._wp / MAX( SQRT(ztaui*ztaui+ztauj*ztauj), 1.e-12 ) * tmask(ji,jj,1)
307                  zWlc2(ji,jj) = 0.5_wp * z1_norm * ( MAX( ut0sd(ji,jj)*ztaui + vt0sd(ji,jj)*ztauj, 0._wp ) )**2
308            END_2D
309         CALL lbc_lnk      ( 'zdftke', zWlc2, 'T', 1. )
310!
311         ELSE                          ! Surface Stokes drift deduced from surface stress
312            !                                ! Wlc = u_s   with u_s = 0.016*U_10m, the surface stokes drift  (Axell 2002, Eq.44)
313            !                                ! using |tau| = rho_air Cd |U_10m|^2 , it comes:
314            !                                ! Wlc = 0.016 * [|tau|/(rho_air Cdrag) ]^1/2   and thus:
315            !                                ! 1/2 Wlc^2 = 0.5 * 0.016 * 0.016 |tau| /( rho_air Cdrag )
316            zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )      ! to convert stress in 10m wind using a constant drag
317            DO_2D( 1, 1, 1, 1 )
318               zWlc2(ji,jj) = zcof * taum(ji,jj)
319            END_2D
320            !
321         ENDIF
322         !
323         !                       !* Depth of the LC circulation  (Axell 2002, Eq.47)
324         !                             !- LHS of Eq.47
325         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm)
326         DO jk = 2, jpk
327            zpelc(:,:,jk)  = zpelc(:,:,jk-1) +   &
328               &        MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm)
329         END DO
330         !
331         !                             !- compare LHS to RHS of Eq.47
332         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
333         DO_3DS( 1, 1, 1, 1, jpkm1, 2, -1 )
334            IF( zpelc(ji,jj,jk) > zWlc2(ji,jj) )   imlc(ji,jj) = jk
335         END_3D
336         !                               ! finite LC depth
337         DO_2D( 1, 1, 1, 1 )
338            zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm)
339         END_2D
340         !
341         zcof = 0.016 / SQRT( zrhoa * zcdrag )
342         DO_2D( 0, 0, 0, 0 )
343            zus = SQRT( 2. * zWlc2(ji,jj) )             ! Stokes drift
344            zus3(ji,jj) = MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * zus * zus * zus * tmask(ji,jj,1) ! zus > 0. ok
345         END_2D
346         DO_3D( 0, 0, 0, 0, 2, jpkm1 )                  !* TKE Langmuir circulation source term added to en
347            IF ( zus3(ji,jj) /= 0._wp ) THEN               
348               IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN
349                  !                                           ! vertical velocity due to LC
350                  zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) )
351                  !                                           ! TKE Langmuir circulation source term
352                  en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * zus3(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj)
353               ENDIF
354            ENDIF
355         END_3D
356         !
357      ENDIF
358      !
359      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
360      !                     !  Now Turbulent kinetic energy (output in en)
361      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
362      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
363      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
364      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
365      !
366      IF( nn_pdl == 1 ) THEN          !* Prandtl number = F( Ri )
367         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
368            !                             ! local Richardson number
369            IF (rn2b(ji,jj,jk) <= 0.0_wp) then
370                zri = 0.0_wp
371            ELSE
372                zri = rn2b(ji,jj,jk) * p_avm(ji,jj,jk) / ( p_sh2(ji,jj,jk) + rn_bshear )
373            ENDIF
374            !                             ! inverse of Prandtl number
375            apdlr(ji,jj,jk) = MAX(  0.1_wp,  ri_cri / MAX( ri_cri , zri )  )
376         END_3D
377      ENDIF
378      !         
379      DO_3D( 0, 0, 0, 0, 2, jpkm1 )   !* Matrix and right hand side in en
380         zcof   = zfact1 * tmask(ji,jj,jk)
381         !                                   ! A minimum of 2.e-5 m2/s is imposed on TKE vertical
382         !                                   ! eddy coefficient (ensure numerical stability)
383         zzd_up = zcof * MAX(  p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk  ) , 2.e-5_wp  )   &  ! upper diagonal
384            &          /    (  e3t(ji,jj,jk  ,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
385         zzd_lw = zcof * MAX(  p_avm(ji,jj,jk  ) + p_avm(ji,jj,jk-1) , 2.e-5_wp  )   &  ! lower diagonal
386            &          /    (  e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk  ,Kmm)  )
387         !
388         zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
389         zd_lw(ji,jj,jk) = zzd_lw
390         zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * wmask(ji,jj,jk)
391         !
392         !                                   ! right hand side in en
393         en(ji,jj,jk) = en(ji,jj,jk) + rn_Dt * (  p_sh2(ji,jj,jk)                        &   ! shear
394            &                                 - p_avt(ji,jj,jk) * rn2(ji,jj,jk)          &   ! stratification
395            &                                 + zfact3 * dissl(ji,jj,jk) * en(ji,jj,jk)  &   ! dissipation
396            &                                ) * wmask(ji,jj,jk)
397      END_3D
398      !
399      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
400      !                     !  Surface boundary condition on tke if
401      !                     !  coupling with waves
402      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
403      !
404      IF ( cpl_phioc .and. ln_phioc )  THEN
405         SELECT CASE (nn_bc_surf) ! Boundary Condition using surface TKE flux from waves
406
407         CASE ( 0 ) ! Dirichlet BC
408            DO_2D( 0, 0, 0, 0 )    ! en(1)   = rn_ebb taum / rho0  (min value rn_emin0)
409               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
410               en(ji,jj,1) = MAX( rn_emin0, .5 * ( 15.8 * phioc(ji,jj) / rho0 )**(2./3.) )  * tmask(ji,jj,1)
411               zdiag(ji,jj,1) = 1._wp/en(ji,jj,1)  ! choose to keep coherence with former estimation of
412            END_2D
413
414         CASE ( 1 ) ! Neumann BC
415            DO_2D( 0, 0, 0, 0 )
416               IF ( phioc(ji,jj) < 0 )  phioc(ji,jj) = 0._wp
417               en(ji,jj,2)    = en(ji,jj,2) + ( rn_Dt * phioc(ji,jj) / rho0 ) /e3w(ji,jj,2,Kmm)
418               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) )
419               zdiag(ji,jj,2) = zdiag(ji,jj,2) + zd_lw(ji,jj,2)
420               zdiag(ji,jj,1) = 1._wp
421               zd_lw(ji,jj,2) = 0._wp
422            END_2D
423
424         END SELECT
425
426      ENDIF
427      !
428      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
429      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
430         zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
431      END_3D
432!XC : commented to allow for neumann boundary condition
433!      DO_2D( 0, 0, 0, 0 )
434!         zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
435!      END_2D
436      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
437         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)
438      END_3D
439      DO_2D( 0, 0, 0, 0 )                          ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
440         en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
441      END_2D
442      DO_3DS( 0, 0, 0, 0, jpk-2, 2, -1 )
443         en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
444      END_3D
445      DO_3D( 0, 0, 0, 0, 2, jpkm1 )                ! set the minimum value of tke
446         en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
447      END_3D
448      !
449      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
450      !                            !  TKE due to surface and internal wave breaking
451      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
452!!gm BUG : in the exp  remove the depth of ssh !!!
453!!gm       i.e. use gde3w in argument (gdepw(:,:,:,Kmm))
454      !
455      ! penetration is partly switched off below sea-ice if nn_eice/=0
456      !
457      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
458         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
459            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
460               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
461         END_3D
462      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
463         DO_2D( 0, 0, 0, 0 )
464            jk = nmln(ji,jj)
465            en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
466               &                                 * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
467         END_2D
468      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
469         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
470            ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
471            zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
472            ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1)    ! module of the mean stress
473            zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
474            zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
475            en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) )   &
476               &                        * MAX( 0._wp, 1._wp - zice_fra(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
477         END_3D
478      ENDIF
479      !
480   END SUBROUTINE tke_tke
481
482
483   SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt )
484      !!----------------------------------------------------------------------
485      !!                   ***  ROUTINE tke_avn  ***
486      !!
487      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
488      !!
489      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
490      !!              the tke_tke routine). First, the now mixing lenth is
491      !!      computed from en and the strafification (N^2), then the mixings
492      !!      coefficients are computed.
493      !!              - Mixing length : a first evaluation of the mixing lengh
494      !!      scales is:
495      !!                      mxl = sqrt(2*en) / N 
496      !!      where N is the brunt-vaisala frequency, with a minimum value set
497      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
498      !!        The mixing and dissipative length scale are bound as follow :
499      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
500      !!                        zmxld = zmxlm = mxl
501      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
502      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
503      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
504      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
505      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
506      !!                    the surface to obtain ldown. the resulting length
507      !!                    scales are:
508      !!                        zmxld = sqrt( lup * ldown )
509      !!                        zmxlm = min ( lup , ldown )
510      !!              - Vertical eddy viscosity and diffusivity:
511      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
512      !!                      avt = max( avmb, pdlr * avm ) 
513      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
514      !!
515      !! ** Action  : - avt, avm : now vertical eddy diffusivity and viscosity (w-point)
516      !!----------------------------------------------------------------------
517      USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d   ! ocean vertical physics
518      !!
519      INTEGER                   , INTENT(in   ) ::   Kbb, Kmm       ! ocean time level indices
520      REAL(wp), DIMENSION(:,:,:), INTENT(  out) ::   p_avm, p_avt   ! vertical eddy viscosity & diffusivity (w-points)
521      !
522      INTEGER  ::   ji, jj, jk   ! dummy loop indices
523      REAL(wp) ::   zrn2, zraug, zcoef, zav   ! local scalars
524      REAL(wp) ::   zdku,   zdkv, zsqen       !   -      -
525      REAL(wp) ::   zemxl, zemlm, zemlp, zmaxice       !   -      -
526      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmxlm, zmxld   ! 3D workspace
527      !!--------------------------------------------------------------------
528      !
529      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
530      !                     !  Mixing length
531      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
532      !
533      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
534      !
535      ! initialisation of interior minimum value (avoid a 2d loop with mikt)
536      zmxlm(:,:,:)  = rmxl_min   
537      zmxld(:,:,:)  = rmxl_min
538      !
539      IF(ln_sdw .AND. ln_mxhsw) THEN
540         zmxlm(:,:,1)= vkarmn * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
541         ! from terray et al 1999 and mellor and blumberg 2004 it should be 0.85 and not 1.6
542         zcoef       = vkarmn * ( (rn_ediff*rn_ediss)**0.25 ) / rn_ediff
543         zmxlm(:,:,1)= zcoef * MAX ( 1.6 * hsw(:,:) , 0.02 )        ! surface mixing length = F(wave height)
544      ELSE
545      !
546         IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rho0*g)
547         !
548            zraug = vkarmn * 2.e5_wp / ( rho0 * grav )
549#if ! defined key_si3 && ! defined key_cice
550            DO_2D( 0, 0, 0, 0 )                  ! No sea-ice
551               zmxlm(ji,jj,1) =  zraug * taum(ji,jj) * tmask(ji,jj,1)
552            END_2D
553#else
554            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
555            !
556            CASE( 0 )                      ! No scaling under sea-ice
557               DO_2D( 0, 0, 0, 0 )
558                  zmxlm(ji,jj,1) = zraug * taum(ji,jj) * tmask(ji,jj,1)
559               END_2D
560               !
561            CASE( 1 )                      ! scaling with constant sea-ice thickness
562               DO_2D( 0, 0, 0, 0 )
563                  zmxlm(ji,jj,1) =  ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
564                     &                          fr_i(ji,jj)   * rn_mxlice           ) * tmask(ji,jj,1)
565               END_2D
566               !
567            CASE( 2 )                      ! scaling with mean sea-ice thickness
568               DO_2D( 0, 0, 0, 0 )
569#if defined key_si3
570                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
571                     &                         fr_i(ji,jj)   * hm_i(ji,jj) * 2._wp ) * tmask(ji,jj,1)
572#elif defined key_cice
573                  zmaxice = MAXVAL( h_i(ji,jj,:) )
574                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
575                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
576#endif
577               END_2D
578               !
579            CASE( 3 )                      ! scaling with max sea-ice thickness
580               DO_2D( 0, 0, 0, 0 )
581                  zmaxice = MAXVAL( h_i(ji,jj,:) )
582                  zmxlm(ji,jj,1) = ( ( 1._wp - fr_i(ji,jj) ) * zraug * taum(ji,jj) + &
583                     &                         fr_i(ji,jj)   * zmaxice             ) * tmask(ji,jj,1)
584               END_2D
585               !
586            END SELECT
587#endif
588            !
589            DO_2D( 0, 0, 0, 0 )
590               zmxlm(ji,jj,1) = MAX( rn_mxl0, zmxlm(ji,jj,1) )
591            END_2D
592            !
593         ELSE
594            zmxlm(:,:,1) = rn_mxl0
595         ENDIF
596      ENDIF
597      !
598      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
599         zrn2 = MAX( rn2(ji,jj,jk), rsmall )
600         zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
601      END_3D
602      !
603      !                     !* Physical limits for the mixing length
604      !
605      zmxld(:,:, 1 ) = zmxlm(:,:,1)   ! surface set to the minimum value
606      zmxld(:,:,jpk) = rmxl_min       ! last level  set to the minimum value
607      !
608      SELECT CASE ( nn_mxl )
609      !
610 !!gm Not sure of that coding for ISF....
611      ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm)
612      CASE ( 0 )           ! bounded by the distance to surface and bottom
613         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
614            zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk),   &
615            &            gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) )
616            ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
617            zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
618               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
619            zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk)   &
620               &            + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk))
621         END_3D
622         !
623      CASE ( 1 )           ! bounded by the vertical scale factor
624         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
625            zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) )
626            zmxlm(ji,jj,jk) = zemxl
627            zmxld(ji,jj,jk) = zemxl
628         END_3D
629         !
630      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
631         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom :
632            zmxlm(ji,jj,jk) =   &
633               &    MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
634         END_3D
635         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface :
636            zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
637            zmxlm(ji,jj,jk) = zemxl
638            zmxld(ji,jj,jk) = zemxl
639         END_3D
640         !
641      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
642         DO_3D( 0, 0, 0, 0, 2, jpkm1 )        ! from the surface to the bottom : lup
643            zmxld(ji,jj,jk) =    &
644               &    MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) )
645         END_3D
646         DO_3DS( 0, 0, 0, 0, jpkm1, 2, -1 )   ! from the bottom to the surface : ldown
647            zmxlm(ji,jj,jk) =   &
648               &    MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) )
649         END_3D
650         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
651            zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
652            zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
653            zmxlm(ji,jj,jk) = zemlm
654            zmxld(ji,jj,jk) = zemlp
655         END_3D
656         !
657      END SELECT
658      !
659      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
660      !                     !  Vertical eddy viscosity and diffusivity  (avm and avt)
661      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
662      DO_3D( 0, 0, 0, 0, 1, jpkm1 )   !* vertical eddy viscosity & diffivity at w-points
663         zsqen = SQRT( en(ji,jj,jk) )
664         zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
665         p_avm(ji,jj,jk) = MAX( zav,                  avmb(jk) ) * wmask(ji,jj,jk)
666         p_avt(ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
667         dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
668      END_3D
669      !
670      !
671      IF( nn_pdl == 1 ) THEN          !* Prandtl number case: update avt
672         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
673            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)
674         END_3D
675      ENDIF
676      !
677      IF(sn_cfctl%l_prtctl) THEN
678         CALL prt_ctl( tab3d_1=en   , clinfo1=' tke  - e: ', tab3d_2=p_avt, clinfo2=' t: ', kdim=jpk)
679         CALL prt_ctl( tab3d_1=p_avm, clinfo1=' tke  - m: ', kdim=jpk )
680      ENDIF
681      !
682   END SUBROUTINE tke_avn
683
684
685   SUBROUTINE zdf_tke_init( Kmm )
686      !!----------------------------------------------------------------------
687      !!                  ***  ROUTINE zdf_tke_init  ***
688      !!                     
689      !! ** Purpose :   Initialization of the vertical eddy diffivity and
690      !!              viscosity when using a tke turbulent closure scheme
691      !!
692      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
693      !!              called at the first timestep (nit000)
694      !!
695      !! ** input   :   Namlist namzdf_tke
696      !!
697      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
698      !!----------------------------------------------------------------------
699      USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag
700      !!
701      INTEGER, INTENT(in) ::   Kmm          ! time level index
702      INTEGER             ::   ji, jj, jk   ! dummy loop indices
703      INTEGER             ::   ios
704      !!
705      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb   , rn_emin  ,  &
706         &                 rn_emin0, rn_bshear, nn_mxl   , ln_mxl0  ,  &
707         &                 rn_mxl0 , nn_mxlice, rn_mxlice,             &
708         &                 nn_pdl  , ln_lc    , rn_lc    ,             &
709         &                 nn_etau , nn_htau  , rn_efr   , nn_eice  ,  &   
710         &                 nn_bc_surf, nn_bc_bot, ln_mxhsw
711      !!----------------------------------------------------------------------
712      !
713      READ  ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
714901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist' )
715
716      READ  ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
717902   IF( ios >  0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist' )
718      IF(lwm) WRITE ( numond, namzdf_tke )
719      !
720      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
721      !
722      IF(lwp) THEN                    !* Control print
723         WRITE(numout,*)
724         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
725         WRITE(numout,*) '~~~~~~~~~~~~'
726         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
727         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
728         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
729         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
730         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
731         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
732         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
733         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
734         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
735         WRITE(numout,*) '         surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
736         WRITE(numout,*) '         surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
737         IF( ln_mxl0 ) THEN
738            WRITE(numout,*) '      type of scaling under sea-ice               nn_mxlice = ', nn_mxlice
739            IF( nn_mxlice == 1 ) &
740            WRITE(numout,*) '      ice thickness when scaling under sea-ice    rn_mxlice = ', rn_mxlice
741            SELECT CASE( nn_mxlice )             ! Type of scaling under sea-ice
742            CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   No scaling under sea-ice'
743            CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   scaling with constant sea-ice thickness'
744            CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   scaling with mean sea-ice thickness'
745            CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   scaling with max sea-ice thickness'
746            CASE DEFAULT
747               CALL ctl_stop( 'zdf_tke_init: wrong value for nn_mxlice, should be 0,1,2,3 or 4')
748            END SELECT
749         ENDIF
750         WRITE(numout,*) '      Langmuir cells parametrization              ln_lc     = ', ln_lc
751         WRITE(numout,*) '         coef to compute vertical velocity of LC     rn_lc  = ', rn_lc
752         IF ( cpl_phioc .and. ln_phioc )  THEN
753            SELECT CASE( nn_bc_surf)             ! Type of scaling under sea-ice
754            CASE( 0 )   ;   WRITE(numout,*) '  nn_bc_surf=0 ==>>> DIRICHLET SBC using surface TKE flux from waves'
755            CASE( 1 )   ;   WRITE(numout,*) '  nn_bc_surf=1 ==>>> NEUMANN SBC using surface TKE flux from waves'
756            END SELECT
757         ENDIF
758         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
759         WRITE(numout,*) '          type of tke penetration profile            nn_htau   = ', nn_htau
760         WRITE(numout,*) '          fraction of TKE that penetrates            rn_efr    = ', rn_efr
761         WRITE(numout,*) '      langmuir & surface wave breaking under ice  nn_eice = ', nn_eice
762         SELECT CASE( nn_eice ) 
763         CASE( 0 )   ;   WRITE(numout,*) '   ==>>>   no impact of ice cover on langmuir & surface wave breaking'
764         CASE( 1 )   ;   WRITE(numout,*) '   ==>>>   weigthed by 1-TANH( fr_i(:,:) * 10 )'
765         CASE( 2 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-fr_i(:,:)'
766         CASE( 3 )   ;   WRITE(numout,*) '   ==>>>   weighted by 1-MIN( 1, 4 * fr_i(:,:) )'
767         CASE DEFAULT
768            CALL ctl_stop( 'zdf_tke_init: wrong value for nn_eice, should be 0,1,2, or 3')
769         END SELECT     
770         WRITE(numout,*)
771         WRITE(numout,*) '   ==>>>   critical Richardson nb with your parameters  ri_cri = ', ri_cri
772         WRITE(numout,*)
773      ENDIF
774      !
775      IF( ln_zdfiwm ) THEN          ! Internal wave-driven mixing
776         rn_emin  = 1.e-10_wp             ! specific values of rn_emin & rmxl_min are used
777         rmxl_min = 1.e-03_wp             ! associated avt minimum = molecular salt diffusivity (10^-9 m2/s)
778         IF(lwp) WRITE(numout,*) '   ==>>>   Internal wave-driven mixing case:   force   rn_emin = 1.e-10 and rmxl_min = 1.e-3'
779      ELSE                          ! standard case : associated avt minimum = molecular viscosity (10^-6 m2/s)
780         rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
781         IF(lwp) WRITE(numout,*) '   ==>>>   minimum mixing length with your parameters rmxl_min = ', rmxl_min
782      ENDIF
783      !
784      !                              ! allocate tke arrays
785      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
786      !
787      !                               !* Check of some namelist values
788      IF( nn_mxl  < 0   .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
789      IF( nn_pdl  < 0   .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
790      IF( nn_htau < 0   .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
791      IF( nn_etau == 3 .AND. .NOT. ln_cpl )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
792      !
793      IF( ln_mxl0 ) THEN
794         IF(lwp) WRITE(numout,*)
795         IF(lwp) WRITE(numout,*) '   ==>>>   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
796         rn_mxl0 = rmxl_min
797      ENDIF
798     
799      IF( nn_etau == 2  )   CALL zdf_mxl( nit000, Kmm )      ! Initialization of nmln
800
801      !                               !* depth of penetration of surface tke
802      IF( nn_etau /= 0 ) THEN     
803         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
804         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
805            htau(:,:) = 10._wp
806         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
807            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
808         END SELECT
809      ENDIF
810      !                                !* read or initialize all required files
811      CALL tke_rst( nit000, 'READ' )      ! (en, avt_k, avm_k, dissl)
812      !
813   END SUBROUTINE zdf_tke_init
814
815
816   SUBROUTINE tke_rst( kt, cdrw )
817      !!---------------------------------------------------------------------
818      !!                   ***  ROUTINE tke_rst  ***
819      !!                     
820      !! ** Purpose :   Read or write TKE file (en) in restart file
821      !!
822      !! ** Method  :   use of IOM library
823      !!                if the restart does not contain TKE, en is either
824      !!                set to rn_emin or recomputed
825      !!----------------------------------------------------------------------
826      USE zdf_oce , ONLY : en, avt_k, avm_k   ! ocean vertical physics
827      !!
828      INTEGER         , INTENT(in) ::   kt     ! ocean time-step
829      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
830      !
831      INTEGER ::   jit, jk              ! dummy loop indices
832      INTEGER ::   id1, id2, id3, id4   ! local integers
833      !!----------------------------------------------------------------------
834      !
835      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
836         !                                   ! ---------------
837         IF( ln_rstart ) THEN                   !* Read the restart file
838            id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
839            id2 = iom_varid( numror, 'avt_k', ldstop = .FALSE. )
840            id3 = iom_varid( numror, 'avm_k', ldstop = .FALSE. )
841            id4 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
842            !
843            IF( MIN( id1, id2, id3, id4 ) > 0 ) THEN      ! fields exist
844               CALL iom_get( numror, jpdom_auto, 'en'   , en    )
845               CALL iom_get( numror, jpdom_auto, 'avt_k', avt_k )
846               CALL iom_get( numror, jpdom_auto, 'avm_k', avm_k )
847               CALL iom_get( numror, jpdom_auto, 'dissl', dissl )
848            ELSE                                          ! start TKE from rest
849               IF(lwp) WRITE(numout,*)
850               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without TKE scheme, set en to background values'
851               en   (:,:,:) = rn_emin * wmask(:,:,:)
852               dissl(:,:,:) = 1.e-12_wp
853               ! avt_k, avm_k already set to the background value in zdf_phy_init
854            ENDIF
855         ELSE                                   !* Start from rest
856            IF(lwp) WRITE(numout,*)
857            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set en to the background value'
858            en   (:,:,:) = rn_emin * wmask(:,:,:)
859            dissl(:,:,:) = 1.e-12_wp
860            ! avt_k, avm_k already set to the background value in zdf_phy_init
861         ENDIF
862         !
863      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
864         !                                   ! -------------------
865         IF(lwp) WRITE(numout,*) '---- tke_rst ----'
866         CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
867         CALL iom_rstput( kt, nitrst, numrow, 'avt_k', avt_k )
868         CALL iom_rstput( kt, nitrst, numrow, 'avm_k', avm_k )
869         CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
870         !
871      ENDIF
872      !
873   END SUBROUTINE tke_rst
874
875   !!======================================================================
876END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.