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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 2590

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

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

Download in other formats: