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.

zdftmx.F90 in branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90 @ 7525

Last change on this file since 7525 was 7525, checked in by mocavero, 7 years ago
changes after review
Property svn:keywords set to `Id`
File size: 60.5 KB

Line
1	MODULE zdftmx
2	!!========================================================================
3	!! * MODULE zdftmx *
4	!! Ocean physics: vertical tidal mixing coefficient
5	!!========================================================================
6	!! History : 1.0 ! 2004-04 (L. Bessieres, G. Madec) Original code
7	!! - ! 2006-08 (A. Koch-Larrouy) Indonesian strait
8	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
9	!!----------------------------------------------------------------------
10	#if defined key_zdftmx
11	!!----------------------------------------------------------------------
12	!! 'key_zdftmx' Tidal vertical mixing
13	!!----------------------------------------------------------------------
14	!! zdf_tmx : global momentum & tracer Kz with tidal induced Kz
15	!! tmx_itf : Indonesian momentum & tracer Kz with tidal induced Kz
16	!!----------------------------------------------------------------------
17	USE oce ! ocean dynamics and tracers variables
18	USE dom_oce ! ocean space and time domain variables
19	USE zdf_oce ! ocean vertical physics variables
20	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
21	USE eosbn2 ! ocean equation of state
22	USE phycst ! physical constants
23	USE prtctl ! Print control
24	USE in_out_manager ! I/O manager
25	USE iom ! I/O Manager
26	USE lib_mpp ! MPP library
27	USE wrk_nemo ! work arrays
28	USE timing ! Timing
29	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
30
31	IMPLICIT NONE
32	PRIVATE
33
34	PUBLIC zdf_tmx ! called in step module
35	PUBLIC zdf_tmx_init ! called in opa module
36	PUBLIC zdf_tmx_alloc ! called in nemogcm module
37
38	LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .TRUE. !: tidal mixing flag
39
40	! !!* Namelist namzdf_tmx : tidal mixing *
41	REAL(wp) :: rn_htmx ! vertical decay scale for turbulence (meters)
42	REAL(wp) :: rn_n2min ! threshold of the Brunt-Vaisala frequency (s-1)
43	REAL(wp) :: rn_tfe ! tidal dissipation efficiency (St Laurent et al. 2002)
44	REAL(wp) :: rn_me ! mixing efficiency (Osborn 1980)
45	LOGICAL :: ln_tmx_itf ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
46	REAL(wp) :: rn_tfe_itf ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
47
48	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: en_tmx ! energy available for tidal mixing (W/m2)
49	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mask_itf ! mask to use over Indonesian area
50	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: az_tmx ! coefficient used to evaluate the tidal induced Kz
51
52	!! * Substitutions
53	# include "vectopt_loop_substitute.h90"
54	!!----------------------------------------------------------------------
55	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
56	!! $Id$
57	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
58	!!----------------------------------------------------------------------
59	CONTAINS
60
61	INTEGER FUNCTION zdf_tmx_alloc()
62	!!----------------------------------------------------------------------
63	!! * FUNCTION zdf_tmx_alloc *
64	!!----------------------------------------------------------------------
65	ALLOCATE(en_tmx(jpi,jpj), mask_itf(jpi,jpj), az_tmx(jpi,jpj,jpk), STAT=zdf_tmx_alloc )
66	!
67	IF( lk_mpp ) CALL mpp_sum ( zdf_tmx_alloc )
68	IF( zdf_tmx_alloc /= 0 ) CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
69	END FUNCTION zdf_tmx_alloc
70
71
72	SUBROUTINE zdf_tmx( kt )
73	!!----------------------------------------------------------------------
74	!! * ROUTINE zdf_tmx *
75	!!
76	!! ** Purpose : add to the vertical mixing coefficients the effect of
77	!! tidal mixing (Simmons et al 2004).
78	!!
79	!! ** Method : - tidal-induced vertical mixing is given by:
80	!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
81	!! where az_tmx is a coefficient that specified the 3D space
82	!! distribution of the faction of tidal energy taht is used
83	!! for mixing. Its expression is set in zdf_tmx_init routine,
84	!! following Simmons et al. 2004.
85	!! NB: a specific bounding procedure is performed on av_tide
86	!! so that the input tidal energy is actually almost used. The
87	!! basic maximum value is 60 cm2/s, but values of 300 cm2/s
88	!! can be reached in area where bottom stratification is too
89	!! weak.
90	!!
91	!! - update av_tide in the Indonesian Through Flow area
92	!! following Koch-Larrouy et al. (2007) parameterisation
93	!! (see tmx_itf routine).
94	!!
95	!! - update the model vertical eddy viscosity and diffusivity:
96	!! avt = avt + av_tides
97	!! avm = avm + av_tides
98	!! avmu = avmu + mi(av_tides)
99	!! avmv = avmv + mj(av_tides)
100	!!
101	!! ** Action : avt, avm, avmu, avmv increased by tidal mixing
102	!!
103	!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
104	!! Koch-Larrouy et al. 2007, GRL.
105	!!----------------------------------------------------------------------
106	INTEGER, INTENT(in) :: kt ! ocean time-step
107	!
108	INTEGER :: ji, jj, jk ! dummy loop indices
109	REAL(wp) :: ztpc ! scalar workspace
110	REAL(wp), POINTER, DIMENSION(:,:) :: zkz
111	REAL(wp), POINTER, DIMENSION(:,:,:) :: zav_tide
112	!!----------------------------------------------------------------------
113	!
114	IF( nn_timing == 1 ) CALL timing_start('zdf_tmx')
115	!
116	CALL wrk_alloc( jpi,jpj, zkz )
117	CALL wrk_alloc( jpi,jpj,jpk, zav_tide )
118	!
119	! ! ----------------------- !
120	! ! Standard tidal mixing ! (compute zav_tide)
121	! ! ----------------------- !
122	! !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
123	!$OMP PARALLEL
124	!$OMP DO schedule(static) private(jk,jj,ji)
125	DO jk = 1, jpk
126	DO jj = 1, jpj
127	DO ji = 1, jpi
128	zav_tide(ji,jj,jk) = MIN( 60.e-4, az_tmx(ji,jj,jk) / MAX( rn_n2min, rn2(ji,jj,jk) ) )
129	END DO
130	END DO
131	END DO
132	!$OMP END DO NOWAIT
133
134	!$OMP DO schedule(static) private(jj, ji)
135	DO jj = 1, jpj
136	DO ji = 1, jpi
137	zkz(ji,jj) = 0.e0 !* Associated potential energy consummed over the whole water column
138	END DO
139	END DO
140	DO jk = 2, jpkm1
141	!$OMP DO schedule(static) private(jj, ji)
142	DO jj = 1, jpj
143	DO ji = 1, jpi
144	zkz(ji,jj) = zkz(ji,jj) + e3w_n(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
145	END DO
146	END DO
147	END DO
148
149	!$OMP DO schedule(static) private(jj, ji)
150	DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
151	DO ji = 1, jpi
152	IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
153	END DO
154	END DO
155
156	DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
157	!$OMP DO schedule(static) private(jj, ji)
158	DO jj = 1, jpj
159	DO ji = 1, jpi
160	zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s
161	END DO
162	END DO
163	!$OMP END DO NOWAIT
164	END DO
165	!$OMP END PARALLEL
166
167	IF( kt == nit000 ) THEN !* check at first time-step: diagnose the energy consumed by zav_tide
168	ztpc = 0._wp
169	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,ztpc)
170	DO jk= 1, jpk
171	DO jj= 1, jpj
172	DO ji= 1, jpi
173	ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj) &
174	& * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
175	END DO
176	END DO
177	END DO
178	ztpc= rau0 / ( rn_tfe * rn_me ) * ztpc
179	IF( lk_mpp ) CALL mpp_sum( ztpc )
180	IF(lwp) WRITE(numout,*)
181	IF(lwp) WRITE(numout,) ' N Total power consumption by av_tide : ztpc = ', ztpc 1.e-12 ,'TW'
182	ENDIF
183
184	! ! ----------------------- !
185	! ! ITF tidal mixing ! (update zav_tide)
186	! ! ----------------------- !
187	IF( ln_tmx_itf ) CALL tmx_itf( kt, zav_tide )
188
189	! ! ----------------------- !
190	! ! Update mixing coefs !
191	! ! ----------------------- !
192	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
193	DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing
194	DO jj = 1, jpj
195	DO ji = 1, jpi ! vector opt.
196	avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
197	avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
198	END DO
199	END DO
200	DO jj = 2, jpjm1
201	DO ji = fs_2, fs_jpim1 ! vector opt.
202	avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * wumask(ji,jj,jk)
203	avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * wvmask(ji,jj,jk)
204	END DO
205	END DO
206	END DO
207	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! lateral boundary condition
208
209	! !* output tidal mixing coefficient
210	CALL iom_put( "av_tide", zav_tide )
211
212	IF(ln_ctl) CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
213	!
214	CALL wrk_dealloc( jpi,jpj, zkz )
215	CALL wrk_dealloc( jpi,jpj,jpk, zav_tide )
216	!
217	IF( nn_timing == 1 ) CALL timing_stop('zdf_tmx')
218	!
219	END SUBROUTINE zdf_tmx
220
221
222	SUBROUTINE tmx_itf( kt, pav )
223	!!----------------------------------------------------------------------
224	!! * ROUTINE tmx_itf *
225	!!
226	!! ** Purpose : modify the vertical eddy diffusivity coefficients
227	!! (pav) in the Indonesian Through Flow area (ITF).
228	!!
229	!! ** Method : - Following Koch-Larrouy et al. (2007), in the ITF defined
230	!! by msk_itf (read in a file, see tmx_init), the tidal
231	!! mixing coefficient is computed with :
232	!! * q=1 (i.e. all the tidal energy remains trapped in
233	!! the area and thus is used for mixing)
234	!! * the vertical distribution of the tifal energy is a
235	!! proportional to N above the thermocline (d(N^2)/dz > 0)
236	!! and to N^2 below the thermocline (d(N^2)/dz < 0)
237	!!
238	!! ** Action : av_tide updated in the ITF area (msk_itf)
239	!!
240	!! References : Koch-Larrouy et al. 2007, GRL
241	!!----------------------------------------------------------------------
242	INTEGER , INTENT(in ) :: kt ! ocean time-step
243	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pav ! Tidal mixing coef.
244	!!
245	INTEGER :: ji, jj, jk ! dummy loop indices
246	REAL(wp) :: zcoef, ztpc ! temporary scalar
247	REAL(wp), DIMENSION(:,:) , POINTER :: zkz ! 2D workspace
248	REAL(wp), DIMENSION(:,:) , POINTER :: zsum1 , zsum2 , zsum ! - -
249	REAL(wp), DIMENSION(:,:,:), POINTER :: zempba_3d_1, zempba_3d_2 ! 3D workspace
250	REAL(wp), DIMENSION(:,:,:), POINTER :: zempba_3d , zdn2dz ! - -
251	REAL(wp), DIMENSION(:,:,:), POINTER :: zavt_itf ! - -
252	!!----------------------------------------------------------------------
253	!
254	IF( nn_timing == 1 ) CALL timing_start('tmx_itf')
255	!
256	CALL wrk_alloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
257	CALL wrk_alloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
258
259	! ! compute the form function using N2 at each time step
260	!$OMP PARALLEL
261	!$OMP DO schedule(static) private(jj, ji)
262	DO jj = 1, jpj
263	DO ji = 1, jpi
264	zempba_3d_1(ji,jj,jpk) = 0.e0
265	zempba_3d_2(ji,jj,jpk) = 0.e0
266	END DO
267	END DO
268	!$OMP DO schedule(static) private(jk)
269	DO jk = 1, jpkm1
270	zdn2dz (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1) ! Vertical profile of dN2/dz
271	zempba_3d_1(:,:,jk) = SQRT( MAX( 0.e0, rn2(:,:,jk) ) ) ! - - of N
272	zempba_3d_2(:,:,jk) = MAX( 0.e0, rn2(:,:,jk) ) ! - - of N^2
273	END DO
274	!$OMP END DO NOWAIT
275	!
276	!$OMP DO schedule(static) private(jj, ji)
277	DO jj = 1, jpj
278	DO ji = 1, jpi
279	zsum (ji,jj) = 0.e0
280	zsum1(ji,jj) = 0.e0
281	zsum2(ji,jj) = 0.e0
282	END DO
283	END DO
284	DO jk= 2, jpk
285	!$OMP DO schedule(static) private(jj,ji)
286	DO jj= 1, jpj
287	DO ji= 1, jpi
288	zsum1(ji,jj) = zsum1(ji,jj) + zempba_3d_1(ji,jj,jk) * e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
289	zsum2(ji,jj) = zsum2(ji,jj) + zempba_3d_2(ji,jj,jk) * e3w_n(ji,jj,jk) * wmask(ji,jj,jk)
290	END DO
291	END DO
292	END DO
293	!$OMP DO schedule(static) private(jj,ji)
294	DO jj = 1, jpj
295	DO ji = 1, jpi
296	IF( zsum1(ji,jj) /= 0.e0 ) zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
297	IF( zsum2(ji,jj) /= 0.e0 ) zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)
298	END DO
299	END DO
300
301	DO jk= 1, jpk
302	!$OMP DO schedule(static) private(jj,ji,zcoef,ztpc)
303	DO jj = 1, jpj
304	DO ji = 1, jpi
305	zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) ) ! =0 if dN2/dz > 0, =1 otherwise
306	ztpc = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) * zcoef &
307	& + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef )
308	!
309	zempba_3d(ji,jj,jk) = ztpc
310	zsum (ji,jj) = zsum(ji,jj) + ztpc * e3w_n(ji,jj,jk)
311	END DO
312	END DO
313	END DO
314	!$OMP DO schedule(static) private(jj,ji)
315	DO jj = 1, jpj
316	DO ji = 1, jpi
317	IF( zsum(ji,jj) > 0.e0 ) zsum(ji,jj) = 1.e0 / zsum(ji,jj)
318	END DO
319	END DO
320
321	! ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min)
322	zcoef = rn_tfe_itf / ( rn_tfe * rau0 )
323	!$OMP DO schedule(static) private(jk,jj,ji)
324	DO jk = 1, jpk
325	DO jj = 1, jpj
326	DO ji = 1, jpi
327	zavt_itf(ji,jj,jk) = MIN( 10.e-4, zcoef * en_tmx(ji,jj) * zsum(ji,jj) * zempba_3d(ji,jj,jk) &
328	& / MAX( rn_n2min, rn2(ji,jj,jk) ) * tmask(ji,jj,jk) )
329	END DO
330	END DO
331	END DO
332
333	!$OMP DO schedule(static) private(jj, ji)
334	DO jj = 1, jpj
335	DO ji = 1, jpi
336	zkz(ji,jj) = 0.e0 ! Associated potential energy consummed over the whole water column
337	END DO
338	END DO
339	DO jk = 2, jpkm1
340	!$OMP DO schedule(static) private(jj,ji)
341	DO jj = 1, jpj
342	DO ji = 1, jpi
343	zkz(ji,jj) = zkz(ji,jj) + e3w_n(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zavt_itf(ji,jj,jk) * wmask(ji,jj,jk)
344	END DO
345	END DO
346	END DO
347
348	!$OMP DO schedule(static) private(jj,ji)
349	DO jj = 1, jpj ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
350	DO ji = 1, jpi
351	IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / zkz(ji,jj)
352	END DO
353	END DO
354
355	!$OMP DO schedule(static) private(jk,jj,ji)
356	DO jk = 2, jpkm1 ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
357	DO jj = 1, jpj
358	DO ji = 1, jpi
359	zavt_itf(ji,jj,jk) = zavt_itf(ji,jj,jk) * MIN( zkz(ji,jj), 120./10. ) * wmask(ji,jj,jk) ! kz max = 120 cm2/s
360	END DO
361	END DO
362	END DO
363	!$OMP END DO NOWAIT
364	!$OMP END PARALLEL
365
366	IF( kt == nit000 ) THEN ! diagnose the nergy consumed by zavt_itf
367	ztpc = 0.e0
368	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,ztpc)
369	DO jk= 1, jpk
370	DO jj= 1, jpj
371	DO ji= 1, jpi
372	ztpc = ztpc + e1e2t(ji,jj) * e3w_n(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) &
373	& * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
374	END DO
375	END DO
376	END DO
377	IF( lk_mpp ) CALL mpp_sum( ztpc )
378	ztpc= rau0 * ztpc / ( rn_me * rn_tfe_itf )
379	IF(lwp) WRITE(numout,) ' N Total power consumption by zavt_itf: ztpc = ', ztpc 1.e-12 ,'TW'
380	ENDIF
381
382	! ! Update pav with the ITF mixing coefficient
383	!$OMP PARALLEL DO schedule(static) private(jk)
384	DO jk = 2, jpkm1
385	pav(:,:,jk) = pav (:,:,jk) * ( 1.e0 - mask_itf(:,:) ) &
386	& + zavt_itf(:,:,jk) * mask_itf(:,:)
387	END DO
388	!
389	CALL wrk_dealloc( jpi,jpj, zkz, zsum1 , zsum2 , zsum )
390	CALL wrk_dealloc( jpi,jpj,jpk, zempba_3d_1, zempba_3d_2, zempba_3d, zdn2dz, zavt_itf )
391	!
392	IF( nn_timing == 1 ) CALL timing_stop('tmx_itf')
393	!
394	END SUBROUTINE tmx_itf
395
396
397	SUBROUTINE zdf_tmx_init
398	!!----------------------------------------------------------------------
399	!! * ROUTINE zdf_tmx_init *
400	!!
401	!! ** Purpose : Initialization of the vertical tidal mixing, Reading
402	!! of M2 and K1 tidal energy in nc files
403	!!
404	!! ** Method : - Read the namtmx namelist and check the parameters
405	!!
406	!! - Read the input data in NetCDF files :
407	!! M2 and K1 tidal energy. The total tidal energy, en_tmx,
408	!! is the sum of M2, K1 and S2 energy where S2 is assumed
409	!! to be: S2=(1/2)^2 * M2
410	!! mask_itf, a mask array that determine where substituing
411	!! the standard Simmons et al. (2005) formulation with the
412	!! one of Koch_Larrouy et al. (2007).
413	!!
414	!! - Compute az_tmx, a 3D coefficient that allows to compute
415	!! the standard tidal-induced vertical mixing as follows:
416	!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
417	!! with az_tmx a bottom intensified coefficient is given by:
418	!! az_tmx(z) = en_tmx / ( rau0 * rn_htmx ) * EXP( -(H-z)/rn_htmx )
419	!! / ( 1. - EXP( - H /rn_htmx ) )
420	!! where rn_htmx the characteristic length scale of the bottom
421	!! intensification, en_tmx the tidal energy, and H the ocean depth
422	!!
423	!! ** input : - Namlist namtmx
424	!! - NetCDF file : M2_ORCA2.nc, K1_ORCA2.nc, and mask_itf.nc
425	!!
426	!! ** Action : - Increase by 1 the nstop flag is setting problem encounter
427	!! - defined az_tmx used to compute tidal-induced mixing
428	!!
429	!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
430	!! Koch-Larrouy et al. 2007, GRL.
431	!!----------------------------------------------------------------------
432	INTEGER :: ji, jj, jk ! dummy loop indices
433	INTEGER :: inum ! local integer
434	INTEGER :: ios
435	REAL(wp) :: ztpc, ze_z ! local scalars
436	REAL(wp), DIMENSION(:,:) , POINTER :: zem2, zek1 ! read M2 and K1 tidal energy
437	REAL(wp), DIMENSION(:,:) , POINTER :: zkz ! total M2, K1 and S2 tidal energy
438	REAL(wp), DIMENSION(:,:) , POINTER :: zfact ! used for vertical structure function
439	REAL(wp), DIMENSION(:,:) , POINTER :: zhdep ! Ocean depth
440	REAL(wp), DIMENSION(:,:,:), POINTER :: zpc, zav_tide ! power consumption
441	!!
442	NAMELIST/namzdf_tmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
443	!!----------------------------------------------------------------------
444	!
445	IF( nn_timing == 1 ) CALL timing_start('zdf_tmx_init')
446	!
447	CALL wrk_alloc( jpi,jpj, zem2, zek1, zkz, zfact, zhdep )
448	CALL wrk_alloc( jpi,jpj,jpk, zpc, zav_tide )
449	!
450	REWIND( numnam_ref ) ! Namelist namzdf_tmx in reference namelist : Tidal Mixing
451	READ ( numnam_ref, namzdf_tmx, IOSTAT = ios, ERR = 901)
452	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp )
453	!
454	REWIND( numnam_cfg ) ! Namelist namzdf_tmx in configuration namelist : Tidal Mixing
455	READ ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 )
456	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
457	IF(lwm) WRITE ( numond, namzdf_tmx )
458	!
459	IF(lwp) THEN ! Control print
460	WRITE(numout,*)
461	WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
462	WRITE(numout,*) '~~~~~~~~~~~~'
463	WRITE(numout,*) ' Namelist namzdf_tmx : set tidal mixing parameters'
464	WRITE(numout,*) ' Vertical decay scale for turbulence = ', rn_htmx
465	WRITE(numout,*) ' Brunt-Vaisala frequency threshold = ', rn_n2min
466	WRITE(numout,*) ' Tidal dissipation efficiency = ', rn_tfe
467	WRITE(numout,*) ' Mixing efficiency = ', rn_me
468	WRITE(numout,*) ' ITF specific parameterisation = ', ln_tmx_itf
469	WRITE(numout,*) ' ITF tidal dissipation efficiency = ', rn_tfe_itf
470	ENDIF
471	! ! allocate tmx arrays
472	IF( zdf_tmx_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
473
474	IF( ln_tmx_itf ) THEN ! read the Indonesian Through Flow mask
475	CALL iom_open('mask_itf',inum)
476	CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) !
477	CALL iom_close(inum)
478	ENDIF
479	! ! read M2 tidal energy flux : W/m2 ( zem2 < 0 )
480	CALL iom_open('M2rowdrg',inum)
481	CALL iom_get (inum, jpdom_data, 'field',zem2,1) !
482	CALL iom_close(inum)
483	! ! read K1 tidal energy flux : W/m2 ( zek1 < 0 )
484	CALL iom_open('K1rowdrg',inum)
485	CALL iom_get (inum, jpdom_data, 'field',zek1,1) !
486	CALL iom_close(inum)
487	! ! Total tidal energy ( M2, S2 and K1 with S2=(1/2)^2 * M2 )
488	! ! only the energy available for mixing is taken into account,
489	! ! (mixing efficiency tidal dissipation efficiency)
490	!$OMP PARALLEL
491
492	!$OMP DO schedule(static) private(jj, ji)
493	DO jj = 1, jpj
494	DO ji = 1, jpi
495	en_tmx(ji,jj) = - rn_tfe * rn_me * ( zem2(ji,jj) * 1.25 + zek1(ji,jj) ) * ssmask(ji,jj)
496	END DO
497	END DO
498
499	!============
500	!TG: Bug for VVL? Should this section be moved out of _init and be updated at every timestep?
501	!!gm : you are right, but tidal mixing acts in deep ocean (H>500m) where e3 is O(100m)
502	!! the error is thus ~1% which I feel comfortable with, compared to uncertainties in tidal energy dissipation.
503	! ! Vertical structure (az_tmx)
504	!$OMP DO schedule(static) private(jj, ji)
505	DO jj = 1, jpj ! part independent of the level
506	DO ji = 1, jpi
507	zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean
508	zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
509	IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
510	END DO
511	END DO
512	!$OMP DO schedule(static) private(jk, jj, ji)
513	DO jk= 1, jpk ! complete with the level-dependent part
514	DO jj = 1, jpj
515	DO ji = 1, jpi
516	az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-gdepw_0(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
517	END DO
518	END DO
519	END DO
520	!$OMP END DO NOWAIT
521	!$OMP END PARALLEL
522	!===========
523	!
524	IF( nprint == 1 .AND. lwp ) THEN
525	! Control print
526	! Total power consumption due to vertical mixing
527	! zpc = rau0 * 1/rn_me * rn2 * zav_tide
528	ztpc = 0._wp
529	!$OMP PARALLEL
530	!$OMP DO schedule(static) private(jk, jj, ji)
531	DO jk = 1, jpk
532	DO jj = 1, jpj
533	DO ji = 1, jpi
534	zav_tide(ji,jj,jk) = 0.e0
535	END DO
536	END DO
537	END DO
538	!$OMP DO schedule(static) private(jk)
539	DO jk = 2, jpkm1
540	zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
541	END DO
542	!
543	!$OMP DO schedule(static) private(jk, jj, ji)
544	DO jk= 1, jpk
545	DO jj = 1, jpj
546	DO ji = 1, jpi
547	zpc(ji,jj,jk) = MAX(rn_n2min,rn2(ji,jj,jk)) * zav_tide(ji,jj,jk)
548	END DO
549	END DO
550	END DO
551	!$OMP DO schedule(static) private(jk, jj, ji, ztpc)
552	DO jk= 2, jpkm1
553	DO jj = 1, jpj
554	DO ji = 1, jpi
555	ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
556	END DO
557	END DO
558	END DO
559	!$OMP END DO NOWAIT
560	!$OMP END PARALLEL
561	IF( lk_mpp ) CALL mpp_sum( ztpc )
562	ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
563	!
564	WRITE(numout,*)
565	WRITE(numout,) ' Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc 1.e-12 ,'TW'
566	!
567	! control print 2
568	!$OMP PARALLEL
569	!$OMP DO schedule(static) private(jk, jj, ji)
570	DO jk= 1, jpk
571	DO jj = 1, jpj
572	DO ji = 1, jpi
573	zav_tide(ji,jj,jk) = MIN( zav_tide(ji,jj,jk), 60.e-4 )
574	zkz(ji,jj) = 0._wp
575	END DO
576	END DO
577	END DO
578
579	DO jk = 2, jpkm1
580	!$OMP DO schedule(static) private(jj, ji)
581	DO jj = 1, jpj
582	DO ji = 1, jpi
583	zkz(ji,jj) = zkz(ji,jj) + e3w_n(ji,jj,jk) * MAX(0.e0, rn2(ji,jj,jk)) * rau0 * zav_tide(ji,jj,jk) * wmask(ji,jj,jk)
584	END DO
585	END DO
586	END DO
587	! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
588	!$OMP DO schedule(static) private(jj, ji)
589	DO jj = 1, jpj
590	DO ji = 1, jpi
591	IF( zkz(ji,jj) /= 0.e0 ) THEN
592	zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
593	ENDIF
594	END DO
595	END DO
596	ztpc = 1.e50
597	!$OMP DO schedule(static) private(jj, ji, ztpc)
598	DO jj = 1, jpj
599	DO ji = 1, jpi
600	IF( zkz(ji,jj) /= 0.e0 ) THEN
601	ztpc = Min( zkz(ji,jj), ztpc)
602	ENDIF
603	END DO
604	END DO
605	!$OMP END PARALLEL
606	WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) )
607	!$OMP PARALLEL
608	!
609	!$OMP DO schedule(static) private(jk)
610	DO jk = 2, jpkm1
611	zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) * wmask(:,:,jk) !kz max = 300 cm2/s
612	END DO
613	ztpc = 0._wp
614	!$OMP DO schedule(static) private(jk, jj, ji)
615	DO jk= 1, jpk
616	DO jj = 1, jpj
617	DO ji = 1, jpi
618	zpc(ji,jj,jk) = Max(0.e0,rn2(ji,jj,jk)) * zav_tide(ji,jj,jk)
619	END DO
620	END DO
621	END DO
622	!$OMP DO schedule(static) private(jk, jj, ji, ztpc)
623	DO jk= 1, jpk
624	DO jj = 1, jpj
625	DO ji = 1, jpi
626	ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
627	END DO
628	END DO
629	END DO
630	!$OMP END DO NOWAIT
631	!$OMP END PARALLEL
632	IF( lk_mpp ) CALL mpp_sum( ztpc )
633	ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
634	WRITE(numout,) ' 2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc 1.e-12 ,'TW'
635	!!gm bug mpp in these diagnostics
636	DO jk = 1, jpk
637	ze_z = SUM( e1e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) ) &
638	& / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) )
639	ztpc = 1.e50
640	!$OMP PARALLEL DO schedule(static) private(ztpc, jj, ji)
641	DO jj = 1, jpj
642	DO ji = 1, jpi
643	IF( zav_tide(ji,jj,jk) /= 0.e0 ) ztpc = MIN( ztpc, zav_tide(ji,jj,jk) )
644	END DO
645	END DO
646	WRITE(numout,) ' N2 min - jk= ', jk,' ', ze_z 1.e4,' cm2/s min= ',ztpc*1.e4, &
647	& 'max= ', MAXVAL(zav_tide(:,:,jk) )*1.e4, ' cm2/s'
648	END DO
649
650	WRITE(numout,) ' e_tide : ', SUM( e1e2ten_tmx ) / ( rn_tfe * rn_me ) * 1.e-12, 'TW'
651	WRITE(numout,*)
652	WRITE(numout,*) ' Initial profile of tidal vertical mixing'
653
654	DO jk = 1, jpk
655	!$OMP PARALLEL DO schedule(static) private(jj, ji)
656	DO jj = 1,jpj
657	DO ji = 1,jpi
658	zkz(ji,jj) = az_tmx(ji,jj,jk) /MAX( rn_n2min, rn2(ji,jj,jk) )
659	END DO
660	END DO
661
662	ze_z = SUM( e1e2t(:,:) * zkz (:,:) * tmask_i(:,:) ) &
663	& / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask(:,:,jk) * tmask_i(:,:) ) )
664	WRITE(numout,) ' jk= ', jk,' ', ze_z 1.e4,' cm2/s'
665
666	END DO
667	DO jk = 1, jpk
668	!$OMP PARALLEL DO schedule(static) private(jj, ji)
669	DO jj = 1,jpj
670	DO ji = 1,jpi
671	zkz(ji,jj) = az_tmx(ji,jj,jk) /rn_n2min
672	END DO
673	END DO
674
675	ze_z = SUM( e1e2t(:,:) * zkz (:,:) * tmask_i(:,:) ) &
676	& / MAX( 1.e-20, SUM( e1e2t(:,:) * wmask(:,:,jk) * tmask_i(:,:) ) )
677	WRITE(numout,*)
678	WRITE(numout,) ' N2 min - jk= ', jk,' ', ze_z 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4, &
679	& 'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
680	END DO
681	!!gm end bug mpp
682	!
683	ENDIF
684	!
685	CALL wrk_dealloc( jpi,jpj, zem2, zek1, zkz, zfact, zhdep )
686	CALL wrk_dealloc( jpi,jpj,jpk, zpc )
687	!
688	IF( nn_timing == 1 ) CALL timing_stop('zdf_tmx_init')
689	!
690	END SUBROUTINE zdf_tmx_init
691
692	#elif defined key_zdftmx_new
693	!!----------------------------------------------------------------------
694	!! 'key_zdftmx_new' Internal wave-driven vertical mixing
695	!!----------------------------------------------------------------------
696	!! zdf_tmx : global momentum & tracer Kz with wave induced Kz
697	!! zdf_tmx_init : global momentum & tracer Kz with wave induced Kz
698	!!----------------------------------------------------------------------
699	USE oce ! ocean dynamics and tracers variables
700	USE dom_oce ! ocean space and time domain variables
701	USE zdf_oce ! ocean vertical physics variables
702	USE zdfddm ! ocean vertical physics: double diffusive mixing
703	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
704	USE eosbn2 ! ocean equation of state
705	USE phycst ! physical constants
706	USE prtctl ! Print control
707	USE in_out_manager ! I/O manager
708	USE iom ! I/O Manager
709	USE lib_mpp ! MPP library
710	USE wrk_nemo ! work arrays
711	USE timing ! Timing
712	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
713
714	IMPLICIT NONE
715	PRIVATE
716
717	PUBLIC zdf_tmx ! called in step module
718	PUBLIC zdf_tmx_init ! called in nemogcm module
719	PUBLIC zdf_tmx_alloc ! called in nemogcm module
720
721	LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .TRUE. !: wave-driven mixing flag
722
723	! !!* Namelist namzdf_tmx : internal wave-driven mixing *
724	INTEGER :: nn_zpyc ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
725	LOGICAL :: ln_mevar ! variable (=T) or constant (=F) mixing efficiency
726	LOGICAL :: ln_tsdiff ! account for differential T/S wave-driven mixing (=T) or not (=F)
727
728	REAL(wp) :: r1_6 = 1._wp / 6._wp
729
730	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ebot_tmx ! power available from high-mode wave breaking (W/m2)
731	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: epyc_tmx ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
732	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ecri_tmx ! power available from low-mode, critical slope wave breaking (W/m2)
733	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbot_tmx ! WKB decay scale for high-mode energy dissipation (m)
734	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcri_tmx ! decay scale for low-mode critical slope dissipation (m)
735	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emix_tmx ! local energy density available for mixing (W/kg)
736	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bflx_tmx ! buoyancy flux Kz * N^2 (W/kg)
737	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pcmap_tmx ! vertically integrated buoyancy flux (W/m2)
738	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T)
739	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zav_wave ! Internal wave-induced diffusivity
740
741	!! * Substitutions
742	# include "zdfddm_substitute.h90"
743	# include "vectopt_loop_substitute.h90"
744	!!----------------------------------------------------------------------
745	!! NEMO/OPA 4.0 , NEMO Consortium (2016)
746	!! $Id$
747	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
748	!!----------------------------------------------------------------------
749	CONTAINS
750
751	INTEGER FUNCTION zdf_tmx_alloc()
752	!!----------------------------------------------------------------------
753	!! * FUNCTION zdf_tmx_alloc *
754	!!----------------------------------------------------------------------
755	ALLOCATE( ebot_tmx(jpi,jpj), epyc_tmx(jpi,jpj), ecri_tmx(jpi,jpj) , &
756	& hbot_tmx(jpi,jpj), hcri_tmx(jpi,jpj), emix_tmx(jpi,jpj,jpk), &
757	& bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk), &
758	& zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc )
759	!
760	IF( lk_mpp ) CALL mpp_sum ( zdf_tmx_alloc )
761	IF( zdf_tmx_alloc /= 0 ) CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
762	END FUNCTION zdf_tmx_alloc
763
764
765	SUBROUTINE zdf_tmx( kt )
766	!!----------------------------------------------------------------------
767	!! * ROUTINE zdf_tmx *
768	!!
769	!! ** Purpose : add to the vertical mixing coefficients the effect of
770	!! breaking internal waves.
771	!!
772	!! ** Method : - internal wave-driven vertical mixing is given by:
773	!! Kz_wave = min( 100 cm2/s, f( Reb = emix_tmx /( Nu * N^2 ) )
774	!! where emix_tmx is the 3D space distribution of the wave-breaking
775	!! energy and Nu the molecular kinematic viscosity.
776	!! The function f(Reb) is linear (constant mixing efficiency)
777	!! if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T.
778	!!
779	!! - Compute emix_tmx, the 3D power density that allows to compute
780	!! Reb and therefrom the wave-induced vertical diffusivity.
781	!! This is divided into three components:
782	!! 1. Bottom-intensified low-mode dissipation at critical slopes
783	!! emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx )
784	!! / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx
785	!! where hcri_tmx is the characteristic length scale of the bottom
786	!! intensification, ecri_tmx a map of available power, and H the ocean depth.
787	!! 2. Pycnocline-intensified low-mode dissipation
788	!! emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
789	!! / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
790	!! where epyc_tmx is a map of available power, and nn_zpyc
791	!! is the chosen stratification-dependence of the internal wave
792	!! energy dissipation.
793	!! 3. WKB-height dependent high mode dissipation
794	!! emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx)
795	!! / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) )
796	!! where hbot_tmx is the characteristic length scale of the WKB bottom
797	!! intensification, ebot_tmx is a map of available power, and z_wkb is the
798	!! WKB-stretched height above bottom defined as
799	!! z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) )
800	!! / SUM( sqrt(rn2(z')) * e3w(z') )
801	!!
802	!! - update the model vertical eddy viscosity and diffusivity:
803	!! avt = avt + av_wave
804	!! avm = avm + av_wave
805	!! avmu = avmu + mi(av_wave)
806	!! avmv = avmv + mj(av_wave)
807	!!
808	!! - if namelist parameter ln_tsdiff = T, account for differential mixing:
809	!! avs = avt + av_wave * diffusivity_ratio(Reb)
810	!!
811	!! ** Action : - Define emix_tmx used to compute internal wave-induced mixing
812	!! - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing
813	!!
814	!! References : de Lavergne et al. 2015, JPO; 2016, in prep.
815	!!----------------------------------------------------------------------
816	INTEGER, INTENT(in) :: kt ! ocean time-step
817	!
818	INTEGER :: ji, jj, jk ! dummy loop indices
819	REAL(wp) :: ztpc ! scalar workspace
820	REAL(wp), DIMENSION(:,:) , POINTER :: zfact ! Used for vertical structure
821	REAL(wp), DIMENSION(:,:) , POINTER :: zhdep ! Ocean depth
822	REAL(wp), DIMENSION(:,:,:), POINTER :: zwkb ! WKB-stretched height above bottom
823	REAL(wp), DIMENSION(:,:,:), POINTER :: zweight ! Weight for high mode vertical distribution
824	REAL(wp), DIMENSION(:,:,:), POINTER :: znu_t ! Molecular kinematic viscosity (T grid)
825	REAL(wp), DIMENSION(:,:,:), POINTER :: znu_w ! Molecular kinematic viscosity (W grid)
826	REAL(wp), DIMENSION(:,:,:), POINTER :: zReb ! Turbulence intensity parameter
827	!!----------------------------------------------------------------------
828	!
829	IF( nn_timing == 1 ) CALL timing_start('zdf_tmx')
830	!
831	CALL wrk_alloc( jpi,jpj, zfact, zhdep )
832	CALL wrk_alloc( jpi,jpj,jpk, zwkb, zweight, znu_t, znu_w, zReb )
833
834	! ! ----------------------------- !
835	! ! Internal wave-driven mixing ! (compute zav_wave)
836	! ! ----------------------------- !
837	!
838	! !* Critical slope mixing: distribute energy over the time-varying ocean depth,
839	! using an exponential decay from the seafloor.
840	!$OMP PARALLEL
841	!$OMP DO schedule(static) private(jj,ji)
842	DO jj = 1, jpj ! part independent of the level
843	DO ji = 1, jpi
844	zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean
845	zfact(ji,jj) = rau0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) ) )
846	IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj)
847	END DO
848	END DO
849
850	!$OMP DO schedule(static) private(jk)
851	DO jk = 2, jpkm1 ! complete with the level-dependent part
852	emix_tmx(:,:,jk) = zfact(:,:) * ( EXP( ( gde3w_n(:,:,jk ) - zhdep(:,:) ) / hcri_tmx(:,:) ) &
853	& - EXP( ( gde3w_n(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) ) ) * wmask(:,:,jk) &
854	& / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
855	END DO
856	!$OMP END DO NOWAIT
857	!$OMP END PARALLEL
858
859	! !* Pycnocline-intensified mixing: distribute energy over the time-varying
860	! !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
861
862	SELECT CASE ( nn_zpyc )
863
864	CASE ( 1 ) ! Dissipation scales as N (recommended)
865
866	!$OMP PARALLEL
867	!$OMP DO schedule(static) private(jj, ji)
868	DO jj = 1, jpj
869	DO ji = 1, jpi
870	zfact(ji,jj) = 0._wp
871	END DO
872	END DO
873	DO jk = 2, jpkm1 ! part independent of the level
874	!$OMP DO schedule(static) private(jj,ji)
875	DO jj = 1, jpj ! part independent of the level
876	DO ji = 1, jpi
877	zfact(ji,jj) = zfact(ji,jj) + e3w_n(ji,jj,jk) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk)
878	END DO
879	END DO
880	END DO
881
882	!$OMP DO schedule(static) private(jj,ji)
883	DO jj = 1, jpj
884	DO ji = 1, jpi
885	IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
886	END DO
887	END DO
888
889	!$OMP DO schedule(static) private(jk)
890	DO jk = 2, jpkm1 ! complete with the level-dependent part
891	emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk)
892	END DO
893	!$OMP END DO NOWAIT
894	!$OMP END PARALLEL
895
896	CASE ( 2 ) ! Dissipation scales as N^2
897
898	!$OMP PARALLEL
899	!$OMP DO schedule(static) private(jj, ji)
900	DO jj = 1, jpj
901	DO ji = 1, jpi
902	zfact(ji,jj) = 0._wp
903	END DO
904	END DO
905
906	DO jk = 2, jpkm1
907	!$OMP DO schedule(static) private(jj,ji)
908	DO jj = 1, jpj
909	DO ji = 1, jpi
910	zfact(ji,jj) = zfact(ji,jj) + e3w_n(ji,jj,jk) * MAX( 0._wp, rn2(ji,jj,jk) ) * wmask(ji,jj,jk)
911	END DO
912	END DO
913	END DO
914
915	!$OMP DO schedule(static) private(jj,ji)
916	DO jj= 1, jpj
917	DO ji = 1, jpi
918	IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
919	END DO
920	END DO
921
922	!$OMP DO schedule(static) private(jk)
923	DO jk = 2, jpkm1 ! complete with the level-dependent part
924	emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
925	END DO
926	!$OMP END DO NOWAIT
927	!$OMP END PARALLEL
928
929	END SELECT
930
931	! !* WKB-height dependent mixing: distribute energy over the time-varying
932	! !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
933
934	!$OMP PARALLEL
935	!$OMP DO schedule(static) private(jk,jj,ji)
936	DO jk = 1, jpk
937	DO jj = 1, jpj
938	DO ji = 1, jpi
939	zwkb(ji,jj,jk) = 0._wp
940	END DO
941	END DO
942	END DO
943	!$OMP DO schedule(static) private(jj,ji)
944	DO jj = 1, jpj
945	DO ji = 1, jpi
946	zfact(ji,jj) = 0._wp
947	END DO
948	END DO
949	DO jk = 2, jpkm1
950	!$OMP DO schedule(static) private(jj,ji)
951	DO jj = 1, jpj
952	DO ji = 1, jpi
953	zfact(ji,jj) = zfact(ji,jj) + e3w_n(ji,jj,jk) * SQRT( MAX( 0._wp, rn2(ji,jj,jk) ) ) * wmask(ji,jj,jk)
954	zwkb(ji,jj,jk) = zfact(ji,jj)
955	END DO
956	END DO
957	END DO
958
959	!$OMP DO schedule(static) private(jk,jj,ji)
960	DO jk = 2, jpkm1
961	DO jj = 1, jpj
962	DO ji = 1, jpi
963	IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) &
964	& * tmask(ji,jj,jk) / zfact(ji,jj)
965	END DO
966	END DO
967	END DO
968
969	!$OMP DO schedule(static) private(jj, ji)
970	DO jj = 1, jpj
971	DO ji = 1, jpi
972	zwkb(ji,jj,1) = zhdep(ji,jj) * tmask(ji,jj,1)
973	END DO
974	END DO
975	!$OMP END DO NOWAIT
976	!$OMP DO schedule(static) private(jk,jj,ji)
977	DO jk = 1, jpk
978	DO jj = 1, jpj
979	DO ji = 1, jpi
980	zweight(ji,jj,jk) = 0._wp
981	END DO
982	END DO
983	END DO
984
985	!$OMP DO schedule(static) private(jk)
986	DO jk = 2, jpkm1
987	zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk) &
988	& * ( EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) ) )
989	END DO
990
991	!$OMP DO schedule(static) private(jj, ji)
992	DO jj = 1, jpj
993	DO ji = 1, jpi
994	zfact(ji,jj) = 0._wp
995	END DO
996	END DO
997
998	DO jk = 2, jpkm1 ! part independent of the level
999	!$OMP DO schedule(static) private(jj,ji)
1000	DO jj = 1, jpj
1001	DO ji = 1, jpi
1002	zfact(ji,jj) = zfact(ji,jj) + zweight(ji,jj,jk)
1003	END DO
1004	END DO
1005	END DO
1006
1007	!$OMP DO schedule(static) private(jj,ji)
1008	DO jj = 1, jpj
1009	DO ji = 1, jpi
1010	IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
1011	END DO
1012	END DO
1013
1014	!$OMP DO schedule(static) private(jk)
1015	DO jk = 2, jpkm1 ! complete with the level-dependent part
1016	emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk) &
1017	& / ( gde3w_n(:,:,jk) - gde3w_n(:,:,jk-1) )
1018	END DO
1019	!$OMP END DO NOWAIT
1020
1021
1022	! Calculate molecular kinematic viscosity
1023	!$OMP DO schedule(static) private(jj, ji)
1024	DO jj = 1, jpj
1025	DO ji = 1, jpi
1026	znu_t(ji,jj,jk) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * tsn(ji,jj,jk,jp_tem) &
1027	& + 0.00694_wp * tsn(ji,jj,jk,jp_tem) * tsn(ji,jj,jk,jp_tem) &
1028	& + 0.02305_wp * tsn(ji,jj,jk,jp_sal) ) * tmask(ji,jj,jk) * r1_rau0
1029	END DO
1030	END DO
1031	!$OMP DO schedule(static) private(jk)
1032	DO jk = 2, jpkm1
1033	znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
1034	END DO
1035
1036	! Calculate turbulence intensity parameter Reb
1037	!$OMP DO schedule(static) private(jk)
1038	DO jk = 2, jpkm1
1039	zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
1040	END DO
1041
1042	! Define internal wave-induced diffusivity
1043	!$OMP DO schedule(static) private(jk)
1044	DO jk = 2, jpkm1
1045	zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6 ! This corresponds to a constant mixing efficiency of 1/6
1046	END DO
1047	!$OMP END DO NOWAIT
1048	!$OMP END PARALLEL
1049
1050	IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the
1051	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
1052	DO jk = 2, jpkm1 ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
1053	DO jj = 1, jpj
1054	DO ji = 1, jpi
1055	IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
1056	zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
1057	ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
1058	zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
1059	ENDIF
1060	END DO
1061	END DO
1062	END DO
1063	ENDIF
1064
1065	!$OMP PARALLEL DO schedule(static) private(jk)
1066	DO jk = 2, jpkm1 ! Bound diffusivity by molecular value and 100 cm2/s
1067	zav_wave(:,:,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp ) * wmask(:,:,jk)
1068	END DO
1069
1070	IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave
1071	ztpc = 0._wp
1072	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji,ztpc)
1073	DO jk = 2, jpkm1
1074	DO jj = 1, jpj
1075	DO ji = 1, jpi
1076	ztpc = ztpc + e3w_n(ji,jj,jk) * e1e2t(ji,jj) &
1077	& * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
1078	END DO
1079	END DO
1080	END DO
1081	IF( lk_mpp ) CALL mpp_sum( ztpc )
1082	ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing
1083
1084	IF(lwp) THEN
1085	WRITE(numout,*)
1086	WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)'
1087	WRITE(numout,*) '~~~~~~~ '
1088	WRITE(numout,*)
1089	WRITE(numout,) ' Total power consumption by av_wave: ztpc = ', ztpc 1.e-12_wp, 'TW'
1090	ENDIF
1091	ENDIF
1092
1093	! ! ----------------------- !
1094	! ! Update mixing coefs !
1095	! ! ----------------------- !
1096	!
1097	IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature
1098	!$OMP PARALLEL
1099	!$OMP DO schedule(static) private(jk,jj,ji)
1100	DO jk = 2, jpkm1 ! Calculate S/T diffusivity ratio as a function of Reb
1101	DO jj = 1, jpj
1102	DO ji = 1, jpi
1103	zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp * &
1104	& TANH( 0.92_wp * ( LOG10( MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 ) ) - 0.60_wp ) ) &
1105	& ) * wmask(ji,jj,jk)
1106	END DO
1107	END DO
1108	END DO
1109	!$OMP DO schedule(static) private(jk)
1110	DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with wave-driven mixing
1111	fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)
1112	avt (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
1113	avm (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
1114	END DO
1115	!$OMP END DO NOWAIT
1116	!$OMP END PARALLEL
1117	CALL iom_put( "av_ratio", zav_ratio )
1118	!
1119	ELSE !* update momentum & tracer diffusivity with wave-driven mixing
1120	!$OMP PARALLEL DO schedule(static) private(jk)
1121	DO jk = 2, jpkm1
1122	fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
1123	avt (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
1124	avm (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
1125	END DO
1126	ENDIF
1127
1128	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
1129	DO jk = 2, jpkm1 !* update momentum diffusivity at wu and wv points
1130	DO jj = 2, jpjm1
1131	DO ji = fs_2, fs_jpim1 ! vector opt.
1132	avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj ,jk) ) * wumask(ji,jj,jk)
1133	avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji ,jj+1,jk) ) * wvmask(ji,jj,jk)
1134	END DO
1135	END DO
1136	END DO
1137	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! lateral boundary condition
1138
1139	! !* output internal wave-driven mixing coefficient
1140	CALL iom_put( "av_wave", zav_wave )
1141	!* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx),
1142	! vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx)
1143	IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN
1144	!$OMP PARALLEL
1145	!$OMP DO schedule(static) private(jk,jj,ji)
1146	DO jk = 1, jpk
1147	DO jj = 1, jpj
1148	DO ji = 1, jpi
1149	bflx_tmx(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk)
1150	END DO
1151	END DO
1152	END DO
1153	!$OMP END DO NOWAIT
1154	!$OMP DO schedule(static) private(jj, ji)
1155	DO jj = 1, jpj
1156	DO ji = 1, jpi
1157	pcmap_tmx(ji,jj) = 0._wp
1158	END DO
1159	END DO
1160	DO jk = 2, jpkm1
1161	!$OMP DO schedule(static) private(jj, ji)
1162	DO jj = 1, jpj
1163	DO ji = 1, jpi
1164	pcmap_tmx(ji,jj) = pcmap_tmx(ji,jj) + e3w_n(ji,jj,jk) * bflx_tmx(ji,jj,jk) * wmask(ji,jj,jk)
1165	END DO
1166	END DO
1167	END DO
1168	!$OMP DO schedule(static) private(jj, ji)
1169	DO jj = 1, jpj
1170	DO ji = 1, jpi
1171	pcmap_tmx(ji,jj) = rau0 * pcmap_tmx(ji,jj)
1172	END DO
1173	END DO
1174	!$OMP END PARALLEL
1175	CALL iom_put( "bflx_tmx", bflx_tmx )
1176	CALL iom_put( "pcmap_tmx", pcmap_tmx )
1177	ENDIF
1178	CALL iom_put( "bn2", rn2 )
1179	CALL iom_put( "emix_tmx", emix_tmx )
1180
1181	CALL wrk_dealloc( jpi,jpj, zfact, zhdep )
1182	CALL wrk_dealloc( jpi,jpj,jpk, zwkb, zweight, znu_t, znu_w, zReb )
1183
1184	IF(ln_ctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
1185	!
1186	IF( nn_timing == 1 ) CALL timing_stop('zdf_tmx')
1187	!
1188	END SUBROUTINE zdf_tmx
1189
1190
1191	SUBROUTINE zdf_tmx_init
1192	!!----------------------------------------------------------------------
1193	!! * ROUTINE zdf_tmx_init *
1194	!!
1195	!! ** Purpose : Initialization of the wave-driven vertical mixing, reading
1196	!! of input power maps and decay length scales in netcdf files.
1197	!!
1198	!! ** Method : - Read the namzdf_tmx namelist and check the parameters
1199	!!
1200	!! - Read the input data in NetCDF files :
1201	!! power available from high-mode wave breaking (mixing_power_bot.nc)
1202	!! power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc)
1203	!! power available from critical slope wave-breaking (mixing_power_cri.nc)
1204	!! WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc)
1205	!! decay scale for critical slope wave-breaking (decay_scale_cri.nc)
1206	!!
1207	!! ** input : - Namlist namzdf_tmx
1208	!! - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc,
1209	!! decay_scale_bot.nc decay_scale_cri.nc
1210	!!
1211	!! ** Action : - Increase by 1 the nstop flag is setting problem encounter
1212	!! - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx
1213	!!
1214	!! References : de Lavergne et al. 2015, JPO; 2016, in prep.
1215	!!
1216	!!----------------------------------------------------------------------
1217	INTEGER :: ji, jj, jk ! dummy loop indices
1218	INTEGER :: inum ! local integer
1219	INTEGER :: ios
1220	REAL(wp) :: zbot, zpyc, zcri ! local scalars
1221	!!
1222	NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff
1223	!!----------------------------------------------------------------------
1224	!
1225	IF( nn_timing == 1 ) CALL timing_start('zdf_tmx_init')
1226	!
1227	REWIND( numnam_ref ) ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing
1228	READ ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901)
1229	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp )
1230	!
1231	REWIND( numnam_cfg ) ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing
1232	READ ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 )
1233	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
1234	IF(lwm) WRITE ( numond, namzdf_tmx_new )
1235	!
1236	IF(lwp) THEN ! Control print
1237	WRITE(numout,*)
1238	WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing'
1239	WRITE(numout,*) '~~~~~~~~~~~~'
1240	WRITE(numout,*) ' Namelist namzdf_tmx_new : set wave-driven mixing parameters'
1241	WRITE(numout,*) ' Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc
1242	WRITE(numout,*) ' Variable (T) or constant (F) mixing efficiency = ', ln_mevar
1243	WRITE(numout,*) ' Differential internal wave-driven mixing (T) or not (F) = ', ln_tsdiff
1244	ENDIF
1245
1246	! The new wave-driven mixing parameterization elevates avt and avm in the interior, and
1247	! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should
1248	! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
1249	avmb(:) = 1.4e-6_wp ! viscous molecular value
1250	avtb(:) = 1.e-10_wp ! very small diffusive minimum (background avt is specified in zdf_tmx)
1251
1252	!$OMP PARALLEL DO schedule(static) private(jj, ji)
1253	DO jj = 1, jpj
1254	DO ji = 1, jpi
1255	avtb_2d(ji,jj) = 1.e0_wp ! uniform
1256	END DO
1257	END DO
1258
1259	IF(lwp) THEN ! Control print
1260	WRITE(numout,*)
1261	WRITE(numout,*) ' Force the background value applied to avm & avt in TKE to be everywhere ', &
1262	& 'the viscous molecular value & a very small diffusive value, resp.'
1263	ENDIF
1264
1265	IF( .NOT.lk_zdfddm ) CALL ctl_stop( 'STOP', 'zdf_tmx_init_new : key_zdftmx_new requires key_zdfddm' )
1266
1267	! ! allocate tmx arrays
1268	IF( zdf_tmx_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
1269	!
1270	! ! read necessary fields
1271	CALL iom_open('mixing_power_bot',inum) ! energy flux for high-mode wave breaking [W/m2]
1272	CALL iom_get (inum, jpdom_data, 'field', ebot_tmx, 1 )
1273	CALL iom_close(inum)
1274	!
1275	CALL iom_open('mixing_power_pyc',inum) ! energy flux for pynocline-intensified wave breaking [W/m2]
1276	CALL iom_get (inum, jpdom_data, 'field', epyc_tmx, 1 )
1277	CALL iom_close(inum)
1278	!
1279	CALL iom_open('mixing_power_cri',inum) ! energy flux for critical slope wave breaking [W/m2]
1280	CALL iom_get (inum, jpdom_data, 'field', ecri_tmx, 1 )
1281	CALL iom_close(inum)
1282	!
1283	CALL iom_open('decay_scale_bot',inum) ! spatially variable decay scale for high-mode wave breaking [m]
1284	CALL iom_get (inum, jpdom_data, 'field', hbot_tmx, 1 )
1285	CALL iom_close(inum)
1286	!
1287	CALL iom_open('decay_scale_cri',inum) ! spatially variable decay scale for critical slope wave breaking [m]
1288	CALL iom_get (inum, jpdom_data, 'field', hcri_tmx, 1 )
1289	CALL iom_close(inum)
1290
1291	!$OMP PARALLEL DO schedule(static) private(jj, ji)
1292	DO jj = 1, jpj
1293	DO ji = 1, jpi
1294	ebot_tmx(ji,jj) = ebot_tmx(ji,jj) * ssmask(ji,jj)
1295	epyc_tmx(ji,jj) = epyc_tmx(ji,jj) * ssmask(ji,jj)
1296	ecri_tmx(ji,jj) = ecri_tmx(ji,jj) * ssmask(ji,jj)
1297
1298	! Set once for all to zero the first and last vertical levels of appropriate variables
1299	emix_tmx (ji,jj, 1 ) = 0._wp
1300	emix_tmx (ji,jj,jpk) = 0._wp
1301	zav_ratio(ji,jj, 1 ) = 0._wp
1302	zav_ratio(ji,jj,jpk) = 0._wp
1303	zav_wave (ji,jj, 1 ) = 0._wp
1304	zav_wave (ji,jj,jpk) = 0._wp
1305	END DO
1306	END DO
1307
1308	zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) )
1309	zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) )
1310	zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) )
1311	IF(lwp) THEN
1312	WRITE(numout,) ' High-mode wave-breaking energy: ', zbot 1.e-12_wp, 'TW'
1313	WRITE(numout,) ' Pycnocline-intensifed wave-breaking energy: ', zpyc 1.e-12_wp, 'TW'
1314	WRITE(numout,) ' Critical slope wave-breaking energy: ', zcri 1.e-12_wp, 'TW'
1315	ENDIF
1316	!
1317	IF( nn_timing == 1 ) CALL timing_stop('zdf_tmx_init')
1318	!
1319	END SUBROUTINE zdf_tmx_init
1320
1321	#else
1322	!!----------------------------------------------------------------------
1323	!! Default option Dummy module NO Tidal MiXing
1324	!!----------------------------------------------------------------------
1325	LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .FALSE. !: tidal mixing flag
1326	CONTAINS
1327	SUBROUTINE zdf_tmx_init ! Dummy routine
1328	WRITE(,) 'zdf_tmx: You should not have seen this print! error?'
1329	END SUBROUTINE zdf_tmx_init
1330	SUBROUTINE zdf_tmx( kt ) ! Dummy routine
1331	WRITE(,) 'zdf_tmx: You should not have seen this print! error?', kt
1332	END SUBROUTINE zdf_tmx
1333	#endif
1334
1335	!!======================================================================
1336	END MODULE zdftmx

Note: See TracBrowser for help on using the repository browser.

Download in other formats: