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

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

Move from jpk to jpkf - trim sub-domains in z
  • Property svn:keywords set to Id
File size: 51.7 KB
Line 
1MODULE zdftke
2   !!======================================================================
3   !!                       ***  MODULE  zdftke  ***
4   !! Ocean physics:  vertical mixing coefficient computed from the tke
5   !!                 turbulent closure parameterization
6   !!=====================================================================
7   !! History :  OPA  !  1991-03  (b. blanke)  Original code
8   !!            7.0  !  1991-11  (G. Madec)   bug fix
9   !!            7.1  !  1992-10  (G. Madec)   new mixing length and eav
10   !!            7.2  !  1993-03  (M. Guyon)   symetrical conditions
11   !!            7.3  !  1994-08  (G. Madec, M. Imbard)  nn_pdl flag
12   !!            7.5  !  1996-01  (G. Madec)   s-coordinates
13   !!            8.0  !  1997-07  (G. Madec)   lbc
14   !!            8.1  !  1999-01  (E. Stretta) new option for the mixing length
15   !!  NEMO      1.0  !  2002-06  (G. Madec) add tke_init routine
16   !!             -   !  2004-10  (C. Ethe )  1D configuration
17   !!            2.0  !  2006-07  (S. Masson)  distributed restart using iom
18   !!            3.0  !  2008-05  (C. Ethe,  G.Madec) : update TKE physics:
19   !!                 !           - tke penetration (wind steering)
20   !!                 !           - suface condition for tke & mixing length
21   !!                 !           - Langmuir cells
22   !!             -   !  2008-05  (J.-M. Molines, G. Madec)  2D form of avtb
23   !!             -   !  2008-06  (G. Madec)  style + DOCTOR name for namelist parameters
24   !!             -   !  2008-12  (G. Reffray) stable discretization of the production term
25   !!            3.2  !  2009-06  (G. Madec, S. Masson) TKE restart compatible with key_cpl
26   !!                 !                                + cleaning of the parameters + bugs correction
27   !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
28   !!----------------------------------------------------------------------
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      ALLOCATE(                                                                    &
121#if defined key_c1d
122         &      e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) ,                          &
123         &      e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) ,                          &
124#endif
125         &      en   (jpi,jpj,jpk) , htau (jpi,jpj)     , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc )
126         !
127      IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
128      IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
129      !
130   END FUNCTION zdf_tke_alloc
131
132
133   SUBROUTINE zdf_tke( kt )
134      !!----------------------------------------------------------------------
135      !!                   ***  ROUTINE zdf_tke  ***
136      !!
137      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
138      !!              coefficients using a turbulent closure scheme (TKE).
139      !!
140      !! ** Method  :   The time evolution of the turbulent kinetic energy (tke)
141      !!              is computed from a prognostic equation :
142      !!         d(en)/dt = avm (d(u)/dz)**2             ! shear production
143      !!                  + d( avm d(en)/dz )/dz         ! diffusion of tke
144      !!                  + avt N^2                      ! stratif. destruc.
145      !!                  - rn_ediss / emxl en**(2/3)    ! Kolmogoroff dissipation
146      !!      with the boundary conditions:
147      !!         surface: en = max( rn_emin0, rn_ebb * taum )
148      !!         bottom : en = rn_emin
149      !!      The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
150      !!
151      !!        The now Turbulent kinetic energy is computed using the following
152      !!      time stepping: implicit for vertical diffusion term, linearized semi
153      !!      implicit for kolmogoroff dissipation term, and explicit forward for
154      !!      both buoyancy and shear production terms. Therefore a tridiagonal
155      !!      linear system is solved. Note that buoyancy and shear terms are
156      !!      discretized in a energy conserving form (Bruchard 2002).
157      !!
158      !!        The dissipative and mixing length scale are computed from en and
159      !!      the stratification (see tke_avn)
160      !!
161      !!        The now vertical eddy vicosity and diffusivity coefficients are
162      !!      given by:
163      !!              avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
164      !!              avt = max( avmb, pdl * avm                 ) 
165      !!              eav = max( avmb, avm )
166      !!      where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
167      !!      given by an empirical funtion of the localRichardson number if nn_pdl=1
168      !!
169      !! ** Action  :   compute en (now turbulent kinetic energy)
170      !!                update avt, avmu, avmv (before vertical eddy coef.)
171      !!
172      !! References : Gaspar et al., JGR, 1990,
173      !!              Blanke and Delecluse, JPO, 1991
174      !!              Mellor and Blumberg, JPO 2004
175      !!              Axell, JGR, 2002
176      !!              Bruchard OM 2002
177      !!----------------------------------------------------------------------
178      INTEGER, INTENT(in) ::   kt   ! ocean time step
179      !!----------------------------------------------------------------------
180      !
181      CALL tke_tke      ! now tke (en)
182      !
183      CALL tke_avn      ! now avt, avm, avmu, avmv
184      !
185   END SUBROUTINE zdf_tke
186
187
188   SUBROUTINE tke_tke
189      !!----------------------------------------------------------------------
190      !!                   ***  ROUTINE tke_tke  ***
191      !!
192      !! ** Purpose :   Compute the now Turbulente Kinetic Energy (TKE)
193      !!
194      !! ** Method  : - TKE surface boundary condition
195      !!              - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
196      !!              - source term due to shear (saved in avmu, avmv arrays)
197      !!              - Now TKE : resolution of the TKE equation by inverting
198      !!                a tridiagonal linear system by a "methode de chasse"
199      !!              - increase TKE due to surface and internal wave breaking
200      !!
201      !! ** Action  : - en : now turbulent kinetic energy)
202      !!              - avmu, avmv : production of TKE by shear at u and v-points
203      !!                (= Kz dz[Ub] * dz[Un] )
204      !! ---------------------------------------------------------------------
205      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released, iwrk_in_use, iwrk_not_released
206      USE oce     , ONLY:   zdiag => ua , zd_up => va , zd_lw => ta   ! (ua,va,ta) used as workspace
207      USE wrk_nemo, ONLY:   imlc  => iwrk_2d_1   ! 2D INTEGER workspace
208      USE wrk_nemo, ONLY:   zhlc  =>  wrk_2d_1   ! 2D REAL workspace
209      USE wrk_nemo, ONLY:   zpelc =>  wrk_3d_1   ! 3D REAL workspace
210
211      !! DCSE_NEMO: need additional directives for renamed module variables
212!FTRANS zdiag zd_up zd_lw :I :I :z
213!FTRANS zpelc :I :I :z
214      !
215      INTEGER  ::   ji, jj, jk                      ! dummy loop arguments
216!!bfr      INTEGER  ::   ikbu, ikbv, ikbum1, ikbvm1      ! temporary scalar
217!!bfr      INTEGER  ::   ikbt, ikbumm1, ikbvmm1          ! temporary scalar
218      REAL(wp) ::   zrhoa  = 1.22                   ! Air density kg/m3
219      REAL(wp) ::   zcdrag = 1.5e-3                 ! drag coefficient
220      REAL(wp) ::   zbbrau, zesh2                   ! temporary scalars
221      REAL(wp) ::   zfact1, zfact2, zfact3          !    -         -
222      REAL(wp) ::   ztx2  , zty2  , zcof            !    -         -
223      REAL(wp) ::   ztau  , zdif                    !    -         -
224      REAL(wp) ::   zus   , zwlc  , zind            !    -         -
225      REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
226!!bfr      REAL(wp) ::   zebot                           !    -         -
227      !!--------------------------------------------------------------------
228      !
229      IF( iwrk_in_use(2, 1) .OR.   &
230           wrk_in_use(2, 1) .OR.   &
231           wrk_in_use(3, 1)   ) THEN
232         CALL ctl_stop('tke_tke: requested workspace arrays unavailable')   ;   RETURN
233      END IF
234
235      zbbrau = rn_ebb / rau0       ! Local constant initialisation
236      zfact1 = -.5_wp * rdt 
237      zfact2 = 1.5_wp * rdt * rn_ediss
238      zfact3 = 0.5_wp       * rn_ediss
239      !
240      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
241      !                     !  Surface boundary condition on tke
242      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
243      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
244         DO ji = fs_2, fs_jpim1   ! vector opt.
245#if defined key_z_first
246            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask_1(ji,jj)
247#else
248            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
249#endif
250         END DO
251      END DO
252     
253!!bfr   - start commented area
254      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
255      !                     !  Bottom boundary condition on tke
256      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
257      !
258      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
259      ! Tests to date have found the bottom boundary condition on tke to have very little effect.
260      ! The condition is coded here for completion but commented out until there is proof that the
261      ! computational cost is justified
262      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
263      !                     en(bot)   = (rn_ebb0/rau0)*0.5*sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
264!CDIR NOVERRCHK
265!!    DO jj = 2, jpjm1
266!CDIR NOVERRCHK
267!!       DO ji = fs_2, fs_jpim1   ! vector opt.
268!!          ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
269!!                 bfrua(ji  ,jj) * ub(ji  ,jj,mbku(ji  ,jj) )
270!!          zty2 = bfrva(ji,jj  ) * vb(ji,jj  ,mbkv(ji,jj  )) + &
271!!                 bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
272!!          zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )   !  where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
273!!          en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
274!!       END DO
275!!    END DO
276!!bfr   - end commented area
277      !
278      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
279      IF( ln_lc ) THEN      !  Langmuir circulation source term added to tke       (Axell JGR 2002)
280         !                  !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
281         !
282         !                        !* total energy produce by LC : cumulative sum over jk
283#if defined key_z_first
284         DO jj = 1, jpj
285            DO ji = 1, jpi
286               zpelc(ji,jj,1) =  MAX( rn2b(ji,jj,1), 0._wp ) * fsdepw(ji,jj,1) * fse3w(ji,jj,1)
287               DO jk = 2, jpkf
288                  zpelc(ji,jj,jk)  = zpelc(ji,jj,jk-1) + MAX( rn2b(ji,jj,jk), 0._wp ) * fsdepw(ji,jj,jk) * fse3w(ji,jj,jk)
289               END DO
290            END DO
291         END DO
292#else
293         zpelc(:,:,1) =  MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
294         DO jk = 2, jpkf
295            zpelc(:,:,jk)  = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
296         END DO
297#endif
298         !                        !* finite Langmuir Circulation depth
299         zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
300         imlc(:,:) = mbkt(:,:) + 1       ! Initialization to the number of w ocean point (=2 over land)
301#if defined key_z_first
302         DO jj = 1, jpj                  ! Last w-level at which zpelc>=0.5*us*us
303            DO ji = 1, jpi               !      with us=0.016*wind(starting from jpk-1)
304               zus  = zcof * taum(ji,jj)
305               DO jk = jpkfm1, 2, -1
306#else
307         DO jk = jpkfm1, 2, -1
308            DO jj = 1, jpj               ! Last w-level at which zpelc>=0.5*us*us
309               DO ji = 1, jpi            !      with us=0.016*wind(starting from jpk-1)
310                  zus  = zcof * taum(ji,jj)
311#endif
312                  IF( zpelc(ji,jj,jk) > zus )   imlc(ji,jj) = jk
313               END DO
314            END DO
315         END DO
316         !                               ! finite LC depth
317# if defined key_vectopt_loop
318         DO jj = 1, 1
319            DO ji = 1, jpij   ! vector opt. (forced unrolling)
320# else
321         DO jj = 1, jpj 
322            DO ji = 1, jpi
323# endif
324               zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
325            END DO
326         END DO
327         zcof = 0.016 / SQRT( zrhoa * zcdrag )
328#if defined key_z_first
329         DO jj = 2, jpjm1         !* TKE Langmuir circulation source term added to en
330            DO ji = 2, jpim1
331               zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
332               DO jk = 2, jpkfm1
333#else
334!CDIR NOVERRCHK
335         DO jk = 2, jpkfm1         !* TKE Langmuir circulation source term added to en
336!CDIR NOVERRCHK
337            DO jj = 2, jpjm1
338!CDIR NOVERRCHK
339               DO ji = fs_2, fs_jpim1   ! vector opt.
340                  zus  = zcof * SQRT( taum(ji,jj) )           ! Stokes drift
341#endif
342                  !                                           ! vertical velocity due to LC
343                  zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
344                  zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
345                  !                                           ! TKE Langmuir circulation source term
346                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
347               END DO
348            END DO
349         END DO
350         !
351      ENDIF
352      !
353      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
354      !                     !  Now Turbulent kinetic energy (output in en)
355      !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
356      !                     ! Resolution of a tridiagonal linear system by a "methode de chasse"
357      !                     ! computation from level 2 to jpkm1  (e(1) already computed and e(jpk)=0 ).
358      !                     ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
359      !
360#if defined key_z_first
361      !* Shear production at uw- and vw-points (energy conserving form)
362             ! here avmu, avmv used as workspace
363      DO jj = 1, jpj
364         DO ji = 1, jpi
365            DO jk = 2, jpkfm1
366#else
367      DO jk = 2, jpkfm1           !* Shear production at uw- and vw-points (energy conserving form)
368         DO jj = 1, jpj                 ! here avmu, avmv used as workspace
369            DO ji = 1, jpi
370#endif
371               avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) )   &
372                  &                            * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) )   & 
373                  &           / (  fse3uw_n(ji,jj,jk)         &
374                  &              * fse3uw_b(ji,jj,jk) )
375               avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) )   &
376                  &                            * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) )   &
377                  &                            / (  fse3vw_n(ji,jj,jk)               &
378                  &                              *  fse3vw_b(ji,jj,jk)  )
379            END DO
380         END DO
381      END DO
382
383      !
384#if defined key_z_first
385      DO jj = 2, jpjm1
386         DO ji = 2, jpim1
387            DO jk = 2, jpkfm1     !* Matrix and right hand side in en
388#else
389      DO jk = 2, jpkfm1           !* Matrix and right hand side in en
390         DO jj = 2, jpjm1
391            DO ji = fs_2, fs_jpim1   ! vector opt.
392#endif
393               zcof   = zfact1 * tmask(ji,jj,jk)
394               zzd_up = zcof * ( avm  (ji,jj,jk+1) + avm  (ji,jj,jk  ) )   &  ! upper diagonal
395                  &          / ( fse3t(ji,jj,jk  ) * fse3w(ji,jj,jk  ) )
396               zzd_lw = zcof * ( avm  (ji,jj,jk  ) + avm  (ji,jj,jk-1) )   &  ! lower diagonal
397                  &          / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk  ) )
398                  !                                                           ! shear prod. at w-point weightened by mask
399               zesh2  =  ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
400                  &    + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )   
401                  !
402               zd_up(ji,jj,jk) = zzd_up            ! Matrix (zdiag, zd_up, zd_lw)
403               zd_lw(ji,jj,jk) = zzd_lw
404               zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
405               !
406               !                                   ! right hand side in en
407               en(ji,jj,jk) = en(ji,jj,jk) + rdt * (  zesh2  -   avt(ji,jj,jk) * rn2(ji,jj,jk)    &
408                  &                                 + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk)  ) * tmask(ji,jj,jk)
409            END DO
410         END DO
411      END DO
412      !                          !* Matrix inversion from level 2 (tke prescribed at level 1)
413#if defined key_z_first
414      DO jj = 2, jpjm1
415         DO ji = 2, jpim1
416            ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
417            DO jk = 3, jpkfm1
418               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
419            END DO
420            ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
421            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
422            DO jk = 3, jpkfm1
423               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)
424            END DO
425            ! Third recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
426            en(ji,jj,jpkfm1) = zd_lw(ji,jj,jpkfm1) / zdiag(ji,jj,jpkfm1)
427            DO jk = jpkf-2, 2, -1
428               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
429            END DO
430            DO jk = 2, jpkfm1                       ! set the minimum value of tke
431               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
432            END DO
433         END DO
434      END DO
435#else
436      DO jk = 3, jpkfm1                             ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
437         DO jj = 2, jpjm1
438            DO ji = fs_2, fs_jpim1    ! vector opt.
439               zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
440            END DO
441         END DO
442      END DO
443      DO jj = 2, jpjm1                             ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
444         DO ji = fs_2, fs_jpim1    ! vector opt.
445            zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1)    ! Surface boudary conditions on tke
446         END DO
447      END DO
448      DO jk = 3, jpkfm1
449         DO jj = 2, jpjm1
450            DO ji = fs_2, fs_jpim1    ! vector opt.
451               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)
452            END DO
453         END DO
454      END DO
455      DO jj = 2, jpjm1                             ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
456         DO ji = fs_2, fs_jpim1    ! vector opt.
457            en(ji,jj,jpkfm1) = zd_lw(ji,jj,jpkfm1) / zdiag(ji,jj,jpkfm1)
458         END DO
459      END DO
460      DO jk = jpkf-2, 2, -1
461         DO jj = 2, jpjm1
462            DO ji = fs_2, fs_jpim1    ! vector opt.
463               en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
464            END DO
465         END DO
466      END DO
467      DO jk = 2, jpkfm1                             ! set the minimum value of tke
468         DO jj = 2, jpjm1
469            DO ji = fs_2, fs_jpim1   ! vector opt.
470               en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
471            END DO
472         END DO
473      END DO
474#endif
475
476      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
477      !                            !  TKE due to surface and internal wave breaking
478      !                            !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
479      IF( nn_etau == 1 ) THEN           !* penetration below the mixed layer (rn_efr fraction)
480#if defined key_z_first
481         DO jj = 2, jpjm1
482            DO ji = 2, jpim1
483               DO jk = 2, jpkfm1
484#else
485         DO jk = 2, jpkfm1
486            DO jj = 2, jpjm1
487               DO ji = fs_2, fs_jpim1   ! vector opt.
488#endif
489                  en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
490                     &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
491               END DO
492            END DO
493         END DO
494      ELSEIF( nn_etau == 2 ) THEN       !* act only at the base of the mixed layer (jk=nmln)  (rn_efr fraction)
495         DO jj = 2, jpjm1
496            DO ji = fs_2, fs_jpim1   ! vector opt.
497               jk = nmln(ji,jj)
498               en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
499                  &                                               * ( 1._wp - fr_i(ji,jj) )  * tmask(ji,jj,jk)
500            END DO
501         END DO
502      ELSEIF( nn_etau == 3 ) THEN       !* penetration belox the mixed layer (HF variability)
503
504         !! DCSE_NEMO: its probably not worth changing the order of these loops for level first indexing,
505         !!            unless we also make zdif a 2-d (jpi,jpj) array
506!CDIR NOVERRCHK
507         DO jk = 2, jpkfm1
508!CDIR NOVERRCHK
509            DO jj = 2, jpjm1
510!CDIR NOVERRCHK
511               DO ji = fs_2, fs_jpim1   ! vector opt.
512                  ztx2 = utau(ji-1,jj  ) + utau(ji,jj)
513                  zty2 = vtau(ji  ,jj-1) + vtau(ji,jj)
514                  ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 )    ! module of the mean stress
515                  zdif = taum(ji,jj) - ztau                            ! mean of modulus - modulus of the mean
516                  zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add )  ! apply some modifications...
517                  en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) )   &
518                     &                                        * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
519               END DO
520            END DO
521         END DO
522      ENDIF
523      CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
524      !
525      IF( iwrk_not_released(2 ,1) .OR.   &
526           wrk_not_released(2, 1) .OR.   &
527           wrk_not_released(3, 1)  )   CALL ctl_stop( 'tke_tke: failed to release workspace arrays' )
528      !
529   END SUBROUTINE tke_tke
530
531   !! * Reset control of array index permutation
532#  include "zdftke_ftrans.h90"
533#  include "oce_ftrans.h90"
534#  include "dom_oce_ftrans.h90"
535#  include "domvvl_ftrans.h90"
536#  include "sbc_oce_ftrans.h90"
537#  include "zdf_oce_ftrans.h90"
538!FTRANS dissl e_dis e_mix e_pdl e_ric :I :I :z
539
540   SUBROUTINE tke_avn
541      !!----------------------------------------------------------------------
542      !!                   ***  ROUTINE tke_avn  ***
543      !!
544      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
545      !!
546      !! ** Method  :   At this stage, en, the now TKE, is known (computed in
547      !!              the tke_tke routine). First, the now mixing lenth is
548      !!      computed from en and the strafification (N^2), then the mixings
549      !!      coefficients are computed.
550      !!              - Mixing length : a first evaluation of the mixing lengh
551      !!      scales is:
552      !!                      mxl = sqrt(2*en) / N 
553      !!      where N is the brunt-vaisala frequency, with a minimum value set
554      !!      to rmxl_min (rn_mxl0) in the interior (surface) ocean.
555      !!        The mixing and dissipative length scale are bound as follow :
556      !!         nn_mxl=0 : mxl bounded by the distance to surface and bottom.
557      !!                        zmxld = zmxlm = mxl
558      !!         nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
559      !!         nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
560      !!                    less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl
561      !!         nn_mxl=3 : mxl is bounded from the surface to the bottom usings
562      !!                    |d/dz(xml)|<1 to obtain lup, and from the bottom to
563      !!                    the surface to obtain ldown. the resulting length
564      !!                    scales are:
565      !!                        zmxld = sqrt( lup * ldown )
566      !!                        zmxlm = min ( lup , ldown )
567      !!              - Vertical eddy viscosity and diffusivity:
568      !!                      avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
569      !!                      avt = max( avmb, pdlr * avm ) 
570      !!      with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
571      !!
572      !! ** Action  : - avt : now vertical eddy diffusivity (w-point)
573      !!              - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
574      !!----------------------------------------------------------------------
575      USE oce, ONLY:   zmpdl => ua , zmxlm => va , zmxld => ta   ! (ua,va,ta) used as workspace
576      USE arpdebugging, ONLY: dump_array
577      !! DCSE_NEMO: need additional directives for renamed module variables
578!FTRANS zmpdl zmxlm zmxld :I :I :z
579      !!
580      INTEGER  ::   ji, jj, jk   ! dummy loop indices
581      REAL(wp) ::   zrn2, zraug, zcoef, zav     ! local scalars
582      REAL(wp) ::   zdku, zpdlr, zri, zsqen     !   -      -
583      REAL(wp) ::   zdkv, zemxl, zemlm, zemlp   !   -      -
584      INTEGER, PARAMETER :: DUMP_LEVEL = 26 ! ARPDBG - level to dump to disk
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      CALL dump_array(1, 'avmu_tke',avmu(:,:,DUMP_LEVEL), &
794                      withHalos=.TRUE.,atStep=1)
795
796      IF( nn_pdl == 1 ) THEN      !* Prandtl number case: update avt
797#if defined key_z_first
798         DO jj = 2, jpjm1
799            DO ji = 2, jpim1
800               DO jk = 2, jpkfm1
801#else
802         DO jk = 2, jpkfm1
803            DO jj = 2, jpjm1
804               DO ji = fs_2, fs_jpim1   ! vector opt.
805#endif
806                  zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
807                  !                                          ! shear
808                  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) )   &
809                    &  + avmu(ji  ,jj,jk) * ( un(ji  ,jj,jk-1) - un(ji  ,jj,jk) ) * ( ub(ji  ,jj,jk-1) - ub(ji  ,jj,jk) )
810                  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) )   &
811                    &  + avmv(ji,jj  ,jk) * ( vn(ji,jj  ,jk-1) - vn(ji,jj  ,jk) ) * ( vb(ji,jj  ,jk-1) - vb(ji,jj  ,jk) )
812                  !                                          ! local Richardson number
813                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
814                  zpdlr = MAX(  0.1_wp,  0.2 / MAX( 0.2 , zri )  )
815!!gm and even better with the use of the "true" ri_crit=0.22222...  (this change the results!)
816!!gm              zpdlr = MAX(  0.1_wp,  ri_crit / MAX( ri_crit , zri )  )
817                  avt(ji,jj,jk)   = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
818# if defined key_c1d
819                  e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk)    ! c1d configuration : save masked Prandlt number
820                  e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)                            ! c1d config. : save Ri
821# endif
822              END DO
823            END DO
824         END DO
825      ENDIF
826      CALL lbc_lnk( avt, 'W', 1. )                      ! Lateral boundary conditions on avt  (sign unchanged)
827
828      IF(ln_ctl) THEN
829         CALL prt_ctl( tab3d_1=en  , clinfo1=' tke  - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
830         CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke  - u: ', mask1=umask,                   &
831            &          tab3d_2=avmv, clinfo2=       ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
832      ENDIF
833      !
834   END SUBROUTINE tke_avn
835
836   !! * Reset control of array index permutation
837#  include "zdftke_ftrans.h90"
838#  include "oce_ftrans.h90"
839#  include "dom_oce_ftrans.h90"
840#  include "domvvl_ftrans.h90"
841#  include "sbc_oce_ftrans.h90"
842#  include "zdf_oce_ftrans.h90"
843!FTRANS dissl e_dis e_mix e_pdl e_ric :I :I :z
844
845   SUBROUTINE zdf_tke_init
846      !!----------------------------------------------------------------------
847      !!                  ***  ROUTINE zdf_tke_init  ***
848      !!                     
849      !! ** Purpose :   Initialization of the vertical eddy diffivity and
850      !!              viscosity when using a tke turbulent closure scheme
851      !!
852      !! ** Method  :   Read the namzdf_tke namelist and check the parameters
853      !!              called at the first timestep (nit000)
854      !!
855      !! ** input   :   Namlist namzdf_tke
856      !!
857      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
858      !!----------------------------------------------------------------------
859      INTEGER ::   ji, jj, jk   ! dummy loop indices
860      !!
861      NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin  ,   &
862         &                 rn_emin0, rn_bshear, nn_mxl , ln_mxl0  ,   &
863         &                 rn_mxl0 , nn_pdl   , ln_lc  , rn_lc    ,   &
864         &                 nn_etau , nn_htau  , rn_efr   
865      !!----------------------------------------------------------------------
866      !
867      REWIND ( numnam )               !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
868      READ   ( numnam, namzdf_tke )
869      !
870      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
871      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
872      !
873      IF(lwp) THEN                    !* Control print
874         WRITE(numout,*)
875         WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
876         WRITE(numout,*) '~~~~~~~~~~~~'
877         WRITE(numout,*) '   Namelist namzdf_tke : set tke mixing parameters'
878         WRITE(numout,*) '      coef. to compute avt                        rn_ediff  = ', rn_ediff
879         WRITE(numout,*) '      Kolmogoroff dissipation coef.               rn_ediss  = ', rn_ediss
880         WRITE(numout,*) '      tke surface input coef.                     rn_ebb    = ', rn_ebb
881         WRITE(numout,*) '      minimum value of tke                        rn_emin   = ', rn_emin
882         WRITE(numout,*) '      surface minimum value of tke                rn_emin0  = ', rn_emin0
883         WRITE(numout,*) '      background shear (>0)                       rn_bshear = ', rn_bshear
884         WRITE(numout,*) '      mixing length type                          nn_mxl    = ', nn_mxl
885         WRITE(numout,*) '      prandl number flag                          nn_pdl    = ', nn_pdl
886         WRITE(numout,*) '      surface mixing length = F(stress) or not    ln_mxl0   = ', ln_mxl0
887         WRITE(numout,*) '      surface  mixing length minimum value        rn_mxl0   = ', rn_mxl0
888         WRITE(numout,*) '      flag to take into acc.  Langmuir circ.      ln_lc     = ', ln_lc
889         WRITE(numout,*) '      coef to compute verticla velocity of LC     rn_lc     = ', rn_lc
890         WRITE(numout,*) '      test param. to add tke induced by wind      nn_etau   = ', nn_etau
891         WRITE(numout,*) '      flag for computation of exp. tke profile    nn_htau   = ', nn_htau
892         WRITE(numout,*) '      fraction of en which pene. the thermocline  rn_efr    = ', rn_efr
893         WRITE(numout,*)
894         WRITE(numout,*) '      critical Richardson nb with your parameters  ri_cri = ', ri_cri
895      ENDIF
896      !
897      !                              ! allocate tke arrays
898      IF( zdf_tke_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
899      !
900      !                               !* Check of some namelist values
901      IF( nn_mxl  < 0  .OR.  nn_mxl  > 3 )   CALL ctl_stop( 'bad flag: nn_mxl is  0, 1 or 2 ' )
902      IF( nn_pdl  < 0  .OR.  nn_pdl  > 1 )   CALL ctl_stop( 'bad flag: nn_pdl is  0 or 1    ' )
903      IF( nn_htau < 0  .OR.  nn_htau > 1 )   CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
904#if ! key_coupled
905      IF( nn_etau == 3 )   CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
906#endif
907
908      IF( ln_mxl0 ) THEN
909         IF(lwp) WRITE(numout,*) '   use a surface mixing length = F(stress) :   set rn_mxl0 = rmxl_min'
910         rn_mxl0 = rmxl_min
911      ENDIF
912     
913      IF( nn_etau == 2  )   CALL zdf_mxl( nit000 )      ! Initialization of nmln
914
915      !                               !* depth of penetration of surface tke
916      IF( nn_etau /= 0 ) THEN     
917         SELECT CASE( nn_htau )             ! Choice of the depth of penetration
918         CASE( 0 )                                 ! constant depth penetration (here 10 meters)
919            htau(:,:) = 10._wp
920         CASE( 1 )                                 ! F(latitude) : 0.5m to 30m poleward of 40 degrees
921            htau(:,:) = MAX(  0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) )   )           
922         END SELECT
923      ENDIF
924      !                               !* set vertical eddy coef. to the background value
925#if defined key_z_first
926      DO jj = 1, jpj
927         DO ji = 1, jpi
928             avt (ji,jj,:jpk) = avtb(:jpk) * tmask(ji,jj,:jpk)
929             avm (ji,jj,:jpk) = avmb(:jpk) * tmask(ji,jj,:jpk)
930             avmu(ji,jj,:jpk) = avmb(:jpk) * umask(ji,jj,:jpk)
931             avmv(ji,jj,:jpk) = avmb(:jpk) * vmask(ji,jj,:jpk)
932         END DO
933      END DO
934#else
935      DO jk = 1, jpk
936         avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
937         avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
938         avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
939         avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
940      END DO
941#endif
942      dissl(:,:,:) = 1.e-12_wp
943      !                             
944      CALL tke_rst( nit000, 'READ' )  !* read or initialize all required files
945      !
946   END SUBROUTINE zdf_tke_init
947
948
949   SUBROUTINE tke_rst( kt, cdrw )
950     !!---------------------------------------------------------------------
951     !!                   ***  ROUTINE tke_rst  ***
952     !!                     
953     !! ** Purpose :   Read or write TKE file (en) in restart file
954     !!
955     !! ** Method  :   use of IOM library
956     !!                if the restart does not contain TKE, en is either
957     !!                set to rn_emin or recomputed
958     !!----------------------------------------------------------------------
959     INTEGER         , INTENT(in) ::   kt     ! ocean time-step
960     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
961     !
962     INTEGER ::   jit, ji, jj, jk   ! dummy loop indices
963     INTEGER ::   id1, id2, id3, id4, id5, id6   ! local integers
964     !!----------------------------------------------------------------------
965     !
966     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise
967        !                                   ! ---------------
968        IF( ln_rstart ) THEN                   !* Read the restart file
969           id1 = iom_varid( numror, 'en'   , ldstop = .FALSE. )
970           id2 = iom_varid( numror, 'avt'  , ldstop = .FALSE. )
971           id3 = iom_varid( numror, 'avm'  , ldstop = .FALSE. )
972           id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
973           id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
974           id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
975           !
976           IF( id1 > 0 ) THEN                       ! 'en' exists
977              CALL iom_get( numror, jpdom_autoglo, 'en', en )
978              IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN        ! all required arrays exist
979                 CALL iom_get( numror, jpdom_autoglo, 'avt'  , avt   )
980                 CALL iom_get( numror, jpdom_autoglo, 'avm'  , avm   )
981                 CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu  )
982                 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv  )
983                 CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
984              ELSE                                                 ! one at least array is missing
985                 CALL tke_avn                                          ! compute avt, avm, avmu, avmv and dissl (approximation)
986              ENDIF
987           ELSE                                     ! No TKE array found: initialisation
988              IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
989              en (:,:,:) = rn_emin * tmask(:,:,:)
990              CALL tke_avn                               ! recompute avt, avm, avmu, avmv and dissl (approximation)
991              DO jit = nit000 + 1, nit000 + 10   ;   CALL zdf_tke( jit )   ;   END DO
992           ENDIF
993        ELSE                                   !* Start from rest
994           en(:,:,:) = rn_emin * tmask(:,:,:)
995#if defined key_z_first
996           DO jj = 1, jpj                           ! set the Kz to the background value
997              DO ji = 1, jpi
998                  avt (ji,jj,:) = avtb(:) * tmask(ji,jj,:)
999                  avm (ji,jj,:) = avmb(:) * tmask(ji,jj,:)
1000                  avmu(ji,jj,:) = avmb(:) * umask(ji,jj,:)
1001                  avmv(ji,jj,:) = avmb(:) * vmask(ji,jj,:)
1002              END DO
1003           END DO
1004#else
1005           DO jk = 1, jpk                           ! set the Kz to the background value
1006              avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
1007              avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
1008              avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
1009              avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
1010           END DO
1011#endif
1012
1013        ENDIF
1014        !
1015     ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
1016        !                                   ! -------------------
1017        IF(lwp) WRITE(numout,*) '---- tke-rst ----'
1018        CALL iom_rstput( kt, nitrst, numrow, 'en'   , en    )
1019        CALL iom_rstput( kt, nitrst, numrow, 'avt'  , avt   )
1020        CALL iom_rstput( kt, nitrst, numrow, 'avm'  , avm   )
1021        CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu  )
1022        CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv  )
1023        CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
1024        !
1025     ENDIF
1026     !
1027   END SUBROUTINE tke_rst
1028
1029#else
1030   !!----------------------------------------------------------------------
1031   !!   Dummy module :                                        NO TKE scheme
1032   !!----------------------------------------------------------------------
1033   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftke = .FALSE.   !: TKE flag
1034CONTAINS
1035   SUBROUTINE zdf_tke_init           ! Dummy routine
1036   END SUBROUTINE zdf_tke_init
1037   SUBROUTINE zdf_tke( kt )          ! Dummy routine
1038      WRITE(*,*) 'zdf_tke: You should not have seen this print! error?', kt
1039   END SUBROUTINE zdf_tke
1040   SUBROUTINE tke_rst( kt, cdrw )
1041     CHARACTER(len=*) ::   cdrw
1042     WRITE(*,*) 'tke_rst: You should not have seen this print! error?', kt, cdwr
1043   END SUBROUTINE tke_rst
1044#endif
1045
1046   !!======================================================================
1047END MODULE zdftke
Note: See TracBrowser for help on using the repository browser.