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

Context Navigation

source: trunk/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 558

Last change on this file since 558 was 516, checked in by opalod, 18 years ago

nemo_v1_update_072 : CT : - lights modifications to ensure good control prints for the restartability & reproductibility tests

build all ocean.output_XXXX files only if logical namelist ln_ctl is TRUE, in all case keep the initial ocean.output file (require a new namelist namctl)
conditionning some prints using nprint parameter

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 74.2 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 : 8.1 ! 00-03 (W.G. Large, J. Chanut) Original code
8	!! 8.1 ! 02-06 (J.M. Molines) for real case CLIPPER
9	!! 8.2 ! 03-10 (Chanut J.) re-writting
10	!! 9.0 ! 05-01 (C. Ethe) Free form, F90
11	!!----------------------------------------------------------------------
12	#if defined key_zdfkpp \|\| defined key_esopa
13	!!----------------------------------------------------------------------
14	!! 'key_zdfkpp' KPP scheme
15	!!----------------------------------------------------------------------
16	!!----------------------------------------------------------------------
17	!! zdf_kpp : update momentum and tracer Kz from a kpp scheme
18	!! zdf_kpp_init : initialization, namelist read, and parameters control
19	!!----------------------------------------------------------------------
20	USE oce ! ocean dynamics and active tracers
21	USE dom_oce ! ocean space and time domain
22	USE zdf_oce ! ocean vertical physics
23	USE phycst ! physical constants
24	USE taumod ! surface stress
25	USE eosbn2 ! equation of state
26	USE ocesbc ! thermohaline fluxes
27	USE zdfddm ! double diffusion mixing
28	USE in_out_manager ! I/O manager
29	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
30	USE prtctl ! Print control
31
32	IMPLICIT NONE
33	PRIVATE
34
35	PUBLIC zdf_kpp ! routine called by step.F90
36	PUBLIC tra_kpp ! routine called by step.F90
37
38	LOGICAL, PUBLIC, PARAMETER :: lk_zdfkpp = .TRUE. !: KPP vertical mixing flag
39	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: ghats !: non-local scalar mixing term (gamma/<ws>o)
40	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: &
41	wt0 , & !: surface temperature flux for non local flux
42	ws0 , & !: surface salinity flux for non local flux
43	hkpp !: boundary layer depht
44
45	INTEGER :: & !!! kpp namelist (namkpp)
46	nave = 1 , & ! = 0/1 flag for horizontal average on avt, avmu, avmv
47	navb = 0 ! = 0/1 flag for constant or profile background avt
48
49	REAL(wp) :: & !!! ** Interior Mixing
50	difmiw = 1.2e-04_wp , & ! constant internal wave viscosity (m2/s)
51	difsiw = 1.2e-05_wp , & ! constant internal wave diffusivity (m2/s)
52	Riinfty = 0.8_wp , & ! local Richardson Number limit for shear instability
53	difri = 5.e-03_wp , & ! maximum shear mixing at Rig = 0 (m2/s)
54	bvsqcon = -1.e-09_wp , & ! Brunt-Vaisala squared (1/s**2) for maximum convection
55	difcon = 1._wp ! maximum mixing in interior convection (m2/s)
56	#if defined key_zdfddm
57	REAL(wp) :: & !!! ** Double diffusion Mixing
58	difssf = 1.e-03_wp , & ! maximum salt fingering mixing
59	Rrho0 = 1.9_wp , & ! limit for salt fingering mixing
60	difsdc = 1.5e-06_wp ! maximum diffusive convection mixing
61	#endif
62	LOGICAL :: &
63	ln_kpprimix = .TRUE. ! Shear instability mixing
64
65	REAL(wp) :: & !!! General constants
66	epsln = 1.0e-20_wp , & ! a small positive number
67	pthird = 1._wp/3._wp , & ! 1/3
68	pfourth = 1._wp/4._wp ! 1/4
69
70	REAL(wp) :: & !!! Boundary Layer Turbulence Parameters
71	vonk = 0.4_wp , & ! von Karman's constant
72	epsilon = 0.1_wp , & ! nondimensional extent of the surface layer
73	rconc1 = 5.0_wp , & ! standard flux profile function parmaeters
74	rconc2 = 16.0_wp , & ! " "
75	rconcm = 8.38_wp , & ! momentum flux profile fit
76	rconam = 1.26_wp , & ! " "
77	rzetam = -.20_wp , & ! " "
78	rconcs = 98.96_wp , & ! scalar flux profile fit
79	rconas = -28.86_wp , & ! " "
80	rzetas = -1.0_wp ! " "
81	REAL(wp) :: & !!! Boundary Layer Depth Diagnostic
82	Ricr = 0.3_wp , & ! critical bulk Richardson Number
83	rcekman = 0.7_wp , & ! coefficient for ekman depth
84	rcmonob = 1.0_wp , & ! coefficient for Monin-Obukhov depth
85	rconcv = 1.7_wp , & ! ratio of interior buoyancy frequency to buoyancy frequency at entrainment depth
86	hbf = 1.0_wp , & ! fraction of bound. layer depth to which absorbed solar
87	! ! rad. and contributes to surf. buo. forcing
88	Vtc ! function of rconcv,rconcs,epsilon,vonk,Ricr
89	REAL(wp) :: & !!! Nonlocal Boundary Layer Mixing
90	rcstar = 5.0_wp , & ! coefficient for convective nonlocal transport
91	rcs = 1.0e-3_wp , & ! conversion: mm/s ==> m/s
92	rcg ! non-dimensional coefficient for nonlocal transport
93
94	#if ! defined key_kppcustom
95	REAL(wp), DIMENSION(jpk,jpk) :: &
96	del ! array for reference mean values of vertical integration
97	#endif
98
99	#if defined key_kpplktb
100	INTEGER, PARAMETER :: & !!! Parameters for lookup table for turbulent velocity scales
101	nilktb = 892 , & ! number of values for zehat in KPP lookup table
102	njlktb = 482 , & ! number of values for ustar in KPP lookup table
103	nilktbm1 = nilktb - 1 , & !
104	njlktbm1 = njlktb - 1 !
105
106	REAL(wp), DIMENSION(nilktb,njlktb) :: &
107	wmlktb , & ! lookup table for the turbulent vertical velocity scale for momentum
108	wslktb ! lookup table for the turbulent vertical velocity scale for tracers
109
110	REAL(wp) :: &
111	dehatmin = -4.e-7_wp , & ! minimum limit for zhat in lookup table (m3/s3)
112	dehatmax = 0._wp , & ! maximum limit for zhat in lookup table (m3/s3)
113	ustmin = 0._wp , & ! minimum limit for ustar in lookup table (m/s)
114	ustmax = 0.04_wp , & ! maximum limit for ustar in lookup table (m/s)
115	dezehat , & ! delta zhat in lookup table
116	deustar ! delta ustar in lookup table
117	#endif
118	REAL(wp), DIMENSION(jpk) :: & !!! attenuation coef
119	ratt
120	!! already defines in module traqsr, but only if the solar radiation penetration is considered
121	REAL(wp) :: & !!! * penetrative solar radiation coefficient *
122	rabs = 0.58_wp , & ! fraction associated with xsi1
123	xsi1 = 0.35_wp , & ! first depth of extinction
124	xsi2 = 23.0_wp ! second depth of extinction
125	! ! (default values: water type Ib)
126
127	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
128	etmean , & ! coefficient used for horizontal smoothing
129	eumean , & ! at t-, u- and v-points
130	evmean
131
132	#if defined key_cfg_1d
133	REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: &
134	rig , & ! gradient Richardson number
135	rib , & ! bulk Richardson number
136	buof , & ! buoyancy forcing
137	mols ! moning-Obukhov length scale
138	REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: &
139	ekdp ! Ekman depth
140	#endif
141
142	INTEGER :: & !
143	jip = 62 , jjp = 111
144
145	!! * Substitutions
146	# include "domzgr_substitute.h90"
147	# include "vectopt_loop_substitute.h90"
148	!!----------------------------------------------------------------------
149	!! OPA 9.0 , LOCEAN-IPSL (2005)
150	!! $Header$
151	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
152	!!----------------------------------------------------------------------
153
154	CONTAINS
155
156	SUBROUTINE zdf_kpp ( kt )
157	!!----------------------------------------------------------------------
158	!! * ROUTINE zdf_kpp *
159	!!
160	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
161	!! coefficients and non local mixing using K-profile parameterization
162	!!
163	!! ** Method : The boundary layer depth hkpp is diagnosed at tracer points
164	!! from profiles of buoyancy, and shear, and the surface forcing.
165	!! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
166	!!
167	!! Kx = hkpp Wx(sigma) G(sigma)
168	!!
169	!! and the non local term ghat = Cs / Ws(sigma) / hkpp
170	!! Below hkpp the coefficients are the sum of mixing due to internal waves
171	!! shear instability and double diffusion.
172	!!
173	!! -1- Compute the now interior vertical mixing coefficients at all depths.
174	!! -2- Diagnose the boundary layer depth.
175	!! -3- Compute the now boundary layer vertical mixing coefficients.
176	!! -4- Compute the now vertical eddy vicosity and diffusivity.
177	!! -5- Smoothing
178	!!
179	!! N.B. The computation is done from jk=2 to jpkm1
180	!! Surface value of avt avmu avmv are set once a time to zero
181	!! in routine zdf_kpp_init.
182	!!
183	!! ** Action : update the non-local terms ghats
184	!! update avt, avmu, avmv (before vertical eddy coef.)
185	!!
186	!! References : Large W.G., Mc Williams J.C. and Doney S.C.
187	!! Reviews of Geophysics, 32, 4, November 1994
188	!! Comments in the code refer to this paper, particularly
189	!! the equation number. (LMD94, here after)
190	!!----------------------------------------------------------------------
191	#if defined key_zdfddm
192	USE oce , zviscos => ua, & ! temp. array for viscosities use ua as workspace
193	& zdiffut => ta, & ! temp. array for diffusivities use sa as workspace
194	& zdiffus => sa ! temp. array for diffusivities use sa as workspace
195	#else
196	USE oce , zviscos => ua, & ! temp. array for viscosities use ua as workspace
197	& zdiffut => ta ! temp. array for diffusivities use sa as workspace
198	#endif
199	!!
200	INTEGER, INTENT( in ) :: kt ! ocean time step
201	!!
202	INTEGER :: ji, jj, jk ! dummy loop indices
203	INTEGER :: ikbot, jkmax, jkm1, jkp2 !
204
205	REAL(wp), DIMENSION(jpi,jpj) :: & !!! Surface buoyancy forcing, friction velocity
206	zBo, zBosol, zustar !
207	!
208	REAL(wp) :: & !
209	ztx, zty, ztau, zflageos, & !
210	zstabl, zbuofdep,zucube, & !
211	zrhos, zalbet, zbeta, & !
212	zthermal, zhalin, zatt1 !
213
214	REAL(wp) :: & !!! Bulk richardson number
215	zref, zt, zs, zh, & !
216	zu, zv, zrh, & !
217	zrib, zrinum, & !
218	zdVsq, zVtsq !
219
220	REAL(wp) :: & !!! Velocity scales
221	zehat, zeta, zhrib, zsig, & !
222	zscale, zwst, zws, zwm
223
224	#if defined key_kpplktb
225	INTEGER :: & !!! Lookup table or Analytical functions
226	il, jl !
227	REAL(wp) :: & !
228	ud, zfrac, ufrac, & !
229	zwam, zwbm, zwas, zwbs !
230	#else
231	REAL(wp) :: & !
232	zwsun, zwmun, &
233	zcons, zconm, zwcons, zwconm !
234	#endif
235
236	REAL(wp) :: & !!! In situ density
237	zsr, zbw, ze, & !
238	zb, zd, zc, zaw, za, & !
239	zb1, za1, zkw, zk0, & !
240	zcomp , zrhd, zrhdr,zbvzed !
241
242	#if ! defined key_kppcustom
243	!! * local declarations
244	INTEGER :: &
245	jm ! dummy loop indices
246	REAL(wp) :: & !!! Compression terms
247	zr1, zr2, zr3, zr4, & !
248	zrhop !
249	#endif
250
251	REAL(wp) :: & !
252	zflag, ztemp, zrn2, & !
253	zdep21, zdep32, zdep43
254
255	REAL(wp) :: & !!! Interior richardson mixing
256	zdku2, zdkv2, ze3sqr, & !
257	zsh2, zri, zfri !
258
259	REAL(wp), DIMENSION(jpi,0:2) :: & !!! Moning-Obukov limitation
260	zmoek
261	REAL(wp), DIMENSION(jpi) :: &
262	zmoa, zekman
263	REAL(wp) :: &
264	zmob, zek
265
266	REAL(wp), DIMENSION(jpi,4) :: & !!! The pipe
267	zdepw, zdift, zvisc
268	REAL(wp), DIMENSION(jpi,3) :: &
269	zdept
270	REAL(wp), DIMENSION(jpi,2) :: &
271	zriblk
272	REAL(wp), DIMENSION(jpi,jpk) :: & !
273	zmask
274	REAL(wp), DIMENSION(jpi) :: & !
275	zhmax, zria, zhbl
276	REAL(wp) :: & !
277	zflagri, zflagek, & !
278	zflagmo, zflagh, zflagkb !
279	REAL(wp), DIMENSION(jpi) :: & !!! Shape function (G)
280	za2m, za3m, zkmpm, &
281	za2t, za3t, zkmpt
282	REAL(wp) :: & !
283	zdelta, zdelta2, & !
284	zdzup, zdzdn, zdzh, & !
285	zvath, zgat1, zdat1, & !
286	zkm1m, zkm1t
287	REAL(wp), DIMENSION(jpi,jpk) :: & !!! Boundary layer diffusivities/viscosities
288	zblcm, zblct
289	#if defined key_zdfddm
290	REAL(wp) :: & !!! double diffusion mixing
291	zrrau, zds, &
292	zavdds, zavddt,zinr
293	REAL(wp), DIMENSION(jpi,4) :: &
294	zdifs
295	REAL(wp), DIMENSION(jpi) :: &
296	za2s, za3s, zkmps
297	REAL(wp) :: &
298	zkm1s
299	REAL(wp), DIMENSION(jpi,jpk) :: &
300	zblcs
301	#endif
302	!!--------------------------------------------------------------------
303
304
305	! Initialization (first time-step only)
306	! --------------
307	IF( kt == nit000 ) CALL zdf_kpp_init
308
309	zviscos(:,:,:) = 0.
310	zblcm (:,: ) = 0.
311	zdiffut(:,:,:) = 0.
312	zblct (:,: ) = 0.
313	#if defined key_zdfddm
314	zdiffus(:,:,:) = 0.
315	zblcs (:,: ) = 0.
316	#endif
317	ghats(:,:,:) = 0.
318
319	zBo (:,:) = 0.
320	zBosol(:,:) = 0.
321	zustar(:,:) = 0.
322
323
324	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
325	! I. Interior diffusivity and viscosity at w points ( T interfaces)
326	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
327	DO jk = 2, jpkm1
328	DO jj = 2, jpjm1
329	DO ji = fs_2, fs_jpim1
330	! Mixing due to internal waves breaking
331	! -------------------------------------
332	avmu(ji,jj,jk) = difmiw
333	avt (ji,jj,jk) = difsiw
334	! Mixing due to vertical shear instability
335	! -------------------------------------
336	IF( ln_kpprimix ) THEN
337	! Compute the gradient Richardson number at interfaces (zri):
338	! LMD94, eq. 27 (is vertical smoothing needed : Rig=N^2 / (dz(u))^2 + (dz(v))^2
339	zdku2 = ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
340	& * ( un(ji - 1,jj,jk - 1) - un(ji - 1,jj,jk) ) &
341	& + ( un(ji ,jj,jk - 1) - un(ji ,jj,jk) ) &
342	& * ( un(ji ,jj,jk - 1) - un(ji ,jj,jk) )
343
344	zdkv2 = ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
345	& * ( vn(ji,jj - 1,jk - 1) - vn(ji,jj - 1,jk) ) &
346	& + ( vn(ji, jj,jk - 1) - vn(ji, jj,jk) ) &
347	& * ( vn(ji, jj,jk - 1) - vn(ji, jj,jk) )
348
349	ze3sqr = 1. / ( fse3w(ji,jj,jk) * fse3w(ji,jj,jk) )
350	! Square of vertical shear at interfaces
351	zsh2 = 0.5 * 0.5 * ( zdku2 + zdkv2 ) * ze3sqr
352	zri = MAX( rn2(ji,jj,jk), 0. ) / ( zsh2 + epsln )
353	#if defined key_cfg_1d
354	! save the gradient richardson number
355	rig(ji,jj,jk) = zri * tmask(ji,jj,jk)
356	#endif
357	! Evaluate f of Ri (zri) for shear instability store in zfri
358	! LMD94, eq. 28a,b,c, figure 3 ; Rem: p1 is 3, hard coded
359	zfri = MAX( zri , 0. )
360	zfri = MIN( zfri / Riinfty , 1.0 )
361	zfri = ( 1.0 - zfri * zfri )
362	zfri = zfri * zfri * zfri
363	! add shear contribution to mixing coef.
364	avmu(ji,jj,jk) = avmu(ji,jj,jk) + difri * zfri
365	avt (ji,jj,jk) = avt (ji,jj,jk) + difri * zfri
366	ENDIF
367	#if defined key_zdfddm
368	avs (ji,jj,jk) = avt (ji,jj,jk)
369	! Double diffusion mixing ; NOT IN ROUTINE ZDFDDM.F90
370	! ------------------------------------------------------------------
371	! only retains positive value of rrau
372	zrrau = MAX( rrau(ji,jj,jk), epsln )
373	zds = sn(ji,jj,jk-1) - sn(ji,jj,jk)
374	IF( zrrau > 1. .AND. zds > 0.) THEN
375	!
376	! Salt fingering case.
377	!---------------------
378	! Compute interior diffusivity for double diffusive mixing of
379	! salinity. Upper bound "zrrau" by "Rrho0"; (Rrho0=1.9, difcoefnuf=0.001).
380	! After that set interior diffusivity for double diffusive mixing
381	! of temperature
382	zavdds = MIN( zrrau, Rrho0 )
383	zavdds = ( zavdds - 1.0 ) / ( Rrho0 - 1.0 )
384	zavdds = 1.0 - zavdds * zavdds
385	zavdds = zavdds * zavdds * zavdds
386	zavdds = difssf * zavdds
387	zavddt = 0.7 * zavdds
388	ELSEIF( zrrau < 1. .AND. zrrau > 0. .AND. zds < 0.) THEN
389	!
390	! Diffusive convection case.
391	!---------------------------
392	! Compute interior diffusivity for double diffusive mixing of
393	! temperature (Marmorino and Caldwell, 1976);
394	! Compute interior diffusivity for double diffusive mixing of salinity
395	zinr = 1. / zrrau
396	zavddt = 0.909 * EXP( 4.6 * EXP( -0.54* ( zinr - 1. ) ) )
397	zavddt = difsdc * zavddt
398	IF( zrrau < 0.5) THEN
399	zavdds = zavddt * 0.15 * zrrau
400	ELSE
401	zavdds = zavddt * (1.85 * zrrau - 0.85 )
402	ENDIF
403	ELSE
404	zavddt = 0.
405	zavdds = 0.
406	ENDIF
407	! Add double diffusion contribution to temperature and salinity mixing coefficients.
408	avt (ji,jj,jk) = avt (ji,jj,jk) + zavddt
409	avs (ji,jj,jk) = avs (ji,jj,jk) + zavdds
410	#endif
411	END DO
412	END DO
413	END DO
414
415
416	! Radiative (zBosol) and non radiative (zBo) surface buoyancy
417	!JMM at the time zdfkpp is called, q still holds the sum q + qsr
418	!---------------------------------------------------------------------
419	DO jj = 2, jpjm1
420	DO ji = fs_2, fs_jpim1
421	IF( neos < 1) THEN
422	zt = tn(ji,jj,1)
423	zs = sn(ji,jj,1) - 35.0
424	zh = fsdept(ji,jj,1)
425	! potential volumic mass
426	zrhos = rhop(ji,jj,1)
427	zalbet = ( ( ( - 0.255019e-07 * zt + 0.298357e-05 ) * zt & ! ratio alpha/beta
428	& - 0.203814e-03 ) * zt &
429	& + 0.170907e-01 ) * zt &
430	& + 0.665157e-01 &
431	& + ( - 0.678662e-05 * zs &
432	& - 0.846960e-04 * zt + 0.378110e-02 ) * zs &
433	& + ( ( - 0.302285e-13 * zh &
434	& - 0.251520e-11 * zs &
435	& + 0.512857e-12 * zt * zt ) * zh &
436	& - 0.164759e-06 * zs &
437	& +( 0.791325e-08 * zt - 0.933746e-06 ) * zt &
438	& + 0.380374e-04 ) * zh
439
440	zbeta = ( ( -0.415613e-09 * zt + 0.555579e-07 ) * zt & ! beta
441	& - 0.301985e-05 ) * zt &
442	& + 0.785567e-03 &
443	& + ( 0.515032e-08 * zs &
444	& + 0.788212e-08 * zt - 0.356603e-06 ) * zs &
445	& +( ( 0.121551e-17 * zh &
446	& - 0.602281e-15 * zs &
447	& - 0.175379e-14 * zt + 0.176621e-12 ) * zh &
448	& + 0.408195e-10 * zs &
449	& + ( - 0.213127e-11 * zt + 0.192867e-09 ) * zt &
450	& - 0.121555e-07 ) * zh
451
452	zthermal = zbeta * zalbet / ( rcp * zrhos + epsln )
453	zhalin = zbeta * sn(ji,jj,1) * rcs
454	ELSE
455	zrhos = rhop(ji,jj,1) + rau0 * ( 1. - tmask(ji,jj,1) )
456	zthermal = ralpha / ( rcp * zrhos + epsln )
457	zhalin = rbeta * sn(ji,jj,1) * rcs
458	ENDIF
459	! Radiative surface buoyancy force
460	zBosol(ji,jj) = grav * zthermal * qsr(ji,jj)
461	! Non radiative surface buoyancy force
462	zBo (ji,jj) = grav * zthermal * ( qt(ji,jj) - qsr(ji,jj) ) - grav * zhalin * emp(ji,jj)
463	! Surface Temperature flux for non-local term
464	wt0(ji,jj) = - qt(ji,jj) * ro0cpr * tmask(ji,jj,1)
465	! Surface salinity flux for non-local term
466	ws0(ji,jj) = - ( emp(ji,jj) * sn(ji,jj,1) * rcs ) * tmask(ji,jj,1)
467	ENDDO
468	ENDDO
469
470	zflageos = 0.5 + SIGN( 0.5, neos - 1. )
471	! Compute surface buoyancy forcing, Monin Obukhov and Ekman depths
472	!------------------------------------------------------------------
473	DO jj = 2, jpjm1
474	DO ji = fs_2, fs_jpim1
475	! Reference surface density = density at first T point level
476	zrhos = rhop(ji,jj,1) + zflageos * rau0 * ( 1. - tmask(ji,jj,1) )
477	! Friction velocity (zustar), at T-point : LMD94 eq. 2
478	ztx = 0.5 * ( taux(ji,jj) + taux(ji - 1, jj ) )
479	zty = 0.5 * ( tauy(ji,jj) + tauy(ji , jj - 1) )
480	ztau = SQRT( ztx * ztx + zty * zty )
481	zustar(ji,jj) = SQRT( ztau / ( 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 = mbathy(ji,jj) - 1 ! 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_cfg_1d
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 ( neos < 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_cfg_1d
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 ( nave )
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	END SUBROUTINE zdf_kpp
1242
1243
1244	SUBROUTINE tra_kpp
1245	!!----------------------------------------------------------------------
1246	!! * ROUTINE tra_kpp *
1247	!!
1248	!! ** Purpose : compute and add to the tracer trend the non-local
1249	!! tracer flux
1250	!!
1251	!! ** Method : ???
1252	!!
1253	!! history :
1254	!! 9.0 ! 05-11 (G. Madec) Original code
1255	!!----------------------------------------------------------------------
1256	!! * Modules used
1257	USE oce, ONLY : ztrdt => ua, & ! use ua as 3D workspace
1258	ztrds => va ! use va as 3D workspace
1259	!!----------------------------------------------------------------------
1260
1261	IF( kt == nit000 ) THEN
1262	IF(lwp) WRITE(numout,*)
1263	IF(lwp) WRITE(numout,*) 'tra_kpp : KPP non-local tracer fluxes'
1264	IF(lwp) WRITE(numout,*) '~~~~~~~ '
1265	ENDIF
1266
1267
1268	! Save ta and sa trends
1269	IF( l_trdtra ) THEN
1270	ztrdt(:,:,:) = ta(:,:,:)
1271	ztrds(:,:,:) = sa(:,:,:)
1272	ENDIF
1273
1274	! add non-local temperature and salinity flux ( in convective case only)
1275	DO jk = 1, jpkm1
1276	DO jj = 2, jpjm1
1277	DO ji = fs_2, fs_jpim1
1278	ta(ji,jj,jk) = ta(ji,jj,jk) &
1279	& - ( ghats(ji,jj,jk ) * avt(ji,jj,jk ) &
1280	& - ghats(ji,jj,jk+1) * avt(ji,jj,jk+1) ) * wt0(ji,jj) / fse3t(ji,jj,jk)
1281	sa(ji,jj,jk) = sa(ji,jj,jk) &
1282	& - ( ghats(ji,jj,jk ) * fsavs(ji,jj,jk ) &
1283	& - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) * ws0(ji,jj) / fse3t(ji,jj,jk)
1284	END DO
1285	END DO
1286	END DO
1287
1288	! save the non-local tracer flux trends for diagnostic
1289	IF( l_trdtra ) THEN
1290	ztrdt(:,:,:) = ta(:,:,:) - ztrdt(:,:,:)
1291	ztrds(:,:,:) = sa(:,:,:) - ztrds(:,:,:)
1292	!!bug gm jpttdzdf ==> jpttkpp
1293	CALL trd_mod(ztrdt, ztrds, jptra_trd_zdf, 'TRA', kt)
1294	ENDIF
1295
1296	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' kpp - Ta: ', mask1=tmask, &
1297	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
1298	ENDIF
1299
1300	END SUBROUTINE tra_kpp
1301
1302
1303	SUBROUTINE zdf_kpp_init
1304	!!----------------------------------------------------------------------
1305	!! * ROUTINE zdf_kpp_init *
1306	!!
1307	!! ** Purpose : Initialization of the vertical eddy diffivity and
1308	!! viscosity when using a kpp turbulent closure scheme
1309	!!
1310	!! ** Method : Read the namkpp namelist and check the parameters
1311	!! called at the first timestep (nit000)
1312	!!
1313	!! ** input : Namlist namkpp
1314	!!
1315	!!
1316	!! history :
1317	!! 8.1 ! 00-02 (J. Chanut) KPP Mixing
1318	!! 9.0 ! 05-01 (C. Ethe) F90 : free form
1319	!!----------------------------------------------------------------------
1320	!! * local declarations
1321
1322	INTEGER :: &
1323	ji, jj, jk ! dummy loop indices
1324
1325	#if ! defined key_kppcustom
1326	INTEGER :: &
1327	jm ! dummy loop indices
1328	REAL(wp) :: & !!! tempory scalars
1329	zref, zdist
1330	#endif
1331
1332	#if defined key_kpplktb
1333	REAL(wp) :: & !!! tempory scalars
1334	zustar, &
1335	zucube, zustvk, &
1336	zeta, zehat
1337	#endif
1338	REAL(wp) :: & !!! tempory scalars
1339	zhbf
1340	LOGICAL :: &
1341	ll_kppcustom, & ! 1st ocean level taken as surface layer
1342	ll_kpplktb ! Lookup table for turbul. velocity scales
1343
1344	NAMELIST/namkpp/ ln_kpprimix, difmiw, difsiw, Riinfty, difri, bvsqcon, difcon, nave, navb
1345
1346	!!----------------------------------------------------------------------
1347
1348	! Read Namelist namtke : Turbulente Kinetic Energy
1349	! --------------------
1350	REWIND ( numnam )
1351	READ ( numnam, namkpp )
1352
1353	! Parameter control and print
1354	! ---------------------------
1355	REWIND( numnam )
1356	READ ( numnam, namkpp )
1357
1358	! Control print
1359	IF(lwp) THEN
1360	WRITE(numout,*)
1361	WRITE(numout,*) 'zdf_kpp_init : kpp turbulent closure scheme'
1362	WRITE(numout,*) '~~~~~~~~~~~~'
1363	WRITE(numout,*) ' Namelist namkpp : set tke mixing parameters'
1364	WRITE(numout,*) ' Shear instability mixing ln_kpprimix = ', ln_kpprimix
1365	WRITE(numout,*) ' max. internal wave viscosity difmiw = ', difmiw
1366	WRITE(numout,*) ' max. internal wave diffusivity difsiw = ', difsiw
1367	WRITE(numout,*) ' Richardson Number limit for shear instability Riinfty = ', Riinfty
1368	WRITE(numout,*) ' max. shear mixing at Rig = 0 difri = ', difri
1369	WRITE(numout,*) ' Brunt-Vaisala squared for max. convection bvsqcon = ', bvsqcon
1370	WRITE(numout,*) ' max. mix. in interior convec. difcon = ', difcon
1371	WRITE(numout,*) ' horizontal average flag nave = ', nave
1372	WRITE(numout,*) ' constant background or profile navb = ', navb
1373	ENDIF
1374
1375	ll_kppcustom = .FALSE.
1376	ll_kpplktb = .FALSE.
1377
1378	#if defined key_kppcustom
1379	ll_kppcustom = .TRUE.
1380	#endif
1381	#if defined key_kpplktb
1382	ll_kpplktb = .TRUE.
1383	#endif
1384	IF(lwp) THEN
1385	WRITE(numout,*) ' Lookup table for turbul. velocity scales ll_kpplktb = ', ll_kpplktb
1386	WRITE(numout,*) ' 1st ocean level taken as surface layer ll_kppcustom = ', ll_kppcustom
1387	WRITE(numout,*) ' '
1388	ENDIF
1389
1390	IF( lk_zdfddm) THEN
1391	IF(lwp) THEN
1392	WRITE(numout,*)
1393	WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
1394	WRITE(numout,*) ' CAUTION : done in routine zdfkpp, not in routine zdfddm '
1395	WRITE(numout,*) ' '
1396	ENDIF
1397	ENDIF
1398
1399
1400
1401	!set constants not in namelist
1402	!-----------------------------
1403	Vtc = rconcv * SQRT( 0.2 / ( rconcs * epsilon ) ) / ( vonk * vonk * Ricr )
1404	rcg = rcstar * vonk * ( rconcs * vonk * epsilon )**pthird
1405
1406	IF(lwp) THEN
1407	WRITE(numout,*) ' Constant value for unreso. turbul. velocity shear Vtc = ', Vtc
1408	WRITE(numout,*) ' Non-dimensional coef. for nonlocal transport rcg = ', rcg
1409	WRITE(numout,*) ' '
1410	ENDIF
1411
1412	! ratt is the attenuation coefficient for solar flux
1413	! Should be different is s_coordinate
1414	DO jk = 1, jpk
1415	zhbf = - fsdept(1,1,jk) * hbf
1416	ratt(jk) = 1.0 - ( rabs * EXP( zhbf / xsi1 ) + ( 1.0 - rabs ) * EXP( zhbf / xsi2 ) )
1417	ENDDO
1418
1419	! Horizontal average : initialization of weighting arrays
1420	! -------------------
1421
1422	SELECT CASE ( nave )
1423
1424	CASE ( 0 ) ! no horizontal average
1425	IF(lwp) WRITE(numout,*) ' no horizontal average on avt, avmu, avmv'
1426	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
1427	! weighting mean arrays etmean, eumean and evmean
1428	! ( 1 1 ) ( 1 )
1429	! avt = 1/4 ( 1 1 ) avmu = 1/2 ( 1 1 ) avmv= 1/2 ( 1 )
1430	!
1431	etmean(:,:,:) = 0.e0
1432	eumean(:,:,:) = 0.e0
1433	evmean(:,:,:) = 0.e0
1434
1435	DO jk = 1, jpkm1
1436	DO jj = 2, jpjm1
1437	DO ji = 2, jpim1 ! vector opt.
1438	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
1439	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
1440	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
1441
1442	eumean(ji,jj,jk) = umask(ji,jj,jk) &
1443	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji+1,jj ,jk) )
1444
1445	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
1446	& / MAX( 1., tmask(ji,jj,jk) + tmask(ji ,jj+1,jk) )
1447	END DO
1448	END DO
1449	END DO
1450
1451	CASE ( 1 ) ! horizontal average
1452	IF(lwp) WRITE(numout,*) ' horizontal average on avt, avmu, avmv'
1453	! weighting mean arrays etmean, eumean and evmean
1454	! ( 1/2 1 1/2 ) ( 1/2 1/2 ) ( 1/2 1 1/2 )
1455	! avt = 1/8 ( 1 2 1 ) avmu = 1/4 ( 1 1 ) avmv= 1/4 ( 1/2 1 1/2 )
1456	! ( 1/2 1 1/2 ) ( 1/2 1/2 )
1457	etmean(:,:,:) = 0.e0
1458	eumean(:,:,:) = 0.e0
1459	evmean(:,:,:) = 0.e0
1460
1461	DO jk = 1, jpkm1
1462	DO jj = 2, jpjm1
1463	DO ji = fs_2, fs_jpim1 ! vector opt.
1464	etmean(ji,jj,jk) = tmask(ji, jj,jk) &
1465	& / MAX( 1., 2.* tmask(ji,jj,jk) &
1466	& +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) &
1467	& +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
1468	& +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) &
1469	& +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) )
1470
1471	eumean(ji,jj,jk) = umask(ji,jj,jk) &
1472	& / MAX( 1., tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) &
1473	& +.5 * ( tmask(ji,jj-1,jk) + tmask(ji+1,jj-1,jk) &
1474	& +tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) ) )
1475
1476	evmean(ji,jj,jk) = vmask(ji,jj,jk) &
1477	& / MAX( 1., tmask(ji ,jj,jk) + tmask(ji ,jj+1,jk) &
1478	& +.5 * ( tmask(ji-1,jj,jk) + tmask(ji-1,jj+1,jk) &
1479	& +tmask(ji+1,jj,jk) + tmask(ji+1,jj+1,jk) ) )
1480	END DO
1481	END DO
1482	END DO
1483
1484	CASE DEFAULT
1485	WRITE(ctmp1,*) ' bad flag value for nave = ', nave
1486	CALL ctl_opa( ctmp1 )
1487
1488	END SELECT
1489
1490
1491
1492	! Background eddy viscosity and diffusivity profil
1493	! ------------------------------------------------
1494	IF( navb == 0 ) THEN
1495	! Define avmb, avtb from namelist parameter
1496	avmb(:) = avm0
1497	avtb(:) = avt0
1498	ELSE
1499	! Background profile of avt (fit a theoretical/observational profile (Krauss 1990)
1500	avmb(:) = avm0
1501	avtb(:) = 1.e-5 + 2.8e-8 * gdepw(:) ! m2/s
1502	ENDIF
1503
1504	! Increase the background in the surface layers
1505	!! avmb(1) = 1. * avmb(1) ; avtb(1) = 1. * avtb(1)
1506	!! avmb(2) = 1. * avmb(2) ; avtb(2) = 1. * avtb(2)
1507	!! avmb(3) = 1. * avmb(3) ; avtb(3) = 1. * avtb(3)
1508	!! avmb(4) = 1. * avmb(4) ; avtb(4) = 1. * avtb(4)
1509
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	END SUBROUTINE zdf_kpp_init
1575
1576	#else
1577	!!----------------------------------------------------------------------
1578	!! Dummy module : NO KPP scheme
1579	!!----------------------------------------------------------------------
1580	LOGICAL, PUBLIC, PARAMETER :: lk_zdfkpp = .FALSE. !: KPP flag
1581	CONTAINS
1582	SUBROUTINE zdf_kpp( kt ) ! Empty routine
1583	WRITE(,) 'zdf_kpp: You should not have seen this print! error?', kt
1584	END SUBROUTINE zdf_kpp
1585	SUBROUTINE tra_kpp( kt ) ! Empty routine
1586	WRITE(,) 'tra_kpp: You should not have seen this print! error?', kt
1587	END SUBROUTINE tra_kpp
1588	#endif
1589
1590	!!======================================================================
1591	END MODULE zdfkpp

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

Download in other formats: