Context Navigation

zdfkpp.F90 @ 9383

Last change on this file since 9383 was 9383, checked in by andmirek, 6 years ago
#2050 fixes and changes
File size: 79.0 KB

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

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

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

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfkpp.F90 @ 9383

Download in other formats: