Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

zdfkpp.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 82.0 KB

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 11738

Download in other formats: