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

source: branches/2016/dev_v3_6_STABLE_OMP/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 6712

Last change on this file since 6712 was 6712, checked in by mocavero, 8 years ago

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