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

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 @ 3229

Last change on this file since 3229 was 3229, checked in by charris, 12 years ago
Added timing calls to most significant routines in LDF, SBC and ZDF.
Property svn:keywords set to `Id`
File size: 44.5 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
[1239]	28	!!----------------------------------------------------------------------
[1531]	29	#if defined key_zdftke \|\| defined key_esopa
[1239]	30	!!----------------------------------------------------------------------
[1531]	31	!! 'key_zdftke' TKE vertical physics
[1239]	32	!!----------------------------------------------------------------------
[2528]	33	!! zdf_tke : update momentum and tracer Kz from a tke scheme
	34	!! tke_tke : tke time stepping: update tke at now time step (en)
	35	!! tke_avn : compute mixing length scale and deduce avm and avt
	36	!! zdf_tke_init : initialization, namelist read, and parameters control
	37	!! tke_rst : read/write tke restart in ocean restart file
[1239]	38	!!----------------------------------------------------------------------
[2528]	39	USE oce ! ocean: dynamics and active tracers variables
	40	USE phycst ! physical constants
	41	USE dom_oce ! domain: ocean
	42	USE domvvl ! domain: variable volume layer
[1492]	43	USE sbc_oce ! surface boundary condition: ocean
[2528]	44	USE zdf_oce ! vertical physics: ocean variables
	45	USE zdfmxl ! vertical physics: mixed layer
	46	USE restart ! ocean restart
[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
[3186]	52	USE wrk_nemo ! work arrays
[3229]	53	USE timing ! Timing
[1239]	54
	55	IMPLICIT NONE
	56	PRIVATE
	57
[2528]	58	PUBLIC zdf_tke ! routine called in step module
	59	PUBLIC zdf_tke_init ! routine called in opa module
	60	PUBLIC tke_rst ! routine called in step module
[1239]	61
[2715]	62	LOGICAL , PUBLIC, PARAMETER :: lk_zdftke = .TRUE. !: TKE vertical mixing flag
[1239]	63
[2528]	64	! !! Namelist namzdf_tke
	65	LOGICAL :: ln_mxl0 = .FALSE. ! mixing length scale surface value as function of wind stress or not
	66	INTEGER :: nn_mxl = 2 ! type of mixing length (=0/1/2/3)
	67	REAL(wp) :: rn_mxl0 = 0.04_wp ! surface min value of mixing length (kappaz_o=0.40.1 m) [m]
	68	INTEGER :: nn_pdl = 1 ! Prandtl number or not (ratio avt/avm) (=0/1)
	69	REAL(wp) :: rn_ediff = 0.1_wp ! coefficient for avt: avt=rn_ediffmxlsqrt(e)
	70	REAL(wp) :: rn_ediss = 0.7_wp ! coefficient of the Kolmogoroff dissipation
	71	REAL(wp) :: rn_ebb = 3.75_wp ! coefficient of the surface input of tke
	72	REAL(wp) :: rn_emin = 0.7071e-6_wp ! minimum value of tke [m2/s2]
	73	REAL(wp) :: rn_emin0 = 1.e-4_wp ! surface minimum value of tke [m2/s2]
	74	REAL(wp) :: rn_bshear = 1.e-20_wp ! background shear (>0) currently a numerical threshold (do not change it)
	75	INTEGER :: nn_etau = 0 ! type of depth penetration of surface tke (=0/1/2/3)
	76	INTEGER :: nn_htau = 0 ! type of tke profile of penetration (=0/1)
	77	REAL(wp) :: rn_efr = 1.0_wp ! fraction of TKE surface value which penetrates in the ocean
	78	LOGICAL :: ln_lc = .FALSE. ! Langmuir cells (LC) as a source term of TKE or not
	79	REAL(wp) :: rn_lc = 0.15_wp ! coef to compute vertical velocity of Langmuir cells
[1239]	80
[2528]	81	REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
	82	REAL(wp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m]
	83	REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3)
	84	REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3)
[1239]	85
[2715]	86	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy [m2/s2]
	87	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau)
	88	REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation
	89	#if defined key_c1d
	90	! !! 1D cfg only ('key_c1d')
	91	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_dis, e_mix !: dissipation and mixing turbulent lengh scales
	92	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers
	93	#endif
[1492]	94
[1239]	95	!! * Substitutions
	96	# include "domzgr_substitute.h90"
	97	# include "vectopt_loop_substitute.h90"
	98	!!----------------------------------------------------------------------
[2715]	99	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[2528]	100	!! $Id$
	101	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1239]	102	!!----------------------------------------------------------------------
	103	CONTAINS
	104
[2715]	105	INTEGER FUNCTION zdf_tke_alloc()
	106	!!----------------------------------------------------------------------
	107	!! * FUNCTION zdf_tke_alloc *
	108	!!----------------------------------------------------------------------
	109	ALLOCATE( &
	110	#if defined key_c1d
	111	& e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , &
	112	& e_pdl(jpi,jpj,jpk) , e_ric(jpi,jpj,jpk) , &
	113	#endif
	114	& en (jpi,jpj,jpk) , htau (jpi,jpj) , dissl(jpi,jpj,jpk) , STAT= zdf_tke_alloc )
	115	!
	116	IF( lk_mpp ) CALL mpp_sum ( zdf_tke_alloc )
	117	IF( zdf_tke_alloc /= 0 ) CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
	118	!
	119	END FUNCTION zdf_tke_alloc
	120
	121
[1531]	122	SUBROUTINE zdf_tke( kt )
[1239]	123	!!----------------------------------------------------------------------
[1531]	124	!! * ROUTINE zdf_tke *
[1239]	125	!!
	126	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
[1492]	127	!! coefficients using a turbulent closure scheme (TKE).
[1239]	128	!!
[1492]	129	!! ** Method : The time evolution of the turbulent kinetic energy (tke)
	130	!! is computed from a prognostic equation :
	131	!! d(en)/dt = avm (d(u)/dz)**2 ! shear production
	132	!! + d( avm d(en)/dz )/dz ! diffusion of tke
	133	!! + avt N^2 ! stratif. destruc.
	134	!! - rn_ediss / emxl en**(2/3) ! Kolmogoroff dissipation
[1239]	135	!! with the boundary conditions:
[1695]	136	!! surface: en = max( rn_emin0, rn_ebb * taum )
[1239]	137	!! bottom : en = rn_emin
[1492]	138	!! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff)
	139	!!
	140	!! The now Turbulent kinetic energy is computed using the following
	141	!! time stepping: implicit for vertical diffusion term, linearized semi
	142	!! implicit for kolmogoroff dissipation term, and explicit forward for
	143	!! both buoyancy and shear production terms. Therefore a tridiagonal
	144	!! linear system is solved. Note that buoyancy and shear terms are
	145	!! discretized in a energy conserving form (Bruchard 2002).
	146	!!
	147	!! The dissipative and mixing length scale are computed from en and
	148	!! the stratification (see tke_avn)
	149	!!
	150	!! The now vertical eddy vicosity and diffusivity coefficients are
	151	!! given by:
	152	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
	153	!! avt = max( avmb, pdl * avm )
[1239]	154	!! eav = max( avmb, avm )
[1492]	155	!! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and
	156	!! given by an empirical funtion of the localRichardson number if nn_pdl=1
[1239]	157	!!
	158	!! ** Action : compute en (now turbulent kinetic energy)
	159	!! update avt, avmu, avmv (before vertical eddy coef.)
	160	!!
	161	!! References : Gaspar et al., JGR, 1990,
	162	!! Blanke and Delecluse, JPO, 1991
	163	!! Mellor and Blumberg, JPO 2004
	164	!! Axell, JGR, 2002
[1492]	165	!! Bruchard OM 2002
[1239]	166	!!----------------------------------------------------------------------
[1492]	167	INTEGER, INTENT(in) :: kt ! ocean time step
	168	!!----------------------------------------------------------------------
[1481]	169	!
[2528]	170	CALL tke_tke ! now tke (en)
[1492]	171	!
[2528]	172	CALL tke_avn ! now avt, avm, avmu, avmv
	173	!
[1531]	174	END SUBROUTINE zdf_tke
[1239]	175
[1492]	176
[1481]	177	SUBROUTINE tke_tke
[1239]	178	!!----------------------------------------------------------------------
[1492]	179	!! * ROUTINE tke_tke *
	180	!!
	181	!! ** Purpose : Compute the now Turbulente Kinetic Energy (TKE)
	182	!!
	183	!! ** Method : - TKE surface boundary condition
[2528]	184	!! - source term due to Langmuir cells (Axell JGR 2002) (ln_lc=T)
[1492]	185	!! - source term due to shear (saved in avmu, avmv arrays)
	186	!! - Now TKE : resolution of the TKE equation by inverting
	187	!! a tridiagonal linear system by a "methode de chasse"
	188	!! - increase TKE due to surface and internal wave breaking
	189	!!
	190	!! ** Action : - en : now turbulent kinetic energy)
	191	!! - avmu, avmv : production of TKE by shear at u and v-points
	192	!! (= Kz dz[Ub] * dz[Un] )
[1239]	193	!! ---------------------------------------------------------------------
[2977]	194	USE oce , ONLY: zdiag => ua ! (ua,va) used as workspace
	195	USE oce , ONLY: tsa ! (tsa) used as workspace
	196	!!
[1705]	197	INTEGER :: ji, jj, jk ! dummy loop arguments
[2528]	198	!!bfr INTEGER :: ikbu, ikbv, ikbum1, ikbvm1 ! temporary scalar
	199	!!bfr INTEGER :: ikbt, ikbumm1, ikbvmm1 ! temporary scalar
[1705]	200	REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3
	201	REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient
	202	REAL(wp) :: zbbrau, zesh2 ! temporary scalars
	203	REAL(wp) :: zfact1, zfact2, zfact3 ! - -
	204	REAL(wp) :: ztx2 , zty2 , zcof ! - -
	205	REAL(wp) :: ztau , zdif ! - -
	206	REAL(wp) :: zus , zwlc , zind ! - -
	207	REAL(wp) :: zzd_up, zzd_lw ! - -
[2528]	208	!!bfr REAL(wp) :: zebot ! - -
[2977]	209	REAL(wp), POINTER, DIMENSION(:,:,:) :: zd_up, zd_lw
[3155]	210	INTEGER , POINTER, DIMENSION(:,: ) :: imlc
	211	REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc
	212	REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc
[1239]	213	!!--------------------------------------------------------------------
[1492]	214	!
[3229]	215	IF( nn_timing == 1 ) CALL timing_start('tke_tke')
	216	!
[3155]	217	CALL wrk_alloc( jpi,jpj, imlc ) ! integer
	218	CALL wrk_alloc( jpi,jpj, zhlc )
	219	CALL wrk_alloc( jpi,jpj,jpk, zpelc )
[2977]	220	!
	221	zd_up => tsa(:,:,:,1)
	222	zd_lw => tsa(:,:,:,2)
[2715]	223
[1695]	224	zbbrau = rn_ebb / rau0 ! Local constant initialisation
[2528]	225	zfact1 = -.5_wp * rdt
	226	zfact2 = 1.5_wp * rdt * rn_ediss
	227	zfact3 = 0.5_wp * rn_ediss
[1492]	228	!
	229	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	230	! ! Surface boundary condition on tke
	231	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1695]	232	DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0)
[1481]	233	DO ji = fs_2, fs_jpim1 ! vector opt.
[1695]	234	en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
[1481]	235	END DO
	236	END DO
[2528]	237
	238	!!bfr - start commented area
[1492]	239	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	240	! ! Bottom boundary condition on tke
	241	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1719]	242	!
	243	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	244	! Tests to date have found the bottom boundary condition on tke to have very little effect.
	245	! The condition is coded here for completion but commented out until there is proof that the
	246	! computational cost is justified
	247	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	248	! en(bot) = (rn_ebb0/rau0)0.5sqrt(u_botfr^2+v_botfr^2) (min value rn_emin)
[1662]	249	!CDIR NOVERRCHK
[1719]	250	!! DO jj = 2, jpjm1
[1662]	251	!CDIR NOVERRCHK
[1719]	252	!! DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	253	!! ztx2 = bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) + &
	254	!! bfrua(ji ,jj) * ub(ji ,jj,mbku(ji ,jj) )
	255	!! zty2 = bfrva(ji,jj ) * vb(ji,jj ,mbkv(ji,jj )) + &
	256	!! bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1) )
[1719]	257	!! zebot = 0.001875_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000.
[2528]	258	!! en (ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * tmask(ji,jj,1)
[1719]	259	!! END DO
	260	!! END DO
[2528]	261	!!bfr - end commented area
[1492]	262	!
	263	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[2528]	264	IF( ln_lc ) THEN ! Langmuir circulation source term added to tke (Axell JGR 2002)
[1492]	265	! !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[1239]	266	!
[1492]	267	! !* total energy produce by LC : cumulative sum over jk
[2528]	268	zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * fsdepw(:,:,1) * fse3w(:,:,1)
[1239]	269	DO jk = 2, jpk
[2528]	270	zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * fsdepw(:,:,jk) * fse3w(:,:,jk)
[1239]	271	END DO
[1492]	272	! !* finite Langmuir Circulation depth
[1705]	273	zcof = 0.5 * 0.016 * 0.016 / ( zrhoa * zcdrag )
[2528]	274	imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land)
[1239]	275	DO jk = jpkm1, 2, -1
[1492]	276	DO jj = 1, jpj ! Last w-level at which zpelc>=0.5usus
	277	DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1)
[1705]	278	zus = zcof * taum(ji,jj)
[1239]	279	IF( zpelc(ji,jj,jk) > zus ) imlc(ji,jj) = jk
	280	END DO
	281	END DO
	282	END DO
[1492]	283	! ! finite LC depth
	284	# if defined key_vectopt_loop
	285	DO jj = 1, 1
	286	DO ji = 1, jpij ! vector opt. (forced unrolling)
	287	# else
	288	DO jj = 1, jpj
[1239]	289	DO ji = 1, jpi
[1492]	290	# endif
[1239]	291	zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj))
	292	END DO
	293	END DO
[1705]	294	zcof = 0.016 / SQRT( zrhoa * zcdrag )
[1239]	295	!CDIR NOVERRCHK
[1492]	296	DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en
[1239]	297	!CDIR NOVERRCHK
	298	DO jj = 2, jpjm1
	299	!CDIR NOVERRCHK
	300	DO ji = fs_2, fs_jpim1 ! vector opt.
[1705]	301	zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift
[1492]	302	! ! vertical velocity due to LC
[1239]	303	zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) )
	304	zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )
[1492]	305	! ! TKE Langmuir circulation source term
	306	en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)
[1239]	307	END DO
	308	END DO
	309	END DO
	310	!
	311	ENDIF
[1492]	312	!
	313	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	314	! ! Now Turbulent kinetic energy (output in en)
	315	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	316	! ! Resolution of a tridiagonal linear system by a "methode de chasse"
	317	! ! computation from level 2 to jpkm1 (e(1) already computed and e(jpk)=0 ).
	318	! ! zdiag : diagonal zd_up : upper diagonal zd_lw : lower diagonal
	319	!
	320	DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form)
	321	DO jj = 1, jpj ! here avmu, avmv used as workspace
	322	DO ji = 1, jpi
	323	avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) &
	324	& * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) &
	325	& / ( fse3uw_n(ji,jj,jk) &
	326	& * fse3uw_b(ji,jj,jk) )
	327	avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) &
	328	& * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) &
	329	& / ( fse3vw_n(ji,jj,jk) &
	330	& * fse3vw_b(ji,jj,jk) )
	331	END DO
	332	END DO
	333	END DO
	334	!
	335	DO jk = 2, jpkm1 !* Matrix and right hand side in en
[1239]	336	DO jj = 2, jpjm1
	337	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	338	zcof = zfact1 * tmask(ji,jj,jk)
	339	zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal
	340	& / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk ) )
	341	zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal
	342	& / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )
	343	! ! shear prod. at w-point weightened by mask
[2528]	344	zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
	345	& + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
[1492]	346	!
	347	zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw)
	348	zd_lw(ji,jj,jk) = zzd_lw
[2528]	349	zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)
[1239]	350	!
[1492]	351	! ! right hand side in en
[1481]	352	en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) &
	353	& + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) * tmask(ji,jj,jk)
[1239]	354	END DO
	355	END DO
	356	END DO
[1492]	357	! !* Matrix inversion from level 2 (tke prescribed at level 1)
[1239]	358	DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1
	359	DO jj = 2, jpjm1
	360	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	361	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]	362	END DO
	363	END DO
	364	END DO
	365	DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1
	366	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	367	zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke
[1239]	368	END DO
	369	END DO
	370	DO jk = 3, jpkm1
	371	DO jj = 2, jpjm1
	372	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	373	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]	374	END DO
	375	END DO
	376	END DO
	377	DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk
	378	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	379	en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1)
[1239]	380	END DO
	381	END DO
	382	DO jk = jpk-2, 2, -1
	383	DO jj = 2, jpjm1
	384	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	385	en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk)
[1239]	386	END DO
	387	END DO
	388	END DO
	389	DO jk = 2, jpkm1 ! set the minimum value of tke
	390	DO jj = 2, jpjm1
	391	DO ji = fs_2, fs_jpim1 ! vector opt.
	392	en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk)
	393	END DO
	394	END DO
	395	END DO
	396
[1492]	397	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	398	! ! TKE due to surface and internal wave breaking
	399	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[2528]	400	IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction)
[1492]	401	DO jk = 2, jpkm1
[1239]	402	DO jj = 2, jpjm1
	403	DO ji = fs_2, fs_jpim1 ! vector opt.
[1492]	404	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &
[2528]	405	& * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
[1239]	406	END DO
	407	END DO
[1492]	408	END DO
[2528]	409	ELSEIF( nn_etau == 2 ) THEN !* act only at the base of the mixed layer (jk=nmln) (rn_efr fraction)
[1492]	410	DO jj = 2, jpjm1
	411	DO ji = fs_2, fs_jpim1 ! vector opt.
	412	jk = nmln(ji,jj)
	413	en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &
[2528]	414	& * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
[1239]	415	END DO
[1492]	416	END DO
[2528]	417	ELSEIF( nn_etau == 3 ) THEN !* penetration belox the mixed layer (HF variability)
[1705]	418	!CDIR NOVERRCHK
	419	DO jk = 2, jpkm1
	420	!CDIR NOVERRCHK
	421	DO jj = 2, jpjm1
	422	!CDIR NOVERRCHK
	423	DO ji = fs_2, fs_jpim1 ! vector opt.
	424	ztx2 = utau(ji-1,jj ) + utau(ji,jj)
	425	zty2 = vtau(ji ,jj-1) + vtau(ji,jj)
[2528]	426	ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! module of the mean stress
	427	zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean
	428	zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications...
[1705]	429	en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &
[2528]	430	& * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk)
[1705]	431	END DO
	432	END DO
	433	END DO
[1239]	434	ENDIF
[1492]	435	CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)
	436	!
[3155]	437	CALL wrk_dealloc( jpi,jpj, imlc ) ! integer
	438	CALL wrk_dealloc( jpi,jpj, zhlc )
	439	CALL wrk_dealloc( jpi,jpj,jpk, zpelc )
[2715]	440	!
[3229]	441	IF( nn_timing == 1 ) CALL timing_stop('tke_tke')
	442	!
[1239]	443	END SUBROUTINE tke_tke
	444
[1492]	445
	446	SUBROUTINE tke_avn
[1239]	447	!!----------------------------------------------------------------------
[1492]	448	!! * ROUTINE tke_avn *
[1239]	449	!!
[1492]	450	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
	451	!!
	452	!! ** Method : At this stage, en, the now TKE, is known (computed in
	453	!! the tke_tke routine). First, the now mixing lenth is
	454	!! computed from en and the strafification (N^2), then the mixings
	455	!! coefficients are computed.
	456	!! - Mixing length : a first evaluation of the mixing lengh
	457	!! scales is:
	458	!! mxl = sqrt(2*en) / N
	459	!! where N is the brunt-vaisala frequency, with a minimum value set
[2528]	460	!! to rmxl_min (rn_mxl0) in the interior (surface) ocean.
[1492]	461	!! The mixing and dissipative length scale are bound as follow :
	462	!! nn_mxl=0 : mxl bounded by the distance to surface and bottom.
	463	!! zmxld = zmxlm = mxl
	464	!! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl
	465	!! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is
	466	!! less than 1 (\|d/dz(mxl)\|<1) and zmxld = zmxlm = mxl
	467	!! nn_mxl=3 : mxl is bounded from the surface to the bottom usings
	468	!! \|d/dz(xml)\|<1 to obtain lup, and from the bottom to
	469	!! the surface to obtain ldown. the resulting length
	470	!! scales are:
	471	!! zmxld = sqrt( lup * ldown )
	472	!! zmxlm = min ( lup , ldown )
	473	!! - Vertical eddy viscosity and diffusivity:
	474	!! avm = max( avtb, rn_ediff * zmxlm * en^1/2 )
	475	!! avt = max( avmb, pdlr * avm )
	476	!! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise.
	477	!!
	478	!! ** Action : - avt : now vertical eddy diffusivity (w-point)
	479	!! - avmu, avmv : now vertical eddy viscosity at uw- and vw-points
[1239]	480	!!----------------------------------------------------------------------
[2977]	481	USE oce, ONLY: zmpdl => ua ! ua used as workspace
	482	USE oce, ONLY: tsa ! use tsa as workspace
[1492]	483	!!
[2715]	484	INTEGER :: ji, jj, jk ! dummy loop indices
	485	REAL(wp) :: zrn2, zraug, zcoef, zav ! local scalars
	486	REAL(wp) :: zdku, zpdlr, zri, zsqen ! - -
	487	REAL(wp) :: zdkv, zemxl, zemlm, zemlp ! - -
[2977]	488	REAL(wp), POINTER, DIMENSION(:,:,:) :: zmxlm, zmxld
[1239]	489	!!--------------------------------------------------------------------
[2977]	490	!
[3229]	491	IF( nn_timing == 1 ) CALL timing_start('tke_avn')
	492	!
[2977]	493	zmxlm => tsa(:,:,:,1)
	494	zmxld => tsa(:,:,:,2)
[1239]	495
[1492]	496	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	497	! ! Mixing length
	498	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	499	!
	500	! !* Buoyancy length scale: l=sqrt(2e/n*2)
	501	!
[2528]	502	IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn2.e5taum/(rau0*g)
	503	zraug = vkarmn * 2.e5_wp / ( rau0 * grav )
	504	zmxlm(:,:,1) = MAX( rn_mxl0, zraug * taum(:,:) )
	505	ELSE ! surface set to the minimum value
	506	zmxlm(:,:,1) = rn_mxl0
[1239]	507	ENDIF
[2528]	508	zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value
[1239]	509	!
	510	!CDIR NOVERRCHK
[2528]	511	DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2)
[1239]	512	!CDIR NOVERRCHK
	513	DO jj = 2, jpjm1
	514	!CDIR NOVERRCHK
	515	DO ji = fs_2, fs_jpim1 ! vector opt.
	516	zrn2 = MAX( rn2(ji,jj,jk), rsmall )
[2528]	517	zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) )
[1239]	518	END DO
	519	END DO
	520	END DO
[1492]	521	!
	522	! !* Physical limits for the mixing length
	523	!
[2528]	524	zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value
	525	zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value
[1492]	526	!
[1239]	527	SELECT CASE ( nn_mxl )
	528	!
	529	CASE ( 0 ) ! bounded by the distance to surface and bottom
	530	DO jk = 2, jpkm1
	531	DO jj = 2, jpjm1
	532	DO ji = fs_2, fs_jpim1 ! vector opt.
	533	zemxl = MIN( fsdepw(ji,jj,jk), zmxlm(ji,jj,jk), &
[2528]	534	& fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) )
[1239]	535	zmxlm(ji,jj,jk) = zemxl
	536	zmxld(ji,jj,jk) = zemxl
	537	END DO
	538	END DO
	539	END DO
	540	!
	541	CASE ( 1 ) ! bounded by the vertical scale factor
	542	DO jk = 2, jpkm1
	543	DO jj = 2, jpjm1
	544	DO ji = fs_2, fs_jpim1 ! vector opt.
	545	zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) )
	546	zmxlm(ji,jj,jk) = zemxl
	547	zmxld(ji,jj,jk) = zemxl
	548	END DO
	549	END DO
	550	END DO
	551	!
	552	CASE ( 2 ) ! \|dk[xml]\| bounded by e3t :
	553	DO jk = 2, jpkm1 ! from the surface to the bottom :
	554	DO jj = 2, jpjm1
	555	DO ji = fs_2, fs_jpim1 ! vector opt.
	556	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
	557	END DO
	558	END DO
	559	END DO
	560	DO jk = jpkm1, 2, -1 ! from the bottom to the surface :
	561	DO jj = 2, jpjm1
	562	DO ji = fs_2, fs_jpim1 ! vector opt.
	563	zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
	564	zmxlm(ji,jj,jk) = zemxl
	565	zmxld(ji,jj,jk) = zemxl
	566	END DO
	567	END DO
	568	END DO
	569	!
	570	CASE ( 3 ) ! lup and ldown, \|dk[xml]\| bounded by e3t :
	571	DO jk = 2, jpkm1 ! from the surface to the bottom : lup
	572	DO jj = 2, jpjm1
	573	DO ji = fs_2, fs_jpim1 ! vector opt.
	574	zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )
	575	END DO
	576	END DO
	577	END DO
	578	DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown
	579	DO jj = 2, jpjm1
	580	DO ji = fs_2, fs_jpim1 ! vector opt.
	581	zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )
	582	END DO
	583	END DO
	584	END DO
	585	!CDIR NOVERRCHK
	586	DO jk = 2, jpkm1
	587	!CDIR NOVERRCHK
	588	DO jj = 2, jpjm1
	589	!CDIR NOVERRCHK
	590	DO ji = fs_2, fs_jpim1 ! vector opt.
	591	zemlm = MIN ( zmxld(ji,jj,jk), zmxlm(ji,jj,jk) )
	592	zemlp = SQRT( zmxld(ji,jj,jk) * zmxlm(ji,jj,jk) )
	593	zmxlm(ji,jj,jk) = zemlm
	594	zmxld(ji,jj,jk) = zemlp
	595	END DO
	596	END DO
	597	END DO
	598	!
	599	END SELECT
[1492]	600	!
[1239]	601	# if defined key_c1d
[1492]	602	e_dis(:,:,:) = zmxld(:,:,:) ! c1d configuration : save mixing and dissipation turbulent length scales
[1239]	603	e_mix(:,:,:) = zmxlm(:,:,:)
	604	# endif
	605
[1492]	606	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	607	! ! Vertical eddy viscosity and diffusivity (avmu, avmv, avt)
	608	! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[1239]	609	!CDIR NOVERRCHK
[1492]	610	DO jk = 1, jpkm1 !* vertical eddy viscosity & diffivity at w-points
[1239]	611	!CDIR NOVERRCHK
	612	DO jj = 2, jpjm1
	613	!CDIR NOVERRCHK
	614	DO ji = fs_2, fs_jpim1 ! vector opt.
	615	zsqen = SQRT( en(ji,jj,jk) )
	616	zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen
[1492]	617	avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)
	618	avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
[1239]	619	dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk)
	620	END DO
	621	END DO
	622	END DO
[1492]	623	CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged)
	624	!
	625	DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points
[1239]	626	DO jj = 2, jpjm1
	627	DO ji = fs_2, fs_jpim1 ! vector opt.
[1481]	628	avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk)
	629	avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk)
[1239]	630	END DO
	631	END DO
	632	END DO
	633	CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions
[1492]	634	!
	635	IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt
[1239]	636	DO jk = 2, jpkm1
	637	DO jj = 2, jpjm1
	638	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	639	zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk)
[1492]	640	! ! shear
	641	zdku = avmu(ji-1,jj,jk) * ( un(ji-1,jj,jk-1) - un(ji-1,jj,jk) ) * ( ub(ji-1,jj,jk-1) - ub(ji-1,jj,jk) ) &
	642	& + avmu(ji ,jj,jk) * ( un(ji ,jj,jk-1) - un(ji ,jj,jk) ) * ( ub(ji ,jj,jk-1) - ub(ji ,jj,jk) )
	643	zdkv = avmv(ji,jj-1,jk) * ( vn(ji,jj-1,jk-1) - vn(ji,jj-1,jk) ) * ( vb(ji,jj-1,jk-1) - vb(ji,jj-1,jk) ) &
	644	& + avmv(ji,jj ,jk) * ( vn(ji,jj ,jk-1) - vn(ji,jj ,jk) ) * ( vb(ji,jj ,jk-1) - vb(ji,jj ,jk) )
	645	! ! local Richardson number
[2528]	646	zri = MAX( rn2b(ji,jj,jk), 0._wp ) * zcoef / (zdku + zdkv + rn_bshear )
	647	zpdlr = MAX( 0.1_wp, 0.2 / MAX( 0.2 , zri ) )
[1492]	648	!!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!)
[2528]	649	!!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) )
[1492]	650	avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)
	651	# if defined key_c1d
	652	e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number
[1239]	653	e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri
	654	# endif
	655	END DO
	656	END DO
	657	END DO
	658	ENDIF
	659	CALL lbc_lnk( avt, 'W', 1. ) ! Lateral boundary conditions on avt (sign unchanged)
	660
	661	IF(ln_ctl) THEN
	662	CALL prt_ctl( tab3d_1=en , clinfo1=' tke - e: ', tab3d_2=avt, clinfo2=' t: ', ovlap=1, kdim=jpk)
	663	CALL prt_ctl( tab3d_1=avmu, clinfo1=' tke - u: ', mask1=umask, &
	664	& tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk )
	665	ENDIF
	666	!
[3229]	667	IF( nn_timing == 1 ) CALL timing_stop('tke_avn')
	668	!
[1492]	669	END SUBROUTINE tke_avn
[1239]	670
[1492]	671
[2528]	672	SUBROUTINE zdf_tke_init
[1239]	673	!!----------------------------------------------------------------------
[2528]	674	!! * ROUTINE zdf_tke_init *
[1239]	675	!!
	676	!! ** Purpose : Initialization of the vertical eddy diffivity and
[1492]	677	!! viscosity when using a tke turbulent closure scheme
[1239]	678	!!
[1601]	679	!! ** Method : Read the namzdf_tke namelist and check the parameters
[1492]	680	!! called at the first timestep (nit000)
[1239]	681	!!
[1601]	682	!! ** input : Namlist namzdf_tke
[1239]	683	!!
	684	!! ** Action : Increase by 1 the nstop flag is setting problem encounter
	685	!!----------------------------------------------------------------------
	686	INTEGER :: ji, jj, jk ! dummy loop indices
	687	!!
[2528]	688	NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , &
	689	& rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , &
	690	& rn_mxl0 , nn_pdl , ln_lc , rn_lc , &
	691	& nn_etau , nn_htau , rn_efr
[1239]	692	!!----------------------------------------------------------------------
[2715]	693	!
[1601]	694	REWIND ( numnam ) !* Read Namelist namzdf_tke : Turbulente Kinetic Energy
	695	READ ( numnam, namzdf_tke )
[2715]	696	!
[2528]	697	ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number
	698	rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity
[2715]	699	!
[1492]	700	IF(lwp) THEN !* Control print
[1239]	701	WRITE(numout,*)
[2528]	702	WRITE(numout,*) 'zdf_tke_init : tke turbulent closure scheme - initialisation'
	703	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	704	WRITE(numout,*) ' Namelist namzdf_tke : set tke mixing parameters'
[1705]	705	WRITE(numout,*) ' coef. to compute avt rn_ediff = ', rn_ediff
	706	WRITE(numout,*) ' Kolmogoroff dissipation coef. rn_ediss = ', rn_ediss
	707	WRITE(numout,*) ' tke surface input coef. rn_ebb = ', rn_ebb
	708	WRITE(numout,*) ' minimum value of tke rn_emin = ', rn_emin
	709	WRITE(numout,*) ' surface minimum value of tke rn_emin0 = ', rn_emin0
	710	WRITE(numout,*) ' background shear (>0) rn_bshear = ', rn_bshear
	711	WRITE(numout,*) ' mixing length type nn_mxl = ', nn_mxl
	712	WRITE(numout,*) ' prandl number flag nn_pdl = ', nn_pdl
	713	WRITE(numout,*) ' surface mixing length = F(stress) or not ln_mxl0 = ', ln_mxl0
[2528]	714	WRITE(numout,*) ' surface mixing length minimum value rn_mxl0 = ', rn_mxl0
	715	WRITE(numout,*) ' flag to take into acc. Langmuir circ. ln_lc = ', ln_lc
	716	WRITE(numout,*) ' coef to compute verticla velocity of LC rn_lc = ', rn_lc
[1705]	717	WRITE(numout,*) ' test param. to add tke induced by wind nn_etau = ', nn_etau
	718	WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau
	719	WRITE(numout,*) ' fraction of en which pene. the thermocline rn_efr = ', rn_efr
[1239]	720	WRITE(numout,*)
[1601]	721	WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri
[1239]	722	ENDIF
[2715]	723	!
	724	! ! allocate tke arrays
	725	IF( zdf_tke_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tke_init : unable to allocate arrays' )
	726	!
[1492]	727	! !* Check of some namelist values
[2528]	728	IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' )
	729	IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' )
	730	IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' )
	731	#if ! key_coupled
	732	IF( nn_etau == 3 ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )
	733	#endif
[1239]	734
[2528]	735	IF( ln_mxl0 ) THEN
	736	IF(lwp) WRITE(numout,*) ' use a surface mixing length = F(stress) : set rn_mxl0 = rmxl_min'
	737	rn_mxl0 = rmxl_min
	738	ENDIF
	739
[1492]	740	IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln
[1239]	741
[1492]	742	! !* depth of penetration of surface tke
	743	IF( nn_etau /= 0 ) THEN
[1601]	744	SELECT CASE( nn_htau ) ! Choice of the depth of penetration
[2528]	745	CASE( 0 ) ! constant depth penetration (here 10 meters)
	746	htau(:,:) = 10._wp
	747	CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees
	748	htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) )
[1492]	749	END SELECT
	750	ENDIF
	751	! !* set vertical eddy coef. to the background value
[1239]	752	DO jk = 1, jpk
	753	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
[1481]	754	avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
[1239]	755	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
	756	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
	757	END DO
[2528]	758	dissl(:,:,:) = 1.e-12_wp
[2715]	759	!
	760	CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files
[1239]	761	!
[2528]	762	END SUBROUTINE zdf_tke_init
[1239]	763
	764
[1531]	765	SUBROUTINE tke_rst( kt, cdrw )
[1239]	766	!!---------------------------------------------------------------------
[1531]	767	!! * ROUTINE tke_rst *
[1239]	768	!!
	769	!! ** Purpose : Read or write TKE file (en) in restart file
	770	!!
	771	!! ** Method : use of IOM library
	772	!! if the restart does not contain TKE, en is either
[1537]	773	!! set to rn_emin or recomputed
[1239]	774	!!----------------------------------------------------------------------
[2715]	775	INTEGER , INTENT(in) :: kt ! ocean time-step
	776	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
[1239]	777	!
[1481]	778	INTEGER :: jit, jk ! dummy loop indices
[2715]	779	INTEGER :: id1, id2, id3, id4, id5, id6 ! local integers
[1239]	780	!!----------------------------------------------------------------------
	781	!
[1481]	782	IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
	783	! ! ---------------
	784	IF( ln_rstart ) THEN !* Read the restart file
	785	id1 = iom_varid( numror, 'en' , ldstop = .FALSE. )
	786	id2 = iom_varid( numror, 'avt' , ldstop = .FALSE. )
	787	id3 = iom_varid( numror, 'avm' , ldstop = .FALSE. )
	788	id4 = iom_varid( numror, 'avmu' , ldstop = .FALSE. )
	789	id5 = iom_varid( numror, 'avmv' , ldstop = .FALSE. )
	790	id6 = iom_varid( numror, 'dissl', ldstop = .FALSE. )
	791	!
	792	IF( id1 > 0 ) THEN ! 'en' exists
[1239]	793	CALL iom_get( numror, jpdom_autoglo, 'en', en )
[1481]	794	IF( MIN( id2, id3, id4, id5, id6 ) > 0 ) THEN ! all required arrays exist
	795	CALL iom_get( numror, jpdom_autoglo, 'avt' , avt )
	796	CALL iom_get( numror, jpdom_autoglo, 'avm' , avm )
	797	CALL iom_get( numror, jpdom_autoglo, 'avmu' , avmu )
	798	CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv )
	799	CALL iom_get( numror, jpdom_autoglo, 'dissl', dissl )
[1492]	800	ELSE ! one at least array is missing
	801	CALL tke_avn ! compute avt, avm, avmu, avmv and dissl (approximation)
[1481]	802	ENDIF
	803	ELSE ! No TKE array found: initialisation
	804	IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without tke scheme, en computed by iterative loop'
[1239]	805	en (:,:,:) = rn_emin * tmask(:,:,:)
[1492]	806	CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation)
[1531]	807	DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO
[1239]	808	ENDIF
[1481]	809	ELSE !* Start from rest
	810	en(:,:,:) = rn_emin * tmask(:,:,:)
	811	DO jk = 1, jpk ! set the Kz to the background value
	812	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
	813	avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)
	814	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
	815	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
	816	END DO
[1239]	817	ENDIF
[1481]	818	!
	819	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
	820	! ! -------------------
[1531]	821	IF(lwp) WRITE(numout,*) '---- tke-rst ----'
[1601]	822	CALL iom_rstput( kt, nitrst, numrow, 'en' , en )
	823	CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt )
	824	CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm )
	825	CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu )
	826	CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv )
	827	CALL iom_rstput( kt, nitrst, numrow, 'dissl', dissl )
[1481]	828	!
[1239]	829	ENDIF
	830	!
[1531]	831	END SUBROUTINE tke_rst
[1239]	832
	833	#else
	834	!!----------------------------------------------------------------------
	835	!! Dummy module : NO TKE scheme
	836	!!----------------------------------------------------------------------
[1531]	837	LOGICAL, PUBLIC, PARAMETER :: lk_zdftke = .FALSE. !: TKE flag
[1239]	838	CONTAINS
[2528]	839	SUBROUTINE zdf_tke_init ! Dummy routine
	840	END SUBROUTINE zdf_tke_init
	841	SUBROUTINE zdf_tke( kt ) ! Dummy routine
[1531]	842	WRITE(,) 'zdf_tke: You should not have seen this print! error?', kt
	843	END SUBROUTINE zdf_tke
	844	SUBROUTINE tke_rst( kt, cdrw )
[1492]	845	CHARACTER(len=*) :: cdrw
[1531]	846	WRITE(,) 'tke_rst: You should not have seen this print! error?', kt, cdwr
	847	END SUBROUTINE tke_rst
[1239]	848	#endif
	849
	850	!!======================================================================
[1531]	851	END MODULE zdftke

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

Download in other formats: