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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago

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