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.

zdfosm.F90 in NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90 @ 14412

Last change on this file since 14412 was 14412, checked in by agn, 3 years ago
revert to simple hmle=zmld
Property svn:keywords set to `Id`
File size: 163.2 KB

Line
1	MODULE zdfosm
2	!!======================================================================
3	!! * MODULE zdfosm *
4	!! Ocean physics: vertical mixing coefficient compute from the OSMOSIS
5	!! turbulent closure parameterization
6	!!=====================================================================
7	!! History : NEMO 4.0 ! A. Grant, G. Nurser
8	!! 15/03/2017 Changed calculation of pycnocline thickness in unstable conditions and stable conditions AG
9	!! 15/03/2017 Calculation of pycnocline gradients for stable conditions changed. Pycnocline gradients now depend on stability of the OSBL. A.G
10	!! 06/06/2017 (1) Checks on sign of buoyancy jump in calculation of OSBL depth. A.G.
11	!! (2) Removed variable zbrad0, zbradh and zbradav since they are not used.
12	!! (3) Approximate treatment for shear turbulence.
13	!! Minimum values for zustar and zustke.
14	!! Add velocity scale, zvstr, that tends to zustar for large Langmuir numbers.
15	!! Limit maximum value for Langmuir number.
16	!! Use zvstr in definition of stability parameter zhol.
17	!! (4) Modified parametrization of entrainment flux, changing original coefficient 0.0485 for Langmuir contribution to 0.135 * zla
18	!! (5) For stable boundary layer add factor that depends on length of timestep to 'slow' collapse and growth. Make sure buoyancy jump not negative.
19	!! (6) For unstable conditions when growth is over multiple levels, limit change to maximum of one level per cycle through loop.
20	!! (7) Change lower limits for loops that calculate OSBL averages from 1 to 2. Large gradients between levels 1 and 2 can cause problems.
21	!! (8) Change upper limits from ibld-1 to ibld.
22	!! (9) Calculation of pycnocline thickness in unstable conditions. Check added to ensure that buoyancy jump is positive before calculating Ri.
23	!! (10) Thickness of interface layer at base of the stable OSBL set by Richardson number. Gives continuity in transition from unstable OSBL.
24	!! (11) Checks that buoyancy jump is poitive when calculating pycnocline profiles.
25	!! (12) Replace zwstrl with zvstr in calculation of eddy viscosity.
26	!! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information
27	!! (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence.
28	!! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added.
29	!! (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out)
30	!! (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out)
31	!! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made,
32	!! (a) Pycnocline temperature and salinity profies changed for unstable layers
33	!! (b) The stable OSBL depth parametrization changed.
34	!! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code.
35	!! 23/05/19 (20) Old code where key_osmldpth1` is not set removed, together with the key key_osmldpth1
36	!!----------------------------------------------------------------------
37
38	!!----------------------------------------------------------------------
39	!! 'ln_zdfosm' OSMOSIS scheme
40	!!----------------------------------------------------------------------
41	!! zdf_osm : update momentum and tracer Kz from osm scheme
42	!! zdf_osm_init : initialization, namelist read, and parameters control
43	!! osm_rst : read (or initialize) and write osmosis restart fields
44	!! tra_osm : compute and add to the T & S trend the non-local flux
45	!! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD)
46	!! dyn_osm : compute and add to u & v trensd the non-local flux
47	!!
48	!! Subroutines in revised code.
49	!!----------------------------------------------------------------------
50	USE oce ! ocean dynamics and active tracers
51	! uses wn from previous time step (which is now wb) to calculate hbl
52	USE dom_oce ! ocean space and time domain
53	USE zdf_oce ! ocean vertical physics
54	USE sbc_oce ! surface boundary condition: ocean
55	USE sbcwave ! surface wave parameters
56	USE phycst ! physical constants
57	USE eosbn2 ! equation of state
58	USE traqsr ! details of solar radiation absorption
59	USE zdfddm ! double diffusion mixing (avs array)
60	! USE zdfmxl ! mixed layer depth
61	USE iom ! I/O library
62	USE lib_mpp ! MPP library
63	USE trd_oce ! ocean trends definition
64	USE trdtra ! tracers trends
65	!
66	USE in_out_manager ! I/O manager
67	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
68	USE prtctl ! Print control
69	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
70
71	IMPLICIT NONE
72	PRIVATE
73
74	PUBLIC zdf_osm ! routine called by step.F90
75	PUBLIC zdf_osm_init ! routine called by nemogcm.F90
76	PUBLIC osm_rst ! routine called by step.F90
77	PUBLIC tra_osm ! routine called by step.F90
78	PUBLIC trc_osm ! routine called by trcstp.F90
79	PUBLIC dyn_osm ! routine called by step.F90
80
81	PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90
82
83	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux
84	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamv !: non-local v-momentum flux
85	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamt !: non-local temperature flux (gamma/<ws>o)
86	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghams !: non-local salinity flux (gamma/<ws>o)
87	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt
88	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth
89	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline
90	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth
91	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift.
92
93	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts
94	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmle ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization
95	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdx_mle ! zonal buoyancy gradient in ML
96	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdy_mle ! meridional buoyancy gradient in ML
97	INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_prof ! level of base of MLE layer.
98
99	! !! Namelist namzdf_osm
100	LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la
101
102	LOGICAL :: ln_osm_mle !: flag to activate the Mixed Layer Eddy (MLE) parameterisation
103
104	REAL(wp) :: rn_osm_la ! Turbulent Langmuir number
105	REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift
106	REAL(wp) :: rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by
107	REAL(wp) :: rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl
108	LOGICAL :: ln_zdfosm_ice_shelter ! flag to activate ice sheltering
109	REAL(wp) :: rn_osm_hbl0 = 10._wp ! Initial value of hbl for 1D runs
110	INTEGER :: nn_ave ! = 0/1 flag for horizontal average on avt
111	INTEGER :: nn_osm_wave = 0 ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
112	INTEGER :: nn_osm_SD_reduce ! = 0/1/2 flag for getting effective stokes drift from surface value
113	LOGICAL :: ln_dia_osm ! Use namelist rn_osm_la
114
115	LOGICAL :: ln_kpprimix = .true. ! Shear instability mixing
116	REAL(wp) :: rn_riinfty = 0.7 ! local Richardson Number limit for shear instability
117	REAL(wp) :: rn_difri = 0.005 ! maximum shear mixing at Rig = 0 (m2/s)
118	LOGICAL :: ln_convmix = .true. ! Convective instability mixing
119	REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s)
120
121	! OSMOSIS mixed layer eddy parametrization constants
122	INTEGER :: nn_osm_mle ! = 0/1 flag for horizontal average on avt
123	REAL(wp) :: rn_osm_mle_ce ! MLE coefficient
124	! ! parameters used in nn_osm_mle = 0 case
125	REAL(wp) :: rn_osm_mle_lf ! typical scale of mixed layer front
126	REAL(wp) :: rn_osm_mle_time ! time scale for mixing momentum across the mixed layer
127	! ! parameters used in nn_osm_mle = 1 case
128	REAL(wp) :: rn_osm_mle_lat ! reference latitude for a 5 km scale of ML front
129	LOGICAL :: ln_osm_hmle_limit ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
130	REAL(wp) :: rn_osm_hmle_limit ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
131	REAL(wp) :: rn_osm_mle_rho_c ! Density criterion for definition of MLD used by FK
132	REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation
133	REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld
134	REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case
135	REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base.
136	REAL(wp) :: rn_osm_bl_thresh ! Threshold buoyancy for deepening of OSBL base.
137	REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE.
138
139
140	! !!! General constants
141	REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero
142	REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw
143	REAL(wp) :: pthird = 1._wp/3._wp ! 1/3
144	REAL(wp) :: p2third = 2._wp/3._wp ! 2/3
145
146	INTEGER :: idebug = 236
147	INTEGER :: jdebug = 228
148	!!----------------------------------------------------------------------
149	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
150	!! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $
151	!! Software governed by the CeCILL license (see ./LICENSE)
152	!!----------------------------------------------------------------------
153	CONTAINS
154
155	INTEGER FUNCTION zdf_osm_alloc()
156	!!----------------------------------------------------------------------
157	!! * FUNCTION zdf_osm_alloc *
158	!!----------------------------------------------------------------------
159	ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &
160	& hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
161	& etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
162
163	ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), &
164	& mld_prof(jpi,jpj), STAT= zdf_osm_alloc )
165
166	! ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & ! would ths be better ?
167	! & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
168	! & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
169	! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
170	! IF ( ln_osm_mle ) THEN
171	! Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc )
172	! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays')
173	! ENDIF
174
175	IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
176	CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
177	END FUNCTION zdf_osm_alloc
178
179
180	SUBROUTINE zdf_osm( kt, p_avm, p_avt )
181	!!----------------------------------------------------------------------
182	!! * ROUTINE zdf_osm *
183	!!
184	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
185	!! coefficients and non local mixing using the OSMOSIS scheme
186	!!
187	!! ** Method : The boundary layer depth hosm is diagnosed at tracer points
188	!! from profiles of buoyancy, and shear, and the surface forcing.
189	!! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
190	!!
191	!! Kx = hosm Wx(sigma) G(sigma)
192	!!
193	!! and the non local term ghamt = Cs / Ws(sigma) / hosm
194	!! Below hosm the coefficients are the sum of mixing due to internal waves
195	!! shear instability and double diffusion.
196	!!
197	!! -1- Compute the now interior vertical mixing coefficients at all depths.
198	!! -2- Diagnose the boundary layer depth.
199	!! -3- Compute the now boundary layer vertical mixing coefficients.
200	!! -4- Compute the now vertical eddy vicosity and diffusivity.
201	!! -5- Smoothing
202	!!
203	!! N.B. The computation is done from jk=2 to jpkm1
204	!! Surface value of avt are set once a time to zero
205	!! in routine zdf_osm_init.
206	!!
207	!! ** Action : update the non-local terms ghamts
208	!! update avt (before vertical eddy coef.)
209	!!
210	!! References : Large W.G., Mc Williams J.C. and Doney S.C.
211	!! Reviews of Geophysics, 32, 4, November 1994
212	!! Comments in the code refer to this paper, particularly
213	!! the equation number. (LMD94, here after)
214	!!----------------------------------------------------------------------
215	INTEGER , INTENT(in ) :: kt ! ocean time step
216	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
217	!!
218	INTEGER :: ji, jj, jk ! dummy loop indices
219
220	INTEGER :: jl ! dummy loop indices
221
222	INTEGER :: ikbot, jkmax, jkm1, jkp2 !
223
224	REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube !
225	REAL(wp) :: zbeta, zthermal !
226	REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
227	REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm !
228
229	REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density
230	INTEGER :: jm ! dummy loop indices
231	REAL(wp) :: zr1, zr2, zr3, zr4, zrhop ! Compression terms
232	REAL(wp) :: zflag, zrn2, zdep21, zdep32, zdep43
233	REAL(wp) :: zesh2, zri, zfri ! Interior richardson mixing
234	REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
235	REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages
236	! Scales
237	REAL(wp), DIMENSION(jpi,jpj) :: zrad0 ! Surface solar temperature flux (deg m/s)
238	REAL(wp), DIMENSION(jpi,jpj) :: zradh ! Radiative flux at bl base (Buoyancy units)
239	REAL(wp), DIMENSION(jpi,jpj) :: zradav ! Radiative flux, bl average (Buoyancy Units)
240	REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity
241	REAL(wp), DIMENSION(jpi,jpj) :: zwstrl ! Langmuir velocity scale
242	REAL(wp), DIMENSION(jpi,jpj) :: zvstr ! Velocity scale that ends to zustar for large Langmuir number.
243	REAL(wp), DIMENSION(jpi,jpj) :: zwstrc ! Convective velocity scale
244	REAL(wp), DIMENSION(jpi,jpj) :: zuw0 ! Surface u-momentum flux
245	REAL(wp), DIMENSION(jpi,jpj) :: zvw0 ! Surface v-momentum flux
246	REAL(wp), DIMENSION(jpi,jpj) :: zwth0 ! Surface heat flux (Kinematic)
247	REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux
248	REAL(wp), DIMENSION(jpi,jpj) :: zwb0 ! Surface buoyancy flux
249	REAL(wp), DIMENSION(jpi,jpj) :: zwb0tot ! Total surface buoyancy flux including insolation
250	REAL(wp), DIMENSION(jpi,jpj) :: zwthav ! Heat flux - bl average
251	REAL(wp), DIMENSION(jpi,jpj) :: zwsav ! freshwater flux - bl average
252	REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average
253	REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux
254	REAL(wp), DIMENSION(jpi,jpj) :: zwb_min
255
256
257	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL
258	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux
259	REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff
260	REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK
261
262	REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift
263	REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number
264	REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
265	REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
266	REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer
267	LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl
268	LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers
269	LOGICAL, DIMENSION(jpi,jpj) :: lpyc ! OSBL pycnocline present
270	LOGICAL, DIMENSION(jpi,jpj) :: lflux ! surface flux extends below OSBL into MLE layer.
271	LOGICAL, DIMENSION(jpi,jpj) :: lmle ! MLE layer increases in hickness.
272
273	! mixed-layer variables
274
275	INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
276	INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
277	INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level
278	INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer
279
280	REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
281	REAL(wp) :: zugrad,zvgrad ! temporary variables for calculating pycnocline shear
282
283	REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid
284	REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid
285
286	REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid
287	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid
288
289	REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid
290	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
291	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients
292	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext ! external temperature/salinity and buoyancy gradients
293	REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization.
294
295	REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl ! averages over the depth of the blayer
296	REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml ! averages over the depth of the mixed layer
297	REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle ! averages over the depth of the MLE layer
298	REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer
299	REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
300	REAL(wp), DIMENSION(jpi,jpj) :: zdt_mle,zds_mle,zdu_mle,zdv_mle,zdb_mle ! difference between MLE layer average and parameter at base of blayer
301	! REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
302	REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
303	REAL(wp) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline
304	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline
305	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline
306	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline
307	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc ! u-shear across the pycnocline
308	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc ! v-shear across the pycnocline
309	REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient.
310	! Flux-gradient relationship variables
311	REAL(wp), DIMENSION(jpi, jpj) :: zshear ! Shear production.
312
313	REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale.
314
315	REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline.
316	REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
317	REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/
318	REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms.
319	REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
320
321	! For calculating Ri#-dependent mixing
322	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du ! u-shear^2
323	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv ! v-shear^2
324	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
325
326	! Temporary variables
327	INTEGER :: inhml
328	REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
329	REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables
330	REAL(wp) :: zthick, zz0, zz1 ! temporary variables
331	REAL(wp) :: zvel_max, zhbl_s ! temporary variables
332	REAL(wp) :: zfac, ztmp ! temporary variable
333	REAL(wp) :: zus_x, zus_y ! temporary Stokes drift
334	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
335	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
336	REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc
337	REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML.
338	REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc
339	REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max
340	INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme
341	REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau
342	REAL(wp) :: zzeta_s = 0._wp
343	REAL(wp) :: zzeta_v = 0.46
344	REAL(wp) :: zabsstke
345	REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness
346	REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc
347
348	! For debugging
349	INTEGER :: ikt
350	!!--------------------------------------------------------------------
351	!
352	ibld(:,:) = 0 ; imld(:,:) = 0
353	zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp
354	zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:) = 0._wp ; zuw0(:,:) = 0._wp
355	zvw0(:,:) = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:) = 0._wp ; zwb0(:,:) = 0._wp
356	zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:) = 0._wp ; zwb_ent(:,:) = 0._wp
357	zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
358	zhol(:,:) = 0._wp ; zwb0tot(:,:) = 0._wp
359	lconv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE.
360	! mixed layer
361	! no initialization of zhbl or zhml (or zdh?)
362	zhbl(:,:) = 1._wp ; zhml(:,:) = 1._wp ; zdh(:,:) = 1._wp ; zdhdt(:,:) = 0._wp
363	zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp
364	zv_bl(:,:) = 0._wp ; zb_bl(:,:) = 0._wp
365	zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp
366	zt_mle(:,:) = 0._wp ; zs_mle(:,:) = 0._wp ; zu_mle(:,:) = 0._wp
367	zb_mle(:,:) = 0._wp
368	zv_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp
369	zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp
370	zdt_ml(:,:) = 0._wp ; zds_ml(:,:) = 0._wp ; zdu_ml(:,:) = 0._wp ; zdv_ml(:,:) = 0._wp
371	zdb_ml(:,:) = 0._wp
372	zdt_mle(:,:) = 0._wp ; zds_mle(:,:) = 0._wp ; zdu_mle(:,:) = 0._wp
373	zdv_mle(:,:) = 0._wp ; zdb_mle(:,:) = 0._wp
374	zwth_ent = 0._wp ; zws_ent = 0._wp
375	!
376	zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
377	zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
378	!
379	zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp
380
381	IF ( ln_osm_mle ) THEN ! only initialise arrays if needed
382	zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp
383	zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp
384	zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp
385	zhmle(:,:) = 0._wp ; zmld(:,:) = 0._wp
386	ENDIF
387	zwb_fk_b(:,:) = 0._wp ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt
388
389	! Flux-Gradient arrays.
390	zsc_wth_1(:,:) = 0._wp ; zsc_ws_1(:,:) = 0._wp ; zsc_uw_1(:,:) = 0._wp
391	zsc_uw_2(:,:) = 0._wp ; zsc_vw_1(:,:) = 0._wp ; zsc_vw_2(:,:) = 0._wp
392	zhbl_t(:,:) = 0._wp ; zdhdt(:,:) = 0._wp
393
394	zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
395	ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp
396
397	! hbl = MAX(hbl,epsln)
398	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
399	! Calculate boundary layer scales
400	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
401
402	! Assume two-band radiation model for depth of OSBL
403	zz0 = rn_abs ! surface equi-partition in 2-bands
404	zz1 = 1. - rn_abs
405	DO jj = 2, jpjm1
406	DO ji = 2, jpim1
407	! Surface downward irradiance (so always +ve)
408	zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp
409	! Downwards irradiance at base of boundary layer
410	zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
411	! Downwards irradiance averaged over depth of the OSBL
412	zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
413	& + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
414	END DO
415	END DO
416	! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
417	DO jj = 2, jpjm1
418	DO ji = 2, jpim1
419	zthermal = rab_n(ji,jj,1,jp_tem)
420	zbeta = rab_n(ji,jj,1,jp_sal)
421	! Upwards surface Temperature flux for non-local term
422	zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
423	! Upwards surface salinity flux for non-local term
424	zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
425	! Non radiative upwards surface buoyancy flux
426	zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj)
427	! Total upwards surface buoyancy flux
428	zwb0tot(ji,jj) = zwb0(ji,jj) - grav * zthermal * zrad0(ji,jj)
429	! turbulent heat flux averaged over depth of OSBL
430	zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
431	! turbulent salinity flux averaged over depth of the OBSL
432	zwsav(ji,jj) = 0.5 * zws0(ji,jj)
433	! turbulent buoyancy flux averaged over the depth of the OBSBL
434	zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj)
435	! Surface upward velocity fluxes
436	zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rau0 * tmask(ji,jj,1)
437	zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rau0 * tmask(ji,jj,1)
438	! Friction velocity (zustar), at T-point : LMD94 eq. 2
439	zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
440	zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
441	zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
442	END DO
443	END DO
444	! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
445	SELECT CASE (nn_osm_wave)
446	! Assume constant La#=0.3
447	CASE(0)
448	DO jj = 2, jpjm1
449	DO ji = 2, jpim1
450	zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
451	zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
452	! Linearly
453	zustke(ji,jj) = MAX ( SQRT( zus_xzus_x + zus_yzus_y), 1.0e-8 )
454	dstokes(ji,jj) = rn_osm_dstokes
455	END DO
456	END DO
457	! Assume Pierson-Moskovitz wind-wave spectrum
458	CASE(1)
459	DO jj = 2, jpjm1
460	DO ji = 2, jpim1
461	! Use wind speed wndm included in sbc_oce module
462	zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
463	dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
464	END DO
465	END DO
466	! Use ECMWF wave fields as output from SBCWAVE
467	CASE(2)
468	zfac = 2.0_wp * rpi / 16.0_wp
469
470	DO jj = 2, jpjm1
471	DO ji = 2, jpim1
472	IF (hsw(ji,jj) > 1.e-4) THEN
473	! Use wave fields
474	zabsstke = SQRT(ut0sd(ji,jj)2 + vt0sd(ji,jj)2)
475	zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8)
476	dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)hsw(ji,jj) / ( MAX(zabsstke wmp(ji,jj), 1.0e-7 ) ), 5.0e-1)
477	ELSE
478	! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run)
479	! .. so default to Pierson-Moskowitz
480	zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
481	dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
482	END IF
483	END DO
484	END DO
485	END SELECT
486
487	IF (ln_zdfosm_ice_shelter) THEN
488	! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
489	DO jj = 2, jpjm1
490	DO ji = 2, jpim1
491	zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj))
492	dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj))
493	END DO
494	END DO
495	END IF
496
497	SELECT CASE (nn_osm_SD_reduce)
498	! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van Roekel (2012) or Grant (2020).
499	CASE(0)
500	! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas.
501	! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation.
502	! It could represent the effects of the spread of wave directions
503	! around the mean wind. The effect of this adjustment needs to be tested.
504	IF(nn_osm_wave > 0) THEN
505	zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1)
506	END IF
507	CASE(1)
508	! van Roekel (2012): consider average SD over top 10% of boundary layer
509	! assumes approximate depth profile of SD from Breivik (2016)
510	zsqrtpi = SQRT(rpi)
511	z_two_thirds = 2.0_wp / 3.0_wp
512
513	DO jj = 2, jpjm1
514	DO ji = 2, jpim1
515	zthickness = rn_osm_hblfrac*hbl(ji,jj)
516	z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
517	zsqrt_depth = SQRT(z2k_times_thickness)
518	zexp_depth = EXP(-z2k_times_thickness)
519	zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth &
520	& - z_two_thirds * ( zsqrtpizsqrt_depthz2k_times_thickness * ERFC(zsqrt_depth) &
521	& + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness
522
523	END DO
524	END DO
525	CASE(2)
526	! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer
527	! assumes approximate depth profile of SD from Breivik (2016)
528	zsqrtpi = SQRT(rpi)
529
530	DO jj = 2, jpjm1
531	DO ji = 2, jpim1
532	zthickness = rn_osm_hblfrac*hbl(ji,jj)
533	z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
534
535	IF(z2k_times_thickness < 50._wp) THEN
536	zsqrt_depth = SQRT(z2k_times_thickness)
537	zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness)
538	ELSE
539	! asymptotic expansion of sqrt(pi)zsqrt_depthEXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness
540	! See Abramowitz and Stegun, Eq. 7.1.23
541	! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness2) - (15/8)/(z2k_times_thickness3)
542	zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp
543	END IF
544	zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp)
545	dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj)
546	zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc)
547	END DO
548	END DO
549	END SELECT
550
551	! Langmuir velocity scale (zwstrl), La # (zla)
552	! mixed scale (zvstr), convective velocity scale (zwstrc)
553	DO jj = 2, jpjm1
554	DO ji = 2, jpim1
555	! Langmuir velocity scale (zwstrl), at T-point
556	zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
557	zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)
558	IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))
559	! Velocity scale that tends to zustar for large Langmuir numbers
560	zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + &
561	& ( 1.0 - EXP( -0.5 * zla(ji,jj)*2 ) ) zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
562
563	! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
564	! Note zustke and zwstrl are not amended.
565	!
566	! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
567	IF ( zwbav(ji,jj) > 0.0) THEN
568	zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
569	zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
570	ELSE
571	zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln )
572	ENDIF
573	END DO
574	END DO
575
576	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
577	! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
578	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
579	! BL must be always 4 levels deep.
580	! For calculation of lateral buoyancy gradients for FK in
581	! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must
582	! previously exist for hbl also.
583
584	! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway
585	! ##########################################################################
586	hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) )
587	ibld(:,:) = 4
588	DO jk = 5, jpkm1
589	DO jj = 1, jpj
590	DO ji = 1, jpi
591	IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
592	ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
593	ENDIF
594	END DO
595	END DO
596	END DO
597	! ##########################################################################
598
599	DO jj = 2, jpjm1
600	DO ji = 2, jpim1
601	zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
602	imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) )) , 1 ))
603	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
604	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
605	END DO
606	END DO
607	! Averages over well-mixed and boundary layer
608	jp_ext(:,:) = 2
609	CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl)
610	! jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1
611	CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
612	! Velocity components in frame aligned with surface stress.
613	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
614	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
615	! Determine the state of the OSBL, stable/unstable, shear/no shear
616	CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear )
617
618	IF ( ln_osm_mle ) THEN
619	! Fox-Kemper Scheme
620	mld_prof = 4
621	DO jk = 5, jpkm1
622	DO jj = 2, jpjm1
623	DO ji = 2, jpim1
624	IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
625	END DO
626	END DO
627	END DO
628	jp_ext_mle(:,:) = 2
629	CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle)
630
631	DO jj = 2, jpjm1
632	DO ji = 2, jpim1
633	zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
634	END DO
635	END DO
636
637	!! External gradient
638	CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
639	CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
640	CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext )
641	CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
642	CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
643	ELSE ! ln_osm_mle
644	! FK not selected, Boundary Layer only.
645	lpyc(:,:) = .TRUE.
646	lflux(:,:) = .FALSE.
647	lmle(:,:) = .FALSE.
648	DO jj = 2, jpjm1
649	DO ji = 2, jpim1
650	IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
651	END DO
652	END DO
653	ENDIF ! ln_osm_mle
654
655	! Test if pycnocline well resolved
656	DO jj = 2, jpjm1
657	DO ji = 2,jpim1
658	IF (lconv(ji,jj) ) THEN
659	ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj))
660	IF ( ztmp > 6 ) THEN
661	! pycnocline well resolved
662	jp_ext(ji,jj) = 1
663	ELSE
664	! pycnocline poorly resolved
665	jp_ext(ji,jj) = 0
666	ENDIF
667	ELSE
668	! Stable conditions
669	jp_ext(ji,jj) = 0
670	ENDIF
671	END DO
672	END DO
673
674	CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )
675	! jp_ext = ibld-imld+1
676	CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
677	! Rate of change of hbl
678	CALL zdf_osm_calculate_dhdt( zdhdt )
679	DO jj = 2, jpjm1
680	DO ji = 2, jpim1
681	zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it
682	! adjustment to represent limiting by ocean bottom
683	IF ( zhbl_t(ji,jj) >= gdepw_n(ji, jj, mbkt(ji,jj) + 1 ) ) THEN
684	zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:))
685	lpyc(ji,jj) = .FALSE.
686	ENDIF
687	END DO
688	END DO
689
690	imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index
691	ibld(:,:) = 4
692
693	DO jk = 4, jpkm1
694	DO jj = 2, jpjm1
695	DO ji = 2, jpim1
696	IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
697	ibld(ji,jj) = jk
698	ENDIF
699	END DO
700	END DO
701	END DO
702
703	!
704	! Step through model levels taking account of buoyancy change to determine the effect on dhdt
705	!
706	CALL zdf_osm_timestep_hbl( zdhdt )
707	! is external level in bounds?
708
709	CALL zdf_osm_vertical_average( ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )
710	!
711	!
712	! Check to see if lpyc needs to be changed
713
714	CALL zdf_osm_pycnocline_thickness( dh, zdh )
715
716	DO jj = 2, jpjm1
717	DO ji = 2, jpim1
718	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(ji,jj) == 1 ) lpyc(ji,jj) = .FALSE.
719	END DO
720	END DO
721
722	dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms.
723	!
724	! Average over the depth of the mixed layer in the convective boundary layer
725	! jp_ext = ibld - imld +1
726	CALL zdf_osm_vertical_average( imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )
727	! rotate mean currents and changes onto wind align co-ordinates
728	!
729	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
730	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
731	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
732	! Pycnocline gradients for scalars and velocity
733	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
734
735	CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
736	CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc )
737	CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )
738	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
739	! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
740	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
741	CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
742
743	!
744	! calculate non-gradient components of the flux-gradient relationships
745	!
746	! Stokes term in scalar flux, flux-gradient relationship
747	WHERE ( lconv )
748	zsc_wth_1 = zwstrl*3 zwth0 / ( zvstr*3 + 0.5 zwstrc**3 + epsln)
749	!
750	zsc_ws_1 = zwstrl*3 zws0 / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
751	ELSEWHERE
752	zsc_wth_1 = 2.0 * zwthav
753	!
754	zsc_ws_1 = 2.0 * zwsav
755	ENDWHERE
756
757
758	DO jj = 2, jpjm1
759	DO ji = 2, jpim1
760	IF ( lconv(ji,jj) ) THEN
761	DO jk = 2, imld(ji,jj)
762	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
763	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
764	!
765	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
766	END DO ! end jk loop
767	ELSE ! else for if (lconv)
768	! Stable conditions
769	DO jk = 2, ibld(ji,jj)
770	zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
771	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) &
772	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
773	!
774	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) &
775	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
776	END DO
777	ENDIF ! endif for check on lconv
778
779	END DO ! end of ji loop
780	END DO ! end of jj loop
781
782	! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero)
783	WHERE ( lconv )
784	zsc_uw_1 = ( zwstrl*3 + 0.5 zwstrc3 )pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 )
785	zsc_uw_2 = ( zwstrl*3 + 0.5 zwstrc3 )pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 )
786	zsc_vw_1 = ff_t * zhml * zustke*3 MIN( zla(8.0/3.0), 0.12 ) / ( ( zvstr3 + 0.5 * zwstrc3 )(2.0/3.0) + epsln )
787	ELSEWHERE
788	zsc_uw_1 = zustar**2
789	zsc_vw_1 = ff_t * zhbl * zustke*3 MIN( zla(8.0/3.0), 0.12 ) / (zvstr2 + epsln)
790	ENDWHERE
791	IF(ln_dia_osm) THEN
792	IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
793	IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
794	END IF
795	DO jj = 2, jpjm1
796	DO ji = 2, jpim1
797	IF ( lconv(ji,jj) ) THEN
798	DO jk = 2, imld(ji,jj)
799	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
800	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) &
801	& + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) &
802	& * ( 1.0 - EXP ( -2.0 * zznd_d ) )
803	!
804	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( - zznd_d ) &
805	& * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
806	END DO ! end jk loop
807	ELSE
808	! Stable conditions
809	DO jk = 2, ibld(ji,jj) ! corrected to ibld
810	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
811	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) &
812	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
813	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
814	END DO ! end jk loop
815	ENDIF
816	END DO ! ji loop
817	END DO ! jj loo
818
819	! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
820
821	WHERE ( lconv )
822	zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
823	zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) * zhml / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
824	ELSEWHERE
825	zsc_wth_1 = 0._wp
826	zsc_ws_1 = 0._wp
827	ENDWHERE
828
829	DO jj = 2, jpjm1
830	DO ji = 2, jpim1
831	IF (lconv(ji,jj) ) THEN
832	DO jk = 2, imld(ji,jj)
833	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
834	! calculate turbulent time scale
835	zl_c = 0.9 * ( 1.0 - EXP ( - 5.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) &
836	& * ( 1.0 - EXP ( -15.0 * ( 1.2 - zznd_ml ) ) )
837	zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) &
838	& * ( 1.0 - EXP ( - 8.0 * ( 1.15 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
839	zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0)
840	! non-gradient buoyancy terms
841	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.4 * zsc_wth_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml )
842	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.4 * zsc_ws_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml )
843	END DO
844
845	IF ( lpyc(ji,jj) ) THEN
846	ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)3 + zwstrc(ji,jj)3 )**pthird
847	ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 )
848	zwth_ent = -0.003 * ( 0.15 * zvstr(ji,jj)3 + zwstrc(ji,jj)3 )*pthird ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj)
849	zws_ent = -0.003 * ( 0.15 * zvstr(ji,jj)3 + zwstrc(ji,jj)3 )*pthird ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj)
850	! Cubic profile used for buoyancy term
851	DO jk = 2, ibld(ji,jj)
852	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
853	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)*2 ) MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc*2 - 0.2 zznd_pyc**3 ), 0.0 )
854
855	ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)*2 ) MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc*2 - 0.2 zznd_pyc**3 ), 0.0 )
856	END DO
857	!
858	IF ( dh(ji,jj) < 0.2*hbl(ji,jj) ) THEN
859	zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj)
860	zdelta_pyc = ( zvstr(ji,jj)3 + zwstrc(ji,jj)3 )pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)3 + zwstrc(ji,jj)3 )p2third / zdh(ji,jj)**2 ) )
861	!
862	zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
863	!
864	zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
865	!
866	zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
867	DO jk = 2, ibld(ji,jj)
868	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
869	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )*2 ) zdh(ji,jj) / ( zvstr(ji,jj)3 + zwstrc(ji,jj)3 )**pthird
870	!
871	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )*2 ) zdh(ji,jj) / ( zvstr(ji,jj)3 + zwstrc(ji,jj)3 )**pthird
872	END DO
873	END IF
874	ENDIF ! End of pycnocline
875	ELSE ! lconv test - stable conditions
876	DO jk = 2, ibld(ji,jj)
877	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
878	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj)
879	END DO
880	ENDIF
881	END DO ! ji loop
882	END DO ! jj loop
883
884	WHERE ( lconv )
885	zsc_uw_1 = -zwb0 * zustar*2 zhml / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
886	zsc_uw_2 = zwb0 * zustke * zhml / ( zvstr*3 + 0.5 zwstrc3 + epsln )(2.0/3.0)
887	zsc_vw_1 = 0._wp
888	ELSEWHERE
889	zsc_uw_1 = 0._wp
890	zsc_vw_1 = 0._wp
891	ENDWHERE
892
893	DO jj = 2, jpjm1
894	DO ji = 2, jpim1
895	IF ( lconv(ji,jj) ) THEN
896	DO jk = 2 , imld(ji,jj)
897	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
898	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) &
899	& * ( 1.0 - EXP( -0.5 * zznd_d ) ) &
900	& * zsc_uw_2(ji,jj) )
901	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
902	END DO ! jk loop
903	ELSE
904	! stable conditions
905	DO jk = 2, ibld(ji,jj)
906	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
907	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
908	END DO
909	ENDIF
910	END DO ! ji loop
911	END DO ! jj loop
912
913	DO jj = 2, jpjm1
914	DO ji = 2, jpim1
915	IF( lconv(ji,jj) ) THEN
916	IF ( lpyc(ji,jj) ) THEN
917	IF ( j_ddh(ji,jj) == 0 ) THEN
918	! Place holding code. Parametrization needs checking for these conditions.
919	zomega = ( 0.15 * zwstrl(ji,jj)3 + zwstrc(ji,jj)3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird
920	zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
921	zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
922	ELSE
923	zomega = ( 0.15 * zwstrl(ji,jj)3 + zwstrc(ji,jj)3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird
924	zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
925	zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
926	ENDIF
927	zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse
928	zc_cubic = zuw_bse - zd_cubic
929	! need ztau_sc_u to be available. Change to array.
930	DO jk = imld(ji,jj), ibld(ji,jj)
931	zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
932	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)*2 ) zuw_bse * &
933	& ( zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc*3 ) ( 0.75 + 0.25 * zznd_pyc )*2 zdbdz_pyc(ji,jj,jk)
934	END DO
935	zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) )
936	zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse
937	zc_cubic = zvw_bse - zd_cubic
938	DO jk = imld(ji,jj), ibld(ji,jj)
939	zznd_pyc = -( gdepw_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj)
940	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)*2 ) zvw_bse * &
941	& ( zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc*3 ) ( 0.75 + 0.25 * zznd_pyc )*2 zdbdz_pyc(ji,jj,jk)
942	END DO
943	ENDIF ! lpyc
944	ENDIF ! lconv
945	END DO ! ji loop
946	END DO ! jj loop
947
948	IF(ln_dia_osm) THEN
949	IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
950	IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
951	END IF
952	! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
953
954	DO jj = 1, jpjm1
955	DO ji = 1, jpim1
956
957	IF ( lconv(ji,jj) ) THEN
958	zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) )
959	zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) )
960	IF ( lpyc(ji,jj) ) THEN
961	! Pycnocline scales
962	zsc_wth_pyc(ji,jj) = -0.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj)
963	zsc_ws_pyc(ji,jj) = -0.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj)
964	ENDIF
965	ELSE
966	zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj)
967	zsc_ws_1(ji,jj) = zws0(ji,jj)
968	ENDIF
969	END DO
970	END DO
971
972	DO jj = 2, jpjm1
973	DO ji = 2, jpim1
974	IF ( lconv(ji,jj) ) THEN
975	DO jk = 2, imld(ji,jj)
976	zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
977	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) &
978	& * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) &
979	& - EXP( - 6.0 * zznd_ml ) ) ) &
980	& * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) )
981	!
982	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) &
983	& * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) &
984	& - EXP( - 6.0 * zznd_ml ) ) ) &
985	& * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
986	END DO
987	!
988	! may need to comment out lpyc block
989	IF ( lpyc(ji,jj) ) THEN
990	! pycnocline
991	DO jk = imld(ji,jj), ibld(ji,jj)
992	zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
993	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )
994	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )
995	END DO
996	ENDIF
997	ELSE
998	IF( zdhdt(ji,jj) > 0. ) THEN
999	DO jk = 2, ibld(ji,jj)
1000	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1001	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1002	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
1003	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
1004	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
1005	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
1006	END DO
1007	ENDIF
1008	ENDIF
1009	ENDDO ! ji loop
1010	END DO ! jj loop
1011
1012	WHERE ( lconv )
1013	zsc_uw_1 = zustar**2
1014	zsc_vw_1 = ff_t * zustke * zhml
1015	ELSEWHERE
1016	zsc_uw_1 = zustar**2
1017	zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
1018	zsc_vw_1 = ff_t * zustke * zhbl
1019	zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )*2 ) zsc_vw_1
1020	ENDWHERE
1021
1022	DO jj = 2, jpjm1
1023	DO ji = 2, jpim1
1024	IF ( lconv(ji,jj) ) THEN
1025	DO jk = 2, imld(ji,jj)
1026	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
1027	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1028	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
1029	& + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml*4 ) - EXP ( -8.0 zznd_ml ) ) * zsc_uw_1(ji,jj)
1030	!
1031	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1032	& + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
1033	END DO
1034
1035	ELSE
1036	DO jk = 2, ibld(ji,jj)
1037	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1038	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1039	IF ( zznd_d <= 2.0 ) THEN
1040	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
1041	&* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
1042	!
1043	ELSE
1044	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
1045	& + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
1046	!
1047	ENDIF
1048
1049	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1050	& + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d*2 ) zsc_vw_1(ji,jj)
1051	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1052	& + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
1053	END DO
1054	ENDIF
1055	END DO
1056	END DO
1057
1058	! DO jj = 1, jpjm1
1059	! DO ji = 1, jpim1
1060	! IF ( lconv(ji,jj) ) THEN
1061	! IF ( lpyc(ji,jj) ) THEN
1062	! zd_cubic = ( 0.948 - 2.13 * zdh(ji,jj) / zhml(ji,jj) ) * zustar(ji,jj)**2
1063	! zc_cubic = -0.474 * zustar(ji,jj)**2 - zd_cubic
1064	! DO jk = imld(ji,jj), ibld(ji,jj)
1065	! zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1066	! ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc**3 )
1067	! END DO
1068	! zc_cubic= 3.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj)
1069	! zd_cubic = -2.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj)
1070	! DO jk = imld(ji,jj), ibld(ji,jj)
1071	! zznd_pyc = -( gdepw_n(ji,jj,jk)-zhbl(ji,jj) ) / zdh(ji,jj)
1072	! ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc**3 )
1073	! END DO
1074	! ENDIF
1075	! ENDIF
1076	! END DO
1077	! END DO
1078
1079	IF(ln_dSia_osm) THEN
1080	IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
1081	IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
1082	IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
1083	IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
1084	IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
1085	IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
1086	END IF
1087	!
1088	! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
1089
1090
1091	! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer.
1092
1093	DO jj = 2, jpjm1
1094	DO ji = 2, jpim1
1095	IF ( .not. lconv(ji,jj) ) THEN
1096	DO jk = 2, ibld(ji,jj)
1097	znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
1098	IF ( znd >= 0.0 ) THEN
1099	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
1100	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
1101	ELSE
1102	ghamu(ji,jj,jk) = 0._wp
1103	ghamv(ji,jj,jk) = 0._wp
1104	ENDIF
1105	END DO
1106	ENDIF
1107	END DO
1108	END DO
1109
1110	! pynocline contributions
1111	DO jj = 2, jpjm1
1112	DO ji = 2, jpim1
1113	IF ( .not. lconv(ji,jj) ) THEN
1114	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
1115	DO jk= 2, ibld(ji,jj)
1116	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
1117	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
1118	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
1119	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
1120	END DO
1121	END IF
1122	END IF
1123	END DO
1124	END DO
1125	IF(ln_dia_osm) THEN
1126	IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
1127	IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
1128	END IF
1129
1130	DO jj=2, jpjm1
1131	DO ji = 2, jpim1
1132	ghamt(ji,jj,ibld(ji,jj)) = 0._wp
1133	ghams(ji,jj,ibld(ji,jj)) = 0._wp
1134	ghamu(ji,jj,ibld(ji,jj)) = 0._wp
1135	ghamv(ji,jj,ibld(ji,jj)) = 0._wp
1136	END DO ! ji loop
1137	END DO ! jj loop
1138
1139	IF(ln_dia_osm) THEN
1140	IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
1141	IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
1142	IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
1143	IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
1144	IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
1145	END IF
1146	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1147	! Need to put in code for contributions that are applied explicitly to
1148	! the prognostic variables
1149	! 1. Entrainment flux
1150	!
1151	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1152
1153
1154
1155	! rotate non-gradient velocity terms back to model reference frame
1156
1157	DO jj = 2, jpjm1
1158	DO ji = 2, jpim1
1159	DO jk = 2, ibld(ji,jj)
1160	ztemp = ghamu(ji,jj,jk)
1161	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
1162	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
1163	END DO
1164	END DO
1165	END DO
1166
1167	IF(ln_dia_osm) THEN
1168	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1169	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
1170	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
1171	END IF
1172
1173	! KPP-style Ri# mixing
1174	IF( ln_kpprimix) THEN
1175	DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form)
1176	DO jj = 1, jpjm1
1177	DO ji = 1, jpim1 ! vector opt.
1178	z3du(ji,jj,jk) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) &
1179	& * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) &
1180	& / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
1181	z3dv(ji,jj,jk) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) &
1182	& * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) &
1183	& / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
1184	END DO
1185	END DO
1186	END DO
1187	!
1188	DO jk = 2, jpkm1
1189	DO jj = 2, jpjm1
1190	DO ji = 2, jpim1 ! vector opt.
1191	! ! shear prod. at w-point weightened by mask
1192	zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
1193	& + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
1194	! ! local Richardson number
1195	zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
1196	zfri = MIN( zri / rn_riinfty , 1.0_wp )
1197	zfri = ( 1.0_wp - zfri * zfri )
1198	zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk)
1199	END DO
1200	END DO
1201	END DO
1202
1203	DO jj = 2, jpjm1
1204	DO ji = 2, jpim1
1205	DO jk = ibld(ji,jj) + 1, jpkm1
1206	zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
1207	zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
1208	END DO
1209	END DO
1210	END DO
1211
1212	END IF ! ln_kpprimix = .true.
1213
1214	! KPP-style set diffusivity large if unstable below BL
1215	IF( ln_convmix) THEN
1216	DO jj = 2, jpjm1
1217	DO ji = 2, jpim1
1218	DO jk = ibld(ji,jj) + 1, jpkm1
1219	IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
1220	END DO
1221	END DO
1222	END DO
1223	END IF ! ln_convmix = .true.
1224
1225
1226
1227	IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing
1228	DO jj = 2 , jpjm1
1229	DO ji = 2, jpim1
1230	IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer
1231	! Calculate MLE flux contribution from surface fluxes
1232	DO jk = 1, ibld(ji,jj)
1233	znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln)
1234	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( zwth0(ji,jj) - zrad0(ji,jj) ) * ( 1.0 - znd )
1235	ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd )
1236	END DO
1237	DO jk = 1, mld_prof(ji,jj)
1238	znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1239	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( zwth0(ji,jj) - zrad0(ji,jj) ) * ( 1.0 - znd )
1240	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd )
1241	END DO
1242	! Viscosity for MLEs
1243	DO jk = 1, mld_prof(ji,jj)
1244	znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1245	zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )*2 ) ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 )
1246	END DO
1247	ELSE
1248	! Surface transports limited to OSBL.
1249	! Viscosity for MLEs
1250	DO jk = 1, mld_prof(ji,jj)
1251	znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1252	zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )*2 ) ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 )
1253	END DO
1254	ENDIF
1255	END DO
1256	END DO
1257	ENDIF
1258
1259	IF(ln_dia_osm) THEN
1260	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1261	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
1262	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
1263	END IF
1264
1265
1266	! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
1267	!CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
1268
1269	! GN 25/8: need to change tmask --> wmask
1270
1271	DO jk = 2, jpkm1
1272	DO jj = 2, jpjm1
1273	DO ji = 2, jpim1
1274	p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1275	p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
1276	END DO
1277	END DO
1278	END DO
1279	! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids
1280	CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1., &
1281	& ghamu, 'W', 1. , ghamv, 'W', 1. )
1282	DO jk = 2, jpkm1
1283	DO jj = 2, jpjm1
1284	DO ji = 2, jpim1
1285	ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
1286	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1287
1288	ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
1289	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1290
1291	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk)
1292	ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk)
1293	END DO
1294	END DO
1295	END DO
1296	! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)
1297	CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
1298	! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged)
1299	! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged)
1300	CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., &
1301	& ghamu, 'U', -1. , ghamv, 'V', -1. )
1302
1303	IF(ln_dia_osm) THEN
1304	SELECT CASE (nn_osm_wave)
1305	! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1).
1306	CASE(0:1)
1307	IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)zustkezcos_wind ) ! x surface Stokes drift
1308	IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)zustkezsin_wind ) ! y surface Stokes drift
1309	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2zustke )
1310	! Stokes drift read in from sbcwave (=2).
1311	CASE(2:3)
1312	IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift
1313	IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift
1314	IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period
1315	IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height
1316	IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.rpi1.026/(0.877grav) )wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum
1317	IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)wndm2tmask(:,:,1) ) ! significant wave height from NP spectrum
1318	IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10
1319	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2 &
1320	& SQRT(ut0sd2 + vt0sd2 ) )
1321	END SELECT
1322	IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL>
1323	IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL>
1324	IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL>
1325	IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL>
1326	IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! <Tw_0>
1327	IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0>
1328	IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*zwb0 ) ! <Sw_0>
1329	IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb buoyancy flux
1330	IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth
1331	IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k
1332	IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base
1333	IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base
1334	IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base
1335	IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base
1336	IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base
1337	IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth
1338	IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth
1339	IF ( iom_use("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*zdt_ml ) ! dt at ml base
1340	IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*zds_ml ) ! ds at ml base
1341	IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml ) ! db at ml base
1342	IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth
1343	IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points
1344	IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc ) ! convective velocity scale
1345	IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale
1346	IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale
1347	IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale
1348	IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir #
1349	IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.rau0tmask(:,:,1)zustar*3 ) ! BL depth internal to zdf_osm routine
1350	IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.rau0tmask(:,:,1)zustar2zustke )
1351	IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine
1352	IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine
1353	IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine
1354	IF ( iom_use("jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext ) ! =1 if pycnocline resolved internal to zdf_osm routine
1355	IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*j_ddh ) ! index forpyc thicknessh internal to zdf_osm routine
1356	IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear ) ! shear production of TKE internal to zdf_osm routine
1357	IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine
1358	IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine
1359	IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux
1360	IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux
1361	IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux
1362	IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML
1363
1364	IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth
1365	IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth
1366	IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux
1367	IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML
1368	IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k
1369	IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt
1370	IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt
1371	IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt
1372	IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt
1373	IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt
1374	IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt
1375	IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1376	IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1377
1378	END IF
1379
1380	CONTAINS
1381	! subroutine code changed, needs syntax checking.
1382	SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
1383
1384	!!---------------------------------------------------------------------
1385	!! * ROUTINE zdf_osm_diffusivity_viscosity *
1386	!!
1387	!! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline.
1388	!!
1389	!! ** Method :
1390	!!
1391	!! !!----------------------------------------------------------------------
1392	REAL(wp), DIMENSION(:,:,:) :: zdiffut
1393	REAL(wp), DIMENSION(:,:,:) :: zviscos
1394	! local
1395
1396	! Scales used to calculate eddy diffusivity and viscosity profiles
1397	REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc
1398	REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr
1399	REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr
1400	REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc
1401	!
1402	REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac
1403
1404	REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375
1405	REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142
1406	REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15
1407
1408	DO jj = 2, jpjm1
1409	DO ji = 2, jpim1
1410	IF ( lconv(ji,jj) ) THEN
1411
1412	zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)3 + zwstrc(ji,jj)3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird
1413	zvel_sc_ml = ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
1414	zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )1.25 ) )2
1415
1416	zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml
1417	zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj)
1418
1419	IF ( lpyc(ji,jj) ) THEN
1420	zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj)
1421
1422	IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN
1423	zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )*pthird zhbl(ji,jj)
1424	ENDIF
1425
1426	zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) )
1427	zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac
1428	zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) )
1429
1430	zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )*2 ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )*2 + 0.0075 ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj)
1431	zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac
1432	IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN
1433	zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )*pthird zhbl(ji,jj)
1434	ENDIF
1435
1436	zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) )
1437	zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac
1438	zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_n_sc(ji,jj) )
1439
1440	zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third
1441	zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln )
1442	ELSE
1443	zbeta_d_sc(ji,jj) = 1.0
1444	zbeta_v_sc(ji,jj) = 1.0
1445	ENDIF
1446	ELSE
1447	zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
1448	zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
1449	END IF
1450	END DO
1451	END DO
1452	!
1453	DO jj = 2, jpjm1
1454	DO ji = 2, jpim1
1455	IF ( lconv(ji,jj) ) THEN
1456	DO jk = 2, imld(ji,jj) ! mixed layer diffusivity
1457	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
1458	!
1459	zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
1460	!
1461	zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) &
1462	& * ( 1.0 - 0.5 * zznd_ml**2 )
1463	END DO
1464	! pycnocline
1465	IF ( lpyc(ji,jj) ) THEN
1466	! Diffusivity profile in the pycnocline given by cubic polynomial.
1467	za_cubic = 0.5
1468	zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj)
1469	zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) &
1470	& - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8)
1471	zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic )
1472	zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
1473	DO jk = imld(ji,jj) , ibld(ji,jj)
1474	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
1475	!
1476	zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc**3 )
1477
1478	zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc*2 - 0.2 zznd_pyc**3 )
1479	END DO
1480	! viscosity profiles.
1481	za_cubic = 0.5
1482	zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)
1483	zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj) ) / MAX(zvispyc_n_sc(ji,jj), 1.e-8)
1484	zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic )
1485	zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
1486	DO jk = imld(ji,jj) , ibld(ji,jj)
1487	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
1488	zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc**3 )
1489	zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc*2 -0.2 zznd_pyc**3 )
1490	END DO
1491	IF ( zdhdt(ji,jj) > 0._wp ) THEN
1492	zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 )
1493	zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 )
1494	ELSE
1495	zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
1496	zviscos(ji,jj,ibld(ji,jj)) = 0._wp
1497	ENDIF
1498	ENDIF
1499	ELSE
1500	! stable conditions
1501	DO jk = 2, ibld(ji,jj)
1502	zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1503	zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
1504	zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
1505	END DO
1506
1507	IF ( zdhdt(ji,jj) > 0._wp ) THEN
1508	zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj))
1509	zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj))
1510	ENDIF
1511	ENDIF ! end if ( lconv )
1512	!
1513	END DO ! end of ji loop
1514	END DO ! end of jj loop
1515
1516	END SUBROUTINE zdf_osm_diffusivity_viscosity
1517
1518	SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear )
1519
1520	!!---------------------------------------------------------------------
1521	!! * ROUTINE zdf_osm_osbl_state *
1522	!!
1523	!! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number
1524	!!
1525	!! ** Method :
1526	!!
1527	!! !!----------------------------------------------------------------------
1528
1529	INTEGER, DIMENSION(jpi,jpj) :: j_ddh ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low.
1530
1531	LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear
1532
1533	REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer.
1534	REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline
1535
1536	! Local Variables
1537
1538	INTEGER :: jj, ji
1539
1540	REAL(wp), DIMENSION(jpi,jpj) :: zekman
1541	REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b ! Richardson numbers
1542	REAL(wp) :: zshear_u, zshear_v, zwb_shr
1543	REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
1544
1545	REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.8
1546	REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03
1547	REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15
1548	REAL, PARAMETER :: rn_ri_p_thresh = 27.0
1549	REAL, PARAMETER :: zri_c = 0.25
1550	REAL, PARAMETER :: zek = 4.0
1551	REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress.
1552
1553	! Determins stability and set flag lconv
1554	DO jj = 2, jpjm1
1555	DO ji = 2, jpim1
1556	IF ( zhol(ji,jj) < 0._wp ) THEN
1557	lconv(ji,jj) = .TRUE.
1558	ELSE
1559	lconv(ji,jj) = .FALSE.
1560	ENDIF
1561	END DO
1562	END DO
1563
1564	zekman(:,:) = EXP( - zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) )
1565
1566	zshear(:,:) = 0._wp
1567	j_ddh(:,:) = 1
1568
1569	DO jj = 2, jpjm1
1570	DO ji = 2, jpim1
1571	IF ( lconv(ji,jj) ) THEN
1572	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1573	zri_p(ji,jj) = MAX ( SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)2 + zdv_bl(ji,jj)2, 1.e-8) ) * ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 &
1574	& / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp )
1575
1576	IF ( ff_t(ji,jj) >= 0._wp ) THEN
1577	! Northern Hemisphere
1578	zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1.e-5 )2 + MAX( -zdv_ml(ji,jj), 1.e-5)2 )
1579	ELSE
1580	! Southern Hemisphere
1581	zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1.e-5 )2 + MAX( zdv_ml(ji,jj), 1.e-5)2 )
1582	ENDIF
1583	zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)*2 zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) )
1584	! Stability Dependence
1585	zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75 * MAX(0._wp,( zri_b(ji,jj) - zri_c ) / zri_c ) )
1586	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1587	! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when !
1588	! full code available !
1589	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1590	IF ( zshear(ji,jj) > 1.e-10 ) THEN
1591	IF ( zri_p(ji,jj) < rn_ri_p_thresh ) THEN
1592	! Growing shear layer
1593	j_ddh(ji,jj) = 0
1594	lshear(ji,jj) = .TRUE.
1595	ELSE
1596	j_ddh(ji,jj) = 1
1597	! IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN
1598	! shear production large enough to determine layer charcteristics, but can't maintain a shear layer.
1599	lshear(ji,jj) = .TRUE.
1600	! ELSE
1601	ENDIF
1602	ELSE
1603	j_ddh(ji,jj) = 2
1604	lshear(ji,jj) = .FALSE.
1605	ENDIF
1606	! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline.
1607	! zshear(ji,jj) = 0.5 * zshear(ji,jj)
1608	! lshear(ji,jj) = .FALSE.
1609	! ENDIF
1610	ELSE ! zdb_bl test, note zshear set to zero
1611	j_ddh(ji,jj) = 2
1612	lshear(ji,jj) = .FALSE.
1613	ENDIF
1614	ENDIF
1615	END DO
1616	END DO
1617
1618	! Calculate entrainment buoyancy flux due to surface fluxes.
1619
1620	DO jj = 2, jpjm1
1621	DO ji = 2, jpim1
1622	IF ( lconv(ji,jj) ) THEN
1623	zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
1624	zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
1625	zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
1626	zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
1627	IF (nn_osm_SD_reduce > 0 ) THEN
1628	! Effective Stokes drift already reduced from surface value
1629	zr_stokes = 1.0_wp
1630	ELSE
1631	! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
1632	! requires further reduction where BL is deep
1633	zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
1634	& * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
1635	END IF
1636	zwb_ent(ji,jj) = - 2.0 * zalpha_c * zrf_conv * zwbav(ji,jj) &
1637	& - zalpha_s * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
1638	& + zr_stokes * ( zalpha_s * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
1639	& - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1640	!
1641	ENDIF
1642	END DO ! ji loop
1643	END DO ! jj loop
1644
1645	zwb_min(:,:) = 0._wp
1646
1647	DO jj = 2, jpjm1
1648	DO ji = 2, jpim1
1649	IF ( lshear(ji,jj) ) THEN
1650	IF ( lconv(ji,jj) ) THEN
1651	! Unstable OSBL
1652	zwb_shr = -za_wb_s * zri_b(ji,jj) * zshear(ji,jj)
1653	IF ( j_ddh(ji,jj) == 0 ) THEN
1654
1655	! ! Developing shear layer, additional shear production possible.
1656
1657	! zshear_u = MAX( zustar(ji,jj)*2 MAX( zdu_ml(ji,jj), 0._wp ) / zhbl(ji,jj), 0._wp )
1658	! zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 )
1659	! zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )
1660
1661	! zwb_shr = zwb_shr - 0.25 * MAX ( zshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 )
1662	! zwb_shr = MAX( zwb_shr, -0.25 * zshear_u )
1663
1664	ENDIF
1665	zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
1666	! zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
1667	ELSE ! IF ( lconv ) THEN - ENDIF
1668	! Stable OSBL - shear production not coded for first attempt.
1669	ENDIF ! lconv
1670	ENDIF ! lshear
1671	IF ( lconv(ji,jj) ) THEN
1672	! Unstable OSBL
1673	zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0tot(ji,jj)
1674	ENDIF ! lconv
1675	END DO ! ji
1676	END DO ! jj
1677	END SUBROUTINE zdf_osm_osbl_state
1678
1679
1680	SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv )
1681	!!---------------------------------------------------------------------
1682	!! * ROUTINE zdf_vertical_average *
1683	!!
1684	!! ** Purpose : Determines vertical averages from surface to jnlev.
1685	!!
1686	!! ** Method : Averages are calculated from the surface to jnlev.
1687	!! The external level used to calculate differences is ibld+ibld_ext
1688	!!
1689	!!----------------------------------------------------------------------
1690
1691	INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over.
1692	INTEGER, DIMENSION(jpi,jpj) :: jp_ext
1693
1694	! Alan: do we need zb?
1695	REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb ! Average temperature and salinity
1696	REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components
1697	REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
1698	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components.
1699
1700	INTEGER :: jk, ji, jj, ibld_ext
1701	REAL(wp) :: zthick, zthermal, zbeta
1702
1703
1704	zt = 0._wp
1705	zs = 0._wp
1706	zu = 0._wp
1707	zv = 0._wp
1708	DO jj = 2, jpjm1 ! Vertical slab
1709	DO ji = 2, jpim1
1710	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1711	zbeta = rab_n(ji,jj,1,jp_sal)
1712	! average over depth of boundary layer
1713	zthick = epsln
1714	DO jk = 2, jnlev_av(ji,jj)
1715	zthick = zthick + e3t_n(ji,jj,jk)
1716	zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
1717	zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
1718	zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) &
1719	& * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
1720	& / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
1721	zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) &
1722	& * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
1723	& / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
1724	END DO
1725	zt(ji,jj) = zt(ji,jj) / zthick
1726	zs(ji,jj) = zs(ji,jj) / zthick
1727	zu(ji,jj) = zu(ji,jj) / zthick
1728	zv(ji,jj) = zv(ji,jj) / zthick
1729	zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj)
1730	ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj)
1731	IF ( ibld_ext < mbkt(ji,jj) ) THEN
1732	zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem)
1733	zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal)
1734	zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) &
1735	& / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) )
1736	zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) &
1737	& / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) )
1738	zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
1739	ELSE
1740	zdt(ji,jj) = 0._wp
1741	zds(ji,jj) = 0._wp
1742	zdu(ji,jj) = 0._wp
1743	zdv(ji,jj) = 0._wp
1744	zdb(ji,jj) = 0._wp
1745	ENDIF
1746	END DO
1747	END DO
1748	END SUBROUTINE zdf_osm_vertical_average
1749
1750	SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
1751	!!---------------------------------------------------------------------
1752	!! * ROUTINE zdf_velocity_rotation *
1753	!!
1754	!! ** Purpose : Rotates frame of reference of averaged velocity components.
1755	!!
1756	!! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w
1757	!!
1758	!!----------------------------------------------------------------------
1759
1760	REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle
1761	REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current
1762	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline
1763
1764	INTEGER :: ji, jj
1765	REAL(wp) :: ztemp
1766
1767	DO jj = 2, jpjm1
1768	DO ji = 2, jpim1
1769	ztemp = zu(ji,jj)
1770	zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
1771	zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
1772	ztemp = zdu(ji,jj)
1773	zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
1774	zdv(ji,jj) = zdv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
1775	END DO
1776	END DO
1777	END SUBROUTINE zdf_osm_velocity_rotation
1778
1779	SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
1780	!!---------------------------------------------------------------------
1781	!! * ROUTINE zdf_osm_osbl_state_fk *
1782	!!
1783	!! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme.
1784	!! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined
1785	!! lflux :: determines whether effects of surface flux extend below the base of the OSBL
1786	!! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl.
1787	!!
1788	!! ** Method :
1789	!!
1790	!!
1791	!!----------------------------------------------------------------------
1792
1793	! Outputs
1794	LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle
1795	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk
1796	!
1797	REAL(wp), DIMENSION(jpi,jpj) :: znd_param
1798	REAL(wp) :: zbuoy, ztmp, zpe_mle_layer
1799	REAL(wp) :: zpe_mle_ref, zdbdz_mle_int
1800
1801	znd_param(:,:) = 0._wp
1802
1803	DO jj = 2, jpjm1
1804	DO ji = 2, jpim1
1805	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
1806	zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj)
1807	END DO
1808	END DO
1809	DO jj = 2, jpjm1
1810	DO ji = 2, jpim1
1811	!
1812	IF ( lconv(ji,jj) ) THEN
1813	IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
1814	zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
1815	zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
1816	zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
1817	zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
1818	! Calculate potential energies of actual profile and reference profile.
1819	zpe_mle_layer = 0._wp
1820	zpe_mle_ref = 0._wp
1821	zthermal = rab_n(ji,jj,1,jp_tem)
1822	zbeta = rab_n(ji,jj,1,jp_sal)
1823
1824	DO jk = ibld(ji,jj), mld_prof(ji,jj)
1825	zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) )
1826	zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
1827	zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) ) * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
1828	END DO
1829	! Non-dimensional parameter to diagnose the presence of thermocline
1830
1831	znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) )
1832	ENDIF
1833	ENDIF
1834	END DO
1835	END DO
1836
1837	! Diagnosis
1838	DO jj = 2, jpjm1
1839	DO ji = 2, jpim1
1840	IF ( lconv(ji,jj) ) THEN
1841	IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN
1842	IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
1843	! MLE layer growing
1844	IF ( znd_param (ji,jj) > 100. ) THEN
1845	! Thermocline present
1846	lflux(ji,jj) = .FALSE.
1847	lmle(ji,jj) =.FALSE.
1848	ELSE
1849	! Thermocline not present
1850	lflux(ji,jj) = .TRUE.
1851	lmle(ji,jj) = .TRUE.
1852	ENDIF ! znd_param > 100
1853	!
1854	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
1855	lpyc(ji,jj) = .FALSE.
1856	ELSE
1857	lpyc(ji,jj) = .TRUE.
1858	ENDIF
1859	ELSE
1860	! MLE layer restricted to OSBL or just below.
1861	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
1862	! Weak stratification MLE layer can grow.
1863	lpyc(ji,jj) = .FALSE.
1864	lflux(ji,jj) = .TRUE.
1865	lmle(ji,jj) = .TRUE.
1866	ELSE
1867	! Strong stratification
1868	lpyc(ji,jj) = .TRUE.
1869	lflux(ji,jj) = .FALSE.
1870	lmle(ji,jj) = .FALSE.
1871	ENDIF ! zdb_bl < rn_mle_thresh_bl and
1872	ENDIF ! zhmle > 1.2 zhbl
1873	ELSE
1874	lpyc(ji,jj) = .TRUE.
1875	lflux(ji,jj) = .FALSE.
1876	lmle(ji,jj) = .FALSE.
1877	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
1878	ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5
1879	ELSE
1880	! Stable Boundary Layer
1881	lpyc(ji,jj) = .FALSE.
1882	lflux(ji,jj) = .FALSE.
1883	lmle(ji,jj) = .FALSE.
1884	ENDIF ! lconv
1885	END DO
1886	END DO
1887	END SUBROUTINE zdf_osm_osbl_state_fk
1888
1889	SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz )
1890	!!---------------------------------------------------------------------
1891	!! * ROUTINE zdf_osm_external_gradients *
1892	!!
1893	!! ** Purpose : Calculates the gradients below the OSBL
1894	!!
1895	!! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient.
1896	!!
1897	!!----------------------------------------------------------------------
1898
1899	INTEGER, DIMENSION(jpi,jpj) :: jbase
1900	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy.
1901
1902	INTEGER :: jj, ji, jkb, jkb1
1903	REAL(wp) :: zthermal, zbeta
1904
1905
1906	DO jj = 2, jpjm1
1907	DO ji = 2, jpim1
1908	IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
1909	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1910	zbeta = rab_n(ji,jj,1,jp_sal)
1911	jkb = jbase(ji,jj)
1912	jkb1 = MIN(jkb + 1, mbkt(ji,jj))
1913	zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
1914	& / e3t_n(ji,jj,jkb)
1915	zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
1916	& / e3t_n(ji,jj,jkb)
1917	zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
1918	ELSE
1919	zdtdz(ji,jj) = 0._wp
1920	zdsdz(ji,jj) = 0._wp
1921	zdbdz(ji,jj) = 0._wp
1922	END IF
1923	END DO
1924	END DO
1925	END SUBROUTINE zdf_osm_external_gradients
1926
1927	SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha )
1928
1929	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline
1930	REAL(wp), DIMENSION(jpi,jpj) :: zalpha
1931
1932	INTEGER :: jk, jj, ji
1933	REAL(wp) :: ztgrad, zsgrad, zbgrad
1934	REAL(wp) :: zgamma_b_nd, znd
1935	REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc
1936	REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
1937
1938	DO jj = 2, jpjm1
1939	DO ji = 2, jpim1
1940	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
1941	IF ( lconv(ji,jj) ) THEN ! convective conditions
1942	IF ( lpyc(ji,jj) ) THEN
1943	zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
1944	zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) )
1945	zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp )
1946
1947	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
1948	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1949	! Commented lines in this section are not needed in new code, once tested !
1950	! can be removed !
1951	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1952	! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
1953	! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
1954	zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)
1955	zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
1956	DO jk = 2, ibld(ji,jj)
1957	znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp
1958	IF ( znd <= zzeta_m ) THEN
1959	! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &
1960	! & EXP( -6.0 * ( znd -zzeta_m )**2 )
1961	! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &
1962	! & EXP( -6.0 * ( znd -zzeta_m )**2 )
1963	zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &
1964	& EXP( -6.0 * ( znd -zzeta_m )**2 )
1965	ELSE
1966	! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
1967	! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
1968	zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
1969	ENDIF
1970	END DO
1971	ENDIF ! if no pycnocline pycnocline gradients set to zero
1972	ELSE
1973	! stable conditions
1974	! if pycnocline profile only defined when depth steady of increasing.
1975	IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady.
1976	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1977	IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline
1978	ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
1979	ztgrad = zdt_bl(ji,jj) * ztmp
1980	zsgrad = zds_bl(ji,jj) * ztmp
1981	zbgrad = zdb_bl(ji,jj) * ztmp
1982	DO jk = 2, ibld(ji,jj)
1983	znd = gdepw_n(ji,jj,jk) * ztmp
1984	zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1985	zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1986	zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1987	END DO
1988	ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
1989	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
1990	ztgrad = zdt_bl(ji,jj) * ztmp
1991	zsgrad = zds_bl(ji,jj) * ztmp
1992	zbgrad = zdb_bl(ji,jj) * ztmp
1993	DO jk = 2, ibld(ji,jj)
1994	znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp
1995	zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1996	zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1997	zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1998	END DO
1999	ENDIF ! IF (zhol >=0.5)
2000	ENDIF ! IF (zdb_bl> 0.)
2001	ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
2002	ENDIF ! IF (lconv)
2003	ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) )
2004	END DO
2005	END DO
2006
2007	END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
2008
2009	SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
2010	!!---------------------------------------------------------------------
2011	!! * ROUTINE zdf_osm_pycnocline_shear_profiles *
2012	!!
2013	!! ** Purpose : Calculates velocity shear in the pycnocline
2014	!!
2015	!! ** Method :
2016	!!
2017	!!----------------------------------------------------------------------
2018
2019	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
2020
2021	INTEGER :: jk, jj, ji
2022	REAL(wp) :: zugrad, zvgrad, znd
2023	REAL(wp) :: zzeta_v = 0.45
2024	!
2025	DO jj = 2, jpjm1
2026	DO ji = 2, jpim1
2027	!
2028	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
2029	IF ( lconv (ji,jj) ) THEN
2030	! Unstable conditions. Shouldn;t be needed with no pycnocline code.
2031	! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
2032	! & ( ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zhml(ji,jj) ) * &
2033	! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
2034	!Alan is this right?
2035	! zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
2036	! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
2037	! & ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird + epsln ) &
2038	! & )/ (zdh(ji,jj) + epsln )
2039	! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
2040	! znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
2041	! IF ( znd <= 0.0 ) THEN
2042	! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
2043	! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
2044	! ELSE
2045	! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
2046	! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
2047	! ENDIF
2048	! END DO
2049	ELSE
2050	! stable conditions
2051	zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
2052	zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
2053	DO jk = 2, ibld(ji,jj)
2054	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
2055	IF ( znd < 1.0 ) THEN
2056	zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
2057	ELSE
2058	zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
2059	ENDIF
2060	zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
2061	END DO
2062	ENDIF
2063	!
2064	END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
2065	END DO
2066	END DO
2067	END SUBROUTINE zdf_osm_pycnocline_shear_profiles
2068
2069	SUBROUTINE zdf_osm_calculate_dhdt( zdhdt )
2070	!!---------------------------------------------------------------------
2071	!! * ROUTINE zdf_osm_calculate_dhdt *
2072	!!
2073	!! ** Purpose : Calculates the rate at which hbl changes.
2074	!!
2075	!! ** Method :
2076	!!
2077	!!----------------------------------------------------------------------
2078
2079	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! Rate of change of hbl
2080
2081	INTEGER :: jj, ji
2082	REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi
2083	REAL(wp) :: zvel_max,zddhdt
2084	REAL(wp), PARAMETER :: zzeta_m = 0.3
2085	REAL(wp), PARAMETER :: zgamma_c = 2.0
2086	REAL(wp), PARAMETER :: zdhoh = 0.1
2087	REAL(wp), PARAMETER :: zalpha_b = 0.3
2088	REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
2089
2090	DO jj = 2, jpjm1
2091	DO ji = 2, jpim1
2092
2093	IF ( lshear(ji,jj) ) THEN
2094	IF ( lconv(ji,jj) ) THEN ! Convective
2095
2096	IF ( ln_osm_mle ) THEN
2097
2098	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
2099	! Fox-Kemper buoyancy flux average over OSBL
2100	zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * &
2101	(1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
2102	ELSE
2103	zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
2104	ENDIF
2105	zvel_max = ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
2106	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
2107	! OSBL is deepening, entrainment > restratification
2108	IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
2109	zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
2110	zpsi = -zalpha_pyc(ji,jj) * ( zwb0tot(ji,jj) - ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) ) * zdh(ji,jj) / zhbl(ji,jj)
2111	zpsi = zpsi - 4.0 * ( 1.0 + zdh(ji,jj) /zhbl(ji,jj) ) * zgamma_b_nd * ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) )
2112	zpsi = zalpha_b * MAX ( zpsi, 0._wp )
2113	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) + zpsi / MAX( zdb_bl(ji,jj), 1.e-15 )
2114	IF ( j_ddh(ji,jj) == 1 ) THEN
2115	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
2116	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2117	ELSE
2118	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2119	ENDIF
2120	! Relaxation to dh_ref = zari * hbl
2121	zddhdt = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
2122
2123	ELSE IF ( j_ddh(ji,jj) == 0 ) THEN
2124	! Growing shear layer
2125	zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
2126	zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
2127	ELSE
2128	zddhdt = 0._wp
2129	ENDIF ! j_ddh
2130	zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( zalpha_pyc(ji,jj) * zdb_ml(ji,jj) + 2.0 * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) ) * zddhdt / zdb_bl(ji,jj)
2131	ELSE ! zdb_bl >0
2132	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)
2133	ENDIF
2134	ELSE ! zwb_min + 2*zwb_fk_b < 0
2135	! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
2136	zdhdt(ji,jj) = - zvel_mle(ji,jj)
2137
2138
2139	ENDIF
2140
2141	ELSE
2142	! Fox-Kemper not used.
2143
2144	zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
2145	& MAX((zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln)
2146	zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
2147	! added ajgn 23 July as temporay fix
2148
2149	ENDIF ! ln_osm_mle
2150
2151	ELSE ! lconv - Stable
2152	zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
2153	IF ( zdhdt(ji,jj) < 0._wp ) THEN
2154	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
2155	zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
2156	ELSE
2157	zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
2158	ENDIF
2159	zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
2160	ENDIF ! lconv
2161	ELSE ! lshear
2162	IF ( lconv(ji,jj) ) THEN ! Convective
2163
2164	IF ( ln_osm_mle ) THEN
2165
2166	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
2167	! Fox-Kemper buoyancy flux average over OSBL
2168	zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * &
2169	(1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
2170	ELSE
2171	zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
2172	ENDIF
2173	zvel_max = ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
2174	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
2175	! OSBL is deepening, entrainment > restratification
2176	IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
2177	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
2178	ELSE
2179	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)
2180	ENDIF
2181	ELSE
2182	! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
2183	zdhdt(ji,jj) = - zvel_mle(ji,jj)
2184
2185
2186	ENDIF
2187
2188	ELSE
2189	! Fox-Kemper not used.
2190
2191	zvel_max = -zwb_ent(ji,jj) / &
2192	& MAX((zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln)
2193	zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
2194	! added ajgn 23 July as temporay fix
2195
2196	ENDIF ! ln_osm_mle
2197
2198	ELSE ! Stable
2199	zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
2200	IF ( zdhdt(ji,jj) < 0._wp ) THEN
2201	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
2202	zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj)
2203	ELSE
2204	zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
2205	ENDIF
2206	zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
2207	ENDIF ! lconv
2208	ENDIF ! lshear
2209	END DO
2210	END DO
2211	END SUBROUTINE zdf_osm_calculate_dhdt
2212
2213	SUBROUTINE zdf_osm_timestep_hbl( zdhdt )
2214	!!---------------------------------------------------------------------
2215	!! * ROUTINE zdf_osm_timestep_hbl *
2216	!!
2217	!! ** Purpose : Increments hbl.
2218	!!
2219	!! ** Method : If thechange in hbl exceeds one model level the change is
2220	!! is calculated by moving down the grid, changing the buoyancy
2221	!! jump. This is to ensure that the change in hbl does not
2222	!! overshoot a stable layer.
2223	!!
2224	!!----------------------------------------------------------------------
2225
2226
2227	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! rates of change of hbl.
2228
2229	INTEGER :: jk, jj, ji, jm
2230	REAL(wp) :: zhbl_s, zvel_max, zdb
2231	REAL(wp) :: zthermal, zbeta
2232
2233	DO jj = 2, jpjm1
2234	DO ji = 2, jpim1
2235	IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
2236	!
2237	! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
2238	!
2239	zhbl_s = hbl(ji,jj)
2240	jm = imld(ji,jj)
2241	zthermal = rab_n(ji,jj,1,jp_tem)
2242	zbeta = rab_n(ji,jj,1,jp_sal)
2243
2244
2245	IF ( lconv(ji,jj) ) THEN
2246	!unstable
2247
2248	IF( ln_osm_mle ) THEN
2249	zvel_max = ( zwstrl(ji,jj)3 + zwstrc(ji,jj)3 )**p2third / hbl(ji,jj)
2250	ELSE
2251
2252	zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
2253	& ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
2254
2255	ENDIF
2256
2257	DO jk = imld(ji,jj), ibld(ji,jj)
2258	zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
2259	& - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
2260	& 0.0 ) + zvel_max
2261
2262
2263	IF ( ln_osm_mle ) THEN
2264	zhbl_s = zhbl_s + MIN( &
2265	& rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
2266	& e3w_n(ji,jj,jm) )
2267	ELSE
2268	zhbl_s = zhbl_s + MIN( &
2269	& rn_rdt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
2270	& e3w_n(ji,jj,jm) )
2271	ENDIF
2272
2273	! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2274	IF ( zhbl_s >= gdepw_n(ji,jj,mbkt(ji,jj) + 1) ) THEN
2275	zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2276	lpyc(ji,jj) = .FALSE.
2277	ENDIF
2278	IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
2279	END DO
2280	hbl(ji,jj) = zhbl_s
2281	ibld(ji,jj) = jm
2282	ELSE
2283	! stable
2284	DO jk = imld(ji,jj), ibld(ji,jj)
2285	zdb = MAX( &
2286	& grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&
2287	& - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
2288	& 0.0 ) + &
2289	& 2.0 * zvstr(ji,jj)**2 / zhbl_s
2290
2291	! Alan is thuis right? I have simply changed hbli to hbl
2292	zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
2293	zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)*3 / zhbl_s - 0.15 / 2.0 ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * &
2294	& zustar(ji,jj)*3 / zhbl_s ) ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
2295	zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
2296	zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
2297
2298	! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2299	IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN
2300	zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2301	lpyc(ji,jj) = .FALSE.
2302	ENDIF
2303	IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
2304	END DO
2305	ENDIF ! IF ( lconv )
2306	hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
2307	ibld(ji,jj) = MAX(jm, 4 )
2308	ELSE
2309	! change zero or one model level.
2310	hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
2311	ENDIF
2312	zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
2313	END DO
2314	END DO
2315
2316	END SUBROUTINE zdf_osm_timestep_hbl
2317
2318	SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
2319	!!---------------------------------------------------------------------
2320	!! * ROUTINE zdf_osm_pycnocline_thickness *
2321	!!
2322	!! ** Purpose : Calculates thickness of the pycnocline
2323	!!
2324	!! ** Method : The thickness is calculated from a prognostic equation
2325	!! that relaxes the pycnocine thickness to a diagnostic
2326	!! value. The time change is calculated assuming the
2327	!! thickness relaxes exponentially. This is done to deal
2328	!! with large timesteps.
2329	!!
2330	!!----------------------------------------------------------------------
2331
2332	REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness.
2333	!
2334	INTEGER :: jj, ji
2335	INTEGER :: inhml
2336	REAL(wp) :: zari, ztau, zdh_ref, zddhdt
2337	REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
2338
2339	DO jj = 2, jpjm1
2340	DO ji = 2, jpim1
2341
2342	IF ( lshear(ji,jj) ) THEN
2343	IF ( lconv(ji,jj) ) THEN
2344	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
2345	IF ( j_ddh(ji,jj) == 0 ) THEN
2346	! ddhdt for pycnocline determined in osm_calculate_dhdt
2347	! zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
2348	! zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
2349	! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer
2350	! dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) )
2351	ztau = MAX( zdb_bl(ji,jj) * ( 0.625 * hbl(ji,jj) ) / ( a_ddh * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt )
2352	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.625 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
2353	IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.8 * zhbl(ji,jj)
2354	ELSE
2355	! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt
2356	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
2357	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2358	ELSE
2359	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2360	ENDIF
2361	ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt )
2362	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
2363	IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)
2364	ENDIF
2365	ELSE
2366	ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln), 2.0 * rn_rdt )
2367	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.2 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
2368	IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2 * hbl(ji,jj)
2369	ENDIF
2370	ELSE ! lconv
2371	! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL
2372
2373	ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
2374	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
2375	! boundary layer deepening
2376	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
2377	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
2378	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
2379	& / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )
2380	zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
2381	ELSE
2382	zdh_ref = 0.2 * hbl(ji,jj)
2383	ENDIF
2384	ELSE ! IF(dhdt < 0)
2385	zdh_ref = 0.2 * hbl(ji,jj)
2386	ENDIF ! IF (dhdt >= 0)
2387	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
2388	IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse
2389	! Alan: this hml is never defined or used -- do we need it?
2390	ENDIF
2391
2392	ELSE ! lshear
2393	! for lshear = .FALSE. calculate ddhdt here
2394
2395	IF ( lconv(ji,jj) ) THEN
2396
2397	IF( ln_osm_mle ) THEN
2398	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN
2399	! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
2400	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
2401	IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability
2402	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2403	ELSE ! unstable
2404	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2405	ENDIF
2406	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
2407	zdh_ref = zari * hbl(ji,jj)
2408	ELSE
2409	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
2410	zdh_ref = 0.2 * hbl(ji,jj)
2411	ENDIF
2412	ELSE
2413	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
2414	zdh_ref = 0.2 * hbl(ji,jj)
2415	ENDIF
2416	ELSE ! ln_osm_mle
2417	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
2418	IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability
2419	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2420	ELSE ! unstable
2421	zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)*2 / MAX(4.5 zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
2422	ENDIF
2423	ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
2424	zdh_ref = zari * hbl(ji,jj)
2425	ELSE
2426	ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
2427	zdh_ref = 0.2 * hbl(ji,jj)
2428	ENDIF
2429
2430	END IF ! ln_osm_mle
2431
2432	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
2433	! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
2434	IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
2435	! Alan: this hml is never defined or used
2436	ELSE ! IF (lconv)
2437	ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
2438	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
2439	! boundary layer deepening
2440	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
2441	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
2442	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
2443	& / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )
2444	zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
2445	ELSE
2446	zdh_ref = 0.2 * hbl(ji,jj)
2447	ENDIF
2448	ELSE ! IF(dhdt < 0)
2449	zdh_ref = 0.2 * hbl(ji,jj)
2450	ENDIF ! IF (dhdt >= 0)
2451	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
2452	IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref ! can be a problem with dh>hbl for rapid collapse
2453	ENDIF ! IF (lconv)
2454	ENDIF ! lshear
2455
2456	hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
2457	inhml = MAX( INT( dh(ji,jj) / MAX(e3t_n(ji,jj,ibld(ji,jj)), 1.e-3) ) , 1 )
2458	imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
2459	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
2460	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
2461	END DO
2462	END DO
2463
2464	END SUBROUTINE zdf_osm_pycnocline_thickness
2465
2466
2467	SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
2468	!!----------------------------------------------------------------------
2469	!! * ROUTINE zdf_osm_horizontal_gradients *
2470	!!
2471	!! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.
2472	!!
2473	!! ** Method :
2474	!!
2475	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
2476	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
2477
2478
2479	REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points
2480	REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient.
2481	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == !
2482	REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy
2483
2484	INTEGER :: ji, jj, jk ! dummy loop indices
2485	INTEGER :: ii, ij, ik, ikmax ! local integers
2486	REAL(wp) :: zc
2487	REAL(wp) :: zN2_c ! local buoyancy difference from 10m value
2488	REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH
2489	REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv
2490	REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv
2491	!!----------------------------------------------------------------------
2492	!
2493	! !== MLD used for MLE ==!
2494
2495	mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point
2496	zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2
2497	zN2_c = grav * rn_osm_mle_rho_c * r1_rau0 ! convert density criteria into N^2 criteria
2498	DO jk = nlb10, jpkm1
2499	DO jj = 1, jpj ! Mixed layer level: w-level
2500	DO ji = 1, jpi
2501	ikt = mbkt(ji,jj)
2502	zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
2503	IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
2504	END DO
2505	END DO
2506	END DO
2507	DO jj = 1, jpj
2508	DO ji = 1, jpi
2509	mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
2510	zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
2511	END DO
2512	END DO
2513	! ensure mld_prof .ge. ibld
2514	!
2515	ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation
2516	!
2517	ztm(:,:) = 0._wp
2518	zsm(:,:) = 0._wp
2519	DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer
2520	DO jj = 1, jpj
2521	DO ji = 1, jpi
2522	zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points
2523	ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem)
2524	zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal)
2525	END DO
2526	END DO
2527	END DO
2528	! average temperature and salinity.
2529	ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
2530	zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
2531	! calculate horizontal gradients at u & v points
2532
2533	DO jj = 2, jpjm1
2534	DO ji = 1, jpim1
2535	zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
2536	zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
2537	zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))
2538	ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )
2539	ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )
2540	END DO
2541	END DO
2542
2543	DO jj = 1, jpjm1
2544	DO ji = 2, jpim1
2545	zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
2546	zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
2547	zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))
2548	ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )
2549	ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )
2550	END DO
2551	END DO
2552
2553	! ensure salinity > 0 in unset values so EOS doesn't give FP error with fpe0 on
2554	ztsm_midu(:,jpj,jp_sal) = 10.
2555	ztsm_midv(jpi,:,jp_sal) = 10.
2556
2557	CALL eos_rab(ztsm_midu, zmld_midu, zabu)
2558	CALL eos_rab(ztsm_midv, zmld_midv, zabv)
2559
2560	DO jj = 2, jpjm1
2561	DO ji = 1, jpim1
2562	dbdx_mle(ji,jj) = grav(zdtdx(ji,jj)zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
2563	END DO
2564	END DO
2565	DO jj = 1, jpjm1
2566	DO ji = 2, jpim1
2567	dbdy_mle(ji,jj) = grav(zdtdy(ji,jj)zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
2568	END DO
2569	END DO
2570
2571	DO jj = 2, jpjm1
2572	DO ji = 2, jpim1
2573	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
2574	zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &
2575	& + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
2576	END DO
2577	END DO
2578
2579	END SUBROUTINE zdf_osm_zmld_horizontal_gradients
2580	SUBROUTINE zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
2581	!!----------------------------------------------------------------------
2582	!! * ROUTINE zdf_osm_mle_parameters *
2583	!!
2584	!! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.
2585	!!
2586	!! ** Method :
2587	!!
2588	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
2589	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
2590
2591	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == !
2592	INTEGER, DIMENSION(jpi,jpj) :: mld_prof
2593	REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle
2594	INTEGER :: ji, jj, jk ! dummy loop indices
2595	INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers
2596	INTEGER , DIMENSION(jpi,jpj) :: inml_mle
2597	REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle
2598
2599	! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
2600
2601	DO jj = 2, jpjm1
2602	DO ji = 2, jpim1
2603	IF ( lconv(ji,jj) ) THEN
2604	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
2605	! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
2606	zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1)
2607	zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2
2608	ENDIF
2609	END DO
2610	END DO
2611	! Timestep mixed layer eddy depth.
2612	DO jj = 2, jpjm1
2613	DO ji = 2, jpim1
2614	IF ( lmle(ji,jj) ) THEN ! MLE layer growing.
2615	! Buoyancy gradient at base of MLE layer.
2616	zthermal = rab_n(ji,jj,1,jp_tem)
2617	zbeta = rab_n(ji,jj,1,jp_sal)
2618	jkb = mld_prof(ji,jj)
2619	jkb1 = MIN(jkb + 1, mbkt(ji,jj))
2620	!
2621	zbuoy = grav * ( zthermal * tsn(ji,jj,mld_prof(ji,jj)+2,jp_tem) - zbeta * tsn(ji,jj,mld_prof(ji,jj)+2,jp_sal) )
2622	zdb_mle = zb_bl(ji,jj) - zbuoy
2623	! Timestep hmle.
2624	hmle(ji,jj) = hmle(ji,jj) + zwb0tot(ji,jj) * rn_rdt / zdb_mle
2625	ELSE
2626	IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN
2627	hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau
2628	ELSE
2629	hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt /rn_osm_mle_tau
2630	ENDIF
2631	ENDIF
2632	hmle(ji,jj) = MAX(MIN(hmle(ji,jj), ht_n(ji,jj)), gdepw_n(ji,jj,4))
2633	IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), rn_osm_hmle_limit*hbl(ji,jj) )
2634	hmle(ji,jj) = zmld(ji,jj)
2635	END DO
2636	END DO
2637
2638	mld_prof = 4
2639	DO jk = 5, jpkm1
2640	DO jj = 2, jpjm1
2641	DO ji = 2, jpim1
2642	IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
2643	END DO
2644	END DO
2645	END DO
2646	DO jj = 2, jpjm1
2647	DO ji = 2, jpim1
2648	zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj))
2649	END DO
2650	END DO
2651	END SUBROUTINE zdf_osm_mle_parameters
2652
2653	END SUBROUTINE zdf_osm
2654
2655
2656	SUBROUTINE zdf_osm_init
2657	!!----------------------------------------------------------------------
2658	!! * ROUTINE zdf_osm_init *
2659	!!
2660	!! ** Purpose : Initialization of the vertical eddy diffivity and
2661	!! viscosity when using a osm turbulent closure scheme
2662	!!
2663	!! ** Method : Read the namosm namelist and check the parameters
2664	!! called at the first timestep (nit000)
2665	!!
2666	!! ** input : Namlist namosm
2667	!!----------------------------------------------------------------------
2668	INTEGER :: ios ! local integer
2669	INTEGER :: ji, jj, jk ! dummy loop indices
2670	REAL z1_t2
2671	!!
2672	NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
2673	& ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd &
2674	& ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
2675	& ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter
2676	! Namelist for Fox-Kemper parametrization.
2677	NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
2678	& rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit
2679
2680	!!----------------------------------------------------------------------
2681	!
2682	REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model
2683	READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
2684	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
2685
2686	REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
2687	READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
2688	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
2689	IF(lwm) WRITE ( numond, namzdf_osm )
2690
2691	IF(lwp) THEN ! Control print
2692	WRITE(numout,*)
2693	WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
2694	WRITE(numout,*) '~~~~~~~~~~~~'
2695	WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters'
2696	WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la
2697	WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle
2698	WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la
2699	WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd
2700	WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0
2701	WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes
2702	WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave
2703	WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave
2704	SELECT CASE (nn_osm_wave)
2705	CASE(0)
2706	WRITE(numout,*) ' calculated assuming constant La#=0.3'
2707	CASE(1)
2708	WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves'
2709	CASE(2)
2710	WRITE(numout,*) ' calculated from ECMWF wave fields'
2711	END SELECT
2712	WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce', nn_osm_SD_reduce
2713	WRITE(numout,*) ' fraction of hbl to average SD over/fit'
2714	WRITE(numout,*) ' exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac
2715	SELECT CASE (nn_osm_SD_reduce)
2716	CASE(0)
2717	WRITE(numout,*) ' No reduction'
2718	CASE(1)
2719	WRITE(numout,*) ' Average SD over upper rn_osm_hblfrac of BL'
2720	CASE(2)
2721	WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL'
2722	END SELECT
2723	WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter
2724	WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm
2725	WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s'
2726	WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix
2727	WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
2728	WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri
2729	WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix
2730	WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv
2731	ENDIF
2732
2733
2734	! ! Check wave coupling settings !
2735	! ! Further work needed - see ticket #2447 !
2736	IF( nn_osm_wave == 2 ) THEN
2737	IF (.NOT. ( ln_wave .AND. ln_sdw )) &
2738	& CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' )
2739	END IF
2740
2741	! ! allocate zdfosm arrays
2742	IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
2743
2744
2745	IF( ln_osm_mle ) THEN
2746	! Initialise Fox-Kemper parametrization
2747	REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme
2748	READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
2749	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist')
2750
2751	REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme
2752	READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
2753	904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist')
2754	IF(lwm) WRITE ( numond, namosm_mle )
2755
2756	IF(lwp) THEN ! Namelist print
2757	WRITE(numout,*)
2758	WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
2759	WRITE(numout,*) '~~~~~~~~~~~~~'
2760	WRITE(numout,*) ' Namelist namosm_mle : '
2761	WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle
2762	WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce
2763	WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm'
2764	WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's'
2765	WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg'
2766	WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c
2767	WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s'
2768	WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's'
2769	WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit
2770	WRITE(numout,*) ' fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T. rn_osm_hmle_limit = ', rn_osm_hmle_limit
2771	ENDIF !
2772	ENDIF
2773	!
2774	IF(lwp) THEN
2775	WRITE(numout,*)
2776	IF( ln_osm_mle ) THEN
2777	WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
2778	IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation'
2779	IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation'
2780	ELSE
2781	WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
2782	ENDIF
2783	ENDIF
2784	!
2785	IF( ln_osm_mle ) THEN ! MLE initialisation
2786	!
2787	rb_c = grav * rn_osm_mle_rho_c /rau0 ! Mixed Layer buoyancy criteria
2788	IF(lwp) WRITE(numout,*)
2789	IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 '
2790	IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
2791	!
2792	IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation !
2793	!
2794	ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation
2795	rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) )
2796	!
2797	ENDIF
2798	! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
2799	z1_t2 = 2.e-5
2800	do jj=1,jpj
2801	do ji = 1,jpi
2802	r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
2803	end do
2804	end do
2805	! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
2806	! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 )
2807	!
2808	ENDIF
2809
2810	call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle
2811
2812
2813	IF( ln_zdfddm) THEN
2814	IF(lwp) THEN
2815	WRITE(numout,*)
2816	WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
2817	WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm '
2818	ENDIF
2819	ENDIF
2820
2821
2822	!set constants not in namelist
2823	!-----------------------------
2824
2825	IF(lwp) THEN
2826	WRITE(numout,*)
2827	ENDIF
2828
2829	IF (nn_osm_wave == 0) THEN
2830	dstokes(:,:) = rn_osm_dstokes
2831	END IF
2832
2833	! Horizontal average : initialization of weighting arrays
2834	! -------------------
2835
2836	SELECT CASE ( nn_ave )
2837
2838	CASE ( 0 ) ! no horizontal average
2839	IF(lwp) WRITE(numout,*) ' no horizontal average on avt'
2840	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
2841	! weighting mean arrays etmean
2842	! ( 1 1 )
2843	! avt = 1/4 ( 1 1 )
2844	!
2845	etmean(:,:,:) = 0.e0
2846
2847	DO jk = 1, jpkm1
2848	DO jj = 2, jpjm1
2849	DO ji = 2, jpim1 ! vector opt.
2850	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
2851	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
2852	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
2853	END DO
2854	END DO
2855	END DO
2856
2857	CASE ( 1 ) ! horizontal average
2858	IF(lwp) WRITE(numout,*) ' horizontal average on avt'
2859	! weighting mean arrays etmean
2860	! ( 1/2 1 1/2 )
2861	! avt = 1/8 ( 1 2 1 )
2862	! ( 1/2 1 1/2 )
2863	etmean(:,:,:) = 0.e0
2864
2865	DO jk = 1, jpkm1
2866	DO jj = 2, jpjm1
2867	DO ji = 2, jpim1 ! vector opt.
2868	etmean(ji,jj,jk) = tmask(ji, jj,jk) &
2869	& / MAX( 1., 2.* tmask(ji,jj,jk) &
2870	& +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) &
2871	& +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
2872	& +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) &
2873	& +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) )
2874	END DO
2875	END DO
2876	END DO
2877
2878	CASE DEFAULT
2879	WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave
2880	CALL ctl_stop( ctmp1 )
2881
2882	END SELECT
2883
2884	! Initialization of vertical eddy coef. to the background value
2885	! -------------------------------------------------------------
2886	DO jk = 1, jpk
2887	avt (:,:,jk) =