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

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

zdfkpp.F90 in branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/ZDF – NEMO

Context Navigation

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

Last change on this file since 2613 was 2613, checked in by gm, 13 years ago
dynamic mem: #785 ; move the allocation of ice in iceini_2/iceini module + bug fixes (define key_esopa)
Property svn:keywords set to `Id`
File size: 78.6 KB

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

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

Download in other formats: