source: CONFIG/UNIFORM/v6/IPSLCM6/SOURCES/NEMO/zdftke.F90 @ 3295

Last change on this file since 3295 was 3295, checked in by cetlod, 4 years ago

CM6.0.11 : minor correction on SOURCES/NEMO/zdftke.F90

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