Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

zdfkpp.F90 @ 12555

Last change on this file since 12555 was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

svn merge -r 11035:11101 svn+ssh://charris@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package

File size: 78.6 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 12555

Download in other formats: