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

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

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

Download in other formats: