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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90 @ 7698

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

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

Download in other formats: