New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
zdftke.F90 in branches/UKMO/dev_r5518_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

source: branches/UKMO/dev_r5518_nn_etau_revision/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 7294

Last change on this file since 7294 was 7294, checked in by till, 7 years ago

a mod in zdftke.F90 for the dependence on the Prandtl number

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