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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90 @ 3270

Last change on this file since 3270 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 32.2 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	!!----------------------------------------------------------------------
[2715]	10	#if defined key_zdftmx \|\| defined key_esopa
[1418]	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
[1496]	24	USE in_out_manager ! I/O manager
	25	USE iom ! I/O Manager
[2715]	26	USE lib_mpp ! MPP library
	27	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
[1418]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[2528]	32	PUBLIC zdf_tmx ! called in step module
	33	PUBLIC zdf_tmx_init ! called in opa module
[2715]	34	PUBLIC zdf_tmx_alloc ! called in nemogcm module
[1418]	35
	36	LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .TRUE. !: tidal mixing flag
	37
[1601]	38	! !!* Namelist namzdf_tmx : tidal mixing *
[1518]	39	REAL(wp) :: rn_htmx = 500. ! vertical decay scale for turbulence (meters)
	40	REAL(wp) :: rn_n2min = 1.e-8 ! threshold of the Brunt-Vaisala frequency (s-1)
	41	REAL(wp) :: rn_tfe = 1./3. ! tidal dissipation efficiency (St Laurent et al. 2002)
	42	REAL(wp) :: rn_me = 0.2 ! mixing efficiency (Osborn 1980)
	43	LOGICAL :: ln_tmx_itf = .TRUE. ! Indonesian Through Flow (ITF): Koch-Larrouy et al. (2007) parameterization
	44	REAL(wp) :: rn_tfe_itf = 1. ! ITF tidal dissipation efficiency (St Laurent et al. 2002)
[1418]	45
[2715]	46	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: en_tmx ! energy available for tidal mixing (W/m2)
	47	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: mask_itf ! mask to use over Indonesian area
	48	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: az_tmx ! coefficient used to evaluate the tidal induced Kz
[1418]	49
[3211]	50	!! * Control permutation of array indices
	51	# include "oce_ftrans.h90"
	52	# include "dom_oce_ftrans.h90"
	53	# include "zdf_oce_ftrans.h90"
	54	!FTRANS az_tmx :I :I :z
	55
[1418]	56	!! * Substitutions
	57	# include "domzgr_substitute.h90"
	58	# include "vectopt_loop_substitute.h90"
	59	!!----------------------------------------------------------------------
[2715]	60	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2528]	61	!! $Id$
[2715]	62	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1418]	63	!!----------------------------------------------------------------------
	64	CONTAINS
	65
[2715]	66	INTEGER FUNCTION zdf_tmx_alloc()
	67	!!----------------------------------------------------------------------
	68	!! * FUNCTION zdf_tmx_alloc *
	69	!!----------------------------------------------------------------------
	70	ALLOCATE(en_tmx(jpi,jpj), mask_itf(jpi,jpj), az_tmx(jpi,jpj,jpk), STAT=zdf_tmx_alloc )
	71	!
	72	IF( lk_mpp ) CALL mpp_sum ( zdf_tmx_alloc )
	73	IF( zdf_tmx_alloc /= 0 ) CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
	74	END FUNCTION zdf_tmx_alloc
	75
	76
[1418]	77	SUBROUTINE zdf_tmx( kt )
	78	!!----------------------------------------------------------------------
	79	!! * ROUTINE zdf_tmx *
	80	!!
	81	!! ** Purpose : add to the vertical mixing coefficients the effect of
[1496]	82	!! tidal mixing (Simmons et al 2004).
[1418]	83	!!
	84	!! ** Method : - tidal-induced vertical mixing is given by:
[1496]	85	!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
	86	!! where az_tmx is a coefficient that specified the 3D space
	87	!! distribution of the faction of tidal energy taht is used
	88	!! for mixing. Its expression is set in zdf_tmx_init routine,
	89	!! following Simmons et al. 2004.
	90	!! NB: a specific bounding procedure is performed on av_tide
	91	!! so that the input tidal energy is actually almost used. The
	92	!! basic maximum value is 60 cm2/s, but values of 300 cm2/s
	93	!! can be reached in area where bottom stratification is too
	94	!! weak.
[1418]	95	!!
[1496]	96	!! - update av_tide in the Indonesian Through Flow area
	97	!! following Koch-Larrouy et al. (2007) parameterisation
	98	!! (see tmx_itf routine).
[1418]	99	!!
[1496]	100	!! - update the model vertical eddy viscosity and diffusivity:
	101	!! avt = avt + av_tides
[1527]	102	!! avm = avm + av_tides
[1496]	103	!! avmu = avmu + mi(av_tides)
	104	!! avmv = avmv + mj(av_tides)
	105	!!
[1527]	106	!! ** Action : avt, avm, avmu, avmv increased by tidal mixing
[1496]	107	!!
[1418]	108	!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
[1496]	109	!! Koch-Larrouy et al. 2007, GRL.
[1418]	110	!!----------------------------------------------------------------------
[1601]	111	USE oce, zav_tide => ua ! use ua as workspace
[2715]	112	USE wrk_nemo, ONLY: zkz => wrk_2d_1
[1546]	113	!!
[1418]	114	INTEGER, INTENT(in) :: kt ! ocean time-step
	115	!!
	116	INTEGER :: ji, jj, jk ! dummy loop indices
	117	REAL(wp) :: ztpc ! scalar workspace
[3211]	118	#if defined key_z_first
	119	REAL(wp) :: ztpc ! scalar workspace
	120	#endif
[1418]	121	!!----------------------------------------------------------------------
	122
[2715]	123	IF(wrk_in_use(2, 1))THEN
	124	CALL ctl_stop('zdf_tmx : requested workspace array unavailable.') ; RETURN
	125	END IF
[1546]	126	! ! ----------------------- !
	127	! ! Standard tidal mixing ! (compute zav_tide)
	128	! ! ----------------------- !
[1496]	129	! !* First estimation (with n2 bound by rn_n2min) bounded by 60 cm2/s
[1546]	130	zav_tide(:,:,:) = MIN( 60.e-4, az_tmx(:,:,:) / MAX( rn_n2min, rn2(:,:,:) ) )
[1418]	131
[1496]	132	zkz(:,:) = 0.e0 !* Associated potential energy consummed over the whole water column
[1418]	133	DO jk = 2, jpkm1
[1546]	134	zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk)
[1418]	135	END DO
	136
[1496]	137	DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
[1418]	138	DO ji = 1, jpi
	139	IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
	140	END DO
	141	END DO
	142
[3211]	143	#if defined key_z_first
	144	DO jj = 1, jpj
	145	DO ji = 1, jpi
	146	zscal = MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s
	147	DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s
	148	zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * zscal
	149	END DO
	150	END DO
	151	END DO
	152	#else
[1546]	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	zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) !kz max = 300 cm2/s
[1418]	155	END DO
[3211]	156	#endif
[1418]	157
[1546]	158	IF( kt == nit000 ) THEN !* check at first time-step: diagnose the energy consumed by zav_tide
[1418]	159	ztpc = 0.e0
[3211]	160	#if defined key_z_first
	161	DO jj = 1, jpj
	162	DO ji = 1, jpi
	163	DO jk = 1, jpk
	164	#else
[1418]	165	DO jk= 1, jpk
	166	DO jj= 1, jpj
	167	DO ji= 1, jpi
[3211]	168	#endif
[1496]	169	ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) &
[1546]	170	& * MAX( 0.e0, rn2(ji,jj,jk) ) * zav_tide(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
[1418]	171	END DO
	172	END DO
	173	END DO
[1495]	174	ztpc= rau0 / ( rn_tfe * rn_me ) * ztpc
[1418]	175	IF(lwp) WRITE(numout,*)
[1496]	176	IF(lwp) WRITE(numout,) ' N Total power consumption by av_tide : ztpc = ', ztpc 1.e-12 ,'TW'
[1418]	177	ENDIF
[1495]	178
[1546]	179	! ! ----------------------- !
	180	! ! ITF tidal mixing ! (update zav_tide)
	181	! ! ----------------------- !
	182	IF( ln_tmx_itf ) CALL tmx_itf( kt, zav_tide )
[1418]	183
[1546]	184	! ! ----------------------- !
	185	! ! Update mixing coefs !
	186	! ! ----------------------- !
[3211]	187	#if defined key_z_first
	188	!* update momentum & tracer diffusivity with tidal mixing
	189	DO jj = 1, jpj
	190	DO ji = 1, jpi
	191	DO jk = 2, jpkm1
	192	avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk)
	193	avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk)
	194	END DO
	195	END DO
	196	END DO
	197	DO jj = 2, jpjm1
	198	DO ji = 2, fpim1
	199	DO jk = 2, jpkm1
	200	avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * umask(ji,jj,jk)
	201	avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * vmask(ji,jj,jk)
	202	END DO
	203	END DO
	204	END DO
	205	#else
[1495]	206	DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing
[1546]	207	avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk)
	208	avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk)
[1418]	209	DO jj = 2, jpjm1
	210	DO ji = fs_2, fs_jpim1 ! vector opt.
[1546]	211	avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * umask(ji,jj,jk)
	212	avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * vmask(ji,jj,jk)
[1418]	213	END DO
	214	END DO
	215	END DO
[3211]	216	#endif
[1496]	217	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! lateral boundary condition
[1418]	218
[1546]	219	! !* output tidal mixing coefficient
	220	CALL iom_put( "av_tide", zav_tide )
	221
	222	IF(ln_ctl) CALL prt_ctl(tab3d_1=zav_tide , clinfo1=' tmx - av_tide: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
[1418]	223	!
[2715]	224	IF(wrk_not_released(2, 1))THEN
	225	CALL ctl_stop('zdf_tmx : failed to release workspace array.')
	226	END IF
	227	!
[1418]	228	END SUBROUTINE zdf_tmx
	229
	230
[1546]	231	SUBROUTINE tmx_itf( kt, pav )
[1418]	232	!!----------------------------------------------------------------------
	233	!! * ROUTINE tmx_itf *
	234	!!
[1496]	235	!! ** Purpose : modify the vertical eddy diffusivity coefficients
[1546]	236	!! (pav) in the Indonesian Through Flow area (ITF).
[1418]	237	!!
[1496]	238	!! ** Method : - Following Koch-Larrouy et al. (2007), in the ITF defined
	239	!! by msk_itf (read in a file, see tmx_init), the tidal
	240	!! mixing coefficient is computed with :
	241	!! * q=1 (i.e. all the tidal energy remains trapped in
	242	!! the area and thus is used for mixing)
	243	!! * the vertical distribution of the tifal energy is a
	244	!! proportional to N above the thermocline (d(N^2)/dz > 0)
	245	!! and to N^2 below the thermocline (d(N^2)/dz < 0)
[1418]	246	!!
[1496]	247	!! ** Action : av_tide updated in the ITF area (msk_itf)
[1418]	248	!!
	249	!! References : Koch-Larrouy et al. 2007, GRL
	250	!!----------------------------------------------------------------------
[2715]	251	USE wrk_nemo, ONLY: zkz => wrk_2d_5
	252	USE wrk_nemo, ONLY: zsum1 => wrk_2d_2, zsum2 => wrk_2d_3, zsum => wrk_2d_4
	253	USE wrk_nemo, ONLY: zempba_3d_1 => wrk_3d_1, zempba_3d_2 => wrk_3d_2
	254	USE wrk_nemo, ONLY: zempba_3d => wrk_3d_3, zdn2dz => wrk_3d_4
	255	USE wrk_nemo, ONLY: zavt_itf => wrk_3d_5
[3211]	256	!! DCSE_NEMO: need additional directives for renamed module variables
	257	!FTRANS zempba_3d_1 zempba_3d_2 zempba_3d zdn2dz zavt_itf :I :I :z
[2715]	258	!!
[1546]	259	INTEGER , INTENT(in ) :: kt ! ocean time-step
[3211]	260
	261	!! DCSE_NEMO: This style defeats ftrans
	262	! REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pav ! Tidal mixing coef.
	263	!FTRANS pav :I :I :z
	264	REAL(wp), INTENT(inout) :: pav(jpi,jpj,jpk) ! Tidal mixing coef.
[1418]	265	!!
[1495]	266	INTEGER :: ji, jj, jk ! dummy loop indices
	267	REAL(wp) :: zcoef, ztpc ! temporary scalar
[1418]	268	!!----------------------------------------------------------------------
[2715]	269	!
	270	IF( wrk_in_use(2, 2,3,4,5) .OR. wrk_in_use(3, 1,2,3,4,5) )THEN
	271	CALL ctl_stop('tmx_itf : requested workspace arrays unavailable.')
	272	RETURN
	273	END IF
[1418]	274	! ! compute the form function using N2 at each time step
[3211]	275	#if defined key_z_first
	276	DO jj = 1, jpj
	277	DO ji = 1, jpi
	278	DO jk = 1, jpkm1
	279	zdn2dz (ji,jj,jk) = rn2(ji,jj,jk) - rn2(ji,jj,jk+1) ! Vertical profile of dN2/dz
	280	zempba_3d_1(ji,jj,jk) = SQRT( MAX( 0.e0, rn2(ji,jj,jk) ) ) ! - - of N
	281	zempba_3d_2(ji,jj,jk) = MAX( 0.e0, rn2(ji,jj,jk) ) ! - - of N^2
	282	END DO
	283	zempba_3d_1(ji,jj,jpk) = 0.e0
	284	zempba_3d_2(ji,jj,jpk) = 0.e0
	285	END DO
	286	END DO
	287	#else
[1518]	288	zempba_3d_1(:,:,jpk) = 0.e0
	289	zempba_3d_2(:,:,jpk) = 0.e0
	290	DO jk = 1, jpkm1
	291	zdn2dz (:,:,jk) = rn2(:,:,jk) - rn2(:,:,jk+1) ! Vertical profile of dN2/dz
[1495]	292	!CDIR NOVERRCHK
[1518]	293	zempba_3d_1(:,:,jk) = SQRT( MAX( 0.e0, rn2(:,:,jk) ) ) ! - - of N
	294	zempba_3d_2(:,:,jk) = MAX( 0.e0, rn2(:,:,jk) ) ! - - of N^2
[1418]	295	END DO
[3211]	296	#endif
[1518]	297	!
[3211]	298	#if defined key_z_first
	299	DO jj = 1, jpj
	300	DO ji = 1, jpj
	301	zsum1(ji,jj) = 0.e0
	302	zsum2(ji,jj) = 0.e0
	303	DO jk= 2, jpk
	304	zsum1(ji,jj) = zsum1(ji,jj) + zempba_3d_1(ji,jj,jk) * fse3w(ji,jj,jk)
	305	zsum2(ji,jj) = zsum2(ji,jj) + zempba_3d_2(ji,jj,jk) * fse3w(ji,jj,jk)
	306	END DO
	307	IF( zsum1(ji,jj) /= 0.e0 ) zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
	308	IF( zsum2(ji,jj) /= 0.e0 ) zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)
	309	END DO
	310	END DO
	311	#else
[1518]	312	zsum1(:,:) = 0.e0
	313	zsum2(:,:) = 0.e0
[1418]	314	DO jk= 2, jpk
[1518]	315	zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk)
	316	zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk)
[1418]	317	END DO
	318	DO jj = 1, jpj
[1518]	319	DO ji = 1, jpi
	320	IF( zsum1(ji,jj) /= 0.e0 ) zsum1(ji,jj) = 1.e0 / zsum1(ji,jj)
	321	IF( zsum2(ji,jj) /= 0.e0 ) zsum2(ji,jj) = 1.e0 / zsum2(ji,jj)
	322	END DO
[1418]	323	END DO
[3211]	324	#endif
[1418]	325
[3211]	326	zsum (:,:) = 0.e0
	327
	328	#if defined key_z_first
	329	DO jj = 1, jpj
	330	DO ji = 1, jpi
	331	DO jk = 1, jpk
	332	#else
	333	DO jk = 1, jpk
[1418]	334	DO jj = 1, jpj
	335	DO ji = 1, jpi
[3211]	336	#endif
[1518]	337	zcoef = 0.5 - SIGN( 0.5, zdn2dz(ji,jj,jk) ) ! =0 if dN2/dz > 0, =1 otherwise
	338	ztpc = zempba_3d_1(ji,jj,jk) * zsum1(ji,jj) * zcoef &
	339	& + zempba_3d_2(ji,jj,jk) * zsum2(ji,jj) * ( 1. - zcoef )
	340	!
	341	zempba_3d(ji,jj,jk) = ztpc
	342	zsum (ji,jj) = zsum(ji,jj) + ztpc * fse3w(ji,jj,jk)
[1418]	343	END DO
[3211]	344	#if !defined key_z_first
[1418]	345	END DO
	346	END DO
	347	DO jj = 1, jpj
	348	DO ji = 1, jpi
[3211]	349	#endif
[1518]	350	IF( zsum(ji,jj) > 0.e0 ) zsum(ji,jj) = 1.e0 / zsum(ji,jj)
[1418]	351	END DO
	352	END DO
	353
[1495]	354	! ! first estimation bounded by 10 cm2/s (with n2 bounded by rn_n2min)
	355	zcoef = rn_tfe_itf / ( rn_tfe * rau0 )
[3211]	356	#if defined key_z_first
	357	DO jj = 1, jpj
	358	DO ji = 1, jpi
	359	DO jk = 1, jpk
	360	zavt_itf(ji,jj,jk) = MIN( 10.e-4, zcoef * en_tmx(ji,jj) * zsum(ji,jj) * zempba_3d(ji,jj,jk) &
	361	& / MAX( rn_n2min, rn2(ji,jj,jk) ) * tmask(ji,jj,jk) )
	362	END DO
	363	END DO
	364	END DO
	365	#else
[1518]	366	DO jk = 1, jpk
	367	zavt_itf(:,:,jk) = MIN( 10.e-4, zcoef * en_tmx(:,:) * zsum(:,:) * zempba_3d(:,:,jk) &
	368	& / MAX( rn_n2min, rn2(:,:,jk) ) * tmask(:,:,jk) )
[1495]	369	END DO
[3211]	370	#endif
[1418]	371
[3211]	372	#if defined key_z_first
	373	DO jj = 1, jpj
	374	DO ji = 1, jpi
	375	zkz(ji,jj) = 0.e0 ! Associated potential energy consummed over the whole water column
	376	DO jk = 2, jpkm1
	377	zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) &
	378	& * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk)
	379	END DO
	380	END DO
	381	END DO
	382	#else
[1418]	383	zkz(:,:) = 0.e0 ! Associated potential energy consummed over the whole water column
	384	DO jk = 2, jpkm1
[1495]	385	zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk)
[1418]	386	END DO
[3211]	387	#endif
[1418]	388
	389	DO jj = 1, jpj ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx
	390	DO ji = 1, jpi
	391	IF( zkz(ji,jj) /= 0.e0 ) zkz(ji,jj) = en_tmx(ji,jj) * rn_tfe_itf / rn_tfe / zkz(ji,jj)
	392	END DO
	393	END DO
	394
[3211]	395	#if defined key_z_first
	396	DO jj = 1, jpj
	397	DO ji = 1, jpi
	398	zcoef = MIN( zkz(:,:), 120./10. ) ! kz max = 120 cm2/s
	399	DO jk = 2, jpkm1 ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
	400	zavt_itf(ji,jj,jk) = zavt_itf(ji,jj,jk) * zcoef
	401	END DO
	402	END DO
	403	END DO
	404	#else
[1495]	405	DO jk = 2, jpkm1 ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s
	406	zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) ! kz max = 120 cm2/s
[1418]	407	END DO
[3211]	408	#endif
[1418]	409
[3211]	410	IF( kt == nit000 ) THEN ! diagnose the energy consumed by zavt_itf
[1418]	411	ztpc = 0.e0
[3211]	412	#if defined key_z_first
	413	DO jj = 1, jpj
	414	DO ji = 1, jpi
	415	DO jk = 1, jpk
	416	#else
	417	DO jk = 1, jpk
	418	DO jj = 1, jpj
	419	DO ji = 1, jpi
	420	#endif
[1495]	421	ztpc = ztpc + e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) &
	422	& * zavt_itf(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
[1418]	423	END DO
	424	END DO
	425	END DO
[1495]	426	ztpc= rau0 * ztpc / ( rn_me * rn_tfe_itf )
	427	IF(lwp) WRITE(numout,) ' N Total power consumption by zavt_itf: ztpc = ', ztpc 1.e-12 ,'TW'
[1418]	428	ENDIF
	429
[1546]	430	! ! Update pav with the ITF mixing coefficient
[3211]	431	#if defined key_z_first
	432	DO jj = 1, jpj
	433	DO ji = 1, jpi
	434	DO jk = 2, jpkm1
	435	pav(ji,jj,jk) = pav (ji,jj,jk) * ( 1.e0 - mask_itf(ji,jj) ) &
	436	& + zavt_itf(ji,jj,jk) * mask_itf(ji,jj)
	437	END DO
	438	END DO
	439	END DO
	440	#else
[1418]	441	DO jk = 2, jpkm1
[1546]	442	pav(:,:,jk) = pav (:,:,jk) * ( 1.e0 - mask_itf(:,:) ) &
	443	& + zavt_itf(:,:,jk) * mask_itf(:,:)
[1418]	444	END DO
[3211]	445	#endif
[1418]	446	!
[2715]	447	IF( wrk_not_released(2, 2,3,4,5) .OR. &
	448	wrk_not_released(3, 1,2,3,4,5) )THEN
	449	CALL ctl_stop('tmx_itf : failed to release workspace arrays.')
	450	END IF
	451	!
[1418]	452	END SUBROUTINE tmx_itf
	453
[3211]	454	!! * Reset control of array index permutation
	455	# include "oce_ftrans.h90"
	456	# include "dom_oce_ftrans.h90"
	457	# include "zdf_oce_ftrans.h90"
	458	!FTRANS az_tmx :I :I :z
[1418]	459
	460	SUBROUTINE zdf_tmx_init
	461	!!----------------------------------------------------------------------
	462	!! * ROUTINE zdf_tmx_init *
	463	!!
	464	!! ** Purpose : Initialization of the vertical tidal mixing, Reading
[1496]	465	!! of M2 and K1 tidal energy in nc files
[1418]	466	!!
[1496]	467	!! ** Method : - Read the namtmx namelist and check the parameters
	468	!!
	469	!! - Read the input data in NetCDF files :
	470	!! M2 and K1 tidal energy. The total tidal energy, en_tmx,
	471	!! is the sum of M2, K1 and S2 energy where S2 is assumed
	472	!! to be: S2=(1/2)^2 * M2
	473	!! mask_itf, a mask array that determine where substituing
	474	!! the standard Simmons et al. (2005) formulation with the
	475	!! one of Koch_Larrouy et al. (2007).
	476	!!
[1418]	477	!! - Compute az_tmx, a 3D coefficient that allows to compute
[1496]	478	!! the standard tidal-induced vertical mixing as follows:
	479	!! Kz_tides = az_tmx / max( rn_n2min, N^2 )
	480	!! with az_tmx a bottom intensified coefficient is given by:
	481	!! az_tmx(z) = en_tmx / ( rau0 * rn_htmx ) * EXP( -(H-z)/rn_htmx )
	482	!! / ( 1. - EXP( - H /rn_htmx ) )
	483	!! where rn_htmx the characteristic length scale of the bottom
	484	!! intensification, en_tmx the tidal energy, and H the ocean depth
[1418]	485	!!
	486	!! ** input : - Namlist namtmx
[1496]	487	!! - NetCDF file : M2_ORCA2.nc, K1_ORCA2.nc, and mask_itf.nc
[1418]	488	!!
	489	!! ** Action : - Increase by 1 the nstop flag is setting problem encounter
	490	!! - defined az_tmx used to compute tidal-induced mixing
[1496]	491	!!
	492	!! References : Simmons et al. 2004, Ocean Modelling, 6, 3-4, 245-263.
	493	!! Koch-Larrouy et al. 2007, GRL.
[1418]	494	!!----------------------------------------------------------------------
[2715]	495	USE oce , zav_tide => ua ! ua used as workspace
	496	USE wrk_nemo, ONLY: zem2 => wrk_2d_1 ! read M2 and
	497	USE wrk_nemo, ONLY: zek1 => wrk_2d_2 ! K1 tidal energy
	498	USE wrk_nemo, ONLY: zkz => wrk_2d_3 ! total M2, K1 and S2 tidal energy
	499	USE wrk_nemo, ONLY: zfact => wrk_2d_4 ! used for vertical structure function
	500	USE wrk_nemo, ONLY: zhdep => wrk_2d_5 ! Ocean depth
	501	USE wrk_nemo, ONLY: zpc => wrk_3d_1 ! power consumption
[3211]	502
	503	!! DCSE_NEMO: need additional directives for renamed module variables
	504	!FTRANS zpc :I :I :z
	505
[1546]	506	!!
[2715]	507	INTEGER :: ji, jj, jk ! dummy loop indices
	508	INTEGER :: inum ! local integer
	509	REAL(wp) :: ztpc, ze_z ! local scalars
[3211]	510	#if defined key_z_first
	511	REAL(wp) :: zcoef ! local scalar
	512	#endif
	513
[1496]	514	!!
[1601]	515	NAMELIST/namzdf_tmx/ rn_htmx, rn_n2min, rn_tfe, rn_me, ln_tmx_itf, rn_tfe_itf
[1418]	516	!!----------------------------------------------------------------------
	517
[2715]	518	IF( wrk_in_use(2, 1,2,3,4,5) .OR. wrk_in_use(3, 1) ) THEN
	519	CALL ctl_stop('zdf_tmx_init : requested workspace arrays unavailable.') ; RETURN
	520	END IF
	521
[1601]	522	REWIND( numnam ) ! Read Namelist namtmx : Tidal Mixing
	523	READ ( numnam, namzdf_tmx )
[1537]	524
	525	IF(lwp) THEN ! Control print
[1418]	526	WRITE(numout,*)
	527	WRITE(numout,*) 'zdf_tmx_init : tidal mixing'
	528	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	529	WRITE(numout,*) ' Namelist namzdf_tmx : set tidal mixing parameters'
[1537]	530	WRITE(numout,*) ' Vertical decay scale for turbulence = ', rn_htmx
	531	WRITE(numout,*) ' Brunt-Vaisala frequency threshold = ', rn_n2min
	532	WRITE(numout,*) ' Tidal dissipation efficiency = ', rn_tfe
	533	WRITE(numout,*) ' Mixing efficiency = ', rn_me
	534	WRITE(numout,*) ' ITF specific parameterisation = ', ln_tmx_itf
	535	WRITE(numout,*) ' ITF tidal dissipation efficiency = ', rn_tfe_itf
[1418]	536	ENDIF
	537
[2715]	538	! ! allocate tmx arrays
	539	IF( zdf_tmx_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
	540
[1537]	541	IF( ln_tmx_itf ) THEN ! read the Indonesian Through Flow mask
[1518]	542	CALL iom_open('mask_itf',inum)
	543	CALL iom_get (inum, jpdom_data, 'tmaskitf',mask_itf,1) !
	544	CALL iom_close(inum)
	545	ENDIF
[1418]	546
	547	! read M2 tidal energy flux : W/m2 ( zem2 < 0 )
	548	CALL iom_open('M2rowdrg',inum)
	549	CALL iom_get (inum, jpdom_data, 'field',zem2,1) !
	550	CALL iom_close(inum)
	551
	552	! read K1 tidal energy flux : W/m2 ( zek1 < 0 )
	553	CALL iom_open('K1rowdrg',inum)
	554	CALL iom_get (inum, jpdom_data, 'field',zek1,1) !
	555	CALL iom_close(inum)
	556
	557	! Total tidal energy ( M2, S2 and K1 with S2=(1/2)^2 * M2 )
	558	! only the energy available for mixing is taken into account,
	559	! (mixing efficiency tidal dissipation efficiency)
	560	en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1)
	561
[1496]	562	! Vertical structure (az_tmx)
	563	DO jj = 1, jpj ! part independent of the level
[1418]	564	DO ji = 1, jpi
[2528]	565	zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean
[1418]	566	zfact(ji,jj) = rau0 * rn_htmx * ( 1. - EXP( -zhdep(ji,jj) / rn_htmx ) )
	567	IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = en_tmx(ji,jj) / zfact(ji,jj)
	568	END DO
	569	END DO
[3211]	570	#if defined key_z_first
	571	DO jj = 1, jpj
	572	DO ji = 1, jpi
	573	DO jk= 1, jpk ! complete with the level-dependent part
	574	#else
[1418]	575	DO jk= 1, jpk ! complete with the level-dependent part
	576	DO jj = 1, jpj
	577	DO ji = 1, jpi
[3211]	578	#endif
[1418]	579	az_tmx(ji,jj,jk) = zfact(ji,jj) * EXP( -( zhdep(ji,jj)-fsdepw(ji,jj,jk) ) / rn_htmx ) * tmask(ji,jj,jk)
	580	END DO
	581	END DO
	582	END DO
	583
	584	IF( nprint == 1 .AND. lwp ) THEN
	585	! Control print
	586	! Total power consumption due to vertical mixing
[1546]	587	! zpc = rau0 * 1/rn_me * rn2 * zav_tide
[3211]	588	#if defined key_z_first
	589	DO jj = 1, jpj
	590	DO ji = 1, jpi
	591	zav_tide(ji,jj,1) = 0.e0
	592	DO jk = 2, jpkm1
	593	zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
	594	END DO
	595	zav_tide(ji,jj,jpk) = 0.e0
	596	END DO
	597	END DO
	598	#else
[1546]	599	zav_tide(:,:,:) = 0.e0
[1418]	600	DO jk = 2, jpkm1
[1546]	601	zav_tide(:,:,jk) = az_tmx(:,:,jk) / MAX( rn_n2min, rn2(:,:,jk) )
[1418]	602	END DO
[3211]	603	#endif
[1418]	604
	605	ztpc = 0.e0
[1546]	606	zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:)
[3211]	607	#if defined key_z_first
	608	DO jj = 1, jpj
	609	DO ji = 1, jpi
	610	DO jk= 2, jpkm1
	611	#else
[1418]	612	DO jk= 2, jpkm1
	613	DO jj = 1, jpj
	614	DO ji = 1, jpi
[3211]	615	#endif
[1418]	616	ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
	617	END DO
	618	END DO
	619	END DO
	620	ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
	621
	622	WRITE(numout,*)
	623	WRITE(numout,) ' Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc 1.e-12 ,'TW'
	624
	625
	626	! control print 2
[1546]	627	zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 )
[3211]	628	#if defined key_z_first
	629	DO jj = 1, jpj
	630	DO ji = 1, jpi
	631	zkz(ji,jj) = 0.e0
	632	DO jk = 2, jpkm1
	633	#else
[1418]	634	zkz(:,:) = 0.e0
	635	DO jk = 2, jpkm1
[3211]	636	DO jj = 1, jpj
	637	DO ji = 1, jpi
	638	#endif
	639	zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) &
	640	& * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk)
	641	END DO
[1418]	642	END DO
	643	END DO
	644	! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz
	645	DO jj = 1, jpj
	646	DO ji = 1, jpi
	647	IF( zkz(ji,jj) /= 0.e0 ) THEN
	648	zkz(ji,jj) = en_tmx(ji,jj) / zkz(ji,jj)
	649	ENDIF
	650	END DO
	651	END DO
	652	ztpc = 1.e50
	653	DO jj = 1, jpj
	654	DO ji = 1, jpi
	655	IF( zkz(ji,jj) /= 0.e0 ) THEN
[3211]	656	ztpc = MIN( zkz(ji,jj), ztpc)
[1418]	657	ENDIF
	658	END DO
	659	END DO
[3211]	660	WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', MAXVAL(zkz(:,:) )
[1418]	661
[3211]	662	#if defined key_z_first
	663	DO jj = 1, jpj
	664	DO ji = 1, jpi
	665	zcoef = MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s
	666	DO jk = 2, jpkm1
	667	zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * zcoef
	668	END DO
	669	END DO
	670	END DO
	671	#else
[1418]	672	DO jk = 2, jpkm1
[1546]	673	zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) !kz max = 300 cm2/s
[1418]	674	END DO
[3211]	675	#endif
[1418]	676	ztpc = 0.e0
[3211]	677	zpc(:,:,:) = MAX(0.e0,rn2(:,:,:)) * zav_tide(:,:,:)
	678	#if defined key_z_first
	679	DO jj = 1, jpj
	680	DO ji = 1, jpi
	681	DO jk= 1, jpk
	682	#else
[1418]	683	DO jk= 1, jpk
	684	DO jj = 1, jpj
	685	DO ji = 1, jpi
[3211]	686	#endif
[1418]	687	ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)
	688	END DO
	689	END DO
	690	END DO
	691	ztpc= rau0 * 1/(rn_tfe * rn_me) * ztpc
	692	WRITE(numout,) ' 2 Total power consumption of the tidally driven part of Kz : ztpc = ', ztpc 1.e-12 ,'TW'
	693
	694	DO jk = 1, jpk
[1546]	695	ze_z = SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) ) &
[1418]	696	& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
	697	ztpc = 1.E50
	698	DO jj = 1, jpj
	699	DO ji = 1, jpi
[1546]	700	IF( zav_tide(ji,jj,jk) /= 0.e0 ) ztpc =Min( ztpc, zav_tide(ji,jj,jk) )
[1418]	701	END DO
	702	END DO
	703	WRITE(numout,) ' N2 min - jk= ', jk,' ', ze_z 1.e4,' cm2/s min= ',ztpc*1.e4, &
[1546]	704	& 'max= ', MAXVAL(zav_tide(:,:,jk) )*1.e4, ' cm2/s'
[1418]	705	END DO
	706
	707	WRITE(numout,) ' e_tide : ', SUM( e1te2ten_tmx ) / ( rn_tfe rn_me ) * 1.e-12, 'TW'
	708	WRITE(numout,*)
	709	WRITE(numout,*) ' Initial profile of tidal vertical mixing'
	710	DO jk = 1, jpk
	711	DO jj = 1,jpj
	712	DO ji = 1,jpi
	713	zkz(ji,jj) = az_tmx(ji,jj,jk) /MAX( rn_n2min, rn2(ji,jj,jk) )
	714	END DO
	715	END DO
	716	ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) &
	717	& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
	718	WRITE(numout,) ' jk= ', jk,' ', ze_z 1.e4,' cm2/s'
	719	END DO
	720	DO jk = 1, jpk
	721	zkz(:,:) = az_tmx(:,:,jk) /rn_n2min
	722	ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) &
	723	& / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )
	724	WRITE(numout,*)
	725	WRITE(numout,) ' N2 min - jk= ', jk,' ', ze_z 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4, &
	726	& 'max= ', MAXVAL(zkz)*1.e4, ' cm2/s'
	727	END DO
	728	!
	729	ENDIF
[2528]	730	!
[2715]	731	IF(wrk_not_released(2, 1,2,3,4,5) .OR. &
	732	wrk_not_released(3, 1) ) CALL ctl_stop( 'zdf_tmx_init : failed to release workspace arrays' )
	733	!
[1418]	734	END SUBROUTINE zdf_tmx_init
	735
	736	#else
	737	!!----------------------------------------------------------------------
	738	!! Default option Dummy module NO Tidal MiXing
	739	!!----------------------------------------------------------------------
	740	LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .FALSE. !: tidal mixing flag
	741	CONTAINS
[2528]	742	SUBROUTINE zdf_tmx_init ! Dummy routine
	743	WRITE(,) 'zdf_tmx: You should not have seen this print! error?'
	744	END SUBROUTINE zdf_tmx_init
	745	SUBROUTINE zdf_tmx( kt ) ! Dummy routine
[1418]	746	WRITE(,) 'zdf_tmx: You should not have seen this print! error?', kt
	747	END SUBROUTINE zdf_tmx
	748	#endif
	749
	750	!!======================================================================
	751	END MODULE zdftmx

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

Download in other formats: