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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 6748

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

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