Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

zdftmx.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 51.3 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90 @ 7910

Download in other formats: