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 @ 4425

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

Removed dump_array calls
  • Property svn:keywords set to Id
File size: 51.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   !!----------------------------------------------------------------------
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      !! DCSE_NEMO: need additional directives for renamed module variables
579!FTRANS zmpdl zmxlm zmxld :I :I :z
580      !!
581      INTEGER  ::   ji, jj, jk   ! dummy loop indices
582      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
583      REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      -
584      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
585      !!--------------------------------------------------------------------
586
587      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
588      !                     !  Mixing length
589      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
590      !
591      !                     !* Buoyancy length scale: l=sqrt(2*e/n**2)
592      !
593      IF( ln_mxl0 ) THEN            ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g)
594         zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
595         zmxlm(:,:,1) = MAX(  rn_mxl0,  zraug * taum(:,:)  )
596      ELSE                          ! surface set to the minimum value
597         zmxlm(:,:,1) = rn_mxl0
598      ENDIF
599
600#if defined key_z_first
601      DO jj = 2, jpjm1
602         DO ji = 2, jpim1
603            zmxlm(ji,jj,jpkf) = rmxl_min     ! last level set to the interior minium value
604            DO jk = 2, jpkfm1        ! interior value : l=sqrt(2*e/n^2)
605#else
606      zmxlm(:,:,jpkf)  = rmxl_min    ! last level set to the interior minium value
607      !
608!CDIR NOVERRCHK
609      DO jk = 2, jpkfm1              ! interior value : l=sqrt(2*e/n^2)
610!CDIR NOVERRCHK
611         DO jj = 2, jpjm1
612!CDIR NOVERRCHK
613            DO ji = fs_2, fs_jpim1   ! vector opt.
614#endif
615               zrn2 = MAX( rn2(ji,jj,jk), rsmall )
616               zmxlm(ji,jj,jk) = MAX(  rmxl_min,  SQRT( 2._wp * en(ji,jj,jk) / zrn2 )  )
617            END DO
618         END DO
619      END DO
620      !
621      !                     !* Physical limits for the mixing length
622      !
623      zmxld(:,:, 1  ) = zmxlm(:,:,1)   ! surface set to the zmxlm   value
624      zmxld(:,:,jpkf) = rmxl_min       ! last level  set to the minimum value
625      !
626      SELECT CASE ( nn_mxl )
627      !
628      CASE ( 0 )           ! bounded by the distance to surface and bottom
629#if defined key_z_first
630         DO jj = 2, jpjm1
631            DO ji = 2, jpim1
632               DO jk = 2, jpkfm1
633#else
634         DO jk = 2, jpkfm1
635            DO jj = 2, jpjm1
636               DO ji = fs_2, fs_jpim1   ! vector opt.
637#endif
638                  zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk),   &
639                  &            fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
640                  zmxlm(ji,jj,jk) = zemxl
641                  zmxld(ji,jj,jk) = zemxl
642               END DO
643            END DO
644         END DO
645         !
646      CASE ( 1 )           ! bounded by the vertical scale factor
647#if defined key_z_first
648         DO jj = 2, jpjm1
649            DO ji = 2, jpim1
650               DO jk = 2, jpkfm1
651#else
652         DO jk = 2, jpkfm1
653            DO jj = 2, jpjm1
654               DO ji = fs_2, fs_jpim1   ! vector opt.
655#endif
656                  zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
657                  zmxlm(ji,jj,jk) = zemxl
658                  zmxld(ji,jj,jk) = zemxl
659               END DO
660            END DO
661         END DO
662         !
663      CASE ( 2 )           ! |dk[xml]| bounded by e3t :
664#if defined key_z_first
665         DO jj = 2, jpjm1
666            DO ji = 2, jpim1
667               DO jk = 2, jpkfm1   ! from the surface to the bottom :
668                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
669               END DO
670               DO jk = jpkfm1, 2, -1     ! from the bottom to the surface :
671                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
672                  zmxlm(ji,jj,jk) = zemxl
673                  zmxld(ji,jj,jk) = zemxl
674               END DO
675            END DO
676         END DO
677#else
678         DO jk = 2, jpkfm1         ! from the surface to the bottom :
679            DO jj = 2, jpjm1
680               DO ji = fs_2, fs_jpim1   ! vector opt.
681                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
682               END DO
683            END DO
684         END DO
685         DO jk = jpkfm1, 2, -1     ! from the bottom to the surface :
686            DO jj = 2, jpjm1
687               DO ji = fs_2, fs_jpim1   ! vector opt.
688                  zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
689                  zmxlm(ji,jj,jk) = zemxl
690                  zmxld(ji,jj,jk) = zemxl
691               END DO
692            END DO
693         END DO
694#endif
695         !
696      CASE ( 3 )           ! lup and ldown, |dk[xml]| bounded by e3t :
697#if defined key_z_first
698         DO jj = 2, jpjm1
699            DO ji = 2, jpim1
700               DO jk = 2, jpkfm1         ! from the surface to the bottom : lup
701                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
702               END DO
703               DO jk = jpkfm1, 2, -1     ! from the bottom to the surface : ldown
704                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
705               END DO
706               DO jk = 2, jpkfm1
707                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
708                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
709                  zmxlm(ji,jj,jk) = zemlm
710                  zmxld(ji,jj,jk) = zemlp
711               END DO
712            END DO
713         END DO
714#else
715         DO jk = 2, jpkfm1         ! from the surface to the bottom : lup
716            DO jj = 2, jpjm1
717               DO ji = fs_2, fs_jpim1   ! vector opt.
718                  zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
719               END DO
720            END DO
721         END DO
722         DO jk = jpkfm1, 2, -1     ! from the bottom to the surface : ldown
723            DO jj = 2, jpjm1
724               DO ji = fs_2, fs_jpim1   ! vector opt.
725                  zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
726               END DO
727            END DO
728         END DO
729!CDIR NOVERRCHK
730         DO jk = 2, jpkfm1
731!CDIR NOVERRCHK
732            DO jj = 2, jpjm1
733!CDIR NOVERRCHK
734               DO ji = fs_2, fs_jpim1   ! vector opt.
735                  zemlm = MIN ( zmxld(ji,jj,jk),  zmxlm(ji,jj,jk) )
736                  zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
737                  zmxlm(ji,jj,jk) = zemlm
738                  zmxld(ji,jj,jk) = zemlp
739               END DO
740            END DO
741         END DO
742#endif
743         !
744      END SELECT
745      !
746# if defined key_c1d
747      e_dis(:,:,:) = zmxld(:,:,:)      ! c1d configuration : save mixing and dissipation turbulent length scales
748      e_mix(:,:,:) = zmxlm(:,:,:)
749# endif
750
751      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
752      !                     !  Vertical eddy viscosity and diffusivity  (avmu, avmv, avt)
753      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
754#if defined key_z_first
755      DO jj = 2, jpjm1
756         DO ji = 2, jpim1
757            DO jk = 1, jpkfm1      !* vertical eddy viscosity & diffivity at w-points
758#else
759!CDIR NOVERRCHK
760      DO jk = 1, jpkfm1            !* vertical eddy viscosity & diffivity at w-points
761!CDIR NOVERRCHK
762         DO jj = 2, jpjm1
763!CDIR NOVERRCHK
764            DO ji = fs_2, fs_jpim1   ! vector opt.
765#endif
766               zsqen = SQRT( en(ji,jj,jk) )
767               zav   = rn_ediff * zmxlm(ji,jj,jk) * zsqen
768               avm  (ji,jj,jk) = MAX( zav,                  avmb(jk) ) * tmask(ji,jj,jk)
769               avt  (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
770               dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
771            END DO
772         END DO
773      END DO
774      CALL lbc_lnk( avm, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
775      !
776#if defined key_z_first
777      DO jj = 2, jpjm1
778         DO ji = 2, jpim1
779            DO jk = 2, jpkfm1      !* vertical eddy viscosity at u- and v-points
780#else
781      DO jk = 2, jpkfm1            !* vertical eddy viscosity at u- and v-points
782         DO jj = 2, jpjm1
783            DO ji = fs_2, fs_jpim1   ! vector opt.
784#endif
785               avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj  ,jk) ) * umask(ji,jj,jk)
786               avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji  ,jj+1,jk) ) * vmask(ji,jj,jk)
787            END DO
788         END DO
789      END DO
790
791      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! Lateral boundary conditions
792      !
793      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
794#if defined key_z_first
795         DO jj = 2, jpjm1
796            DO ji = 2, jpim1
797               DO jk = 2, jpkfm1
798#else
799         DO jk = 2, jpkfm1
800            DO jj = 2, jpjm1
801               DO ji = fs_2, fs_jpim1   ! vector opt.
802#endif
803                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
804                  !                                          ! shear
805                  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) )   &
806                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
807                  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) )   &
808                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
809                  !                                          ! local Richardson number
810                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
811                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
812!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
813!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
814                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
815# if defined key_c1d
816                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
817                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
818# endif
819              END DO
820            END DO
821         END DO
822      ENDIF
823      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
824
825      IF(ln_ctl) THEN
826         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
827         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
828            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
829      ENDIF
830      !
831   END SUBROUTINE tke_avn
832
833   !! * Reset control of array index permutation
834#  include "zdftke_ftrans.h90"
835#  include "oce_ftrans.h90"
836#  include "dom_oce_ftrans.h90"
837#  include "domvvl_ftrans.h90"
838#  include "sbc_oce_ftrans.h90"
839#  include "zdf_oce_ftrans.h90"
840!FTRANS dissl e_dis e_mix e_pdl e_ric :I :I :z
841
842   SUBROUTINE zdf_tke_init
843      !!----------------------------------------------------------------------
844      !!                  ***  ROUTINE zdf_tke_init  ***
845      !!                     
846      !! ** Purpose :   Initialization of the vertical eddy diffivity and
847      !!              viscosity when using a tke turbulent closure scheme
848      !!
849      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
850      !!              called at the first timestep (nit000)
851      !!
852      !! ** input   :   Namlist namzdf_tke
853      !!
854      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
855      !!----------------------------------------------------------------------
856      INTEGER ::   ji, jj, jk   ! dummy loop indices
857      !!
858      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
859         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
860         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
861         &                 nn_etau , nn_htau  , rn_efr   
862      !!----------------------------------------------------------------------
863      !
864      REWIND ( numnam )               !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
865      READ   ( numnam, namzdf_tke )
866      !
867      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
868      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
869      !
870      IF(lwp) THEN                    !* Control print
871         WRITE(numout,*)
872         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
873         WRITE(numout,*) '~~~~~~~~~~~~'
874         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
875         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
876         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
877         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
878         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
879         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
880         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
881         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
882         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
883         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
884         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
885         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
886         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
887         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
888         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
889         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
890         WRITE(numout,*)
891         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
892      ENDIF
893      !
894      !                              ! allocate tke arrays
895      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
896      !
897      !                               !* Check of some namelist values
898      IF( nn_mxl  < 0  .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
899      IF( nn_pdl  < 0  .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
900      IF( nn_htau < 0  .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
901#if ! key_coupled
902      IF( nn_etau == 3 )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
903#endif
904
905      IF( ln_mxl0 ) THEN
906         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
907         rn_mxl0 = rmxl_min
908      ENDIF
909     
910      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
911
912      !                               !* depth of penetration of surface tke
913      IF( nn_etau /= 0 ) THEN     
914         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
915         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
916            htau(:,:) = 10._wp
917         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
918            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
919         END SELECT
920      ENDIF
921      !                               !* set vertical eddy coef. to the background value
922#if defined key_z_first
923      DO jj = 1, jpj
924         DO ji = 1, jpi
925             avt (ji,jj,:jpk) = avtb(:jpk) * tmask(ji,jj,:jpk)
926             avm (ji,jj,:jpk) = avmb(:jpk) * tmask(ji,jj,:jpk)
927             avmu(ji,jj,:jpk) = avmb(:jpk) * umask(ji,jj,:jpk)
928             avmv(ji,jj,:jpk) = avmb(:jpk) * vmask(ji,jj,:jpk)
929         END DO
930      END DO
931#else
932      DO jk = 1, jpk
933         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
934         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
935         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
936         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
937      END DO
938#endif
939      dissl(:,:,:) = 1.e-12_wp
940      !                             
941      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
942      !
943   END SUBROUTINE zdf_tke_init
944
945
946   SUBROUTINE tke_rst( kt, cdrw )
947     !!---------------------------------------------------------------------
948     !!                   ***  ROUTINE tke_rst  ***
949     !!                     
950     !! ** Purpose :   Read or write TKE file (en) in restart file
951     !!
952     !! ** Method  :   use of IOM library
953     !!                if the restart does not contain TKE, en is either
954     !!                set to rn_emin or recomputed
955     !!----------------------------------------------------------------------
956     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
957     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
958     !
959     INTEGER ::   jit, ji, jj, jk   ! dummy loop indices
960     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
961     !!----------------------------------------------------------------------
962     !
963     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
964        !                                   ! ---------------
965        IF( ln_rstart ) THEN                   !* Read the restart file
966           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
967           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
968           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
969           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
970           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
971           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
972           !
973           IF( id1 > 0 ) THEN                       ! 'en' exists
974              CALL iom_get( numror, jpdom_autoglo, 'en', en )
975              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
976                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
977                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
978                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
979                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
980                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
981              ELSE                                                 ! one at least array is missing
982                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
983              ENDIF
984           ELSE                                     ! No TKE array found: initialisation
985              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
986              en (:,:,:) = rn_emin * tmask(:,:,:)
987              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
988              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
989           ENDIF
990        ELSE                                   !* Start from rest
991           en(:,:,:) = rn_emin * tmask(:,:,:)
992#if defined key_z_first
993           DO jj = 1, jpj                           ! set the Kz to the background value
994              DO ji = 1, jpi
995                  avt (ji,jj,:) = avtb(:) * tmask(ji,jj,:)
996                  avm (ji,jj,:) = avmb(:) * tmask(ji,jj,:)
997                  avmu(ji,jj,:) = avmb(:) * umask(ji,jj,:)
998                  avmv(ji,jj,:) = avmb(:) * vmask(ji,jj,:)
999              END DO
1000           END DO
1001#else
1002           DO jk = 1, jpk                           ! set the Kz to the background value
1003              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1004              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
1005              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1006              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1007           END DO
1008#endif
1009
1010        ENDIF
1011        !
1012     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1013        !                                   ! -------------------
1014        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
1015        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
1016        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt   )
1017        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm   )
1018        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu  )
1019        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  )
1020        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
1021        !
1022     ENDIF
1023     !
1024   END SUBROUTINE tke_rst
1025
1026#else
1027   !!----------------------------------------------------------------------
1028   !!   Dummy module :                                        NO TKE scheme
1029   !!----------------------------------------------------------------------
1030   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
1031CONTAINS
1032   SUBROUTINE zdf_tke_init           ! Dummy routine
1033   END SUBROUTINE zdf_tke_init
1034   SUBROUTINE zdf_tke( kt )          ! Dummy routine
1035      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
1036   END SUBROUTINE zdf_tke
1037   SUBROUTINE tke_rst( kt, cdrw )
1038     CHARACTER(len=*) ::   cdrw
1039     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
1040   END SUBROUTINE tke_rst
1041#endif
1042
1043   !!======================================================================
1044END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.