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

Last change on this file since 3837 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

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