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.

zdfkpp.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/zdfkpp.F90 @ 3281

Last change on this file since 3281 was 3229, checked in by charris, 13 years ago
Added timing calls to most significant routines in LDF, SBC and ZDF.
Property svn:keywords set to `Id`
File size: 78.7 KB

Rev	Line
[255]	1	MODULE zdfkpp
	2	!!======================================================================
	3	!! * MODULE zdfkpp *
	4	!! Ocean physics: vertical mixing coefficient compute from the KPP
	5	!! turbulent closure parameterization
	6	!!=====================================================================
[2528]	7	!! History : OPA ! 2000-03 (W.G. Large, J. Chanut) Original code
	8	!! 8.1 ! 2002-06 (J.M. Molines) for real case CLIPPER
	9	!! 8.2 ! 2003-10 (Chanut J.) re-writting
	10	!! NEMO 1.0 ! 2005-01 (C. Ethe, G. Madec) Free form, F90 + creation of tra_kpp routine
[2715]	11	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase + merge TRC-TRA
[503]	12	!!----------------------------------------------------------------------
[255]	13	#if defined key_zdfkpp \|\| defined key_esopa
	14	!!----------------------------------------------------------------------
	15	!! 'key_zdfkpp' KPP scheme
	16	!!----------------------------------------------------------------------
	17	!! zdf_kpp : update momentum and tracer Kz from a kpp scheme
	18	!! zdf_kpp_init : initialization, namelist read, and parameters control
[2528]	19	!! tra_kpp : compute and add to the T & S trend the non-local flux
	20	!! trc_kpp : compute and add to the passive tracer trend the non-local flux (lk_top=T)
[255]	21	!!----------------------------------------------------------------------
	22	USE oce ! ocean dynamics and active tracers
	23	USE dom_oce ! ocean space and time domain
	24	USE zdf_oce ! ocean vertical physics
[888]	25	USE sbc_oce ! surface boundary condition: ocean
[255]	26	USE phycst ! physical constants
	27	USE eosbn2 ! equation of state
	28	USE zdfddm ! double diffusion mixing
[463]	29	USE in_out_manager ! I/O manager
[2715]	30	USE lib_mpp ! MPP library
[3186]	31	USE wrk_nemo ! work arrays
[463]	32	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	33	USE prtctl ! Print control
[2528]	34	USE trdmod_oce ! ocean trends definition
	35	USE trdtra ! tracers trends
[3229]	36	USE timing ! Timing
[255]	37
	38	IMPLICIT NONE
	39	PRIVATE
	40
[2528]	41	PUBLIC zdf_kpp ! routine called by step.F90
	42	PUBLIC zdf_kpp_init ! routine called by opa.F90
	43	PUBLIC tra_kpp ! routine called by step.F90
	44	PUBLIC trc_kpp ! routine called by trcstp.F90
[255]	45
[1537]	46	LOGICAL , PUBLIC, PARAMETER :: lk_zdfkpp = .TRUE. !: KPP vertical mixing flag
	47
[2715]	48	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghats !: non-local scalar mixing term (gamma/<ws>o)
	49	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wt0 !: surface temperature flux for non local flux
	50	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ws0 !: surface salinity flux for non local flux
	51	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hkpp !: boundary layer depth
[503]	52
[1601]	53	! !!* Namelist namzdf_kpp *
[1537]	54	REAL(wp) :: rn_difmiw = 1.2e-04_wp ! constant internal wave viscosity (m2/s)
	55	REAL(wp) :: rn_difsiw = 1.2e-05_wp ! constant internal wave diffusivity (m2/s)
	56	REAL(wp) :: rn_riinfty = 0.8_wp ! local Richardson Number limit for shear instability
	57	REAL(wp) :: rn_difri = 5.e-03_wp ! maximum shear mixing at Rig = 0 (m2/s)
	58	REAL(wp) :: rn_bvsqcon = -1.e-09_wp ! Brunt-Vaisala squared (1/s**2) for maximum convection
	59	REAL(wp) :: rn_difcon = 1._wp ! maximum mixing in interior convection (m2/s)
	60	INTEGER :: nn_ave = 1 ! = 0/1 flag for horizontal average on avt, avmu, avmv
[255]	61
	62	#if defined key_zdfddm
[2528]	63	! !!! ** Double diffusion Mixing
	64	REAL(wp) :: difssf = 1.e-03_wp ! maximum salt fingering mixing
	65	REAL(wp) :: Rrho0 = 1.9_wp ! limit for salt fingering mixing
	66	REAL(wp) :: difsdc = 1.5e-06_wp ! maximum diffusive convection mixing
[255]	67	#endif
[1601]	68	LOGICAL :: ln_kpprimix = .TRUE. ! Shear instability mixing
[255]	69
[2528]	70	! !!! General constants
	71	REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number
	72	REAL(wp) :: pthird = 1._wp/3._wp ! 1/3
	73	REAL(wp) :: pfourth = 1._wp/4._wp ! 1/4
[255]	74
[2528]	75	! !!! Boundary Layer Turbulence Parameters
	76	REAL(wp) :: vonk = 0.4_wp ! von Karman's constant
	77	REAL(wp) :: epsilon = 0.1_wp ! nondimensional extent of the surface layer
	78	REAL(wp) :: rconc1 = 5.0_wp ! standard flux profile function parmaeters
	79	REAL(wp) :: rconc2 = 16.0_wp ! " "
	80	REAL(wp) :: rconcm = 8.38_wp ! momentum flux profile fit
	81	REAL(wp) :: rconam = 1.26_wp ! " "
	82	REAL(wp) :: rzetam = -.20_wp ! " "
	83	REAL(wp) :: rconcs = 98.96_wp ! scalar flux profile fit
	84	REAL(wp) :: rconas = -28.86_wp ! " "
	85	REAL(wp) :: rzetas = -1.0_wp ! " "
	86
	87	! !!! Boundary Layer Depth Diagnostic
	88	REAL(wp) :: Ricr = 0.3_wp ! critical bulk Richardson Number
	89	REAL(wp) :: rcekman = 0.7_wp ! coefficient for ekman depth
	90	REAL(wp) :: rcmonob = 1.0_wp ! coefficient for Monin-Obukhov depth
	91	REAL(wp) :: rconcv = 1.7_wp ! ratio of interior buoyancy frequency to its value at entrainment depth
	92	REAL(wp) :: hbf = 1.0_wp ! fraction of bound. layer depth to which absorbed solar
	93	! ! rad. and contributes to surf. buo. forcing
	94	REAL(wp) :: Vtc ! function of rconcv,rconcs,epsilon,vonk,Ricr
	95
	96	! !!! Nonlocal Boundary Layer Mixing
	97	REAL(wp) :: rcstar = 5.0_wp ! coefficient for convective nonlocal transport
	98	REAL(wp) :: rcs = 1.0e-3_wp ! conversion: mm/s ==> m/s
	99	REAL(wp) :: rcg ! non-dimensional coefficient for nonlocal transport
[255]	100
	101	#if ! defined key_kppcustom
[2715]	102	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: del ! array for reference mean values of vertical integration
[255]	103	#endif
	104
	105	#if defined key_kpplktb
[2528]	106	! !!! Parameters for lookup table for turbulent velocity scales
	107	INTEGER, PARAMETER :: nilktb = 892 ! number of values for zehat in KPP lookup table
	108	INTEGER, PARAMETER :: njlktb = 482 ! number of values for ustar in KPP lookup table
	109	INTEGER, PARAMETER :: nilktbm1 = nilktb-1 !
	110	INTEGER, PARAMETER :: njlktbm1 = njlktb-1 !
[255]	111
[2528]	112	REAL(wp), DIMENSION(nilktb,njlktb) :: wmlktb ! lookup table for the turbulent vertical velocity scale (momentum)
	113	REAL(wp), DIMENSION(nilktb,njlktb) :: wslktb ! lookup table for the turbulent vertical velocity scale (tracers)
[255]	114
[2528]	115	REAL(wp) :: dehatmin = -4.e-7_wp ! minimum limit for zhat in lookup table (m3/s3)
	116	REAL(wp) :: dehatmax = 0._wp ! maximum limit for zhat in lookup table (m3/s3)
	117	REAL(wp) :: ustmin = 0._wp ! minimum limit for ustar in lookup table (m/s)
	118	REAL(wp) :: ustmax = 0.04_wp ! maximum limit for ustar in lookup table (m/s)
	119	REAL(wp) :: dezehat ! delta zhat in lookup table
	120	REAL(wp) :: deustar ! delta ustar in lookup table
[255]	121	#endif
[2715]	122	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: ratt ! attenuation coef (already defines in module traqsr,
[1601]	123	! ! but only if the solar radiation penetration is considered)
[2528]	124
	125	! !!! * penetrative solar radiation coefficient *
	126	REAL(wp) :: rabs = 0.58_wp ! fraction associated with xsi1
	127	REAL(wp) :: xsi1 = 0.35_wp ! first depth of extinction
	128	REAL(wp) :: xsi2 = 23.0_wp ! second depth of extinction
[255]	129	! ! (default values: water type Ib)
	130
[2715]	131	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean, eumean, evmean ! coeff. used for hor. smoothing at t-, u- & v-points
[2528]	132
[899]	133	#if defined key_c1d
[2715]	134	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rig !: gradient Richardson number
	135	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rib !: bulk Richardson number
	136	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: buof !: buoyancy forcing
	137	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mols !: moning-Obukhov length scale
	138	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ekdp !: Ekman depth
[255]	139	#endif
	140
[1601]	141	INTEGER :: jip = 62 , jjp = 111
[255]	142
	143	!! * Substitutions
	144	# include "domzgr_substitute.h90"
	145	# include "vectopt_loop_substitute.h90"
[896]	146	# include "zdfddm_substitute.h90"
[255]	147	!!----------------------------------------------------------------------
[2715]	148	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
[888]	149	!! $Id$
[2715]	150	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[255]	151	!!----------------------------------------------------------------------
	152	CONTAINS
	153
[2715]	154	INTEGER FUNCTION zdf_kpp_alloc()
	155	!!----------------------------------------------------------------------
	156	!! * FUNCTION zdf_kpp_alloc *
	157	!!----------------------------------------------------------------------
	158	ALLOCATE( ghats(jpi,jpj,jpk), wt0(jpi,jpj), ws0(jpi,jpj), hkpp(jpi,jpj), &
	159	#if ! defined key_kpplktb
	160	& del(jpk,jpk), &
	161	#endif
	162	& ratt(jpk), &
	163	& etmean(jpi,jpj,jpk), eumean(jpi,jpj,jpk), evmean(jpi,jpj,jpk), &
	164	#if defined key_c1d
	165	& rig (jpi,jpj,jpk), rib(jpi,jpj,jpk), buof(jpi,jpj,jpk), &
	166	& mols(jpi,jpj,jpk), ekdp(jpi,jpj), &
	167	#endif
	168	& STAT= zdf_kpp_alloc )
	169	!
	170	IF( lk_mpp ) CALL mpp_sum ( zdf_kpp_alloc )
	171	IF( zdf_kpp_alloc /= 0 ) CALL ctl_warn('zdf_kpp_alloc: failed to allocate arrays')
	172	END FUNCTION zdf_kpp_alloc
	173
	174
[2528]	175	SUBROUTINE zdf_kpp( kt )
[255]	176	!!----------------------------------------------------------------------
	177	!! * ROUTINE zdf_kpp *
	178	!!
	179	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
	180	!! coefficients and non local mixing using K-profile parameterization
	181	!!
	182	!! ** Method : The boundary layer depth hkpp is diagnosed at tracer points
	183	!! from profiles of buoyancy, and shear, and the surface forcing.
	184	!! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
	185	!!
	186	!! Kx = hkpp Wx(sigma) G(sigma)
	187	!!
	188	!! and the non local term ghat = Cs / Ws(sigma) / hkpp
	189	!! Below hkpp the coefficients are the sum of mixing due to internal waves
	190	!! shear instability and double diffusion.
	191	!!
	192	!! -1- Compute the now interior vertical mixing coefficients at all depths.
	193	!! -2- Diagnose the boundary layer depth.
	194	!! -3- Compute the now boundary layer vertical mixing coefficients.
	195	!! -4- Compute the now vertical eddy vicosity and diffusivity.
	196	!! -5- Smoothing
	197	!!
	198	!! N.B. The computation is done from jk=2 to jpkm1
	199	!! Surface value of avt avmu avmv are set once a time to zero
	200	!! in routine zdf_kpp_init.
	201	!!
	202	!! ** Action : update the non-local terms ghats
	203	!! update avt, avmu, avmv (before vertical eddy coef.)
	204	!!
[503]	205	!! References : Large W.G., Mc Williams J.C. and Doney S.C.
[255]	206	!! Reviews of Geophysics, 32, 4, November 1994
	207	!! Comments in the code refer to this paper, particularly
	208	!! the equation number. (LMD94, here after)
	209	!!----------------------------------------------------------------------
[2528]	210	USE oce , zviscos => ua ! temp. array for viscosities use ua as workspace
[2977]	211	USE oce , zdiffut => va ! temp. array for diffusivities use sa as workspace
[503]	212	!!
	213	INTEGER, INTENT( in ) :: kt ! ocean time step
	214	!!
	215	INTEGER :: ji, jj, jk ! dummy loop indices
	216	INTEGER :: ikbot, jkmax, jkm1, jkp2 !
[255]	217
[2528]	218	REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube !
	219	REAL(wp) :: zrhos, zalbet, zbeta, zthermal, zhalin, zatt1 !
	220	REAL(wp) :: zref, zt, zs, zh, zu, zv, zrh ! Bulk richardson number
	221	REAL(wp) :: zrib, zrinum, zdVsq, zVtsq !
	222	REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
[255]	223	#if defined key_kpplktb
[2528]	224	INTEGER :: il, jl ! Lookup table or Analytical functions
	225	REAL(wp) :: ud, zfrac, ufrac, zwam, zwbm, zwas, zwbs !
[255]	226	#else
[2528]	227	REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm !
[255]	228	#endif
[2528]	229	REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density
[255]	230	#if ! defined key_kppcustom
[2528]	231	INTEGER :: jm ! dummy loop indices
	232	REAL(wp) :: zr1, zr2, zr3, zr4, zrhop ! Compression terms
[255]	233	#endif
[2528]	234	REAL(wp) :: zflag, ztemp, zrn2, zdep21, zdep32, zdep43
	235	REAL(wp) :: zdku2, zdkv2, ze3sqr, zsh2, zri, zfri ! Interior richardson mixing
[3155]	236	REAL(wp), POINTER, DIMENSION(:,:) :: zmoek ! Moning-Obukov limitation
[2715]	237	REAL(wp), POINTER, DIMENSION(:) :: zmoa, zekman
	238	REAL(wp) :: zmob, zek
	239	REAL(wp), POINTER, DIMENSION(:,:) :: zdepw, zdift, zvisc ! The pipe
	240	REAL(wp), POINTER, DIMENSION(:,:) :: zdept
	241	REAL(wp), POINTER, DIMENSION(:,:) :: zriblk
	242	REAL(wp), POINTER, DIMENSION(:) :: zhmax, zria, zhbl
[2528]	243	REAL(wp) :: zflagri, zflagek, zflagmo, zflagh, zflagkb !
[2715]	244	REAL(wp), POINTER, DIMENSION(:) :: za2m, za3m, zkmpm, za2t, za3t, zkmpt ! Shape function (G)
[2528]	245	REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
[255]	246	#if defined key_zdfddm
[2528]	247	REAL(wp) :: zrrau, zds, zavdds, zavddt,zinr ! double diffusion mixing
[2715]	248	REAL(wp), POINTER, DIMENSION(:,:) :: zdifs
	249	REAL(wp), POINTER, DIMENSION(:) :: za2s, za3s, zkmps
[2977]	250	REAL(wp) :: zkm1s
[3155]	251	REAL(wp), POINTER, DIMENSION(:,:) :: zblcs
	252	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdiffus
[255]	253	#endif
[3155]	254	REAL(wp), POINTER, DIMENSION(:,:) :: zBo, zBosol, zustar ! Surface buoyancy forcing, friction velocity
	255	REAL(wp), POINTER, DIMENSION(:,:) :: zmask, zblcm, zblct
[255]	256	!!--------------------------------------------------------------------
[3229]	257	!
	258	IF( nn_timing == 1 ) CALL timing_start('zdf_kpp')
	259	!
[3155]	260	CALL wrk_alloc( jpi, zmoa, zekman, zhmax, zria, zhbl )
	261	CALL wrk_alloc( jpi, za2m, za3m, zkmpm, za2t, za3t, zkmpt )
	262	CALL wrk_alloc( jpi,2, zriblk )
	263	CALL wrk_alloc( jpi,3, zmoek, kjstart = 0 )
	264	CALL wrk_alloc( jpi,3, zdept )
	265	CALL wrk_alloc( jpi,4, zdepw, zdift, zvisc )
	266	CALL wrk_alloc( jpi,jpj, zBo, zBosol, zustar )
	267	CALL wrk_alloc( jpi,jpk, zmask, zblcm, zblct, zblcs )
[2715]	268	#if defined key_zdfddm
[3155]	269	CALL wrk_alloc( jpi,4, zdifs )
	270	CALL wrk_alloc( jpi, zmoa, za2s, za3s, zkmps )
	271	CALL wrk_alloc( jpi,jpk, zblcs )
	272	CALL wrk_alloc( jpi,jpi,jpk, zdiffus )
[2715]	273	#endif
	274
[255]	275	zviscos(:,:,:) = 0.
	276	zblcm (:,: ) = 0.
	277	zdiffut(:,:,:) = 0.
	278	zblct (:,: ) = 0.
	279	#if defined key_zdfddm
	280	zdiffus(:,:,:) = 0.
	281	zblcs (:,: ) = 0.
	282	#endif
	283	ghats(:,:,:) = 0.
	284
	285	zBo (:,:) = 0.
	286	zBosol(:,:) = 0.
	287	zustar(:,:) = 0.
	288
	289
	290	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	291	! I. Interior diffusivity and viscosity at w points ( T interfaces)
	292	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	293	DO jk = 2, jpkm1
	294	DO jj = 2, jpjm1
	295	DO ji = fs_2, fs_jpim1
	296	! Mixing due to internal waves breaking
	297	! -------------------------------------
[1537]	298	avmu(ji,jj,jk) = rn_difmiw
	299	avt (ji,jj,jk) = rn_difsiw
[255]	300	! Mixing due to vertical shear instability
	301	! -------------------------------------
	302	IF( ln_kpprimix ) THEN
	303	! Compute the gradient Richardson number at interfaces (zri):
	304	! LMD94, eq. 27 (is vertical smoothing needed : Rig=N^2 / (dz(u))^2 + (dz(v))^2
	305	zdku2 = ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
	306	& * ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
	307	& + ( un(ji ,jj,jk - 1) - un(ji ,jj,jk) ) &
	308	& * ( un(ji ,jj,jk - 1) - un(ji ,jj,jk) )
	309
	310	zdkv2 = ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
	311	& * ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
	312	& + ( vn(ji, jj,jk - 1) - vn(ji, jj,jk) ) &
	313	& * ( vn(ji, jj,jk - 1) - vn(ji, jj,jk) )
	314
	315	ze3sqr = 1. / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
	316	! Square of vertical shear at interfaces
[1082]	317	zsh2 = 0.5 * ( zdku2 + zdkv2 ) * ze3sqr
[255]	318	zri = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + epsln )
[899]	319	#if defined key_c1d
[255]	320	! save the gradient richardson number
	321	rig(ji,jj,jk) = zri * tmask(ji,jj,jk)
	322	#endif
	323	! Evaluate f of Ri (zri) for shear instability store in zfri
	324	! LMD94, eq. 28a,b,c, figure 3 ; Rem: p1 is 3, hard coded
	325	zfri = MAX( zri , 0. )
[1537]	326	zfri = MIN( zfri / rn_riinfty , 1.0 )
[255]	327	zfri = ( 1.0 - zfri * zfri )
	328	zfri = zfri * zfri * zfri
	329	! add shear contribution to mixing coef.
[1537]	330	avmu(ji,jj,jk) = avmu(ji,jj,jk) + rn_difri * zfri
	331	avt (ji,jj,jk) = avt (ji,jj,jk) + rn_difri * zfri
[255]	332	ENDIF
	333	#if defined key_zdfddm
	334	avs (ji,jj,jk) = avt (ji,jj,jk)
	335	! Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90
	336	! ------------------------------------------------------------------
	337	! only retains positive value of rrau
	338	zrrau = MAX( rrau(ji,jj,jk), epsln )
[2977]	339	zds = tsn(ji,jj,jk-1,jp_sal) - tsn(ji,jj,jk,jp_sal)
[255]	340	IF( zrrau > 1. .AND. zds > 0.) THEN
	341	!
	342	! Salt fingering case.
	343	!---------------------
	344	! Compute interior diffusivity for double diffusive mixing of
	345	! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
	346	! After that set interior diffusivity for double diffusive mixing
	347	! of temperature
	348	zavdds = MIN( zrrau, Rrho0 )
	349	zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 )
	350	zavdds = 1.0 - zavdds * zavdds
	351	zavdds = zavdds * zavdds * zavdds
	352	zavdds = difssf * zavdds
	353	zavddt = 0.7 * zavdds
	354	ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN
	355	!
	356	! Diffusive convection case.
	357	!---------------------------
	358	! Compute interior diffusivity for double diffusive mixing of
	359	! temperature (Marmorino and Caldwell, 1976);
	360	! Compute interior diffusivity for double diffusive mixing of salinity
	361	zinr = 1. / zrrau
	362	zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) )
	363	zavddt = difsdc * zavddt
	364	IF( zrrau < 0.5) THEN
	365	zavdds = zavddt * 0.15 * zrrau
	366	ELSE
	367	zavdds = zavddt * (1.85 * zrrau - 0.85 )
	368	ENDIF
	369	ELSE
	370	zavddt = 0.
	371	zavdds = 0.
	372	ENDIF
	373	! Add double diffusion contribution to temperature and salinity mixing coefficients.
	374	avt (ji,jj,jk) = avt (ji,jj,jk) + zavddt
	375	avs (ji,jj,jk) = avs (ji,jj,jk) + zavdds
	376	#endif
	377	END DO
	378	END DO
	379	END DO
	380
	381
	382	! Radiative (zBosol) and non radiative (zBo) surface buoyancy
	383	!JMM at the time zdfkpp is called, q still holds the sum q + qsr
	384	!---------------------------------------------------------------------
	385	DO jj = 2, jpjm1
	386	DO ji = fs_2, fs_jpim1
[1601]	387	IF( nn_eos < 1) THEN
[2977]	388	zt = tsn(ji,jj,1,jp_tem)
	389	zs = tsn(ji,jj,1,jp_sal) - 35.0
[255]	390	zh = fsdept(ji,jj,1)
	391	! potential volumic mass
	392	zrhos = rhop(ji,jj,1)
	393	zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta
	394	& - 0.203814e-03 ) * zt &
	395	& + 0.170907e-01 ) * zt &
	396	& + 0.665157e-01 &
	397	& + ( - 0.678662e-05 * zs &
	398	& - 0.846960e-04 * zt + 0.378110e-02 ) * zs &
	399	& + ( ( - 0.302285e-13 * zh &
	400	& - 0.251520e-11 * zs &
	401	& + 0.512857e-12 * zt * zt ) * zh &
	402	& - 0.164759e-06 * zs &
	403	& +( 0.791325e-08 * zt - 0.933746e-06 ) * zt &
	404	& + 0.380374e-04 ) * zh
	405
	406	zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta
	407	& - 0.301985e-05 ) * zt &
	408	& + 0.785567e-03 &
	409	& + ( 0.515032e-08 * zs &
	410	& + 0.788212e-08 * zt - 0.356603e-06 ) * zs &
	411	& +( ( 0.121551e-17 * zh &
	412	& - 0.602281e-15 * zs &
	413	& - 0.175379e-14 * zt + 0.176621e-12 ) * zh &
	414	& + 0.408195e-10 * zs &
	415	& + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt &
	416	& - 0.121555e-07 ) * zh
	417
	418	zthermal = zbeta * zalbet / ( rcp * zrhos + epsln )
[2977]	419	zhalin = zbeta * tsn(ji,jj,1,jp_sal) * rcs
[255]	420	ELSE
	421	zrhos = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) )
[1601]	422	zthermal = rn_alpha / ( rcp * zrhos + epsln )
[2977]	423	zhalin = rn_beta * tsn(ji,jj,1,jp_sal) * rcs
[255]	424	ENDIF
	425	! Radiative surface buoyancy force
	426	zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
	427	! Non radiative surface buoyancy force
[2528]	428	zBo (ji,jj) = grav * zthermal * qns(ji,jj) - grav * zhalin * ( emps(ji,jj)-rnf(ji,jj) )
[255]	429	! Surface Temperature flux for non-local term
[888]	430	wt0(ji,jj) = - ( qsr(ji,jj) + qns(ji,jj) )* ro0cpr * tmask(ji,jj,1)
[255]	431	! Surface salinity flux for non-local term
[2977]	432	ws0(ji,jj) = - ( ( emps(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) * rcs ) * tmask(ji,jj,1)
[255]	433	ENDDO
	434	ENDDO
	435
[1601]	436	zflageos = 0.5 + SIGN( 0.5, nn_eos - 1. )
[255]	437	! Compute surface buoyancy forcing, Monin Obukhov and Ekman depths
	438	!------------------------------------------------------------------
	439	DO jj = 2, jpjm1
	440	DO ji = fs_2, fs_jpim1
	441	! Reference surface density = density at first T point level
	442	zrhos = rhop(ji,jj,1) + zflageos * rau0 * ( 1. - tmask(ji,jj,1) )
	443	! Friction velocity (zustar), at T-point : LMD94 eq. 2
[1695]	444	zustar(ji,jj) = SQRT( taum(ji,jj) / ( zrhos + epsln ) )
[255]	445	ENDDO
	446	ENDDO
	447
	448	!CDIR NOVERRCHK
	449	! ! ===============
	450	DO jj = 2, jpjm1 ! Vertical slab
	451	! ! ===============
	452
	453	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	454	! II Compute Boundary layer mixing coef. and diagnose the new boundary layer depth
	455	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	456
	457	! Initialization
	458	jkmax = 0
	459	zdept (:,:) = 0.
	460	zdepw (:,:) = 0.
	461	zriblk(:,:) = 0.
	462	zmoek (:,:) = 0.
	463	zvisc (:,:) = 0.
	464	zdift (:,:) = 0.
	465	#if defined key_zdfddm
	466	zdifs (:,:) = 0.
	467	#endif
	468	zmask (:,:) = 0.
	469	DO ji = fs_2, fs_jpim1
	470	zria(ji ) = 0.
	471	! Maximum boundary layer depth
[2528]	472	ikbot = mbkt(ji,jj) ! ikbot is the last T point in the water
[255]	473	zhmax(ji) = fsdept(ji,jj,ikbot) - 0.001
	474	! Compute Monin obukhov length scale at the surface and Ekman depth:
	475	zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * ratt(1)
	476	zekman(ji) = rcekman * zustar(ji,jj) / ( ABS( ff(ji,jj) ) + epsln )
	477	zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)
	478	zmoa(ji) = zucube / ( vonk * ( zbuofdep + epsln ) )
[899]	479	#if defined key_c1d
[255]	480	! store the surface buoyancy forcing
	481	zstabl = 0.5 + SIGN( 0.5, zbuofdep )
	482	buof(ji,jj,1) = zbuofdep * tmask(ji,jj,1)
	483	! store the moning-oboukov length scale at surface
	484	zmob = zstabl * zmoa(ji) + ( 1.0 - zstabl ) * fsdept(ji,jj,1)
	485	mols(ji,jj,1) = MIN( zmob , zhmax(ji) ) * tmask(ji,jj,1)
	486	! store Ekman depth
	487	zek = zstabl * zekman(ji) + ( 1.0 - zstabl ) * fsdept(ji,jj,1)
	488	ekdp(ji,jj ) = MIN( zek , zhmax(ji) ) * tmask(ji,jj,1)
	489	#endif
	490	END DO
	491	! Compute the pipe
	492	! ---------------------
	493	DO jk = 2, jpkm1
	494	DO ji = fs_2, fs_jpim1
	495	! Compute bfsfc = Bo + radiative contribution down to hbf*depht
	496	zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * ratt(jk)
	497	! Flag (zstabl = 1) if positive forcing
	498	zstabl = 0.5 + SIGN( 0.5, zbuofdep)
	499
	500	! Compute bulk richardson number zrib at depht
	501	!-------------------------------------------------------
	502	! [Br - B(d)] * d zrinum
	503	! Rib(z) = ----------------------- = -------------
	504	! \|Vr - V(d)\|^2 + Vt(d)^2 zdVsq + zVtsq
	505	!
	506	! First compute zt,zs,zu,zv = means in the surface layer < epsilon*depht
	507	! Else surface values are taken at the first T level.
	508	! For stability, resolved vertical shear is computed with "before velocities".
	509	zref = epsilon * fsdept(ji,jj,jk)
	510	#if defined key_kppcustom
	511	! zref = gdept(1)
	512	zref = fsdept(ji,jj,1)
[2977]	513	zt = tsn(ji,jj,1,jp_tem)
	514	zs = tsn(ji,jj,1,jp_sal)
[255]	515	zrh = rhop(ji,jj,1)
	516	zu = ( ub(ji,jj,1) + ub(ji - 1,jj ,1) ) / MAX( 1. , umask(ji,jj,1) + umask(ji - 1,jj ,1) )
	517	zv = ( vb(ji,jj,1) + vb(ji ,jj - 1,1) ) / MAX( 1. , vmask(ji,jj,1) + vmask(ji ,jj - 1,1) )
	518	#else
	519	zt = 0.
	520	zs = 0.
	521	zu = 0.
	522	zv = 0.
	523	zrh = 0.
	524	! vertically integration over the upper epsilon*gdept(jk) ; del () array is computed once in zdf_kpp_init
	525	DO jm = 1, jpkm1
[2977]	526	zt = zt + del(jk,jm) * tsn(ji,jj,jm,jp_tem)
	527	zs = zs + del(jk,jm) * tsn(ji,jj,jm,jp_sal)
[255]	528	zu = zu + 0.5 * del(jk,jm) &
	529	& * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
	530	& / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
	531	zv = zv + 0.5 * del(jk,jm) &
	532	& * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
	533	& / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
	534	zrh = zrh + del(jk,jm) * rhop(ji,jj,jm)
	535	END DO
	536	#endif
[2977]	537	zsr = SQRT( ABS( tsn(ji,jj,jk,jp_sal) ) )
[255]	538	! depth
	539	zh = fsdept(ji,jj,jk)
	540	! compute compression terms on density
	541	ze = ( -3.508914e-8zt-1.248266e-8 ) zt-2.595994e-6
	542	zbw = ( 1.296821e-6zt-5.782165e-9 ) zt+1.045941e-4
	543	zb = zbw + ze * zs
	544
	545	zd = -2.042967e-2
	546	zc = (-7.267926e-5zt+2.598241e-3 ) zt+0.1571896
	547	zaw = ( ( 5.939910e-6zt+2.512549e-3 ) zt-0.1028859 ) *zt - 4.721788
	548	za = ( zdzsr + zc ) zs + zaw
	549
	550	zb1 = (-0.1909078zt+7.390729 ) zt-55.87545
	551	za1 = ( ( 2.326469e-3zt+1.553190)zt-65.00517 ) *zt+1044.077
	552	zkw = ( ( (-1.361629e-4zt-1.852732e-2 ) zt-30.41638 ) zt + 2098.925 ) zt+190925.6
	553	zk0 = ( zb1zsr + za1 )zs + zkw
	554	zcomp = 1.0 - zh / ( zk0 - zh * ( za - zh * zb ) )
	555
	556	#if defined key_kppcustom
	557	! potential density of water(zrh = zt,zs at level jk):
	558	zrhdr = zrh / zcomp
	559	#else
	560	! potential density of water(ztref,zsref at level jk):
	561	! compute volumic mass pure water at atm pressure
[1601]	562	IF ( nn_eos < 1 ) THEN
[255]	563	zr1= ( ( ( ( 6.536332e-9zt-1.120083e-6 )zt+1.001685e-4)*zt &
	564	& -9.095290e-3 )zt+6.793952e-2 )zt+999.842594
	565	! seawater volumic mass atm pressure
	566	zr2= ( ( ( 5.3875e-9zt-8.2467e-7 ) zt+7.6438e-5 ) *zt &
	567	& -4.0899e-3 ) *zt+0.824493
	568	zr3= ( -1.6546e-6zt+1.0227e-4 ) zt-5.72466e-3
	569	zr4= 4.8314e-4
	570	! potential volumic mass (reference to the surface)
	571	zrhop= ( zr4zs + zr3zsr + zr2 ) *zs + zr1
	572	zrhdr = zrhop / zcomp
	573	ELSE
	574	zrhdr = zrh / zcomp
	575	ENDIF
	576	#endif
	577
	578	! potential density of ambiant water at level jk :
	579	zrhd = ( rhd(ji,jj,jk) * rau0 + rau0 )
	580
	581	! And now the Rib number numerator .
	582	zrinum = grav * ( zrhd - zrhdr ) / rau0
	583	zrinum = zrinum * ( fsdept(ji,jj,jk) - zref ) * tmask(ji,jj,jk)
	584
	585	! Resolved shear contribution to Rib at depth T-point (zdVsq)
	586	ztx = ( ub( ji , jj ,jk) + ub(ji - 1, jj ,jk) ) &
	587	& / MAX( 1. , umask( ji , jj ,jk) + umask(ji - 1, jj ,jk) )
	588	zty = ( vb( ji , jj ,jk) + vb(ji ,jj - 1,jk) ) &
	589	& / MAX( 1., vmask( ji , jj ,jk) + vmask(ji ,jj - 1,jk) )
	590
	591	zdVsq = ( zu - ztx ) * ( zu - ztx ) + ( zv - zty ) * ( zv - zty )
	592
	593	! Scalar turbulent velocity scale zws for hbl=gdept
	594	zscale = zstabl + ( 1.0 - zstabl ) * epsilon
	595	zehat = vonk * zscale * fsdept(ji,jj,jk) * zbuofdep
	596	zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)
	597	zeta = zehat / ( zucube + epsln )
	598
	599	IF( zehat > 0. ) THEN
	600	! Stable case
	601	zws = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
	602	ELSE
	603	! Unstable case
	604	#if defined key_kpplktb
	605	! use lookup table
	606	zd = zehat - dehatmin
	607	il = INT( zd / dezehat )
	608	il = MIN( il, nilktbm1 )
	609	il = MAX( il, 1 )
	610
	611	ud = zustar(ji,jj) - ustmin
	612	jl = INT( ud / deustar )
	613	jl = MIN( jl, njlktbm1 )
	614	jl = MAX( jl, 1 )
	615
	616	zfrac = zd / dezehat - FLOAT( il )
	617	ufrac = ud / deustar - FLOAT( jl )
	618	zwas = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
	619	zwbs = ( 1. - zfrac ) * wslktb(il,jl ) + zfrac * wslktb(il+1,jl )
	620	!
	621	zws = ( 1. - ufrac ) * zwbs + ufrac * zwas
	622	#else
	623	! use analytical functions:
	624	zcons = 0.5 + SIGN( 0.5 , ( rzetas - zeta ) )
	625	zwcons = vonk * zustar(ji,jj) * ( ( ABS( rconas - rconcs * zeta ) )**pthird )
	626	zwsun = vonk * zustar(ji,jj) * SQRT( ABS ( 1.0 - rconc2 * zeta ) )
	627	!
	628	zws = zcons * zwcons + ( 1.0 - zcons) * zwsun
	629	#endif
	630	ENDIF
	631
	632	! Turbulent shear contribution to Rib (zVtsq) bv frequency at levels ( ie T-point jk)
	633	zrn2 = 0.5 * ( rn2(ji,jj,jk) + rn2(ji,jj,jk+1) )
	634	zbvzed = SQRT( ABS( zrn2 ) )
	635	zVtsq = fsdept(ji,jj,jk) * zws * zbvzed * Vtc
	636
	637	! Finally, the bulk Richardson number at depth fsdept(i,j,k)
	638	zrib = zrinum / ( zdVsq + zVtsq + epsln )
	639
	640	! Find subscripts around the boundary layer depth, build the pipe
	641	! ----------------------------------------------------------------
	642
	643	! Flag (zflagri = 1) if zrib < Ricr
	644	zflagri = 0.5 + SIGN( 0.5, ( Ricr - zrib ) )
	645	! Flag (zflagh = 1) if still within overall boundary layer
	646	zflagh = 0.5 + SIGN( 0.5, ( fsdept(ji,jj,1) - zdept(ji,2) ) )
	647
	648	! Ekman layer depth
	649	zek = zstabl * zekman(ji) + ( 1.0 - zstabl ) * zhmax(ji)
	650	zflag = 0.5 + SIGN( 0.5, ( zek - fsdept(ji,jj,jk-1) ) )
	651	zek = zflag * zek + ( 1.0 - zflag ) * zhmax(ji)
	652	zflagek = 0.5 + SIGN( 0.5, ( zek - fsdept(ji,jj,jk) ) )
	653	! Flag (zflagmo = 1) if still within stable Monin-Obukhov and in water
	654	zmob = zucube / ( vonk * ( zbuofdep + epsln ) )
	655	ztemp = zstabl * zmob + ( 1.0 - zstabl) * zhmax(ji)
	656	ztemp = MIN( ztemp , zhmax(ji) )
	657	zflagmo = 0.5 + SIGN( 0.5, ( ztemp - fsdept(ji,jj,jk) ) )
	658
	659	! No limitation by Monin Obukhov or Ekman depths:
	660	! zflagek = 1.0
	661	! zflagmo = 0.5 + SIGN( 0.5, ( zhmax(ji) - fsdept(ji,jj,jk) ) )
	662
	663	! Load pipe via zflagkb for later calculations
	664	! Flag (zflagkb = 1) if zflagh = 1 and (zflagri = 0 or zflagek = 0 or zflagmo = 0)
	665	zflagkb = zflagh * ( 1.0 - ( zflagri * zflagek * zflagmo ) )
	666
	667	zmask(ji,jk) = zflagh
	668	jkp2 = MIN( jk+2 , ikbot )
	669	jkm1 = MAX( jk-1 , 2 )
	670	jkmax = MAX( jkmax, jk * INT( REAL( zflagh+epsln ) ) )
	671
	672	zdept(ji,1) = zdept(ji,1) + zflagkb * fsdept(ji,jj,jk-1)
	673	zdept(ji,2) = zdept(ji,2) + zflagkb * fsdept(ji,jj,jk )
	674	zdept(ji,3) = zdept(ji,3) + zflagkb * fsdept(ji,jj,jk+1)
	675
	676	zdepw(ji,1) = zdepw(ji,1) + zflagkb * fsdepw(ji,jj,jk-1)
	677	zdepw(ji,2) = zdepw(ji,2) + zflagkb * fsdepw(ji,jj,jk )
	678	zdepw(ji,3) = zdepw(ji,3) + zflagkb * fsdepw(ji,jj,jk+1)
	679	zdepw(ji,4) = zdepw(ji,4) + zflagkb * fsdepw(ji,jj,jkp2)
	680
	681	zriblk(ji,1) = zriblk(ji,1) + zflagkb * zria(ji)
	682	zriblk(ji,2) = zriblk(ji,2) + zflagkb * zrib
	683
	684	zmoek (ji,0) = zmoek (ji,0) + zflagkb * zek
	685	zmoek (ji,1) = zmoek (ji,1) + zflagkb * zmoa(ji)
	686	zmoek (ji,2) = zmoek (ji,2) + zflagkb * ztemp
	687	! Save Monin Obukhov depth
	688	zmoa (ji) = zmob
	689
	690	zvisc(ji,1) = zvisc(ji,1) + zflagkb * avmu(ji,jj,jkm1)
	691	zvisc(ji,2) = zvisc(ji,2) + zflagkb * avmu(ji,jj,jk )
	692	zvisc(ji,3) = zvisc(ji,3) + zflagkb * avmu(ji,jj,jk+1)
	693	zvisc(ji,4) = zvisc(ji,4) + zflagkb * avmu(ji,jj,jkp2)
	694
	695	zdift(ji,1) = zdift(ji,1) + zflagkb * avt (ji,jj,jkm1)
	696	zdift(ji,2) = zdift(ji,2) + zflagkb * avt (ji,jj,jk )
	697	zdift(ji,3) = zdift(ji,3) + zflagkb * avt (ji,jj,jk+1)
	698	zdift(ji,4) = zdift(ji,4) + zflagkb * avt (ji,jj,jkp2)
	699
	700	#if defined key_zdfddm
	701	zdifs(ji,1) = zdifs(ji,1) + zflagkb * avs (ji,jj,jkm1)
	702	zdifs(ji,2) = zdifs(ji,2) + zflagkb * avs (ji,jj,jk )
	703	zdifs(ji,3) = zdifs(ji,3) + zflagkb * avs (ji,jj,jk+1)
	704	zdifs(ji,4) = zdifs(ji,4) + zflagkb * avs (ji,jj,jkp2)
	705	#endif
	706	! Save the Richardson number
	707	zria (ji) = zrib
[899]	708	#if defined key_c1d
[255]	709	! store buoyancy length scale
	710	buof(ji,jj,jk) = zbuofdep * tmask(ji,jj,jk)
	711	! store Monin Obukhov
	712	zmob = zstabl * zmob + ( 1.0 - zstabl) * fsdept(ji,jj,1)
	713	mols(ji,jj,jk) = MIN( zmob , zhmax(ji) ) * tmask(ji,jj,jk)
	714	! Bulk Richardson number
	715	rib(ji,jj,jk) = zrib * tmask(ji,jj,jk)
	716	#endif
	717	END DO
	718	END DO
	719	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	720	! III PROCESS THE PIPE
	721	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	722
	723	DO ji = fs_2, fs_jpim1
	724
	725	! Find the boundary layer depth zhbl
	726	! ----------------------------------------
	727
	728	! Interpolate monin Obukhov and critical Ri mumber depths
	729	ztemp = zdept(ji,2) - zdept(ji,1)
	730	zflag = ( Ricr - zriblk(ji,1) ) / ( zriblk(ji,2) - zriblk(ji,1) + epsln )
	731	zhrib = zdept(ji,1) + zflag * ztemp
	732
	733	IF( zriblk(ji,2) < Ricr ) zhrib = zhmax(ji)
	734
	735	IF( zmoek(ji,2) < zdept(ji,2) ) THEN
	736	IF ( zmoek(ji,1) < 0. ) THEN
	737	zmob = zdept(ji,2) - epsln
	738	ELSE
	739	zmob = ztemp + zmoek(ji,1) - zmoek(ji,2)
	740	zmob = ( zmoek(ji,1) * zdept(ji,2) - zmoek(ji,2) * zdept(ji,1) ) / zmob
	741	zmob = MAX( zmob , zdept(ji,1) + epsln )
	742	ENDIF
	743	ELSE
	744	zmob = zhmax(ji)
	745	ENDIF
	746	ztemp = MIN( zmob , zmoek(ji,0) )
	747
	748	! Finally, the boundary layer depth, zhbl
	749	zhbl(ji) = MAX( fsdept(ji,jj,1) + epsln, MIN( zhrib , ztemp ) )
	750
	751	! Save hkpp for further diagnostics (optional)
	752	hkpp(ji,jj) = zhbl(ji) * tmask(ji,jj,1)
	753
	754	! Correct mask if zhbl < fsdepw(ji,jj,2) for no viscosity/diffusivity enhancement at fsdepw(ji,jj,2)
	755	! zflag = 1 if zhbl(ji) > fsdepw(ji,jj,2)
	756	IF( zhbl(ji) < fsdepw(ji,jj,2) ) zmask(ji,2) = 0.
	757
	758
	759	! Velocity scales at depth zhbl
	760	! -----------------------------------
	761
	762	! Compute bouyancy forcing down to zhbl
	763	ztemp = -hbf * zhbl(ji)
	764	zatt1 = 1.0 - ( rabs * EXP( ztemp / xsi1 ) + ( 1.0 - rabs ) * EXP( ztemp / xsi2 ) )
	765	zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * zatt1
	766	zstabl = 0.5 + SIGN( 0.5 , zbuofdep )
	767
	768	zbuofdep = zbuofdep + zstabl * epsln
	769
	770	zscale = zstabl + ( 1.0 - zstabl ) * epsilon
	771	zehat = vonk * zscale * zhbl(ji) * zbuofdep
	772	zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)
	773	zeta = zehat / ( zucube + epsln )
	774
	775	IF( zehat > 0. ) THEN
	776	! Stable case
	777	zws = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
	778	zwm = zws
	779	ELSE
	780	! Unstable case
	781	#if defined key_kpplktb
	782	! use lookup table
	783	zd = zehat - dehatmin
	784	il = INT( zd / dezehat )
	785	il = MIN( il, nilktbm1 )
	786	il = MAX( il, 1 )
	787
	788	ud = zustar(ji,jj) - ustmin
	789	jl = INT( ud / deustar )
	790	jl = MIN( jl, njlktbm1 )
	791	jl = MAX( jl, 1 )
	792
	793	zfrac = zd / dezehat - FLOAT( il )
	794	ufrac = ud / deustar - FLOAT( jl )
	795	zwas = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
	796	zwbs = ( 1. - zfrac ) * wslktb(il,jl ) + zfrac * wslktb(il+1,jl )
	797	zwam = ( 1. - zfrac ) * wmlktb(il,jl+1) + zfrac * wmlktb(il+1,jl+1)
	798	zwbm = ( 1. - zfrac ) * wmlktb(il,jl ) + zfrac * wmlktb(il+1,jl )
	799	!
	800	zws = ( 1. - ufrac ) * zwbs + ufrac * zwas
	801	zwm = ( 1. - ufrac ) * zwbm + ufrac * zwam
	802	#else
	803	! use analytical functions
	804	zconm = 0.5 + SIGN( 0.5, ( rzetam - zeta) )
	805	zcons = 0.5 + SIGN( 0.5, ( rzetas - zeta) )
	806
	807	! Momentum : zeta < rzetam (zconm = 1)
	808	! Scalars : zeta < rzetas (zcons = 1)
	809	zwconm = zustar(ji,jj) * vonk * ( ( ABS( rconam - rconcm * zeta) )**pthird )
	810	zwcons = zustar(ji,jj) * vonk * ( ( ABS( rconas - rconcs * zeta) )**pthird )
	811
	812	! Momentum : rzetam <= zeta < 0 (zconm = 0)
	813	! Scalars : rzetas <= zeta < 0 (zcons = 0)
	814	zwmun = SQRT( ABS( 1.0 - rconc2 * zeta ) )
	815	zwsun = vonk * zustar(ji,jj) * zwmun
	816	zwmun = vonk * zustar(ji,jj) * SQRT(zwmun)
	817	!
	818	zwm = zconm * zwconm + ( 1.0 - zconm ) * zwmun
	819	zws = zcons * zwcons + ( 1.0 - zcons ) * zwsun
	820
	821	#endif
	822	ENDIF
	823
	824
	825	! Viscosity, diffusivity values and derivatives at h
	826	! --------------------------------------------------------
	827
	828	! check between at which interfaces is located zhbl(ji)
	829	! ztemp = 1, zdepw(ji,2) < zhbl < zdepw(ji,3)
	830	! ztemp = 0, zdepw(ji,1) < zhbl < zdepw(ji,2)
	831	ztemp = 0.5 + SIGN( 0.5, ( zhbl(ji) - zdepw(ji,2) ) )
	832	zdep21 = zdepw(ji,2) - zdepw(ji,1) + epsln
	833	zdep32 = zdepw(ji,3) - zdepw(ji,2) + epsln
	834	zdep43 = zdepw(ji,4) - zdepw(ji,3) + epsln
	835
	836	! Compute R as in LMD94, eq D5b
	837	zdelta = ( zhbl(ji) - zdepw(ji,2) ) * ztemp / zdep32 &
	838	& + ( zhbl(ji) - zdepw(ji,1) ) * ( 1.0 - ztemp ) / zdep21
	839
	840	! Compute the vertical derivative of viscosities (zdzh) at z=zhbl(ji)
	841	zdzup = ( zvisc(ji,2) - zvisc(ji,3) ) * ztemp / zdep32 &
	842	& + ( zvisc(ji,1) - zvisc(ji,2) ) * ( 1.0 - ztemp ) / zdep21
	843
	844	zdzdn = ( zvisc(ji,3) - zvisc(ji,4) ) * ztemp / zdep43 &
	845	& + ( zvisc(ji,2) - zvisc(ji,3) ) * ( 1.0 - ztemp ) / zdep32
	846
	847	! LMD94, eq D5b :
	848	zdzh = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
	849	zdzh = MAX( zdzh , 0. )
	850
	851	! Compute viscosities (zvath) at z=zhbl(ji), LMD94 eq D5a
	852	zvath = ztemp * ( zvisc(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
	853	& + ( 1.0 - ztemp ) * ( zvisc(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
	854
	855	! Compute G (zgat1) and its derivative (zdat1) at z=hbl(ji), LMD94 eq 18
	856
	857	! Vertical derivative of velocity scale divided by velocity scale squared at z=hbl(ji)
	858	! (non zero only in stable conditions)
	859	zflag = -zstabl * rconc1 * zbuofdep / ( zucube * zustar(ji,jj) + epsln )
	860
	861	! G at its derivative at z=hbl:
	862	zgat1 = zvath / ( zhbl(ji) * ( zwm + epsln ) )
	863	zdat1 = -zdzh / ( zwm + epsln ) - zflag * zvath / zhbl(ji)
	864
	865	! G coefficients, LMD94 eq 17
	866	za2m(ji) = -2.0 + 3.0 * zgat1 - zdat1
	867	za3m(ji) = 1.0 - 2.0 * zgat1 + zdat1
	868
	869
	870	! Compute the vertical derivative of temperature diffusivities (zdzh) at z=zhbl(ji)
	871	zdzup = ( zdift(ji,2) - zdift(ji,3) ) * ztemp / zdep32 &
	872	& + ( zdift(ji,1) - zdift(ji,2) ) * ( 1.0 - ztemp ) / zdep21
	873
	874	zdzdn = ( zdift(ji,3) - zdift(ji,4) ) * ztemp / zdep43 &
	875	& + ( zdift(ji,2) - zdift(ji,3) ) * ( 1.0 - ztemp ) / zdep32
	876
	877	! LMD94, eq D5b :
	878	zdzh = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
	879	zdzh = MAX( zdzh , 0. )
	880
	881
	882	! Compute diffusivities (zvath) at z=zhbl(ji), LMD94 eq D5a
	883	zvath = ztemp * ( zdift(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
	884	& + ( 1.0 - ztemp ) * ( zdift(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
	885
	886	! G at its derivative at z=hbl:
	887	zgat1 = zvath / ( zhbl(ji) * ( zws + epsln ) )
	888	zdat1 = -zdzh / ( zws + epsln ) - zflag * zvath / zhbl(ji)
	889
	890	! G coefficients, LMD94 eq 17
	891	za2t(ji) = -2.0 + 3.0 * zgat1 - zdat1
	892	za3t(ji) = 1.0 - 2.0 * zgat1 + zdat1
	893
	894	#if defined key_zdfddm
	895	! Compute the vertical derivative of salinities diffusivities (zdzh) at z=zhbl(ji)
	896	zdzup = ( zdifs(ji,2) - zdifs(ji,3) ) * ztemp / zdep32 &
	897	& + ( zdifs(ji,1) - zdifs(ji,2) ) * ( 1.0 - ztemp ) / zdep21
	898
	899	zdzdn = ( zdifs(ji,3) - zdifs(ji,4) ) * ztemp / zdep43 &
	900	& + ( zdifs(ji,2) - zdifs(ji,3) ) * ( 1.0 - ztemp ) / zdep32
	901
	902	! LMD94, eq D5b :
	903	zdzh = ( 1.0 - zdelta ) * zdzup + zdelta * zdzdn
	904	zdzh = MAX( zdzh , 0. )
	905
	906	! Compute diffusivities (zvath) at z=zhbl(ji), LMD94 eq D5a
	907	zvath = ztemp * ( zdifs(ji,3) + zdzh * ( zdepw(ji,3) - zhbl(ji) ) ) &
	908	& + ( 1.0 - ztemp ) * ( zdifs(ji,2) + zdzh * ( zdepw(ji,2) - zhbl(ji) ) )
	909
	910	! G at its derivative at z=hbl:
	911	zgat1 = zvath / ( zhbl(ji) * ( zws + epsln ) )
	912	zdat1 = -zdzh / ( zws + epsln ) - zflag * zvath / zhbl(ji)
	913
	914	! G coefficients, LMD94 eq 17
	915	za2s(ji) = -2.0 + 3.0 * zgat1 - zdat1
	916	za3s(ji) = 1.0 - 2.0 * zgat1 + zdat1
	917	#endif
	918
	919	!-------------------turn off interior matching here------
	920	! za2(ji,1) = -2.0
	921	! za3(ji,1) = 1.0
	922	! za2(ji,2) = -2.0
	923	! za3(ji,2) = 1.0
	924	!--------------------------------------------------------
	925
	926	! Compute Enhanced Mixing Coefficients (LMD94,eq D6)
	927	! ---------------------------------------------------------------
	928
	929	! Delta
	930	zdelta = ( zhbl(ji) - zdept(ji,1) ) / ( zdept(ji,2) - zdept(ji,1) + epsln )
	931	zdelta2 = zdelta * zdelta
	932
	933	! Mixing coefficients at first level above h (zdept(ji,1))
	934	! and at first interface in the pipe (zdepw(ji,2))
	935
	936	! At first T level above h (zdept(ji,1)) (always in the boundary layer)
	937	zsig = zdept(ji,1) / zhbl(ji)
	938	ztemp = zstabl * zsig + ( 1.0 - zstabl ) * MIN( zsig , epsilon )
	939	zehat = vonk * ztemp * zhbl(ji) * zbuofdep
	940	zeta = zehat / ( zucube + epsln)
	941	zwst = vonk * zustar(ji,jj) / ( ABS( 1.0 + rconc1 * zeta ) + epsln)
	942	zwm = zstabl * zwst + ( 1.0 - zstabl ) * zwm
	943	zws = zstabl * zwst + ( 1.0 - zstabl ) * zws
	944
	945	zkm1m = zhbl(ji) * zwm * zsig * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
	946	zkm1t = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
	947	#if defined key_zdfddm
	948	zkm1s = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
	949	#endif
	950	! At first W level in the pipe (zdepw(ji,2)) (not always in the boundary layer ):
	951	zsig = MIN( zdepw(ji,2) / zhbl(ji) , 1.0 )
	952	ztemp = zstabl * zsig + ( 1.0 - zstabl ) * MIN( zsig , epsilon )
	953	zehat = vonk * ztemp * zhbl(ji) * zbuofdep
	954	zeta = zehat / ( zucube + epsln )
	955	zwst = vonk * zustar(ji,jj) / ( ABS( 1.0 + rconc1 * zeta ) + epsln)
	956	zws = zstabl * zws + ( 1.0 - zstabl ) * zws
	957	zwm = zstabl * zws + ( 1.0 - zstabl ) * zwm
	958
	959	zkmpm(ji) = zhbl(ji) * zwm * zsig * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
	960	zkmpt(ji) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
	961	#if defined key_zdfddm
	962	zkmps(ji) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
	963	#endif
	964
	965	! check if this point is in the boundary layer,else take interior viscosity/diffusivity:
	966	zflag = 0.5 + SIGN( 0.5, ( zhbl(ji) - zdepw(ji,2) ) )
	967	zkmpm(ji) = zkmpm(ji) * zflag + ( 1.0 - zflag ) * zvisc(ji,2)
	968	zkmpt(ji) = zkmpt(ji) * zflag + ( 1.0 - zflag ) * zdift(ji,2)
	969	#if defined key_zdfddm
	970	zkmps(ji) = zkmps(ji) * zflag + ( 1.0 - zflag ) * zdifs(ji,2)
	971	#endif
	972
	973	! Enhanced viscosity/diffusivity at zdepw(ji,2)
	974	ztemp = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1m + zdelta2 * zkmpm(ji)
	975	zkmpm(ji) = ( 1.0 - zdelta ) * zvisc(ji,2) + zdelta * ztemp
	976	ztemp = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1t + zdelta2 * zkmpt(ji)
	977	zkmpt(ji) = ( 1.0 - zdelta ) * zdift(ji,2) + zdelta * ztemp
	978	#if defined key_zdfddm
	979	ztemp = ( 1.0 - 2.0 * zdelta + zdelta2 ) * zkm1s + zdelta2 * zkmps(ji)
	980	zkmps(ji) = ( 1.0 - zdelta ) * zdifs(ji,2) + zdelta * ztemp
	981	#endif
	982
	983	END DO
	984	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	985	! IV. Compute vertical eddy viscosity and diffusivity coefficients
	986	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	987
	988	DO jk = 2, jkmax
	989
	990	! Compute turbulent velocity scales on the interfaces
	991	! --------------------------------------------------------
	992	DO ji = fs_2, fs_jpim1
	993	zbuofdep = zBo(ji,jj) + zBosol(ji,jj) * zatt1
	994	zstabl = 0.5 + SIGN( 0.5 , zbuofdep )
	995	zbuofdep = zbuofdep + zstabl * epsln
	996	zsig = fsdepw(ji,jj,jk) / zhbl(ji)
	997	ztemp = zstabl * zsig + ( 1. - zstabl ) * MIN( zsig , epsilon )
	998	zehat = vonk * ztemp * zhbl(ji) * zbuofdep
	999	zucube = zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj)
	1000	zeta = zehat / ( zucube + epsln )
	1001
	1002	IF( zehat > 0. ) THEN
	1003	! Stable case
	1004	zws = vonk * zustar(ji,jj) / ( 1.0 + rconc1 * zeta )
	1005	zwm = zws
	1006	ELSE
	1007	! Unstable case
	1008	#if defined key_kpplktb
	1009	! use lookup table
	1010	zd = zehat - dehatmin
	1011	il = INT( zd / dezehat )
	1012	il = MIN( il, nilktbm1 )
	1013	il = MAX( il, 1 )
	1014
	1015	ud = zustar(ji,jj) - ustmin
	1016	jl = INT( ud / deustar )
	1017	jl = MIN( jl, njlktbm1 )
	1018	jl = MAX( jl, 1 )
	1019
	1020	zfrac = zd / dezehat - FLOAT( il )
	1021	ufrac = ud / deustar - FLOAT( jl )
	1022	zwas = ( 1. - zfrac ) * wslktb(il,jl+1) + zfrac * wslktb(il+1,jl+1)
	1023	zwbs = ( 1. - zfrac ) * wslktb(il,jl ) + zfrac * wslktb(il+1,jl )
	1024	zwam = ( 1. - zfrac ) * wmlktb(il,jl+1) + zfrac * wmlktb(il+1,jl+1)
	1025	zwbm = ( 1. - zfrac ) * wmlktb(il,jl ) + zfrac * wmlktb(il+1,jl )
	1026	!
	1027	zws = ( 1. - ufrac ) * zwbs + ufrac * zwas
	1028	zwm = ( 1. - ufrac ) * zwbm + ufrac * zwam
	1029	#else
	1030	! use analytical functions
	1031	zconm = 0.5 + SIGN( 0.5, ( rzetam - zeta) )
	1032	zcons = 0.5 + SIGN( 0.5, ( rzetas - zeta) )
	1033
	1034	! Momentum : zeta < rzetam (zconm = 1)
	1035	! Scalars : zeta < rzetas (zcons = 1)
	1036	zwconm = zustar(ji,jj) * vonk * ( ( ABS( rconam - rconcm * zeta) )**pthird )
	1037	zwcons = zustar(ji,jj) * vonk * ( ( ABS( rconas - rconcs * zeta) )**pthird )
	1038
	1039	! Momentum : rzetam <= zeta < 0 (zconm = 0)
	1040	! Scalars : rzetas <= zeta < 0 (zcons = 0)
	1041	zwmun = SQRT( ABS( 1.0 - rconc2 * zeta ) )
	1042	zwsun = vonk * zustar(ji,jj) * zwmun
	1043	zwmun = vonk * zustar(ji,jj) * SQRT(zwmun)
	1044	!
	1045	zwm = zconm * zwconm + ( 1.0 - zconm ) * zwmun
	1046	zws = zcons * zwcons + ( 1.0 - zcons ) * zwsun
	1047
	1048	#endif
	1049	ENDIF
	1050
	1051	zblcm(ji,jk) = zhbl(ji) * zwm * zsig * ( 1.0 + zsig * ( za2m(ji) + zsig * za3m(ji) ) )
	1052	zblct(ji,jk) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2t(ji) + zsig * za3t(ji) ) )
	1053	#if defined key_zdfddm
	1054	zblcs(ji,jk) = zhbl(ji) * zws * zsig * ( 1.0 + zsig * ( za2s(ji) + zsig * za3s(ji) ) )
	1055	#endif
	1056	! Compute Nonlocal transport term = ghats * <ws>o
	1057	! ----------------------------------------------------
	1058	ghats(ji,jj,jk-1) = ( 1. - zstabl ) * rcg / ( zws * zhbl(ji) + epsln ) * tmask(ji,jj,jk)
	1059
	1060	END DO
	1061	END DO
	1062	! Combine interior and boundary layer coefficients and nonlocal term
	1063	! -----------------------------------------------------------------------
	1064	DO jk = 2, jpkm1
	1065	DO ji = fs_2, fs_jpim1
	1066	zflag = zmask(ji,jk) * zmask(ji,jk+1)
	1067	zviscos(ji,jj,jk) = ( 1.0 - zmask(ji,jk) ) * avmu (ji,jj,jk) & ! interior viscosities
	1068	& + zflag * zblcm(ji,jk ) & ! boundary layer viscosities
	1069	& + zmask(ji,jk) * ( 1.0 - zflag ) * zkmpm(ji ) ! viscosity enhancement at W_level near zhbl
	1070
	1071	zviscos(ji,jj,jk) = zviscos(ji,jj,jk) * tmask(ji,jj,jk)
	1072
	1073
	1074	zdiffut(ji,jj,jk) = ( 1.0 - zmask(ji,jk) ) * avt (ji,jj,jk) & ! interior diffusivities
	1075	& + zflag * zblct(ji,jk ) & ! boundary layer diffusivities
	1076	& + zmask(ji,jk) * ( 1.0 - zflag ) * zkmpt(ji ) ! diffusivity enhancement at W_level near zhbl
	1077
	1078	zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) * tmask(ji,jj,jk)
	1079	#if defined key_zdfddm
	1080	zdiffus(ji,jj,jk) = ( 1.0 - zmask(ji,jk) ) * avs (ji,jj,jk) & ! interior diffusivities
	1081	& + zflag * zblcs(ji,jk ) & ! boundary layer diffusivities
	1082	& + zmask(ji,jk) * ( 1.0 - zflag ) * zkmps(ji ) ! diffusivity enhancement at W_level near zhbl
	1083
	1084	zdiffus(ji,jj,jk) = zdiffus(ji,jj,jk) * tmask(ji,jj,jk)
	1085	#endif
	1086	! Non local flux in the boundary layer only
	1087	ghats(ji,jj,jk-1) = zmask(ji,jk) * ghats(ji,jj,jk-1)
	1088
	1089	ENDDO
	1090	END DO
	1091	! ! ===============
	1092	END DO ! End of slab
	1093	! ! ===============
	1094
	1095	! Lateral boundary conditions on zvicos and zdiffus (sign unchanged)
	1096	CALL lbc_lnk( zviscos(:,:,:), 'U', 1. ) ; CALL lbc_lnk( zdiffut(:,:,:), 'W', 1. )
	1097	#if defined key_zdfddm
	1098	CALL lbc_lnk( zdiffus(:,:,:), 'W', 1. )
	1099	#endif
	1100
[1537]	1101	SELECT CASE ( nn_ave )
[255]	1102	!
	1103	CASE ( 0 ) ! no viscosity and diffusivity smoothing
	1104
	1105	DO jk = 2, jpkm1
	1106	DO jj = 2, jpjm1
	1107	DO ji = fs_2, fs_jpim1
	1108	avmu(ji,jj,jk) = ( zviscos(ji,jj,jk) + zviscos(ji+1,jj,jk) ) &
	1109	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
	1110
	1111	avmv(ji,jj,jk) = ( zviscos(ji,jj,jk) + zviscos(ji,jj+1,jk) ) &
	1112	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
	1113
	1114	avt (ji,jj,jk) = zdiffut(ji,jj,jk) * tmask(ji,jj,jk)
	1115	#if defined key_zdfddm
	1116	avs (ji,jj,jk) = zdiffus(ji,jj,jk) * tmask(ji,jj,jk)
	1117	#endif
	1118	END DO
	1119	END DO
	1120	END DO
	1121
	1122	CASE ( 1 ) ! viscosity and diffusivity smoothing
	1123	!
	1124	! ( 1/2 1 1/2 ) ( 1/2 1/2 ) ( 1/2 1 1/2 )
	1125	! avt = 1/8 ( 1 2 1 ) avmu = 1/4 ( 1 1 ) avmv= 1/4 ( 1/2 1 1/2 )
	1126	! ( 1/2 1 1/2 ) ( 1/2 1/2 )
	1127
	1128	DO jk = 2, jpkm1
	1129	DO jj = 2, jpjm1
	1130	DO ji = fs_2, fs_jpim1
	1131
	1132	avmu(ji,jj,jk) = ( zviscos(ji ,jj ,jk) + zviscos(ji+1,jj ,jk) &
	1133	& +.5*( zviscos(ji ,jj-1,jk) + zviscos(ji+1,jj-1,jk) &
	1134	& +zviscos(ji ,jj+1,jk) + zviscos(ji+1,jj+1,jk) ) ) * eumean(ji,jj,jk)
	1135
	1136	avmv(ji,jj,jk) = ( zviscos(ji ,jj ,jk) + zviscos(ji ,jj+1,jk) &
	1137	& +.5*( zviscos(ji-1,jj ,jk) + zviscos(ji-1,jj+1,jk) &
	1138	& +zviscos(ji+1,jj ,jk) + zviscos(ji+1,jj+1,jk) ) ) * evmean(ji,jj,jk)
	1139
	1140	avt (ji,jj,jk) = ( .5*( zdiffut(ji-1,jj+1,jk) + zdiffut(ji-1,jj-1,jk) &
	1141	& +zdiffut(ji+1,jj+1,jk) + zdiffut(ji+1,jj-1,jk) ) &
	1142	& +1.*( zdiffut(ji-1,jj ,jk) + zdiffut(ji ,jj+1,jk) &
	1143	& +zdiffut(ji ,jj-1,jk) + zdiffut(ji+1,jj ,jk) ) &
	1144	& +2.* zdiffut(ji ,jj ,jk) ) * etmean(ji,jj,jk)
	1145	#if defined key_zdfddm
	1146	avs (ji,jj,jk) = ( .5*( zdiffus(ji-1,jj+1,jk) + zdiffus(ji-1,jj-1,jk) &
	1147	& +zdiffus(ji+1,jj+1,jk) + zdiffus(ji+1,jj-1,jk) ) &
	1148	& +1.*( zdiffus(ji-1,jj ,jk) + zdiffus(ji ,jj+1,jk) &
	1149	& +zdiffus(ji ,jj-1,jk) + zdiffus(ji+1,jj ,jk) ) &
	1150	& +2.* zdiffus(ji ,jj ,jk) ) * etmean(ji,jj,jk)
	1151	#endif
	1152	END DO
	1153	END DO
	1154	END DO
	1155
	1156	END SELECT
	1157
	1158	DO jk = 2, jpkm1 ! vertical slab
	1159	!
	1160	! Minimum value on the eddy diffusivity
	1161	! ----------------------------------------
	1162	DO jj = 2, jpjm1
	1163	DO ji = fs_2, fs_jpim1 ! vector opt.
	1164	avt(ji,jj,jk) = MAX( avt(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
	1165	#if defined key_zdfddm
	1166	avs(ji,jj,jk) = MAX( avs(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
	1167	#endif
	1168	END DO
	1169	END DO
	1170
	1171	!
	1172	! Minimum value on the eddy viscosity
	1173	! ----------------------------------------
	1174	DO jj = 1, jpj
	1175	DO ji = 1, jpi
	1176	avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), avmb(jk) ) * umask(ji,jj,jk)
	1177	avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), avmb(jk) ) * vmask(ji,jj,jk)
	1178	END DO
	1179	END DO
	1180	!
	1181	END DO
	1182
	1183	! Lateral boundary conditions on avt (sign unchanged)
	1184	CALL lbc_lnk( hkpp(:,:), 'T', 1. )
	1185
	1186	! Lateral boundary conditions on avt (sign unchanged)
	1187	CALL lbc_lnk( avt(:,:,:), 'W', 1. )
	1188	#if defined key_zdfddm
	1189	CALL lbc_lnk( avs(:,:,:), 'W', 1. )
	1190	#endif
	1191	! Lateral boundary conditions (avmu,avmv) (U- and V- points, sign unchanged)
	1192	CALL lbc_lnk( avmu(:,:,:), 'U', 1. ) ; CALL lbc_lnk( avmv(:,:,:), 'V', 1. )
	1193
[258]	1194	IF(ln_ctl) THEN
	1195	#if defined key_zdfddm
[516]	1196	CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', tab3d_2=avs , clinfo2=' s: ', ovlap=1, kdim=jpk)
[258]	1197	#else
[516]	1198	CALL prt_ctl(tab3d_1=avt , clinfo1=' kpp - t: ', ovlap=1, kdim=jpk)
[258]	1199	#endif
[516]	1200	CALL prt_ctl(tab3d_1=avmu, clinfo1=' kpp - u: ', mask1=umask, &
	1201	& tab3d_2=avmv, clinfo2= ' v: ', mask2=vmask, ovlap=1, kdim=jpk)
[258]	1202	ENDIF
	1203
[3155]	1204	CALL wrk_dealloc( jpi, zmoa, zekman, zhmax, zria, zhbl )
	1205	CALL wrk_dealloc( jpi, za2m, za3m, zkmpm, za2t, za3t, zkmpt )
	1206	CALL wrk_dealloc( jpi,2, zriblk )
	1207	CALL wrk_dealloc( jpi,3, zmoek, kjstart = 0 )
	1208	CALL wrk_dealloc( jpi,3, zdept )
	1209	CALL wrk_dealloc( jpi,4, zdepw, zdift, zvisc )
	1210	CALL wrk_dealloc( jpi,jpj, zBo, zBosol, zustar )
	1211	CALL wrk_dealloc( jpi,jpk, zmask, zblcm, zblct, zblcs )
	1212	#if defined key_zdfddm
	1213	CALL wrk_dealloc( jpi,4, zdifs )
	1214	CALL wrk_dealloc( jpi, zmoa, za2s, za3s, zkmps )
	1215	CALL wrk_dealloc( jpi,jpk, zblcs )
	1216	CALL wrk_dealloc( jpi,jpi,jpk, zdiffus )
	1217	#endif
[2715]	1218	!
[3229]	1219	IF( nn_timing == 1 ) CALL timing_stop('zdf_kpp')
	1220	!
[255]	1221	END SUBROUTINE zdf_kpp
	1222
	1223
[896]	1224	SUBROUTINE tra_kpp( kt )
[463]	1225	!!----------------------------------------------------------------------
	1226	!! * ROUTINE tra_kpp *
	1227	!!
[2528]	1228	!! ** Purpose : compute and add to the tracer trend the non-local tracer flux
[463]	1229	!!
	1230	!! ** Method : ???
	1231	!!----------------------------------------------------------------------
[2528]	1232	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace
[463]	1233	!!----------------------------------------------------------------------
[896]	1234	INTEGER, INTENT(in) :: kt
	1235	INTEGER :: ji, jj, jk
[3229]	1236	!
	1237	IF( nn_timing == 1 ) CALL timing_start('tra_kpp')
	1238	!
[463]	1239	IF( kt == nit000 ) THEN
[2528]	1240	IF(lwp) WRITE(numout,*)
[463]	1241	IF(lwp) WRITE(numout,*) 'tra_kpp : KPP non-local tracer fluxes'
	1242	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	1243	ENDIF
	1244
[2528]	1245	IF( l_trdtra ) THEN !* Save ta and sa trends
	1246	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
	1247	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)
[463]	1248	ENDIF
	1249
	1250	! add non-local temperature and salinity flux ( in convective case only)
	1251	DO jk = 1, jpkm1
[2528]	1252	DO jj = 2, jpjm1
[463]	1253	DO ji = fs_2, fs_jpim1
[2528]	1254	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &
	1255	& - ( ghats(ji,jj,jk ) * avt (ji,jj,jk ) &
	1256	& - ghats(ji,jj,jk+1) * avt (ji,jj,jk+1) ) * wt0(ji,jj) / fse3t(ji,jj,jk)
	1257	tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) &
	1258	& - ( ghats(ji,jj,jk ) * fsavs(ji,jj,jk ) &
	1259	& - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * ws0(ji,jj) / fse3t(ji,jj,jk)
[463]	1260	END DO
	1261	END DO
	1262	END DO
	1263
	1264	! save the non-local tracer flux trends for diagnostic
	1265	IF( l_trdtra ) THEN
[2528]	1266	ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
	1267	ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
[463]	1268	!!bug gm jpttdzdf ==> jpttkpp
[2528]	1269	CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_zdf, ztrdt )
	1270	CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_zdf, ztrds )
	1271	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
[463]	1272	ENDIF
	1273
[2528]	1274	IF(ln_ctl) THEN
	1275	CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' kpp - Ta: ', mask1=tmask, &
	1276	& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
[463]	1277	ENDIF
[3229]	1278	!
	1279	IF( nn_timing == 1 ) CALL timing_stop('tra_kpp')
	1280	!
[463]	1281	END SUBROUTINE tra_kpp
	1282
[2528]	1283	#if defined key_top
	1284	!!----------------------------------------------------------------------
	1285	!! 'key_top' TOP models
	1286	!!----------------------------------------------------------------------
	1287	SUBROUTINE trc_kpp( kt )
	1288	!!----------------------------------------------------------------------
	1289	!! * ROUTINE trc_kpp *
	1290	!!
	1291	!! ** Purpose : compute and add to the tracer trend the non-local
	1292	!! tracer flux
	1293	!!
	1294	!! ** Method : ???
	1295	!!
	1296	!! history :
	1297	!! 9.0 ! 2005-11 (G. Madec) Original code
	1298	!! NEMO 3.3 ! 2010-06 (C. Ethe ) Adapted to passive tracers
	1299	!!----------------------------------------------------------------------
	1300	USE trc
	1301	USE prtctl_trc ! Print control
[2715]	1302	!
	1303	INTEGER, INTENT(in) :: kt ! ocean time-step index
	1304	!
[2528]	1305	INTEGER :: ji, jj, jk, jn ! Dummy loop indices
	1306	REAL(wp) :: ztra, zflx
	1307	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrtrd
	1308	!!----------------------------------------------------------------------
[463]	1309
[2528]	1310	IF( kt == nit000 ) THEN
	1311	IF(lwp) WRITE(numout,*)
	1312	IF(lwp) WRITE(numout,*) 'trc_kpp : KPP non-local tracer fluxes'
	1313	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	1314	ENDIF
	1315
	1316	IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) )
	1317	!
	1318	DO jn = 1, jptra
	1319	!
	1320	IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn)
	1321	! add non-local on passive tracer flux ( in convective case only)
	1322	DO jk = 1, jpkm1
	1323	DO jj = 2, jpjm1
	1324	DO ji = fs_2, fs_jpim1
	1325	! Surface tracer flux for non-local term
	1326	zflx = - ( emps(ji,jj) * tra(ji,jj,1,jn) * rcs ) * tmask(ji,jj,1)
	1327	! compute the trend
	1328	ztra = - ( ghats(ji,jj,jk ) * fsavs(ji,jj,jk ) &
	1329	& - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * zflx / fse3t(ji,jj,jk)
	1330	! add the trend to the general trend
	1331	tra(ji,jj,jk,jn) = tra(ji,jj,jk,jn) + ztra
	1332	END DO
	1333	END DO
	1334	END DO
	1335	! save the non-local tracer flux trends for diagnostic
	1336	IF( l_trdtrc ) ztrtrd(:,:,:) = tra(:,:,:,jn) - ztrtrd(:,:,:)
	1337	CALL trd_tra( kt, 'TRC', jn, jptra_trd_zdf, ztrtrd(:,:,:,jn) )
	1338	!
	1339	END DO
	1340	IF( l_trdtrc ) DEALLOCATE( ztrtrd )
	1341	IF( ln_ctl ) THEN
	1342	WRITE(charout, FMT="(' kpp')") ; CALL prt_ctl_trc_info(charout)
	1343	CALL prt_ctl_trc( tab4d=tra, mask=tmask, clinfo=clname, clinfo2='trd' )
	1344	ENDIF
	1345	!
	1346	END SUBROUTINE trc_kpp
[2715]	1347
	1348	#else
	1349	!!----------------------------------------------------------------------
	1350	!! NO 'key_top' DUMMY routine No TOP models
	1351	!!----------------------------------------------------------------------
	1352	SUBROUTINE trc_kpp( kt ) ! Dummy routine
	1353	INTEGER, INTENT(in) :: kt ! ocean time-step index
	1354	WRITE(,) 'tra_kpp: You should not have seen this print! error?', kt
	1355	END SUBROUTINE trc_kpp
[2528]	1356	#endif
	1357
[255]	1358	SUBROUTINE zdf_kpp_init
	1359	!!----------------------------------------------------------------------
	1360	!! * ROUTINE zdf_kpp_init *
	1361	!!
	1362	!! ** Purpose : Initialization of the vertical eddy diffivity and
	1363	!! viscosity when using a kpp turbulent closure scheme
	1364	!!
	1365	!! ** Method : Read the namkpp namelist and check the parameters
	1366	!! called at the first timestep (nit000)
	1367	!!
	1368	!! ** input : Namlist namkpp
	1369	!!----------------------------------------------------------------------
[2528]	1370	INTEGER :: ji, jj, jk ! dummy loop indices
[255]	1371	#if ! defined key_kppcustom
[2528]	1372	INTEGER :: jm ! dummy loop indices
	1373	REAL(wp) :: zref, zdist ! tempory scalars
[255]	1374	#endif
	1375	#if defined key_kpplktb
[2528]	1376	REAL(wp) :: zustar, zucube, zustvk, zeta, zehat ! tempory scalars
[255]	1377	#endif
[2528]	1378	REAL(wp) :: zhbf ! tempory scalars
	1379	LOGICAL :: ll_kppcustom ! 1st ocean level taken as surface layer
	1380	LOGICAL :: ll_kpplktb ! Lookup table for turbul. velocity scales
[1537]	1381	!!
[1601]	1382	NAMELIST/namzdf_kpp/ ln_kpprimix, rn_difmiw, rn_difsiw, rn_riinfty, rn_difri, rn_bvsqcon, rn_difcon, nn_ave
[255]	1383	!!----------------------------------------------------------------------
[3229]	1384	!
	1385	IF( nn_timing == 1 ) CALL timing_start('zdf_kpp_init')
	1386	!
[1537]	1387	REWIND ( numnam ) ! Read Namelist namkpp : K-Profile Parameterisation
[1601]	1388	READ ( numnam, namzdf_kpp )
[255]	1389
[1537]	1390	IF(lwp) THEN ! Control print
[255]	1391	WRITE(numout,*)
[1537]	1392	WRITE(numout,*) 'zdf_kpp_init : K-Profile Parameterisation'
[255]	1393	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	1394	WRITE(numout,*) ' Namelist namzdf_kpp : set tke mixing parameters'
[1537]	1395	WRITE(numout,*) ' Shear instability mixing ln_kpprimix = ', ln_kpprimix
	1396	WRITE(numout,*) ' max. internal wave viscosity rn_difmiw = ', rn_difmiw
	1397	WRITE(numout,*) ' max. internal wave diffusivity rn_difsiw = ', rn_difsiw
	1398	WRITE(numout,*) ' Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
	1399	WRITE(numout,*) ' max. shear mixing at Rig = 0 rn_difri = ', rn_difri
	1400	WRITE(numout,*) ' Brunt-Vaisala squared for max. convection rn_bvsqcon = ', rn_bvsqcon
	1401	WRITE(numout,*) ' max. mix. in interior convec. rn_difcon = ', rn_difcon
	1402	WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave
[255]	1403	ENDIF
	1404
[2715]	1405	! ! allocate zdfkpp arrays
	1406	IF( zdf_kpp_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_kpp_init : unable to allocate arrays' )
	1407
[255]	1408	ll_kppcustom = .FALSE.
	1409	ll_kpplktb = .FALSE.
	1410
	1411	#if defined key_kppcustom
	1412	ll_kppcustom = .TRUE.
	1413	#endif
	1414	#if defined key_kpplktb
	1415	ll_kpplktb = .TRUE.
	1416	#endif
	1417	IF(lwp) THEN
	1418	WRITE(numout,*) ' Lookup table for turbul. velocity scales ll_kpplktb = ', ll_kpplktb
	1419	WRITE(numout,*) ' 1st ocean level taken as surface layer ll_kppcustom = ', ll_kppcustom
	1420	ENDIF
	1421
	1422	IF( lk_zdfddm) THEN
	1423	IF(lwp) THEN
	1424	WRITE(numout,*)
	1425	WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
	1426	WRITE(numout,*) ' CAUTION : done in routine zdfkpp, not in routine zdfddm '
	1427	ENDIF
	1428	ENDIF
	1429
	1430
	1431	!set constants not in namelist
	1432	!-----------------------------
	1433	Vtc = rconcv * SQRT( 0.2 / ( rconcs * epsilon ) ) / ( vonk * vonk * Ricr )
	1434	rcg = rcstar * vonk * ( rconcs * vonk * epsilon )**pthird
	1435
	1436	IF(lwp) THEN
[1537]	1437	WRITE(numout,*)
[255]	1438	WRITE(numout,*) ' Constant value for unreso. turbul. velocity shear Vtc = ', Vtc
	1439	WRITE(numout,*) ' Non-dimensional coef. for nonlocal transport rcg = ', rcg
	1440	ENDIF
	1441
	1442	! ratt is the attenuation coefficient for solar flux
	1443	! Should be different is s_coordinate
	1444	DO jk = 1, jpk
	1445	zhbf = - fsdept(1,1,jk) * hbf
	1446	ratt(jk) = 1.0 - ( rabs * EXP( zhbf / xsi1 ) + ( 1.0 - rabs ) * EXP( zhbf / xsi2 ) )
	1447	ENDDO
	1448
	1449	! Horizontal average : initialization of weighting arrays
	1450	! -------------------
	1451
[1537]	1452	SELECT CASE ( nn_ave )
[255]	1453
	1454	CASE ( 0 ) ! no horizontal average
	1455	IF(lwp) WRITE(numout,*) ' no horizontal average on avt, avmu, avmv'
	1456	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
	1457	! weighting mean arrays etmean, eumean and evmean
	1458	! ( 1 1 ) ( 1 )
	1459	! avt = 1/4 ( 1 1 ) avmu = 1/2 ( 1 1 ) avmv= 1/2 ( 1 )
	1460	!
	1461	etmean(:,:,:) = 0.e0
	1462	eumean(:,:,:) = 0.e0
	1463	evmean(:,:,:) = 0.e0
	1464
	1465	DO jk = 1, jpkm1
	1466	DO jj = 2, jpjm1
	1467	DO ji = 2, jpim1 ! vector opt.
	1468	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
	1469	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
	1470	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
	1471
	1472	eumean(ji,jj,jk) = umask(ji,jj,jk) &
	1473	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) )
	1474
	1475	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
	1476	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) )
	1477	END DO
	1478	END DO
	1479	END DO
	1480
	1481	CASE ( 1 ) ! horizontal average
	1482	IF(lwp) WRITE(numout,*) ' horizontal average on avt, avmu, avmv'
	1483	! weighting mean arrays etmean, eumean and evmean
	1484	! ( 1/2 1 1/2 ) ( 1/2 1/2 ) ( 1/2 1 1/2 )
	1485	! avt = 1/8 ( 1 2 1 ) avmu = 1/4 ( 1 1 ) avmv= 1/4 ( 1/2 1 1/2 )
	1486	! ( 1/2 1 1/2 ) ( 1/2 1/2 )
	1487	etmean(:,:,:) = 0.e0
	1488	eumean(:,:,:) = 0.e0
	1489	evmean(:,:,:) = 0.e0
	1490
	1491	DO jk = 1, jpkm1
	1492	DO jj = 2, jpjm1
	1493	DO ji = fs_2, fs_jpim1 ! vector opt.
	1494	etmean(ji,jj,jk) = tmask(ji, jj,jk) &
	1495	& / MAX( 1., 2.* tmask(ji,jj,jk) &
	1496	& +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) &
	1497	& +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
	1498	& +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) &
	1499	& +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) )
	1500
	1501	eumean(ji,jj,jk) = umask(ji,jj,jk) &
	1502	& / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) &
	1503	& +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) &
	1504	& +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) )
	1505
	1506	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
	1507	& / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
	1508	& +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) &
	1509	& +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) )
	1510	END DO
	1511	END DO
	1512	END DO
	1513
	1514	CASE DEFAULT
[1537]	1515	WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave
[896]	1516	CALL ctl_stop( ctmp1 )
[255]	1517
	1518	END SELECT
	1519
	1520	! Initialization of vertical eddy coef. to the background value
	1521	! -------------------------------------------------------------
	1522	DO jk = 1, jpk
	1523	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
	1524	avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)
	1525	avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)
	1526	END DO
	1527
	1528	! zero the surface flux for non local term and kpp mixed layer depth
	1529	! ------------------------------------------------------------------
	1530	ghats(:,:,:) = 0.
	1531	wt0 (:,: ) = 0.
	1532	ws0 (:,: ) = 0.
	1533	hkpp (:,: ) = 0. ! just a diagnostic (not essential)
	1534
	1535	#if ! defined key_kppcustom
	1536	! compute arrays (del, wz) for reference mean values
	1537	! (increase speed for vectorization key_kppcustom not defined)
	1538	del(1:jpk, 1:jpk) = 0.
	1539	DO jk = 1, jpk
	1540	zref = epsilon * fsdept(1,1,jk)
	1541	DO jm = 1 , jpk
	1542	zdist = zref - fsdepw(1,1,jm)
	1543	IF( zdist > 0. ) THEN
	1544	del(jk,jm) = MIN( zdist, fse3t(1,1,jm) ) / zref
	1545	ELSE
	1546	del(jk,jm) = 0.
	1547	ENDIF
	1548	ENDDO
	1549	ENDDO
	1550	#endif
	1551
	1552	#if defined key_kpplktb
	1553	! build lookup table for turbulent velocity scales
	1554	dezehat = ( dehatmax - dehatmin ) / nilktbm1
	1555	deustar = ( ustmax - ustmin ) / njlktbm1
	1556
	1557	DO jj = 1, njlktb
	1558	zustar = ( jj - 1) * deustar + ustmin
	1559	zustvk = vonk * zustar
	1560	zucube = zustar * zustar * zustar
	1561	DO ji = 1 , nilktb
	1562	zehat = ( ji - 1 ) * dezehat + dehatmin
	1563	zeta = zehat / ( zucube + epsln )
	1564	IF( zehat >= 0 ) THEN ! Stable case
	1565	wmlktb(ji,jj) = zustvk / ABS( 1.0 + rconc1 * zeta + epsln )
	1566	wslktb(ji,jj) = wmlktb(ji,jj)
	1567	ELSE ! Unstable case
	1568	IF( zeta > rzetam ) THEN
	1569	wmlktb(ji,jj) = zustvk * ABS( 1.0 - rconc2 * zeta )**pfourth
	1570	ELSE
	1571	wmlktb(ji,jj) = zustvk * ABS( rconam - rconcm * zeta )**pthird
	1572	ENDIF
	1573
	1574	IF( zeta > rzetas ) THEN
	1575	wslktb(ji,jj) = zustvk * SQRT( ABS( 1.0 - rconc2 * zeta ) )
	1576	ELSE
	1577	wslktb(ji,jj) = zustvk * ABS( rconas - rconcs * zeta )**pthird
	1578	ENDIF
	1579	ENDIF
	1580	END DO
	1581	END DO
	1582	#endif
[2715]	1583	!
[3229]	1584	IF( nn_timing == 1 ) CALL timing_stop('zdf_kpp_init')
	1585	!
[255]	1586	END SUBROUTINE zdf_kpp_init
	1587
	1588	#else
	1589	!!----------------------------------------------------------------------
	1590	!! Dummy module : NO KPP scheme
	1591	!!----------------------------------------------------------------------
	1592	LOGICAL, PUBLIC, PARAMETER :: lk_zdfkpp = .FALSE. !: KPP flag
	1593	CONTAINS
[2528]	1594	SUBROUTINE zdf_kpp_init ! Dummy routine
	1595	WRITE(,) 'zdf_kpp_init: You should not have seen this print! error?'
	1596	END SUBROUTINE zdf_kpp_init
	1597	SUBROUTINE zdf_kpp( kt ) ! Dummy routine
[255]	1598	WRITE(,) 'zdf_kpp: You should not have seen this print! error?', kt
	1599	END SUBROUTINE zdf_kpp
[2528]	1600	SUBROUTINE tra_kpp( kt ) ! Dummy routine
[463]	1601	WRITE(,) 'tra_kpp: You should not have seen this print! error?', kt
	1602	END SUBROUTINE tra_kpp
[2528]	1603	SUBROUTINE trc_kpp( kt ) ! Dummy routine
	1604	WRITE(,) 'trc_kpp: You should not have seen this print! error?', kt
	1605	END SUBROUTINE trc_kpp
[255]	1606	#endif
	1607
	1608	!!======================================================================
	1609	END MODULE zdfkpp

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

Download in other formats: