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.

zdftke.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 7729

Last change on this file since 7729 was 7698, checked in by mocavero, 7 years ago
update trunk with OpenMP parallelization
Property svn:keywords set to `Id`
File size: 52.3 KB

Rev	Line
[1531]	1	MODULE zdftke
[1239]	2	!!======================================================================
[1531]	3	!! * MODULE zdftke *
[1239]	4	!! Ocean physics: vertical mixing coefficient computed from the tke
	5	!! turbulent closure parameterization
	6	!!=====================================================================
[1492]	7	!! History : OPA ! 1991-03 (b. blanke) Original code
	8	!! 7.0 ! 1991-11 (G. Madec) bug fix
	9	!! 7.1 ! 1992-10 (G. Madec) new mixing length and eav
	10	!! 7.2 ! 1993-03 (M. Guyon) symetrical conditions
	11	!! 7.3 ! 1994-08 (G. Madec, M. Imbard) nn_pdl flag
	12	!! 7.5 ! 1996-01 (G. Madec) s-coordinates
	13	!! 8.0 ! 1997-07 (G. Madec) lbc
	14	!! 8.1 ! 1999-01 (E. Stretta) new option for the mixing length
	15	!! NEMO 1.0 ! 2002-06 (G. Madec) add tke_init routine
	16	!! - ! 2004-10 (C. Ethe ) 1D configuration
	17	!! 2.0 ! 2006-07 (S. Masson) distributed restart using iom
	18	!! 3.0 ! 2008-05 (C. Ethe, G.Madec) : update TKE physics:
	19	!! ! - tke penetration (wind steering)
	20	!! ! - suface condition for tke & mixing length
	21	!! ! - Langmuir cells
	22	!! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb
	23	!! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters
	24	!! - ! 2008-12 (G. Reffray) stable discretization of the production term
	25	!! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl
	26	!! ! + cleaning of the parameters + bugs correction
[2528]	27	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
[5120]	28	!! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability
[1239]	29	!!----------------------------------------------------------------------
[5836]	30	#if defined key_zdftke
[1239]	31	!!----------------------------------------------------------------------
[1531]	32	!! 'key_zdftke' TKE vertical physics
[1239]	33	!!----------------------------------------------------------------------
[3625]	34	!! zdf_tke : update momentum and tracer Kz from a tke scheme
	35	!! tke_tke : tke time stepping: update tke at now time step (en)
	36	!! tke_avn : compute mixing length scale and deduce avm and avt
	37	!! zdf_tke_init : initialization, namelist read, and parameters control
	38	!! tke_rst : read/write tke restart in ocean restart file
[1239]	39	!!----------------------------------------------------------------------
[2528]	40	USE oce ! ocean: dynamics and active tracers variables
	41	USE phycst ! physical constants
	42	USE dom_oce ! domain: ocean
	43	USE domvvl ! domain: variable volume layer
[1492]	44	USE sbc_oce ! surface boundary condition: ocean
[2528]	45	USE zdf_oce ! vertical physics: ocean variables
	46	USE zdfmxl ! vertical physics: mixed layer
[1492]	47	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	48	USE prtctl ! Print control
	49	USE in_out_manager ! I/O manager
	50	USE iom ! I/O manager library
[2715]	51	USE lib_mpp ! MPP library
[3294]	52	USE wrk_nemo ! work arrays
	53	USE timing ! Timing
[3625]	54	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[5656]	55	#if defined key_agrif
	56	USE agrif_opa_interp
	57	USE agrif_opa_update
	58	#endif
[1239]	59
	60	IMPLICIT NONE
	61	PRIVATE
	62
[2528]	63	PUBLIC zdf_tke ! routine called in step module
	64	PUBLIC zdf_tke_init ! routine called in opa module
	65	PUBLIC tke_rst ! routine called in step module
[1239]	66
[2715]	67	LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag
[1239]	68
[4147]	69	! !! Namelist namzdf_tke
	70	LOGICAL :: ln_mxl0 ! mixing length scale surface value as function of wind stress or not
	71	INTEGER :: nn_mxl ! type of mixing length (=0/1/2/3)
	72	REAL(wp) :: rn_mxl0 ! surface min value of mixing length (kappaz_o=0.40.1 m) [m]
	73	INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1)
	74	REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediffmxlsqrt(e)
	75	REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation
	76	REAL(wp) :: rn_ebb ! coefficient of the surface input of tke
	77	REAL(wp) :: rn_emin ! minimum value of tke [m2/s2]
	78	REAL(wp) :: rn_emin0 ! surface minimum value of tke [m2/s2]
	79	REAL(wp) :: rn_bshear ! background shear (>0) currently a numerical threshold (do not change it)
	80	INTEGER :: nn_etau ! type of depth penetration of surface tke (=0/1/2/3)
	81	INTEGER :: nn_htau ! type of tke profile of penetration (=0/1)
	82	REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean
	83	LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not
	84	REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells
[1239]	85
[4147]	86	REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
	87	REAL(wp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m]
[2528]	88	REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3)
	89	REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3)
[1239]	90
[2715]	91	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
	92	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation
[5656]	93	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apdlr ! now mixing lenght of dissipation
[2715]	94	#if defined key_c1d
	95	! !! 1D cfg only ('key_c1d')
	96	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales
	97	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers
	98	#endif
[1492]	99
[1239]	100	!! * Substitutions
	101	# include "vectopt_loop_substitute.h90"
	102	!!----------------------------------------------------------------------
[5836]	103	!! NEMO/OPA 3.7 , NEMO Consortium (2015)
[2528]	104	!! $Id$
	105	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1239]	106	!!----------------------------------------------------------------------
	107	CONTAINS
	108
[2715]	109	INTEGER FUNCTION zdf_tke_alloc()
	110	!!----------------------------------------------------------------------
	111	!! * FUNCTION zdf_tke_alloc *
	112	!!----------------------------------------------------------------------
	113	ALLOCATE( &
	114	#if defined key_c1d
	115	& e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , &
	116	& e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , &
	117	#endif
[5836]	118	& htau (jpi,jpj) , dissl(jpi,jpj,jpk) , &
	119	& apdlr(jpi,jpj,jpk) , STAT= zdf_tke_alloc )
[2715]	120	!
	121	IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc )
	122	IF( zdf_tke_alloc /= 0 ) CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
	123	!
	124	END FUNCTION zdf_tke_alloc
	125
	126
[1531]	127	SUBROUTINE zdf_tke( kt )
[1239]	128	!!----------------------------------------------------------------------
[1531]	129	!! * ROUTINE zdf_tke *
[1239]	130	!!
	131	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
[1492]	132	!! coefficients using a turbulent closure scheme (TKE).
[1239]	133	!!
[1492]	134	!! ** Method : The time evolution of the turbulent kinetic energy (tke)
	135	!! is computed from a prognostic equation :
	136	!! d(en)/dt = avm (d(u)/dz)**2 ! shear production
	137	!! + d( avm d(en)/dz )/dz ! diffusion of tke
	138	!! + avt N^2 ! stratif. destruc.
	139	!! - rn_ediss / emxl en**(2/3) ! Kolmogoroff dissipation
[1239]	140	!! with the boundary conditions:
[1695]	141	!! surface: en = max( rn_emin0, rn_ebb * taum )
[1239]	142	!! bottom : en = rn_emin
[1492]	143	!! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
	144	!!
	145	!! The now Turbulent kinetic energy is computed using the following
	146	!! time stepping: implicit for vertical diffusion term, linearized semi
	147	!! implicit for kolmogoroff dissipation term, and explicit forward for
	148	!! both buoyancy and shear production terms. Therefore a tridiagonal
	149	!! linear system is solved. Note that buoyancy and shear terms are
	150	!! discretized in a energy conserving form (Bruchard 2002).
	151	!!
	152	!! The dissipative and mixing length scale are computed from en and
	153	!! the stratification (see tke_avn)
	154	!!
	155	!! The now vertical eddy vicosity and diffusivity coefficients are
	156	!! given by:
	157	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
	158	!! avt = max( avmb, pdl * avm )
[1239]	159	!! eav = max( avmb, avm )
[1492]	160	!! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
	161	!! given by an empirical funtion of the localRichardson number if nn_pdl=1
[1239]	162	!!
	163	!! ** Action : compute en (now turbulent kinetic energy)
	164	!! update avt, avmu, avmv (before vertical eddy coef.)
	165	!!
	166	!! References : Gaspar et al., JGR, 1990,
	167	!! Blanke and Delecluse, JPO, 1991
	168	!! Mellor and Blumberg, JPO 2004
	169	!! Axell, JGR, 2002
[1492]	170	!! Bruchard OM 2002
[1239]	171	!!----------------------------------------------------------------------
[1492]	172	INTEGER, INTENT(in) :: kt ! ocean time step
[7698]	173	INTEGER :: jk, jj, ji
[1492]	174	!!----------------------------------------------------------------------
[1481]	175	!
[5656]	176	#if defined key_agrif
	177	! interpolation parent grid => child grid for avm_k ( ex : at west border: update column 1 and 2)
	178	IF( .NOT.Agrif_Root() ) CALL Agrif_Tke
	179	#endif
	180	!
[3632]	181	IF( kt /= nit000 ) THEN ! restore before value to compute tke
[7698]	182	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
	183	DO jk = 1, jpk
	184	DO jj = 1, jpj
	185	DO ji = 1, jpi
	186	avt (ji,jj,jk) = avt_k (ji,jj,jk)
	187	avm (ji,jj,jk) = avm_k (ji,jj,jk)
	188	avmu(ji,jj,jk) = avmu_k(ji,jj,jk)
	189	avmv(ji,jj,jk) = avmv_k(ji,jj,jk)
	190	END DO
	191	END DO
	192	END DO
[3632]	193	ENDIF
	194	!
[2528]	195	CALL tke_tke ! now tke (en)
[1492]	196	!
[2528]	197	CALL tke_avn ! now avt, avm, avmu, avmv
	198	!
[7698]	199	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
	200	DO jk = 1, jpk
	201	DO jj = 1, jpj
	202	DO ji = 1, jpi
	203	avt_k (ji,jj,jk) = avt (ji,jj,jk)
	204	avm_k (ji,jj,jk) = avm (ji,jj,jk)
	205	avmu_k(ji,jj,jk) = avmu(ji,jj,jk)
	206	avmv_k(ji,jj,jk) = avmv(ji,jj,jk)
	207	END DO
	208	END DO
	209	END DO
[3632]	210	!
[5656]	211	#if defined key_agrif
	212	! Update child grid f => parent grid
	213	IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only
	214	#endif
	215	!
	216	END SUBROUTINE zdf_tke
[1239]	217
[1492]	218
[1481]	219	SUBROUTINE tke_tke
[1239]	220	!!----------------------------------------------------------------------
[1492]	221	!! * ROUTINE tke_tke *
	222	!!
	223	!! ** Purpose : Compute the now Turbulente Kinetic Energy (TKE)
	224	!!
	225	!! ** Method : - TKE surface boundary condition
[2528]	226	!! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
[1492]	227	!! - source term due to shear (saved in avmu, avmv arrays)
	228	!! - Now TKE : resolution of the TKE equation by inverting
	229	!! a tridiagonal linear system by a "methode de chasse"
	230	!! - increase TKE due to surface and internal wave breaking
	231	!!
	232	!! ** Action : - en : now turbulent kinetic energy)
	233	!! - avmu, avmv : production of TKE by shear at u and v-points
	234	!! (= Kz dz[Ub] * dz[Un] )
[1239]	235	!! ---------------------------------------------------------------------
[1705]	236	INTEGER :: ji, jj, jk ! dummy loop arguments
[2528]	237	!!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar
	238	!!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar
[1705]	239	REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
	240	REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
	241	REAL(wp) :: zbbrau, zesh2 ! temporary scalars
	242	REAL(wp) :: zfact1, zfact2, zfact3 ! - -
	243	REAL(wp) :: ztx2 , zty2 , zcof ! - -
	244	REAL(wp) :: ztau , zdif ! - -
	245	REAL(wp) :: zus , zwlc , zind ! - -
	246	REAL(wp) :: zzd_up, zzd_lw ! - -
[2528]	247	!!bfr REAL(wp) :: zebot ! - -
[5836]	248	INTEGER , POINTER, DIMENSION(:,: ) :: imlc
	249	REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc
	250	REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw, z3du, z3dv
[5656]	251	REAL(wp) :: zri ! local Richardson number
[1239]	252	!!--------------------------------------------------------------------
[1492]	253	!
[3294]	254	IF( nn_timing == 1 ) CALL timing_start('tke_tke')
	255	!
[5836]	256	CALL wrk_alloc( jpi,jpj, imlc ) ! integer
	257	CALL wrk_alloc( jpi,jpj, zhlc )
	258	CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )
[3294]	259	!
[1695]	260	zbbrau = rn_ebb / rau0 ! Local constant initialisation
[2528]	261	zfact1 = -.5_wp * rdt
	262	zfact2 = 1.5_wp * rdt * rn_ediss
	263	zfact3 = 0.5_wp * rn_ediss
[1492]	264	!
[5120]	265	!
[1492]	266	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	267	! ! Surface boundary condition on tke
	268	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[5120]	269	IF ( ln_isfcav ) THEN
[7698]	270	!$OMP PARALLEL DO schedule(static) private(jj, ji)
[5120]	271	DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin
	272	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	273	en(ji,jj,mikt(ji,jj)) = rn_emin * tmask(ji,jj,1)
[5120]	274	END DO
	275	END DO
	276	END IF
[7698]	277	!$OMP PARALLEL DO schedule(static) private(jj, ji)
[1695]	278	DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0)
[1481]	279	DO ji = fs_2, fs_jpim1 ! vector opt.
[5120]	280	en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
[1481]	281	END DO
	282	END DO
[2528]	283
	284	!!bfr - start commented area
[1492]	285	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	286	! ! Bottom boundary condition on tke
	287	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1719]	288	!
	289	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	290	! Tests to date have found the bottom boundary condition on tke to have very little effect.
	291	! The condition is coded here for completion but commented out until there is proof that the
	292	! computational cost is justified
	293	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	294	! en(bot) = (rn_ebb0/rau0)0.5sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
	295	!! DO jj = 2, jpjm1
	296	!! DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	297	!! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
	298	!! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) )
	299	!! zty2 = bfrva(ji,jj ) * vb(ji,jj ,mbkv(ji,jj )) + &
	300	!! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
[1719]	301	!! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
[2528]	302	!! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
[1719]	303	!! END DO
	304	!! END DO
[2528]	305	!!bfr - end commented area
[1492]	306	!
	307	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[2528]	308	IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002)
[1492]	309	! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[1239]	310	!
[1492]	311	! !* total energy produce by LC : cumulative sum over jk
[7698]	312	!$OMP PARALLEL
	313	!$OMP DO schedule(static) private(jj, ji)
	314	DO jj =1, jpj
	315	DO ji=1, jpi
	316	zpelc(ji,jj,1) = MAX( rn2b(ji,jj,1), 0._wp ) * gdepw_n(ji,jj,1) * e3w_n(ji,jj,1)
	317	END DO
	318	END DO
[1239]	319	DO jk = 2, jpk
[7698]	320	!$OMP DO schedule(static) private(jj, ji)
	321	DO jj =1, jpj
	322	DO ji=1, jpi
	323	zpelc(ji,jj,jk) = zpelc(ji,jj,jk-1) + MAX( rn2b(ji,jj,jk), 0._wp ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
	324	END DO
	325	END DO
[1239]	326	END DO
[1492]	327	! !* finite Langmuir Circulation depth
[1705]	328	zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
[7698]	329	!$OMP DO schedule(static) private(jj,ji)
	330	DO jj = 1, jpj
	331	DO ji = 1, jpi
	332	imlc(ji,jj) = mbkt(ji,jj) + 1 ! Initialization to the number of w ocean point (=2 over land)
	333	END DO
	334	END DO
[1239]	335	DO jk = jpkm1, 2, -1
[7698]	336	!$OMP DO schedule(static) private(jj, ji, zus)
[1492]	337	DO jj = 1, jpj ! Last w-level at which zpelc>=0.5usus
	338	DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1)
[1705]	339	zus = zcof * taum(ji,jj)
[1239]	340	IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk
	341	END DO
	342	END DO
	343	END DO
[1492]	344	! ! finite LC depth
[7698]	345	!$OMP DO schedule(static) private(jj, ji)
[1492]	346	DO jj = 1, jpj
[1239]	347	DO ji = 1, jpi
[6140]	348	zhlc(ji,jj) = gdepw_n(ji,jj,imlc(ji,jj))
[1239]	349	END DO
	350	END DO
[1705]	351	zcof = 0.016 / SQRT( zrhoa * zcdrag )
[7698]	352	!$OMP DO schedule(static) private(jk, jj, ji, zus, zind, zwlc)
[1492]	353	DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en
[1239]	354	DO jj = 2, jpjm1
	355	DO ji = fs_2, fs_jpim1 ! vector opt.
[1705]	356	zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift
[1492]	357	! ! vertical velocity due to LC
[6140]	358	zind = 0.5 - SIGN( 0.5, gdepw_n(ji,jj,jk) - zhlc(ji,jj) )
	359	zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )
[1492]	360	! ! TKE Langmuir circulation source term
[6497]	361	en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) &
	362	& / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1239]	363	END DO
	364	END DO
	365	END DO
[7698]	366	!$OMP END PARALLEL
[1239]	367	!
	368	ENDIF
[1492]	369	!
	370	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	371	! ! Now Turbulent kinetic energy (output in en)
	372	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	373	! ! Resolution of a tridiagonal linear system by a "methode de chasse"
	374	! ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
	375	! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
	376	!
[7698]	377	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
[1492]	378	DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form)
[5656]	379	DO jj = 1, jpjm1
	380	DO ji = 1, fs_jpim1 ! vector opt.
	381	z3du(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji+1,jj,jk) ) &
	382	& * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) &
[5803]	383	& * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) &
[6140]	384	& / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
[5656]	385	z3dv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk ) + avm(ji,jj+1,jk) ) &
	386	& * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) &
[5803]	387	& * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) &
[6140]	388	& / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
[1492]	389	END DO
	390	END DO
	391	END DO
	392	!
[5656]	393	IF( nn_pdl == 1 ) THEN !* Prandtl number case: compute apdlr
	394	! Note that zesh2 is also computed in the next loop.
	395	! We decided to compute it twice to keep code readability and avoid an IF case in the DO loops
[7698]	396	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zesh2, zri)
[5656]	397	DO jk = 2, jpkm1
	398	DO jj = 2, jpjm1
	399	DO ji = fs_2, fs_jpim1 ! vector opt.
	400	! ! shear prod. at w-point weightened by mask
	401	zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
	402	& + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
	403	! ! local Richardson number
	404	zri = MAX( rn2b(ji,jj,jk), 0._wp ) * avm(ji,jj,jk) / ( zesh2 + rn_bshear )
	405	apdlr(ji,jj,jk) = MAX( 0.1_wp, ri_cri / MAX( ri_cri , zri ) )
	406
	407	END DO
	408	END DO
	409	END DO
	410	!
	411	ENDIF
[5836]	412	!
[7698]	413	!$OMP PARALLEL
	414	!$OMP DO schedule(static) private(jk, jj, ji, zcof, zzd_up, zzd_lw, zesh2)
[5120]	415	DO jk = 2, jpkm1 !* Matrix and right hand side in en
	416	DO jj = 2, jpjm1
	417	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	418	zcof = zfact1 * tmask(ji,jj,jk)
[6497]	419	# if defined key_zdftmx_new
	420	! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability)
	421	zzd_up = zcof * MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) & ! upper diagonal
	422	& / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk ) )
	423	zzd_lw = zcof * MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) & ! lower diagonal
	424	& / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) )
	425	# else
[1492]	426	zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal
[6140]	427	& / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk ) )
[1492]	428	zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal
[6140]	429	& / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) )
[6497]	430	# endif
[5656]	431	! ! shear prod. at w-point weightened by mask
	432	zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
	433	& + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
	434	!
[1492]	435	zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw)
	436	zd_lw(ji,jj,jk) = zzd_lw
[2528]	437	zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
[1239]	438	!
[1492]	439	! ! right hand side in en
[1481]	440	en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) &
[4990]	441	& + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) &
[5120]	442	& * wmask(ji,jj,jk)
[1239]	443	END DO
[5120]	444	END DO
	445	END DO
	446	! !* Matrix inversion from level 2 (tke prescribed at level 1)
	447	DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
[7698]	448	!$OMP DO schedule(static) private(jj, ji)
[5120]	449	DO jj = 2, jpjm1
	450	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	451	zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1)
[1239]	452	END DO
[5120]	453	END DO
	454	END DO
[7698]	455	!$OMP DO schedule(static) private(jj, ji)
[5836]	456	DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
[5120]	457	DO ji = fs_2, fs_jpim1 ! vector opt.
	458	zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke
	459	END DO
	460	END DO
	461	DO jk = 3, jpkm1
[7698]	462	!$OMP DO schedule(static) private(jj, ji)
[5120]	463	DO jj = 2, jpjm1
	464	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	465	zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1)
[1239]	466	END DO
[5120]	467	END DO
	468	END DO
[7698]	469	!$OMP DO schedule(static) private(jj, ji)
[5836]	470	DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
[5120]	471	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	472	en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
[5120]	473	END DO
	474	END DO
	475	DO jk = jpk-2, 2, -1
[7698]	476	!$OMP DO schedule(static) private(jj, ji)
[5120]	477	DO jj = 2, jpjm1
	478	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	479	en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
[1239]	480	END DO
[5120]	481	END DO
	482	END DO
[7698]	483	!$OMP DO schedule(static) private(jk,jj, ji)
[5120]	484	DO jk = 2, jpkm1 ! set the minimum value of tke
	485	DO jj = 2, jpjm1
	486	DO ji = fs_2, fs_jpim1 ! vector opt.
	487	en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk)
[1239]	488	END DO
	489	END DO
	490	END DO
[7698]	491	!$OMP END PARALLEL
[1239]	492
[1492]	493	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	494	! ! TKE due to surface and internal wave breaking
	495	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[6140]	496	!!gm BUG : in the exp remove the depth of ssh !!!
	497
	498
[2528]	499	IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction)
[7698]	500	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
[1492]	501	DO jk = 2, jpkm1
[1239]	502	DO jj = 2, jpjm1
	503	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	504	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) &
[5120]	505	& * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1239]	506	END DO
	507	END DO
[1492]	508	END DO
[2528]	509	ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction)
[7698]	510	!$OMP PARALLEL DO schedule(static) private(jj, ji, jk)
[1492]	511	DO jj = 2, jpjm1
	512	DO ji = fs_2, fs_jpim1 ! vector opt.
	513	jk = nmln(ji,jj)
[6140]	514	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) &
[5120]	515	& * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1239]	516	END DO
[1492]	517	END DO
[2528]	518	ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability)
[7698]	519	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, ztx2, zty2, ztau, zdif)
[1705]	520	DO jk = 2, jpkm1
	521	DO jj = 2, jpjm1
	522	DO ji = fs_2, fs_jpim1 ! vector opt.
	523	ztx2 = utau(ji-1,jj ) + utau(ji,jj)
	524	zty2 = vtau(ji ,jj-1) + vtau(ji,jj)
[4990]	525	ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress
[2528]	526	zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean
	527	zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications...
[6140]	528	en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw_n(ji,jj,jk) / htau(ji,jj) ) &
[5120]	529	& * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1)
[1705]	530	END DO
	531	END DO
	532	END DO
[1239]	533	ENDIF
[1492]	534	CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)
	535	!
[5836]	536	CALL wrk_dealloc( jpi,jpj, imlc ) ! integer
	537	CALL wrk_dealloc( jpi,jpj, zhlc )
	538	CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw, z3du, z3dv )
[2715]	539	!
[3294]	540	IF( nn_timing == 1 ) CALL timing_stop('tke_tke')
	541	!
[1239]	542	END SUBROUTINE tke_tke
	543
[1492]	544
	545	SUBROUTINE tke_avn
[1239]	546	!!----------------------------------------------------------------------
[1492]	547	!! * ROUTINE tke_avn *
[1239]	548	!!
[1492]	549	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
	550	!!
	551	!! ** Method : At this stage, en, the now TKE, is known (computed in
	552	!! the tke_tke routine). First, the now mixing lenth is
	553	!! computed from en and the strafification (N^2), then the mixings
	554	!! coefficients are computed.
	555	!! - Mixing length : a first evaluation of the mixing lengh
	556	!! scales is:
	557	!! mxl = sqrt(2*en) / N
	558	!! where N is the brunt-vaisala frequency, with a minimum value set
[2528]	559	!! to rmxl_min (rn_mxl0) in the interior (surface) ocean.
[1492]	560	!! The mixing and dissipative length scale are bound as follow :
	561	!! nn_mxl=0 : mxl bounded by the distance to surface and bottom.
	562	!! zmxld = zmxlm = mxl
	563	!! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
	564	!! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
	565	!! less than 1 (\|d/dz(mxl)\|<1) and zmxld = zmxlm = mxl
	566	!! nn_mxl=3 : mxl is bounded from the surface to the bottom usings
	567	!! \|d/dz(xml)\|<1 to obtain lup, and from the bottom to
	568	!! the surface to obtain ldown. the resulting length
	569	!! scales are:
	570	!! zmxld = sqrt( lup * ldown )
	571	!! zmxlm = min ( lup , ldown )
	572	!! - Vertical eddy viscosity and diffusivity:
	573	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
	574	!! avt = max( avmb, pdlr * avm )
	575	!! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
	576	!!
	577	!! ** Action : - avt : now vertical eddy diffusivity (w-point)
	578	!! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
[1239]	579	!!----------------------------------------------------------------------
[2715]	580	INTEGER :: ji, jj, jk ! dummy loop indices
	581	REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars
[5836]	582	REAL(wp) :: zdku, zri, zsqen ! - -
[2715]	583	REAL(wp) :: zdkv, zemxl, zemlm, zemlp ! - -
[3294]	584	REAL(wp), POINTER, DIMENSION(:,:,:) :: zmpdl, zmxlm, zmxld
[1239]	585	!!--------------------------------------------------------------------
[3294]	586	!
	587	IF( nn_timing == 1 ) CALL timing_start('tke_avn')
[1239]	588
[3294]	589	CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
	590
[1492]	591	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	592	! ! Mixing length
	593	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	594	!
	595	! !* Buoyancy length scale: l=sqrt(2e/n*2)
	596	!
[5120]	597	! initialisation of interior minimum value (avoid a 2d loop with mikt)
[7698]	598	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
	599	DO jk = 1, jpk
	600	DO jj = 1, jpj
	601	DO ji = 1, jpi
	602	zmxlm(ji,jj,jk) = rmxl_min
	603	zmxld(ji,jj,jk) = rmxl_min
	604	END DO
	605	END DO
	606	END DO
[5120]	607	!
[2528]	608	IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn2.e5taum/(rau0*g)
[7698]	609	!$OMP PARALLEL DO schedule(static) private(jj, ji, zraug)
[4990]	610	DO jj = 2, jpjm1
	611	DO ji = fs_2, fs_jpim1
[5120]	612	zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
	613	zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) )
[4990]	614	END DO
	615	END DO
	616	ELSE
[7698]	617	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	618	DO jj = 1, jpj
	619	DO ji = 1, jpi
	620	zmxlm(ji,jj,1) = rn_mxl0
	621	END DO
	622	END DO
[1239]	623	ENDIF
	624	!
[7698]	625	!$OMP PARALLEL
	626	!$OMP DO schedule(static) private(jk, jj, ji, zrn2)
[5120]	627	DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2)
	628	DO jj = 2, jpjm1
	629	DO ji = fs_2, fs_jpim1 ! vector opt.
[1239]	630	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
[5836]	631	zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
[1239]	632	END DO
	633	END DO
	634	END DO
[1492]	635	!
	636	! !* Physical limits for the mixing length
	637	!
[7698]	638	!$OMP DO schedule(static) private(jj,ji)
	639	DO jj = 1, jpj
	640	DO ji = 1, jpi
	641	zmxld(ji,jj, 1 ) = zmxlm(ji,jj,1) ! surface set to the minimum value
	642	zmxld(ji,jj,jpk) = rmxl_min ! last level set to the minimum value
	643	END DO
	644	END DO
	645	!$OMP END PARALLEL
[1492]	646	!
[1239]	647	SELECT CASE ( nn_mxl )
	648	!
[5836]	649	!!gm Not sure of that coding for ISF....
[6140]	650	! where wmask = 0 set zmxlm == e3w_n
[1239]	651	CASE ( 0 ) ! bounded by the distance to surface and bottom
[7698]	652	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
[5120]	653	DO jk = 2, jpkm1
	654	DO jj = 2, jpjm1
	655	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	656	zemxl = MIN( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), &
	657	& gdepw_n(ji,jj,mbkt(ji,jj)+1) - gdepw_n(ji,jj,jk) )
[5120]	658	! wmask prevent zmxlm = 0 if jk = mikt(ji,jj)
[6140]	659	zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
	660	zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),e3w_n(ji,jj,jk)) * (1 - wmask(ji,jj,jk))
[1239]	661	END DO
	662	END DO
	663	END DO
	664	!
	665	CASE ( 1 ) ! bounded by the vertical scale factor
[7698]	666	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zemxl)
[5120]	667	DO jk = 2, jpkm1
	668	DO jj = 2, jpjm1
	669	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	670	zemxl = MIN( e3w_n(ji,jj,jk), zmxlm(ji,jj,jk) )
[1239]	671	zmxlm(ji,jj,jk) = zemxl
	672	zmxld(ji,jj,jk) = zemxl
	673	END DO
	674	END DO
	675	END DO
	676	!
	677	CASE ( 2 ) ! \|dk[xml]\| bounded by e3t :
[7698]	678	!$OMP PARALLEL
[5120]	679	DO jk = 2, jpkm1 ! from the surface to the bottom :
[7698]	680	!$OMP DO schedule(static) private(jj, ji)
[5120]	681	DO jj = 2, jpjm1
	682	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	683	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
[1239]	684	END DO
[5120]	685	END DO
	686	END DO
	687	DO jk = jpkm1, 2, -1 ! from the bottom to the surface :
[7698]	688	!$OMP DO schedule(static) private(jj, ji, zemxl)
[5120]	689	DO jj = 2, jpjm1
	690	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	691	zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
[1239]	692	zmxlm(ji,jj,jk) = zemxl
	693	zmxld(ji,jj,jk) = zemxl
	694	END DO
	695	END DO
	696	END DO
[7698]	697	!$OMP END PARALLEL
[1239]	698	!
	699	CASE ( 3 ) ! lup and ldown, \|dk[xml]\| bounded by e3t :
[7698]	700	!$OMP PARALLEL
[5120]	701	DO jk = 2, jpkm1 ! from the surface to the bottom : lup
[7698]	702	!$OMP DO schedule(static) private(jj, ji)
[5120]	703	DO jj = 2, jpjm1
	704	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	705	zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t_n(ji,jj,jk-1), zmxlm(ji,jj,jk) )
[1239]	706	END DO
[5120]	707	END DO
	708	END DO
	709	DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown
[7698]	710	!$OMP DO schedule(static) private(jj, ji)
[5120]	711	DO jj = 2, jpjm1
	712	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	713	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t_n(ji,jj,jk+1), zmxlm(ji,jj,jk) )
[1239]	714	END DO
	715	END DO
	716	END DO
[7698]	717	!$OMP DO schedule(static) private(jk, jj, ji, zemlm, zemlp)
[1239]	718	DO jk = 2, jpkm1
	719	DO jj = 2, jpjm1
	720	DO ji = fs_2, fs_jpim1 ! vector opt.
	721	zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) )
	722	zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
	723	zmxlm(ji,jj,jk) = zemlm
	724	zmxld(ji,jj,jk) = zemlp
	725	END DO
	726	END DO
	727	END DO
[7698]	728	!$OMP END PARALLEL
[1239]	729	!
	730	END SELECT
[1492]	731	!
[1239]	732	# if defined key_c1d
[7698]	733	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
	734	DO jk = 1, jpk
	735	DO jj = 1, jpj
	736	DO ji = 1, jpi
	737	e_dis(ji,jj,jk) = zmxld(ji,jj,jk) ! c1d configuration : save mixing and dissipation turbulent length scales
	738	e_mix(ji,jj,jk) = zmxlm(ji,jj,jk)
	739	END DO
	740	END DO
	741	END DO
[1239]	742	# endif
	743
[1492]	744	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	745	! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt)
	746	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[7698]	747	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji, zsqen, zav)
[1492]	748	DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points
[1239]	749	DO jj = 2, jpjm1
	750	DO ji = fs_2, fs_jpim1 ! vector opt.
	751	zsqen = SQRT( en(ji,jj,jk) )
	752	zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen
[5120]	753	avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk)
	754	avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk)
[1239]	755	dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
	756	END DO
	757	END DO
	758	END DO
[1492]	759	CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)
	760	!
[7698]	761	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
[5120]	762	DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points
	763	DO jj = 2, jpjm1
	764	DO ji = fs_2, fs_jpim1 ! vector opt.
	765	avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk)
	766	avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk)
[4990]	767	END DO
[1239]	768	END DO
	769	END DO
	770	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions
[1492]	771	!
	772	IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt
[7698]	773	!$OMP PARALLEL DO schedule(static) private(jk, jj, ji)
[5120]	774	DO jk = 2, jpkm1
	775	DO jj = 2, jpjm1
	776	DO ji = fs_2, fs_jpim1 ! vector opt.
[5656]	777	avt(ji,jj,jk) = MAX( apdlr(ji,jj,jk) * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
[1492]	778	# if defined key_c1d
[5656]	779	e_pdl(ji,jj,jk) = apdlr(ji,jj,jk) * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number
[5836]	780	!!gm bug NO zri here....
	781	!!gm remove the specific diag for c1d !
[5656]	782	e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri
[1239]	783	# endif
	784	END DO
	785	END DO
	786	END DO
	787	ENDIF
	788	CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged)
	789
	790	IF(ln_ctl) THEN
	791	CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
	792	CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, &
	793	& tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
	794	ENDIF
	795	!
[3294]	796	CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld )
	797	!
	798	IF( nn_timing == 1 ) CALL timing_stop('tke_avn')
	799	!
[1492]	800	END SUBROUTINE tke_avn
[1239]	801
[1492]	802
[2528]	803	SUBROUTINE zdf_tke_init
[1239]	804	!!----------------------------------------------------------------------
[2528]	805	!! * ROUTINE zdf_tke_init *
[1239]	806	!!
	807	!! ** Purpose : Initialization of the vertical eddy diffivity and
[1492]	808	!! viscosity when using a tke turbulent closure scheme
[1239]	809	!!
[1601]	810	!! ** Method : Read the namzdf_tke namelist and check the parameters
[1492]	811	!! called at the first timestep (nit000)
[1239]	812	!!
[1601]	813	!! ** input : Namlist namzdf_tke
[1239]	814	!!
	815	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
	816	!!----------------------------------------------------------------------
	817	INTEGER :: ji, jj, jk ! dummy loop indices
[4147]	818	INTEGER :: ios
[1239]	819	!!
[2528]	820	NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , &
	821	& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
	822	& rn_mxl0 , nn_pdl , ln_lc , rn_lc , &
	823	& nn_etau , nn_htau , rn_efr
[1239]	824	!!----------------------------------------------------------------------
[2715]	825	!
[4147]	826	REWIND( numnam_ref ) ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy
	827	READ ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901)
	828	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in reference namelist', lwp )
	829
	830	REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
	831	READ ( numnam_cfg, namzdf_tke, IOSTAT = ios, ERR = 902 )
	832	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tke in configuration namelist', lwp )
[4624]	833	IF(lwm) WRITE ( numond, namzdf_tke )
[2715]	834	!
[2528]	835	ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number
[6497]	836	# if defined key_zdftmx_new
	837	! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
	838	rn_emin = 1.e-10_wp
	839	rmxl_min = 1.e-03_wp
	840	IF(lwp) THEN ! Control print
	841	WRITE(numout,*)
	842	WRITE(numout,*) 'zdf_tke_init : New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
	843	WRITE(numout,*) '~~~~~~~~~~~~'
	844	ENDIF
	845	# else
[2528]	846	rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity
[6497]	847	# endif
[2715]	848	!
[1492]	849	IF(lwp) THEN !* Control print
[1239]	850	WRITE(numout,*)
[2528]	851	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
	852	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	853	WRITE(numout,*) ' Namelist namzdf_tke : set tke mixing parameters'
[1705]	854	WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff
	855	WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss
	856	WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb
	857	WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin
	858	WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0
	859	WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear
	860	WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl
	861	WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl
	862	WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0
[2528]	863	WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0
	864	WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc
	865	WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc
[1705]	866	WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
	867	WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau
	868	WRITE(numout,*) ' fraction of en which pene. the thermocline rn_efr = ', rn_efr
[1239]	869	WRITE(numout,*)
[1601]	870	WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri
[1239]	871	ENDIF
[2715]	872	!
	873	! ! allocate tke arrays
	874	IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
	875	!
[1492]	876	! !* Check of some namelist values
[4990]	877	IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' )
	878	IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' )
	879	IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
[5407]	880	IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
[1239]	881
[2528]	882	IF( ln_mxl0 ) THEN
	883	IF(lwp) WRITE(numout,*) ' use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min'
	884	rn_mxl0 = rmxl_min
	885	ENDIF
	886
[1492]	887	IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln
[1239]	888
[1492]	889	! !* depth of penetration of surface tke
	890	IF( nn_etau /= 0 ) THEN
[1601]	891	SELECT CASE( nn_htau ) ! Choice of the depth of penetration
[2528]	892	CASE( 0 ) ! constant depth penetration (here 10 meters)
[7698]	893	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	894	DO jj = 1, jpj
	895	DO ji = 1, jpi
	896	htau(ji,jj) = 10._wp
	897	END DO
	898	END DO
[2528]	899	CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
[7698]	900	!$OMP PARALLEL DO schedule(static) private(jj,ji)
	901	DO jj = 1, jpj
	902	DO ji = 1, jpi
	903	htau(ji,jj) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) )
	904	END DO
	905	END DO
[1492]	906	END SELECT
	907	ENDIF
	908	! !* set vertical eddy coef. to the background value
[7698]	909	!$OMP PARALLEL
	910	!$OMP DO schedule(static) private(jk,jj,ji)
[1239]	911	DO jk = 1, jpk
[7698]	912	DO jj = 1, jpj
	913	DO ji = 1, jpi
	914	avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk)
	915	avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk)
	916	avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk)
	917	avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk)
	918	END DO
	919	END DO
[1239]	920	END DO
[7698]	921	!$OMP END DO NOWAIT
	922	!$OMP DO schedule(static) private(jk,jj,ji)
	923	DO jk = 1, jpk
	924	DO jj = 1, jpj
	925	DO ji = 1, jpi
	926	dissl(ji,jj,jk) = 1.e-12_wp
	927	END DO
	928	END DO
	929	END DO
	930	!$OMP END PARALLEL
[2715]	931	!
	932	CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files
[1239]	933	!
[2528]	934	END SUBROUTINE zdf_tke_init
[1239]	935
	936
[1531]	937	SUBROUTINE tke_rst( kt, cdrw )
[1239]	938	!!---------------------------------------------------------------------
[1531]	939	!! * ROUTINE tke_rst *
[1239]	940	!!
	941	!! ** Purpose : Read or write TKE file (en) in restart file
	942	!!
	943	!! ** Method : use of IOM library
	944	!! if the restart does not contain TKE, en is either
[1537]	945	!! set to rn_emin or recomputed
[1239]	946	!!----------------------------------------------------------------------
[2715]	947	INTEGER , INTENT(in) :: kt ! ocean time-step
	948	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[1239]	949	!
[7698]	950	INTEGER :: jit, jk, jj, ji ! dummy loop indices
[2715]	951	INTEGER :: id1, id2, id3, id4, id5, id6 ! local integers
[1239]	952	!!----------------------------------------------------------------------
	953	!
[1481]	954	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	955	! ! ---------------
	956	IF( ln_rstart ) THEN !* Read the restart file
	957	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
	958	id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. )
	959	id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. )
	960	id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
	961	id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
	962	id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
	963	!
	964	IF( id1 > 0 ) THEN ! 'en' exists
[1239]	965	CALL iom_get( numror, jpdom_autoglo, 'en', en )
[1481]	966	IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist
	967	CALL iom_get( numror, jpdom_autoglo, 'avt' , avt )
	968	CALL iom_get( numror, jpdom_autoglo, 'avm' , avm )
	969	CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu )
	970	CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv )
	971	CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
[1492]	972	ELSE ! one at least array is missing
	973	CALL tke_avn ! compute avt, avm, avmu, avmv and dissl (approximation)
[1481]	974	ENDIF
	975	ELSE ! No TKE array found: initialisation
	976	IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
[7698]	977	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
	978	DO jk = 1, jpk
	979	DO jj = 1, jpj
	980	DO ji = 1, jpi
	981	en (ji,jj,jk) = rn_emin * tmask(ji,jj,jk)
	982	END DO
	983	END DO
	984	END DO
[1492]	985	CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation)
[5112]	986	!
[7698]	987	!$OMP PARALLEL DO schedule(static) private(jk,jj,ji)
	988	DO jk = 1, jpk
	989	DO jj = 1, jpj
	990	DO ji = 1, jpi
	991	avt_k (ji,jj,jk) = avt (ji,jj,jk)
	992	avm_k (ji,jj,jk) = avm (ji,jj,jk)
	993	avmu_k(ji,jj,jk) = avmu(ji,jj,jk)
	994	avmv_k(ji,jj,jk) = avmv(ji,jj,jk)
	995	END DO
	996	END DO
	997	END DO
[5112]	998	!
[1531]	999	DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO
[1239]	1000	ENDIF
[1481]	1001	ELSE !* Start from rest
[7698]	1002	!$OMP PARALLEL
	1003	!$OMP DO schedule(static) private(jk,jj,ji)
	1004	DO jk = 1, jpk
	1005	DO jj = 1, jpj
	1006	DO ji = 1, jpi
	1007	en(ji,jj,jk) = rn_emin * tmask(ji,jj,jk)
	1008	END DO
	1009	END DO
	1010	END DO
	1011	!$OMP END DO NOWAIT
	1012	!$OMP DO schedule(static) private(jk)
[1481]	1013	DO jk = 1, jpk ! set the Kz to the background value
[7698]	1014	DO jj = 1, jpj
	1015	DO ji = 1, jpi
	1016	avt (ji,jj,jk) = avtb(jk) * wmask (ji,jj,jk)
	1017	avm (ji,jj,jk) = avmb(jk) * wmask (ji,jj,jk)
	1018	avmu(ji,jj,jk) = avmb(jk) * wumask(ji,jj,jk)
	1019	avmv(ji,jj,jk) = avmb(jk) * wvmask(ji,jj,jk)
	1020	END DO
	1021	END DO
[1481]	1022	END DO
[7698]	1023	!$OMP END PARALLEL
[1239]	1024	ENDIF
[1481]	1025	!
	1026	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	1027	! ! -------------------
[1531]	1028	IF(lwp) WRITE(numout,*) '---- tke-rst ----'
[3632]	1029	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
	1030	CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k )
	1031	CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k )
	1032	CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k )
	1033	CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k )
	1034	CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
[1481]	1035	!
[1239]	1036	ENDIF
	1037	!
[1531]	1038	END SUBROUTINE tke_rst
[1239]	1039
	1040	#else
	1041	!!----------------------------------------------------------------------
	1042	!! Dummy module : NO TKE scheme
	1043	!!----------------------------------------------------------------------
[1531]	1044	LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag
[1239]	1045	CONTAINS
[2528]	1046	SUBROUTINE zdf_tke_init ! Dummy routine
	1047	END SUBROUTINE zdf_tke_init
	1048	SUBROUTINE zdf_tke( kt ) ! Dummy routine
[1531]	1049	WRITE(,) 'zdf_tke: You should not have seen this print! error?', kt
	1050	END SUBROUTINE zdf_tke
	1051	SUBROUTINE tke_rst( kt, cdrw )
[1492]	1052	CHARACTER(len=*) :: cdrw
[1531]	1053	WRITE(,) 'tke_rst: You should not have seen this print! error?', kt, cdwr
	1054	END SUBROUTINE tke_rst
[1239]	1055	#endif
	1056
	1057	!!======================================================================
[1531]	1058	END MODULE zdftke

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

Download in other formats: