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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

zdfkpp.F90 in branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

source: branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 3396

Last change on this file since 3396 was 3396, checked in by acc, 12 years ago
Branch: dev_r3385_NOCS04_HAMF; #665. Stage 1 of 2012 development: porting of changes on old development branch (2011/DEV_r1837_mass_heat_salt_fluxes) into new branch. Corrected a few errors on the way. This branch now compiles but is incomplete. Still missing LIM3 changes which must reside on a certain persons laptop somewhere
Property svn:keywords set to `Id`
File size: 78.7 KB

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

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

Download in other formats: