Context Navigation

zdfosm.F90 @ 14774

Last change on this file since 14774 was 14774, checked in by dancopsey, 3 years ago

Merge in revisions 14705 to 14759 of:

http://forge.ipsl.jussieu.fr/nemo/log/NEMO/branches/NERC/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0

File size: 203.5 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 zdfdrg, ONLY : rCdU_bot ! bottom friction velocity
60	USE zdfddm ! double diffusion mixing (avs array)
61	! USE zdfmxl ! mixed layer depth
62	USE iom ! I/O library
63	USE lib_mpp ! MPP library
64	USE trd_oce ! ocean trends definition
65	USE trdtra ! tracers trends
66	!
67	USE in_out_manager ! I/O manager
68	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
69	USE prtctl ! Print control
70	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
71
72	IMPLICIT NONE
73	PRIVATE
74
75	PUBLIC zdf_osm ! routine called by step.F90
76	PUBLIC zdf_osm_init ! routine called by nemogcm.F90
77	PUBLIC osm_rst ! routine called by step.F90
78	PUBLIC tra_osm ! routine called by step.F90
79	PUBLIC trc_osm ! routine called by trcstp.F90
80	PUBLIC dyn_osm ! routine called by step.F90
81
82	PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90
83
84	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux
85	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamv !: non-local v-momentum flux
86	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamt !: non-local temperature flux (gamma/<ws>o)
87	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghams !: non-local salinity flux (gamma/<ws>o)
88	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt
89	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth
90	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline
91	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth
92	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift.
93
94	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts
95	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmle ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization
96	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdx_mle ! zonal buoyancy gradient in ML
97	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdy_mle ! meridional buoyancy gradient in ML
98	INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_prof ! level of base of MLE layer.
99
100	! !! Namelist namzdf_osm
101	LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la
102
103	LOGICAL :: ln_osm_mle !: flag to activate the Mixed Layer Eddy (MLE) parameterisation
104
105	REAL(wp) :: rn_osm_la ! Turbulent Langmuir number
106	REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift
107	REAL(wp) :: rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by
108	REAL(wp) :: rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl
109	LOGICAL :: ln_zdfosm_ice_shelter ! flag to activate ice sheltering
110	REAL(wp) :: rn_osm_hbl0 = 10._wp ! Initial value of hbl for 1D runs
111	INTEGER :: nn_ave ! = 0/1 flag for horizontal average on avt
112	INTEGER :: nn_osm_wave = 0 ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
113	INTEGER :: nn_osm_SD_reduce ! = 0/1/2 flag for getting effective stokes drift from surface value
114	LOGICAL :: ln_dia_osm ! Use namelist rn_osm_la
115
116	LOGICAL :: ln_kpprimix = .true. ! Shear instability mixing
117	REAL(wp) :: rn_riinfty = 0.7 ! local Richardson Number limit for shear instability
118	REAL(wp) :: rn_difri = 0.005 ! maximum shear mixing at Rig = 0 (m2/s)
119	LOGICAL :: ln_convmix = .true. ! Convective instability mixing
120	REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s)
121
122	#ifdef key_osm_debug
123	INTEGER :: nn_idb = 297, nn_jdb = 193, nn_kdb = 35, nn_narea_db = 109
124	INTEGER :: iloc_db, jloc_db
125	#endif
126	!
127
128	! OSMOSIS mixed layer eddy parametrization constants
129	INTEGER :: nn_osm_mle ! = 0/1 flag for horizontal average on avt
130	REAL(wp) :: rn_osm_mle_ce ! MLE coefficient
131	! ! parameters used in nn_osm_mle = 0 case
132	REAL(wp) :: rn_osm_mle_lf ! typical scale of mixed layer front
133	REAL(wp) :: rn_osm_mle_time ! time scale for mixing momentum across the mixed layer
134	! ! parameters used in nn_osm_mle = 1 case
135	REAL(wp) :: rn_osm_mle_lat ! reference latitude for a 5 km scale of ML front
136	LOGICAL :: ln_osm_hmle_limit ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
137	REAL(wp) :: rn_osm_hmle_limit ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
138	REAL(wp) :: rn_osm_mle_rho_c ! Density criterion for definition of MLD used by FK
139	REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation
140	REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld
141	REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case
142	REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base.
143	REAL(wp) :: rn_osm_bl_thresh ! Threshold buoyancy for deepening of OSBL base.
144	REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE.
145
146
147	! !!! General constants
148	REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero
149	REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw
150	REAL(wp) :: pthird = 1._wp/3._wp ! 1/3
151	REAL(wp) :: p2third = 2._wp/3._wp ! 2/3
152
153	INTEGER :: idebug = 236
154	INTEGER :: jdebug = 228
155	!!----------------------------------------------------------------------
156	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
157	!! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $
158	!! Software governed by the CeCILL license (see ./LICENSE)
159	!!----------------------------------------------------------------------
160	CONTAINS
161
162	INTEGER FUNCTION zdf_osm_alloc()
163	!!----------------------------------------------------------------------
164	!! * FUNCTION zdf_osm_alloc *
165	!!----------------------------------------------------------------------
166	ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &
167	& hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
168	& etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
169
170	ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), &
171	& mld_prof(jpi,jpj), STAT= zdf_osm_alloc )
172
173	! ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & ! would ths be better ?
174	! & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
175	! & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
176	! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
177	! IF ( ln_osm_mle ) THEN
178	! Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc )
179	! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays')
180	! ENDIF
181
182	IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
183	CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
184	END FUNCTION zdf_osm_alloc
185
186
187	SUBROUTINE zdf_osm( kt, p_avm, p_avt )
188	!!----------------------------------------------------------------------
189	!! * ROUTINE zdf_osm *
190	!!
191	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
192	!! coefficients and non local mixing using the OSMOSIS scheme
193	!!
194	!! ** Method : The boundary layer depth hosm is diagnosed at tracer points
195	!! from profiles of buoyancy, and shear, and the surface forcing.
196	!! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
197	!!
198	!! Kx = hosm Wx(sigma) G(sigma)
199	!!
200	!! and the non local term ghamt = Cs / Ws(sigma) / hosm
201	!! Below hosm the coefficients are the sum of mixing due to internal waves
202	!! shear instability and double diffusion.
203	!!
204	!! -1- Compute the now interior vertical mixing coefficients at all depths.
205	!! -2- Diagnose the boundary layer depth.
206	!! -3- Compute the now boundary layer vertical mixing coefficients.
207	!! -4- Compute the now vertical eddy vicosity and diffusivity.
208	!! -5- Smoothing
209	!!
210	!! N.B. The computation is done from jk=2 to jpkm1
211	!! Surface value of avt are set once a time to zero
212	!! in routine zdf_osm_init.
213	!!
214	!! ** Action : update the non-local terms ghamts
215	!! update avt (before vertical eddy coef.)
216	!!
217	!! References : Large W.G., Mc Williams J.C. and Doney S.C.
218	!! Reviews of Geophysics, 32, 4, November 1994
219	!! Comments in the code refer to this paper, particularly
220	!! the equation number. (LMD94, here after)
221	!!----------------------------------------------------------------------
222	INTEGER , INTENT(in ) :: kt ! ocean time step
223	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
224	!!
225	INTEGER :: ji, jj, jk ! dummy loop indices
226
227	INTEGER :: jl ! dummy loop indices
228
229	INTEGER :: ikbot, jkmax, jkm1, jkp2 !
230
231	REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube !
232	REAL(wp) :: zbeta, zthermal !
233	REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
234	REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm !
235
236	REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density
237	INTEGER :: jm ! dummy loop indices
238	REAL(wp) :: zr1, zr2, zr3, zr4, zrhop ! Compression terms
239	REAL(wp) :: zflag, zrn2, zdep21, zdep32, zdep43
240	REAL(wp) :: zesh2, zri, zfri ! Interior richardson mixing
241	REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
242	REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages
243	! Scales
244	REAL(wp), DIMENSION(jpi,jpj) :: zrad0 ! Surface solar temperature flux (deg m/s)
245	REAL(wp), DIMENSION(jpi,jpj) :: zradh ! Radiative flux at bl base (Buoyancy units)
246	REAL(wp), DIMENSION(jpi,jpj) :: zradav ! Radiative flux, bl average (Buoyancy Units)
247	REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity
248	REAL(wp), DIMENSION(jpi,jpj) :: zwstrl ! Langmuir velocity scale
249	REAL(wp), DIMENSION(jpi,jpj) :: zvstr ! Velocity scale that ends to zustar for large Langmuir number.
250	REAL(wp), DIMENSION(jpi,jpj) :: zwstrc ! Convective velocity scale
251	REAL(wp), DIMENSION(jpi,jpj) :: zuw0 ! Surface u-momentum flux
252	REAL(wp), DIMENSION(jpi,jpj) :: zvw0 ! Surface v-momentum flux
253	REAL(wp), DIMENSION(jpi,jpj) :: zwth0 ! Surface heat flux (Kinematic)
254	REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux
255	REAL(wp), DIMENSION(jpi,jpj) :: zwb0 ! Surface buoyancy flux
256	REAL(wp), DIMENSION(jpi,jpj) :: zwb0tot ! Total surface buoyancy flux including insolation
257	REAL(wp), DIMENSION(jpi,jpj) :: zwthav ! Heat flux - bl average
258	REAL(wp), DIMENSION(jpi,jpj) :: zwsav ! freshwater flux - bl average
259	REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average
260	REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux
261	REAL(wp), DIMENSION(jpi,jpj) :: zwb_min
262
263
264	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL
265	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux
266	REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff
267	REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK
268
269	REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift
270	REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number
271	REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
272	REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
273	REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer
274	LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl
275	LOGICAL, DIMENSION(jpi,jpj) :: lshear ! Shear layers
276	LOGICAL, DIMENSION(jpi,jpj) :: lcoup ! Coupling to bottom
277	LOGICAL, DIMENSION(jpi,jpj) :: lpyc ! OSBL pycnocline present
278	LOGICAL, DIMENSION(jpi,jpj) :: lflux ! surface flux extends below OSBL into MLE layer.
279	LOGICAL, DIMENSION(jpi,jpj) :: lmle ! MLE layer increases in hickness.
280
281	! mixed-layer variables
282
283	INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
284	INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
285	INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level
286	INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer
287
288	REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
289	REAL(wp) :: zugrad,zvgrad ! temporary variables for calculating pycnocline shear
290
291	REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid
292	REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid
293
294	REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid
295	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid
296
297	REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid
298	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
299	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext ! external temperature/salinity and buoyancy gradients
300	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext ! external temperature/salinity and buoyancy gradients
301	REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization.
302
303	REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl ! averages over the depth of the blayer
304	REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml ! averages over the depth of the mixed layer
305	REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle ! averages over the depth of the MLE layer
306	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
307	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
308	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
309	! REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
310	REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
311	REAL(wp) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline
312	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline
313	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline
314	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline
315	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc ! u-shear across the pycnocline
316	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc ! v-shear across the pycnocline
317	REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient.
318	! Flux-gradient relationship variables
319	REAL(wp), DIMENSION(jpi, jpj) :: zshear ! Shear production.
320
321	REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale.
322
323	REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline.
324	REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
325	REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/
326	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.
327	REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
328
329	! For calculating Ri#-dependent mixing
330	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du ! u-shear^2
331	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv ! v-shear^2
332	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
333
334	! Temporary variables
335	INTEGER :: inhml
336	REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
337	REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables
338	REAL(wp) :: zthick, zz0, zz1 ! temporary variables
339	REAL(wp) :: zvel_max, zhbl_s ! temporary variables
340	REAL(wp) :: zfac, ztmp ! temporary variable
341	REAL(wp) :: zus_x, zus_y ! temporary Stokes drift
342	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
343	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
344	REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc
345	REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML.
346	REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc
347	REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max
348	INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme
349	REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau
350	REAL(wp) :: zzeta_s = 0._wp
351	REAL(wp) :: zzeta_v = 0.46
352	REAL(wp) :: zabsstke
353	REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness
354	REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc
355
356	! For debugging
357	INTEGER :: ikt
358	REAL(wp) :: zlarge = -1.e10_wp, zero = 0._wp
359	!!--------------------------------------------------------------------
360	!
361	ibld(:,:) = 0 ; imld(:,:) = 0
362	zrad0(:,:) = zlarge ; zradh(:,:) = zlarge ; zradav(:,:) = zlarge ; zustar(:,:) = zlarge
363	zwstrl(:,:) = zlarge ; zvstr(:,:) = zlarge ; zwstrc(:,:) = zlarge ; zuw0(:,:) = zlarge
364	zvw0(:,:) = zlarge ; zwth0(:,:) = zlarge ; zws0(:,:) = zlarge ; zwb0(:,:) = zlarge
365	zwthav(:,:) = zlarge ; zwsav(:,:) = zlarge ; zwbav(:,:) = zlarge ; zwb_ent(:,:) = zlarge
366	zustke(:,:) = zlarge ; zla(:,:) = zlarge ; zcos_wind(:,:) = zlarge ; zsin_wind(:,:) = zlarge
367	zhol(:,:) = zlarge ; zwb0tot(:,:) = zlarge; zalpha_pyc(:,:) = zlarge
368	lconv(:,:) = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ; lmle(:,:) = .FALSE.
369	! mixed layer
370	! no initialization of zhbl or zhml (or zdh?)
371	zhbl(:,:) = zlarge ; zhml(:,:) = zlarge ; zdh(:,:) = zlarge ; zdhdt(:,:) = zlarge
372	zt_bl(:,:) = zlarge ; zs_bl(:,:) = zlarge ; zu_bl(:,:) = zlarge
373	zv_bl(:,:) = zlarge ; zb_bl(:,:) = zlarge
374	zt_ml(:,:) = zlarge ; zs_ml(:,:) = zlarge ; zu_ml(:,:) = zlarge
375	zt_mle(:,:) = zlarge ; zs_mle(:,:) = zlarge ; zu_mle(:,:) = zlarge
376	zb_mle(:,:) = zlarge
377	zv_ml(:,:) = zlarge ; zdt_bl(:,:) = zlarge ; zds_bl(:,:) = zlarge
378	zdu_bl(:,:) = zlarge ; zdv_bl(:,:) = zlarge ; zdb_bl(:,:) = zlarge
379	zdt_ml(:,:) = zlarge ; zds_ml(:,:) = zlarge ; zdu_ml(:,:) = zlarge ; zdv_ml(:,:) = zlarge
380	zdb_ml(:,:) = zlarge
381	zdt_mle(:,:) = zlarge ; zds_mle(:,:) = zlarge ; zdu_mle(:,:) = zlarge
382	zdv_mle(:,:) = zlarge ; zdb_mle(:,:) = zlarge
383	zwth_ent = zlarge ; zws_ent = zlarge
384	!
385	zdtdz_pyc(:,:,:) = zlarge ; zdsdz_pyc(:,:,:) = zlarge ; zdbdz_pyc(:,:,:) = zlarge
386	zdudz_pyc(:,:,:) = zlarge ; zdvdz_pyc(:,:,:) = zlarge
387	zdtdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp ; zdsdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp ; zdbdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp
388	zdudz_pyc(2:jpim1,2:jpjm1,:) = 0._wp ; zdvdz_pyc(2:jpim1,2:jpjm1,:) = 0._wp
389	!
390	zdtdz_bl_ext(:,:) = zlarge ; zdsdz_bl_ext(:,:) = zlarge ; zdbdz_bl_ext(:,:) = zlarge
391
392	IF ( ln_osm_mle ) THEN ! only initialise arrays if needed
393	zdtdx(:,:) = zlarge ; zdtdy(:,:) = zlarge ; zdsdx(:,:) = zlarge
394	zdsdy(:,:) = zlarge ; dbdx_mle(:,:) = zlarge ; dbdy_mle(:,:) = zlarge
395	zwb_fk(:,:) = zlarge ; zvel_mle(:,:) = zlarge; zdiff_mle(:,:) = zlarge
396	zhmle(:,:) = zlarge ; zmld(:,:) = zlarge
397	ENDIF
398	zwb_fk_b(:,:) = zlarge ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt
399
400	! Flux-Gradient arrays.
401	zsc_wth_1(:,:) = zlarge ; zsc_ws_1(:,:) = zlarge ; zsc_uw_1(:,:) = zlarge
402	zsc_uw_2(:,:) = zlarge ; zsc_vw_1(:,:) = zlarge ; zsc_vw_2(:,:) = zlarge
403	zhbl_t(:,:) = zlarge ; zdhdt(:,:) = zlarge
404
405	zdiffut(:,:,:) = zlarge ; zviscos(:,:,:) = zlarge
406	zdiffut(2:jpim1,2:jpjm1,:) = 0._wp ; zviscos(2:jpim1,2:jpjm1,:) = 0._wp
407	ghamt(:,:,:) = zlarge; ghams(:,:,:) = zlarge
408	ghamt(2:jpim1,2:jpjm1,:) = 0._wp; ghams(2:jpim1,2:jpjm1,:) = 0._wp
409	ghamu(:,:,:) = zlarge ; ghamv(:,:,:) = zlarge
410	ghamu(2:jpim1,2:jpjm1,:) = 0._wp ; ghamv(2:jpim1,2:jpjm1,:) = 0._wp
411	zdiff_mle(2:jpim1,2:jpjm1) = 0._wp
412
413
414	#ifdef key_osm_debug
415	IF(mi0(nn_idb)==mi1(nn_idb) .AND. mj0(nn_jdb)==mj1(nn_jdb) .AND. &
416	& mi0(nn_idb) > 1 .AND. mi0(nn_idb) < jpi .AND. mj0(nn_jdb) > 1 .AND. mj0(nn_jdb) < jpj) THEN
417	nn_narea_db = narea
418	iloc_db=mi0(nn_idb); jloc_db=mj0(nn_jdb)
419
420	WRITE(narea+100,*)
421	WRITE(narea+100,'(a,i7)')'timestep=',kt
422	WRITE(narea+100,'(3(a,i7))')'narea=',narea,' nn_idb',nn_idb,' nn_jdb=',nn_jdb
423	WRITE(narea+100,'(4(a,i7))')'iloc_db=',iloc_db,' jloc_db',jloc_db,' jpi=',jpi,' jpj=',jpj
424	ji=iloc_db; jj=jloc_db
425	WRITE(narea+100,'(a,i7,5(a,g10.2))')'mbkt=',mbkt(ji,jj),' ht_n',ht_n(ji,jj),&
426	&' hu_n-',hu_n(ji-1,jj),' hu_n+',hu_n(ji,jj), ' hv_n-',hv_n(ji,jj-1),' hv_n+',hv_n(ji,jj)
427	WRITE(narea+100,*)
428	FLUSH(narea+100)
429	ELSE
430	nn_narea_db = -1000
431	END IF
432	#endif
433
434	! hbl = MAX(hbl,epsln)
435	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
436	! Calculate boundary layer scales
437	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
438
439	! Assume two-band radiation model for depth of OSBL
440	zz0 = rn_abs ! surface equi-partition in 2-bands
441	zz1 = 1. - rn_abs
442	DO jj = 2, jpjm1
443	DO ji = 2, jpim1
444	! Surface downward irradiance (so always +ve)
445	zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp
446	! Downwards irradiance at base of boundary layer
447	zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
448	! Downwards irradiance averaged over depth of the OSBL
449	zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
450	& + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
451	END DO
452	END DO
453	! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
454	DO jj = 2, jpjm1
455	DO ji = 2, jpim1
456	zthermal = rab_n(ji,jj,1,jp_tem)
457	zbeta = rab_n(ji,jj,1,jp_sal)
458	! Upwards surface Temperature flux for non-local term
459	zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
460	! Upwards surface salinity flux for non-local term
461	zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
462	! Non radiative upwards surface buoyancy flux
463	zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj)
464	! Total upwards surface buoyancy flux
465	zwb0tot(ji,jj) = zwb0(ji,jj) - grav * zthermal * ( zrad0(ji,jj) - zradh(ji,jj) )
466	! turbulent heat flux averaged over depth of OSBL
467	zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
468	! turbulent salinity flux averaged over depth of the OBSL
469	zwsav(ji,jj) = 0.5 * zws0(ji,jj)
470	! turbulent buoyancy flux averaged over the depth of the OBSBL
471	zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj)
472	! Surface upward velocity fluxes
473	zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rau0 * tmask(ji,jj,1)
474	zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rau0 * tmask(ji,jj,1)
475	! Friction velocity (zustar), at T-point : LMD94 eq. 2
476	zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
477	zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
478	zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
479	#ifdef key_osm_debug
480	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
481	WRITE(narea+100,'(4(3(a,g11.3),/), 2(a,g11.3),/)') &
482	& 'after calculating fluxes: hbl=', hbl(ji,jj),' zthermal=',zthermal, ' zbeta=', zbeta,&
483	& ' zrad0=', zrad0(ji,jj),' zradh=', zradh(ji,jj), ' zradav=', zradav(ji,jj), &
484	& ' zwth0=', zwth0(ji,jj), ' zwthav=', zwthav(ji,jj), ' zws0=', zws0(ji,jj), &
485	& ' zwb0=', zwb0(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' zwb0tot_in hbl=', zwb0tot(ji,jj) + grav * zthermal * zradh(ji,jj),&
486	& ' zwbav=', zwbav(ji,jj)
487	FLUSH(narea+100)
488	END IF
489	#endif
490	END DO
491	END DO
492	! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
493	SELECT CASE (nn_osm_wave)
494	! Assume constant La#=0.3
495	CASE(0)
496	DO jj = 2, jpjm1
497	DO ji = 2, jpim1
498	zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
499	zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
500	! Linearly
501	zustke(ji,jj) = MAX ( SQRT( zus_xzus_x + zus_yzus_y), 1.0e-8 )
502	dstokes(ji,jj) = rn_osm_dstokes
503	END DO
504	END DO
505	! Assume Pierson-Moskovitz wind-wave spectrum
506	CASE(1)
507	DO jj = 2, jpjm1
508	DO ji = 2, jpim1
509	! Use wind speed wndm included in sbc_oce module
510	zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
511	dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
512	END DO
513	END DO
514	! Use ECMWF wave fields as output from SBCWAVE
515	CASE(2)
516	zfac = 2.0_wp * rpi / 16.0_wp
517
518	DO jj = 2, jpjm1
519	DO ji = 2, jpim1
520	IF (hsw(ji,jj) > 1.e-4) THEN
521	! Use wave fields
522	zabsstke = SQRT(ut0sd(ji,jj)2 + vt0sd(ji,jj)2)
523	zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8)
524	dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)hsw(ji,jj) / ( MAX(zabsstke wmp(ji,jj), 1.0e-7 ) ), 5.0e-1)
525	ELSE
526	! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run)
527	! .. so default to Pierson-Moskowitz
528	zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
529	dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
530	END IF
531	END DO
532	END DO
533	END SELECT
534	#ifdef key_osm_debug
535	IF(narea==nn_narea_db)THEN
536	WRITE(narea+100,'(2(a,g11.3))') &
537	& 'Before reduction: zustke=', zustke(iloc_db,jloc_db),' dstokes =',dstokes(iloc_db,jloc_db)
538	FLUSH(narea+100)
539	END IF
540	#endif
541
542	IF (ln_zdfosm_ice_shelter) THEN
543	! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
544	DO jj = 2, jpjm1
545	DO ji = 2, jpim1
546	zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj))
547	dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj))
548	END DO
549	END DO
550	END IF
551
552	SELECT CASE (nn_osm_SD_reduce)
553	! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van Roekel (2012) or Grant (2020).
554	CASE(0)
555	! The Langmur number from the ECMWF model (or from PM) appears to give La<0.3 for wind-driven seas.
556	! The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3 in this situation.
557	! It could represent the effects of the spread of wave directions
558	! around the mean wind. The effect of this adjustment needs to be tested.
559	IF(nn_osm_wave > 0) THEN
560	zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1)
561	END IF
562	CASE(1)
563	! van Roekel (2012): consider average SD over top 10% of boundary layer
564	! assumes approximate depth profile of SD from Breivik (2016)
565	zsqrtpi = SQRT(rpi)
566	z_two_thirds = 2.0_wp / 3.0_wp
567
568	DO jj = 2, jpjm1
569	DO ji = 2, jpim1
570	zthickness = rn_osm_hblfrac*hbl(ji,jj)
571	z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
572	zsqrt_depth = SQRT(z2k_times_thickness)
573	zexp_depth = EXP(-z2k_times_thickness)
574	zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth &
575	& - z_two_thirds * ( zsqrtpizsqrt_depthz2k_times_thickness * ERFC(zsqrt_depth) &
576	& + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness
577
578	END DO
579	END DO
580	CASE(2)
581	! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer
582	! assumes approximate depth profile of SD from Breivik (2016)
583	zsqrtpi = SQRT(rpi)
584
585	DO jj = 2, jpjm1
586	DO ji = 2, jpim1
587	zthickness = rn_osm_hblfrac*hbl(ji,jj)
588	z2k_times_thickness = zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
589
590	IF(z2k_times_thickness < 50._wp) THEN
591	zsqrt_depth = SQRT(z2k_times_thickness)
592	zexperfc = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP(z2k_times_thickness)
593	ELSE
594	! asymptotic expansion of sqrt(pi)zsqrt_depthEXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness
595	! See Abramowitz and Stegun, Eq. 7.1.23
596	! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness2) - (15/8)/(z2k_times_thickness3)
597	zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp
598	END IF
599	zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp)
600	dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj)
601	zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc)
602	END DO
603	END DO
604	END SELECT
605
606	! Langmuir velocity scale (zwstrl), La # (zla)
607	! mixed scale (zvstr), convective velocity scale (zwstrc)
608	DO jj = 2, jpjm1
609	DO ji = 2, jpim1
610	! Langmuir velocity scale (zwstrl), at T-point
611	zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
612	zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)
613	IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))
614	! Velocity scale that tends to zustar for large Langmuir numbers
615	zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + &
616	& ( 1.0 - EXP( -0.5 * zla(ji,jj)*2 ) ) zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
617
618	! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
619	! Note zustke and zwstrl are not amended.
620	!
621	! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
622	IF ( zwbav(ji,jj) > 0.0) THEN
623	zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
624	zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
625	ELSE
626	zwstrc(ji,jj) = 0.0_wp
627	zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln )
628	ENDIF
629	#ifdef key_osm_debug
630	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
631	WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,3(a,g11.3),/)') &
632	& 'After reduction: zustke=', zustke(ji,jj), ' dstokes=', dstokes(ji,jj), &
633	& ' zustar =', zustar(ji,jj), ' zwstrl=', zwstrl(ji,jj), ' zwstrc=', zwstrc(ji,jj),&
634	& ' zhol=', zhol(ji,jj), ' zla=', zla(ji,jj), ' zvstr=', zvstr(ji,jj)
635	FLUSH(narea+100)
636	END IF
637	#endif
638	END DO
639	END DO
640
641	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
642	! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
643	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
644	! BL must be always 4 levels deep.
645	! For calculation of lateral buoyancy gradients for FK in
646	! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must
647	! previously exist for hbl also.
648
649	! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway
650	! ##########################################################################
651	hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) )
652	ibld(:,:) = 4
653	DO jk = 5, jpkm1
654	DO jj = 1, jpj
655	DO ji = 1, jpi
656	IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
657	ibld(ji,jj) = MIN(mbkt(ji,jj)-2, jk)
658	ENDIF
659	END DO
660	END DO
661	END DO
662	! ##########################################################################
663
664	DO jj = 2, jpjm1
665	DO ji = 2, jpim1
666	zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
667	imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t_n(ji, jj, ibld(ji,jj) - 1 )) , 1 ))
668	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
669	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
670	END DO
671	END DO
672	#ifdef key_osm_debug
673	IF(narea==nn_narea_db) THEN
674	ji=iloc_db; jj=jloc_db
675	WRITE(narea+100,'(2(a,g11.3),/,3(a,g11.3),/,2(a,i7),/)') &
676	& 'Before updating hbl: hbl=', hbl(ji,jj), ' dh=', dh(ji,jj), &
677	&' zhbl =',zhbl(ji,jj) , ' zhml=', zhml(ji,jj), ' zdh=', zdh(ji,jj),&
678	&' imld=', imld(ji,jj), ' ibld=', ibld(ji,jj)
679
680	WRITE(narea+100,'(a,g11.3,a,2g11.3)') 'Physics: ssh ',sshn(ji,jj),' T S surface=',tsn(ji,jj,1,jp_tem),tsn(ji,jj,1,jp_sal)
681	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
682	WRITE(narea+100,'(a,*(g11.3))') ' T[imld-1..ibld+2] =', ( tsn(ji,jj,jk,jp_tem), jk=jl,jm )
683	WRITE(narea+100,'(a,*(g11.3))') ' S[imld-1..ibld+2] =', ( tsn(ji,jj,jk,jp_sal), jk=jl,jm )
684	WRITE(narea+100,'(a,*(g11.3))') ' U+[imld-1..ibld+2] =', ( un(ji,jj,jk), jk=jl,jm )
685	WRITE(narea+100,'(a,*(g11.3))') ' U-[imld-1..ibld+2] =', ( un(ji-1,jj,jk), jk=jl,jm )
686	WRITE(narea+100,'(a,*(g11.3))') ' V+[imld-1..ibld+2] =', ( vn(ji,jj,jk), jk=jl,jm )
687	WRITE(narea+100,'(a,*(g11.3))') ' V-[imld-1..ibld+2] =', ( vn(ji,jj-1,jk), jk=jl,jm )
688	WRITE(narea+100,'(a,*(g11.3))') ' W[imld-1..ibld+2] =', ( wn(ji,jj-1,jk), jk=jl,jm )
689	WRITE(narea+100,*)
690	FLUSH(narea+100)
691	END IF
692	#endif
693
694	! Averages over well-mixed and boundary layer, note BL averages use jp_ext=2 everywhere
695	jp_ext(:,:) = 1 ! ag 19/03
696	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)
697	jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03
698	CALL zdf_osm_vertical_average(imld-1, jp_ext, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
699	#ifdef key_osm_debug
700	IF(narea==nn_narea_db) THEN
701	ji=iloc_db; jj=jloc_db
702	WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') &
703	& 'After averaging, with old hbl (& jp_ext==2), hml: zt_bl=', zt_bl(ji,jj),&
704	& ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),&
705	& 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),&
706	& 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),&
707	& 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),&
708	& 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&
709	& 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)
710	FLUSH(narea+100)
711	END IF
712	#endif
713	! Velocity components in frame aligned with surface stress.
714	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
715	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
716	#ifdef key_osm_debug
717	IF(narea==nn_narea_db) THEN
718	ji=iloc_db; jj=jloc_db
719	WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') &
720	& 'After rotation, with old hbl (& jp_ext==2), hml:', &
721	& 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&
722	& 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)
723	FLUSH(narea+100)
724	END IF
725	#endif
726
727	! Determine the state of the OSBL, stable/unstable, shear/no shear
728	CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear )
729
730	#ifdef key_osm_debug
731	IF(narea==nn_narea_db) THEN
732	ji=iloc_db; jj=jloc_db
733	WRITE(narea+100,'(2(a,l7),a, i7,/,3(a,g11.3),/)') &
734	& 'After zdf_osm_osbl_state: lconv=', lconv(ji,jj), ' lshear=', lshear(ji,jj), ' j_ddh=', j_ddh(ji,jj),&
735	& 'zwb_ent=', zwb_ent(ji,jj), ' zwb_min=', zwb_min(ji,jj), ' zshear=', zshear(ji,jj)
736	FLUSH(narea+100)
737	END IF
738	#endif
739	IF ( ln_osm_mle ) THEN
740	! Fox-Kemper Scheme
741	mld_prof = 4
742	DO jk = 5, jpkm1
743	DO jj = 2, jpjm1
744	DO ji = 2, jpim1
745	IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
746	END DO
747	END DO
748	END DO
749	jp_ext_mle(:,:) = 2
750	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)
751
752	DO jj = 2, jpjm1
753	DO ji = 2, jpim1
754	zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
755	END DO
756	END DO
757	#ifdef key_osm_debug
758	IF(narea==nn_narea_db) THEN
759	ji=iloc_db; jj=jloc_db
760	WRITE(narea+100,'(2(a,g11.3), a, i7,/,2(3(a,g11.3),/),4(a,g11.3),/)') &
761	& 'Before updating hmle: hmle =',hmle(ji,jj) , ' zhmle=', zhmle(ji,jj), ' mld_prof=', mld_prof(ji,jj), &
762	& 'averaging over hmle: zt_mle=', zt_mle(ji,jj), ' zs_mle=', zs_mle(ji,jj), ' zb_mle=', zb_mle(ji,jj),&
763	& 'zdt_mle=', zdt_mle(ji,jj), ' zds_mle=', zds_mle(ji,jj), ' zdb_mle=', zdb_mle(ji,jj),&
764	& 'zu_mle =', zu_mle(ji,jj), ' zv_mle=', zv_mle(ji,jj), ' zdu_mle=', zdu_mle(ji,jj), ' zdv_mle=', zdv_mle(ji,jj)
765	FLUSH(narea+100)
766	END IF
767	#endif
768
769	!! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients
770	CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
771	!! Calculate vertical gradients immediately below zmld
772	CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext )
773	!! calculate max vertical FK flux zwb_fk & set logical descriptors
774	CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
775	!! recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle
776	CALL zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
777	#ifdef key_osm_debug
778	IF(narea==nn_narea_db) THEN
779	ji=iloc_db; jj=jloc_db
780	WRITE(narea+100,'(a,g11.3,a,i7,/, 2(4(a,g11.3),/),2(a,g11.3),/,2(3(a,g11.3),/),a,i7,2(a,g11.3),/,3(a,g11.3),/,/)') &
781	& 'Before updating hmle: zmld =',zmld(ji,jj),' mld_prof=', mld_prof(ji,jj), &
782	& 'zdtdx+=', zdtdx(ji,jj),' zdtdx-=', zdtdx(ji-1,jj),' zdsdx+=', zdsdx(ji,jj),' zdsdx-=',zdsdx(ji-1,jj), &
783	& 'zdtdy+=', zdtdy(ji,jj),' zdtdy-=', zdtdy(ji,jj-1),' zdsdy+=', zdsdy(ji,jj),' zdsdy-=',zdsdy(ji,jj-1), &
784	& 'dbdx_mle+=', dbdx_mle(ji,jj),' dbdx_mle-=', dbdx_mle(ji-1,jj),&
785	& 'dbdy_mle+=', dbdy_mle(ji,jj),' dbdy_mle-=',dbdy_mle(ji,jj-1),' zdbds_mle=',zdbds_mle(ji,jj), &
786	& 'zdtdz_mle_ext=', zdtdz_mle_ext(ji,jj), ' zdsdz_mle_ext=', zdsdz_mle_ext(ji,jj), &
787	& ' zdbdz_mle_ext=', zdbdz_mle_ext(ji,jj), &
788	& 'After updating hmle: mld_prof=', mld_prof(ji,jj),' hmle=', hmle(ji,jj), ' zhmle=', zhmle(ji,jj),&
789	& 'zvel_mle =', zvel_mle(ji,jj), ' zdiff_mle=', zdiff_mle(ji,jj), ' zwb_fk=', zwb_fk(ji,jj)
790	FLUSH(narea+100)
791	END IF
792	#endif
793	ELSE ! ln_osm_mle
794	! FK not selected, Boundary Layer only.
795	lpyc(:,:) = .TRUE.
796	lflux(:,:) = .FALSE.
797	lmle(:,:) = .FALSE.
798	DO jj = 2, jpjm1
799	DO ji = 2, jpim1
800	IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
801	END DO
802	END DO
803	ENDIF ! ln_osm_mle
804
805	!! External gradient below BL needed both with and w/o FK
806	CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext ) ! ag 19/03
807
808	! Test if pycnocline well resolved
809	! DO jj = 2, jpjm1 Removed with ag 19/03 changes. A change in eddy diffusivity/viscosity
810	! DO ji = 2,jpim1 should account for this.
811	! IF (lconv(ji,jj) ) THEN
812	! ztmp = 0.2 * zhbl(ji,jj) / e3w_n(ji,jj,ibld(ji,jj))
813	! IF ( ztmp > 3 ) THEN ! ag 19/03
814	! ! pycnocline well resolved
815	! jp_ext(ji,jj) = 1
816	! ELSE
817	! ! pycnocline poorly resolved
818	! jp_ext(ji,jj) = 0
819	! ENDIF
820	! ELSE
821	! ! Stable conditions
822	! jp_ext(ji,jj) = 0
823	! ENDIF
824	! END DO
825	! END DO
826	#ifdef key_osm_debug
827	IF(narea==nn_narea_db) THEN
828	ji=iloc_db; jj=jloc_db
829	WRITE(narea+100,'(4(a,l7),a,i7,/, 3(a,g11.3),/)') &
830	& 'BL logical descriptors: lconv =',lconv(ji,jj),' lpyc=', lpyc(ji,jj),' lflux=', lflux(ji,jj),' lmle=', lmle(ji,jj),&
831	& ' jp_ext=', jp_ext(ji,jj), &
832	& 'sub-BL strat: zdtdz_bl_ext=', zdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', zdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', zdbdz_bl_ext(ji,jj)
833	FLUSH(narea+100)
834	END IF
835	#endif
836
837	! Recalculate bl averages using jp_ext & ml averages .... note no rotation of u & v here..
838	jp_ext(:,:) = 1 ! ag 19/03
839	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 )
840	jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03
841	CALL zdf_osm_vertical_average(imld-1, jp_ext, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml) ! ag 19/03
842	#ifdef key_osm_debug
843	IF(narea==nn_narea_db) THEN
844	ji=iloc_db; jj=jloc_db
845	WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') &
846	& 'After averaging, with old hbl (&correct jp_ext), hml: zt_bl=', zt_bl(ji,jj),&
847	& ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),&
848	& 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),&
849	& 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),&
850	& 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),&
851	& 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&
852	& 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)
853	FLUSH(narea+100)
854	END IF
855	#endif
856
857
858	! Rate of change of hbl
859	CALL zdf_osm_calculate_dhdt( zdhdt )
860	! Test if surface boundary layer coupled to bottom.
861	lcoup(:,:) = .FALSE. ! ag 19/03
862	DO jj = 2, jpjm1
863	DO ji = 2, jpim1
864	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
865	! adjustment to represent limiting by ocean bottom
866	IF ( mbkt(ji,jj) >2 ) THEN ! to ensure mbkt(ji,jj) - 2 > 0 so no incorrect array access
867	IF( zhbl_t(ji,jj) > gdepw_n(ji, jj, mbkt(ji,jj) - 2 ) ) THEN
868	zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) - 2))! ht_n(:,:))
869	lpyc(ji,jj) = .FALSE.
870	lcoup(ji,jj) = .TRUE. ! ag 19/03
871	ENDIF
872	ENDIF
873	#ifdef key_osm_debug
874	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
875	WRITE(narea+100,'(2(a,g11.3),/,2(a,g11.3),2(a,l7))')'after zdf_osm_calculate_dhdt: zhbl_t=',zhbl_t(ji,jj), 'hbl=', hbl(ji,jj),&
876	& 'delta hbl from dzdhdt', zdhdt(ji,jj)rn_rdt,' delta hbl from w ', wn(ji,jj,ibld(ji,jj))rn_rdt, &
877	& ' lcoup= ',lcoup(ji,jj), ' lpyc= ', lpyc(ji,jj)
878	FLUSH(narea+100)
879	END IF
880	#endif
881	END DO
882	END DO
883
884	imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index
885	ibld(:,:) = 4
886
887	DO jk = 4, jpkm1
888	DO jj = 2, jpjm1
889	DO ji = 2, jpim1
890	IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
891	ibld(ji,jj) = jk
892	ENDIF
893	END DO
894	END DO
895	END DO
896
897	!
898	! Step through model levels taking account of buoyancy change to determine the effect on dhdt
899	!
900	CALL zdf_osm_timestep_hbl( zdhdt )
901	! is external level in bounds?
902
903	! Recalculate BL averages and differences using new BL depth
904	jp_ext(:,:) = 1 ! ag 19/03
905	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 )
906	!
907
908	CALL zdf_osm_pycnocline_thickness( dh, zdh )
909
910	! reset l_pyc before calculating terms in the flux-gradient relationships
911
912	DO jj = 2, jpjm1
913	DO ji = 2, jpim1
914	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) >= mbkt(ji,jj) -2 .or. ibld(ji,jj)-imld(ji,jj) == 1 .or. zdhdt(ji,jj) < 0._wp) THEN ! ag 19/03
915	lpyc(ji,jj) = .FALSE. ! ag 19/03
916	IF ( ibld(ji,jj) >= mbkt(ji,jj) -2 ) THEN
917	imld(ji,jj) = ibld(ji,jj) - 1 ! ag 19/03
918	zdh(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj)) - gdepw_n(ji,jj,imld(ji,jj)) ! ag 19/03
919	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj)) ! ag 19/03
920	dh(ji,jj) = zdh(ji,jj) ! ag 19/03
921	hml(ji,jj) = hbl(ji,jj) - dh(ji,jj) ! ag 19/03
922	#ifdef key_osm_debug
923	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
924	WRITE(narea+100,'(a)')'After setting pycnocline thickness BL running aground: lpyc= F5: ibld(ji,jj) >= mbkt(ji,jj) -2'
925	WRITE(narea+100,'(2(a,i7),2(a,g11.3))')' ibld=',ibld(ji,jj),' imld=',imld(ji,jj), ' zdh=',zdh(ji,jj), ' zhml=',zhml(ji,jj)
926	WRITE(narea+100,'(2(a,g11.3))')'dh=',dh(ji,jj),' hml=',hml(ji,jj)
927	FLUSH(narea+100)
928	END IF
929	#endif
930	ENDIF
931	ENDIF ! ag 19/03
932	END DO
933	END DO
934
935	dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms.
936	!
937	! Average over the depth of the mixed layer in the convective boundary layer
938	! jp_ext = ibld - imld +1
939	! Recalculate ML averages and differences using new ML depth
940	jp_ext(:,:) = ibld(:,:) - imld(:,:) + jp_ext(:,:) + 1 ! ag 19/03
941	CALL zdf_osm_vertical_average( imld-1, jp_ext, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )
942
943	CALL zdf_osm_external_gradients( ibld+1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
944	#ifdef key_osm_debug
945	IF(narea==nn_narea_db) THEN
946	ji=iloc_db; jj=jloc_db
947	WRITE(narea+100,'(4(3(a,g11.3),/), 2(4(a,g11.3),/))') &
948	& 'After averaging, with new hbl (&correct jp_ext), hml: zt_bl=', zt_bl(ji,jj),&
949	& ' zs_bl=', zs_bl(ji,jj), ' zb_bl=', zb_bl(ji,jj),&
950	& 'zdt_bl=', zdt_bl(ji,jj), ' zds_bl=', zds_bl(ji,jj), ' zdb_bl=', zdb_bl(ji,jj),&
951	& 'zt_ml=', zt_ml(ji,jj), ' zs_ml=', zs_ml(ji,jj), ' zb_ml=', zb_ml(ji,jj),&
952	& 'zdt_ml=', zdt_ml(ji,jj), ' zds_ml=', zds_ml(ji,jj), ' zdb_ml=', zdb_ml(ji,jj),&
953	& 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&
954	& 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)
955	FLUSH(narea+100)
956	END IF
957	#endif
958
959	! rotate mean currents and changes onto wind align co-ordinates
960
961	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
962	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
963	#ifdef key_osm_debug
964	IF(narea==nn_narea_db) THEN
965	ji=iloc_db; jj=jloc_db
966	WRITE(narea+100,'(a,/, 2(4(a,g11.3),/))') &
967	& 'After rotation, with new hbl (& correct jp_ext), hml:', &
968	& 'zu_bl =', zu_bl(ji,jj) , ' zv_bl=', zv_bl(ji,jj), ' zdu_bl=', zdu_bl(ji,jj), ' zdv_bl=', zdv_bl(ji,jj),&
969	& 'zu_ml =', zu_ml(ji,jj) , ' zv_ml=', zv_ml(ji,jj), ' zdu_ml=', zdu_ml(ji,jj), ' zdv_ml=', zdv_ml(ji,jj)
970	FLUSH(narea+100)
971	END IF
972	#endif
973	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
974	! Pycnocline gradients for scalars and velocity
975	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
976
977	jp_ext(:,:) = 1 ! ag 19/03
978	CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc )
979	CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )
980	#ifdef key_osm_debug
981	IF(narea==nn_narea_db) THEN
982	ji=iloc_db; jj=jloc_db
983	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
984	WRITE(narea+100,'(a,l7,/,3(a,g11.3),/)') &
985	& 'After pycnocline profiles BL lpyc=', lpyc(ji,jj),&
986	& 'sub-BL strat: zdtdz_bl_ext=', zdtdz_bl_ext(ji,jj),' zdsdz_bl_ext=', zdsdz_bl_ext(ji,jj),' zdbdz_bl_ext=', zdbdz_bl_ext(ji,jj), &
987	& 'Pycnocline: zalpha_pyc=', zalpha_pyc(ji,jj)
988	WRITE(narea+100,'(a,*(g11.3))') ' zdtdz_pyc[imld-1..ibld+2] =', ( zdtdz_pyc(ji,jj,jk), jk=jl,jm )
989	WRITE(narea+100,'(a,*(g11.3))') ' zdsdz_pyc[imld-1..ibld+2] =', ( zdsdz_pyc(ji,jj,jk), jk=jl,jm )
990	WRITE(narea+100,'(a,*(g11.3))') ' zdbdz_pyc[imld-1..ibld+2] =', ( zdbdz_pyc(ji,jj,jk), jk=jl,jm )
991	WRITE(narea+100,'(a,*(g11.3))') ' zdudz_pyc[imld-1..ibld+2] =', ( zdudz_pyc(ji,jj,jk), jk=jl,jm )
992	WRITE(narea+100,'(a,*(g11.3))') ' zdvdz_pyc[imld-1..ibld+2] =', ( zdvdz_pyc(ji,jj,jk), jk=jl,jm )
993	WRITE(narea+100,*)
994	FLUSH(narea+100)
995	END IF
996	#endif
997	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
998	! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
999	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1000	CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
1001	#ifdef key_osm_debug
1002	IF(narea==nn_narea_db) THEN
1003	ji=iloc_db; jj=jloc_db
1004	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1005	WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm )
1006	WRITE(narea+100,'(a,*(g11.3))') ' zviscos[imld-1..ibld+2] =', ( zviscos(ji,jj,jk), jk=jl,jm )
1007	WRITE(narea+100,*)
1008	FLUSH(narea+100)
1009	END IF
1010	#endif
1011
1012	!
1013	! calculate non-gradient components of the flux-gradient relationships
1014	!
1015	! Stokes term in scalar flux, flux-gradient relationship
1016	WHERE ( lconv(2:jpim1,2:jpjm1) )
1017	zsc_wth_1(2:jpim1,2:jpjm1) = zwstrl(2:jpim1,2:jpjm1)*3 zwth0(2:jpim1,2:jpjm1) / &
1018	&( zvstr(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)**3 + epsln)
1019	!
1020	zsc_ws_1(2:jpim1,2:jpjm1) = zwstrl(2:jpim1,2:jpjm1)*3 zws0(2:jpim1,2:jpjm1) / &
1021	&( zvstr(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)**3 + epsln )
1022	ELSEWHERE
1023	zsc_wth_1(2:jpim1,2:jpjm1) = 2.0 * zwthav(2:jpim1,2:jpjm1)
1024	!
1025	zsc_ws_1(2:jpim1,2:jpjm1) = 2.0 * zwsav(2:jpim1,2:jpjm1)
1026	ENDWHERE
1027
1028
1029	DO jj = 2, jpjm1
1030	DO ji = 2, jpim1
1031	IF ( lconv(ji,jj) ) THEN
1032	DO jk = 2, imld(ji,jj)
1033	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1034	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)
1035	!
1036	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)
1037	END DO ! end jk loop
1038	ELSE ! else for if (lconv)
1039	! Stable conditions
1040	DO jk = 2, ibld(ji,jj)
1041	zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1042	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) &
1043	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
1044	!
1045	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 2.15 * EXP ( -0.85 * zznd_d ) &
1046	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
1047	END DO
1048	ENDIF ! endif for check on lconv
1049	END DO ! end of ji loop
1050	END DO ! end of jj loop
1051
1052	! 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)
1053	WHERE ( lconv(2:jpim1,2:jpjm1) )
1054	zsc_uw_1(2:jpim1,2:jpjm1) = ( zwstrl(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)3 )pthird * zustke(2:jpim1,2:jpjm1) / &
1055	& MAX( ( 1.0 - 1.0 * 6.5 * zla(2:jpim1,2:jpjm1)**(8.0/3.0) ), 0.2 )
1056	zsc_uw_2(2:jpim1,2:jpjm1) = ( zwstrl(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)3 )pthird * zustke(2:jpim1,2:jpjm1) / &
1057	& MIN( zla(2:jpim1,2:jpjm1)**(8.0/3.0) + epsln, 0.12 )
1058	zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zhml(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1)*3 MIN( zla(2:jpim1,2:jpjm1)**(8.0/3.0), 0.12 ) / &
1059	& ( ( zvstr(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)3 )(2.0/3.0) + epsln )
1060	ELSEWHERE
1061	zsc_uw_1(2:jpim1,2:jpjm1) = zustar(2:jpim1,2:jpjm1)**2
1062	zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zhbl(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1)*3 &
1063	& MIN( zla(2:jpim1,2:jpjm1)(8.0/3.0), 0.12 ) / (zvstr(2:jpim1,2:jpjm1)2 + epsln)
1064	ENDWHERE
1065	IF(ln_dia_osm) THEN
1066	IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
1067	IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
1068	END IF
1069	DO jj = 2, jpjm1
1070	DO ji = 2, jpim1
1071	IF ( lconv(ji,jj) ) THEN
1072	DO jk = 2, imld(ji,jj)
1073	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1074	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) &
1075	& + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) &
1076	& * ( 1.0 - EXP ( -2.0 * zznd_d ) )
1077	!
1078	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( - zznd_d ) &
1079	& * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
1080	END DO ! end jk loop
1081	ELSE
1082	! Stable conditions
1083	DO jk = 2, ibld(ji,jj) ! corrected to ibld
1084	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1085	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) &
1086	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
1087	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
1088	END DO ! end jk loop
1089	ENDIF
1090	END DO ! ji loop
1091	END DO ! jj loo
1092	#ifdef key_osm_debug
1093	IF(narea==nn_narea_db) THEN
1094	ji=iloc_db; jj=jloc_db
1095	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1096	WRITE(narea+100,'(a,g11.3)')'Stokes contrib to ghamt/s: zsc_wth_1=',zsc_wth_1(ji,jj), ' zsc_ws_1=',zsc_ws_1(ji,jj)
1097	WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm )
1098	WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm )
1099	IF( lconv(ji,jj) ) THEN
1100	WRITE(narea+100,'(3(a,g11.3))')'Stokes contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj), &
1101	&' zsc_uw_2=',zsc_uw_2(ji,jj)
1102	ELSE
1103	WRITE(narea+100,'(2(a,g11.3))')'Stokes contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj)
1104	END IF
1105	WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm )
1106	WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm )
1107	WRITE(narea+100,*)
1108	FLUSH(narea+100)
1109	END IF
1110	#endif
1111
1112	! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
1113
1114	WHERE ( lconv(2:jpim1,2:jpjm1) )
1115	zsc_wth_1(2:jpim1,2:jpjm1) = zwbav(2:jpim1,2:jpjm1) * zwth0(2:jpim1,2:jpjm1) * ( 1.0 + EXP ( 0.2 * zhol(2:jpim1,2:jpjm1) ) ) * &
1116	& zhml(2:jpim1,2:jpjm1) / ( zvstr(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)**3 + epsln )
1117	zsc_ws_1(2:jpim1,2:jpjm1) = zwbav(2:jpim1,2:jpjm1) * zws0(2:jpim1,2:jpjm1) * ( 1.0 + EXP ( 0.2 * zhol(2:jpim1,2:jpjm1) ) ) * &
1118	& zhml(2:jpim1,2:jpjm1) / ( zvstr(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)**3 + epsln )
1119	ELSEWHERE
1120	zsc_wth_1(2:jpim1,2:jpjm1) = 0._wp
1121	zsc_ws_1(2:jpim1,2:jpjm1) = 0._wp
1122	ENDWHERE
1123
1124	DO jj = 2, jpjm1
1125	DO ji = 2, jpim1
1126	IF (lconv(ji,jj) ) THEN
1127	DO jk = 2, imld(ji,jj)
1128	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
1129	! calculate turbulent time scale
1130	zl_c = 0.9 * ( 1.0 - EXP ( - 5.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) &
1131	& * ( 1.0 - EXP ( -15.0 * ( 1.2 - zznd_ml ) ) )
1132	zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml + zznd_ml**3 / 3.0 ) ) ) &
1133	& * ( 1.0 - EXP ( - 8.0 * ( 1.15 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
1134	zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0)
1135	! non-gradient buoyancy terms
1136	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.4 * zsc_wth_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml )
1137	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.4 * zsc_ws_1(ji,jj) * zl_eps / ( 0.15 + zznd_ml )
1138	END DO
1139	IF ( lpyc(ji,jj) ) THEN
1140	ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)3 + zwstrc(ji,jj)3 )**pthird
1141	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 )
1142	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)
1143	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)
1144	! Cubic profile used for buoyancy term
1145	DO jk = 2, ibld(ji,jj)
1146	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1147	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 )
1148
1149	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 )
1150	END DO
1151	#ifdef key_osm_debug
1152	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1153	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1154	WRITE(narea+100,'(3(a,g11.3))')'lpyc= lconv=T: ztau_sc_u=',ztau_sc_u(ji,jj),' zwth_ent=',zwth_ent,' zws_ent=',zws_ent
1155	WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm )
1156	WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm )
1157	WRITE(narea+100,*)
1158	FLUSH(narea+100)
1159	END IF
1160	#endif
1161	!
1162	IF ( dh(ji,jj) < 0.2*hbl(ji,jj) .AND. ibld(ji,jj) - imld(ji,jj) > 3 ) THEN
1163	zbuoy_pyc_sc = 2.0_wp * MAX(zdb_ml(ji,jj), 0._wp) / zdh(ji,jj)
1164	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 ) )
1165	!
1166	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)
1167	!
1168	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)
1169	!
1170	zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
1171	DO jk = 2, ibld(ji,jj)
1172	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1173	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
1174	!
1175	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
1176	END DO
1177	#ifdef key_osm_debug
1178	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1179	WRITE(narea+100,'(2(a,g11.3))')'lpyc= lconv=T,dh<0.2*hbl: zbuoy_pyc_sc=',zbuoy_pyc_sc,' zdelta_pyc=',zdelta_pyc
1180	WRITE(narea+100,'(3(a,g11.3))')'zwt_pyc_sc_1=',zwt_pyc_sc_1,' zws_pyc_sc_1=',zws_pyc_sc_1,' zzeta_pyc=',zzeta_pyc
1181	FLUSH(narea+100)
1182	END IF
1183	#endif
1184
1185	END IF
1186	ENDIF ! End of pycnocline
1187	ELSE ! lconv test - stable conditions
1188	DO jk = 2, ibld(ji,jj)
1189	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
1190	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj)
1191	END DO
1192	ENDIF
1193	END DO ! ji loop
1194	END DO ! jj loop
1195
1196	WHERE ( lconv(2:jpim1,2:jpjm1) )
1197	zsc_uw_1(2:jpim1,2:jpjm1) = -zwb0(2:jpim1,2:jpjm1) * zustar(2:jpim1,2:jpjm1)*2 zhml(2:jpim1,2:jpjm1) / &
1198	& ( zvstr(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)**3 + epsln )
1199	zsc_uw_2(2:jpim1,2:jpjm1) = zwb0(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1) * zhml(2:jpim1,2:jpjm1) / &
1200	& ( zvstr(2:jpim1,2:jpjm1)*3 + 0.5 zwstrc(2:jpim1,2:jpjm1)3 + epsln )(2.0/3.0)
1201	zsc_vw_1(2:jpim1,2:jpjm1) = 0._wp
1202	ELSEWHERE
1203	zsc_uw_1(2:jpim1,2:jpjm1) = 0._wp
1204	zsc_vw_1(2:jpim1,2:jpjm1) = 0._wp
1205	ENDWHERE
1206
1207	DO jj = 2, jpjm1
1208	DO ji = 2, jpim1
1209	IF ( lconv(ji,jj) ) THEN
1210	DO jk = 2 , imld(ji,jj)
1211	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1212	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) &
1213	& * ( 1.0 - EXP( -0.5 * zznd_d ) ) &
1214	& * zsc_uw_2(ji,jj) )
1215	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
1216	END DO ! jk loop
1217	ELSE
1218	! stable conditions
1219	DO jk = 2, ibld(ji,jj)
1220	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
1221	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
1222	END DO
1223	ENDIF
1224	END DO ! ji loop
1225	END DO ! jj loop
1226
1227	DO jj = 2, jpjm1
1228	DO ji = 2, jpim1
1229	IF( lconv(ji,jj) ) THEN
1230	IF ( lpyc(ji,jj) ) THEN
1231	IF ( j_ddh(ji,jj) == 0 ) THEN
1232	! Place holding code. Parametrization needs checking for these conditions.
1233	zomega = ( 0.15 * zwstrl(ji,jj)3 + zwstrc(ji,jj)3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird
1234	zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
1235	zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
1236	ELSE
1237	zomega = ( 0.15 * zwstrl(ji,jj)3 + zwstrc(ji,jj)3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird
1238	zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
1239	zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
1240	ENDIF
1241	zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse
1242	zc_cubic = zuw_bse - zd_cubic
1243	! need ztau_sc_u to be available. Change to array.
1244	DO jk = imld(ji,jj), ibld(ji,jj)
1245	zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1246	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)*2 ) zuw_bse * &
1247	& ( zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc*3 ) ( 0.75 + 0.25 * zznd_pyc )*2 zdbdz_pyc(ji,jj,jk)
1248	END DO
1249	zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) )
1250	zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse
1251	zc_cubic = zvw_bse - zd_cubic
1252	DO jk = imld(ji,jj), ibld(ji,jj)
1253	zznd_pyc = -( gdepw_n(ji,jj,jk) -zhbl(ji,jj) ) / zdh(ji,jj)
1254	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ( ztau_sc_u(ji,jj)*2 ) zvw_bse * &
1255	& ( zc_cubic * zznd_pyc*2 + zd_cubic zznd_pyc*3 ) ( 0.75 + 0.25 * zznd_pyc )*2 zdbdz_pyc(ji,jj,jk)
1256	END DO
1257	ENDIF ! lpyc
1258	ENDIF ! lconv
1259	END DO ! ji loop
1260	END DO ! jj loop
1261
1262	#ifdef key_osm_debug
1263	IF(narea==nn_narea_db) THEN
1264	ji=iloc_db; jj=jloc_db
1265	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1266	WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc contribs to ghamt/s: zsc_wth_1=',zsc_wth_1(ji,jj), ' zsc_ws_1=',zsc_ws_1(ji,jj)
1267	WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm )
1268	WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm )
1269	IF( lconv(ji,jj) ) THEN
1270	WRITE(narea+100,'(3(a,g11.3))')'Stokes + buoy + pyc contribs to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj), &
1271	&' zsc_uw_2=',zsc_uw_2(ji,jj)
1272	ELSE
1273	WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc contribs to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj)
1274	END IF
1275	WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm )
1276	WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm )
1277	WRITE(narea+100,*)
1278	FLUSH(narea+100)
1279	END IF
1280	#endif
1281
1282	IF(ln_dia_osm) THEN
1283	IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
1284	IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
1285	END IF
1286	! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
1287
1288	DO jj = 2, jpjm1
1289	DO ji = 2, jpim1
1290
1291	IF ( lconv(ji,jj) ) THEN
1292	zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) )
1293	zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) )
1294	IF ( lpyc(ji,jj) ) THEN
1295	! Pycnocline scales
1296	zsc_wth_pyc(ji,jj) = -0.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj)
1297	zsc_ws_pyc(ji,jj) = -0.003 * zwstrc(ji,jj) * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj)
1298	ENDIF
1299	ELSE
1300	zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj)
1301	zsc_ws_1(ji,jj) = zws0(ji,jj)
1302	ENDIF
1303	END DO
1304	END DO
1305
1306	DO jj = 2, jpjm1
1307	DO ji = 2, jpim1
1308	IF ( lconv(ji,jj) ) THEN
1309	DO jk = 2, imld(ji,jj)
1310	zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
1311	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) &
1312	& * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) &
1313	& - EXP( - 6.0 * zznd_ml ) ) ) &
1314	& * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) )
1315	!
1316	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) &
1317	& * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) &
1318	& - EXP( - 6.0 * zznd_ml ) ) ) &
1319	& * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
1320	END DO
1321	!
1322	! may need to comment out lpyc block
1323	IF ( lpyc(ji,jj) ) THEN
1324	! pycnocline
1325	DO jk = imld(ji,jj), ibld(ji,jj)
1326	zznd_pyc = - ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
1327	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 ) )
1328	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 ) )
1329	END DO
1330	ENDIF
1331	ELSE
1332	IF( zdhdt(ji,jj) > 0. ) THEN
1333	DO jk = 2, ibld(ji,jj)
1334	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1335	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1336	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
1337	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
1338	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
1339	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
1340	END DO
1341	ENDIF
1342	ENDIF
1343	ENDDO ! ji loop
1344	END DO ! jj loop
1345
1346	WHERE ( lconv(2:jpim1,2:jpjm1) )
1347	zsc_uw_1(2:jpim1,2:jpjm1) = zustar(2:jpim1,2:jpjm1)**2
1348	zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1) * zhml(2:jpim1,2:jpjm1)
1349	ELSEWHERE
1350	zsc_uw_1(2:jpim1,2:jpjm1) = zustar(2:jpim1,2:jpjm1)**2
1351	zsc_uw_2(2:jpim1,2:jpjm1) = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1(2:jpim1,2:jpjm1)
1352	zsc_vw_1(2:jpim1,2:jpjm1) = ff_t(2:jpim1,2:jpjm1) * zustke(2:jpim1,2:jpjm1) * zhbl(2:jpim1,2:jpjm1)
1353	zsc_vw_2(2:jpim1,2:jpjm1) = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )*2 ) zsc_vw_1(2:jpim1,2:jpjm1)
1354	ENDWHERE
1355
1356	DO jj = 2, jpjm1
1357	DO ji = 2, jpim1
1358	IF ( lconv(ji,jj) ) THEN
1359	DO jk = 2, imld(ji,jj)
1360	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
1361	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1362	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
1363	& + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml*4 ) - EXP ( -8.0 zznd_ml ) ) * zsc_uw_1(ji,jj)
1364	!
1365	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1366	& + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
1367	END DO
1368
1369	ELSE
1370	DO jk = 2, ibld(ji,jj)
1371	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1372	zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
1373	IF ( zznd_d <= 2.0 ) THEN
1374	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
1375	&* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
1376	!
1377	ELSE
1378	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
1379	& + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
1380	!
1381	ENDIF
1382
1383	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1384	& + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d*2 ) zsc_vw_1(ji,jj)
1385	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
1386	& + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
1387	END DO
1388	ENDIF
1389	END DO
1390	END DO
1391	#ifdef key_osm_debug
1392	IF(narea==nn_narea_db) THEN
1393	ji=iloc_db; jj=jloc_db
1394	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1395	WRITE(narea+100,'(2(a,g11.3))')'Stokes + buoy + pyc + transport contribs to ghamt/s: zsc_wth_1=',zsc_wth_1(ji,jj), ' zsc_ws_1=',zsc_ws_1(ji,jj)
1396	IF (lpyc(ji,jj)) WRITE(narea+100,'(2(a,g11.3))') 'zsc_wth_pyc=', zsc_wth_pyc(ji,jj), ' zsc_wth_pyc=',zsc_wth_pyc(ji,jj)
1397	WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm )
1398	WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm )
1399	IF( lconv(ji,jj) ) THEN
1400	WRITE(narea+100,'(2(a,g11.3))')'Unstable; transport contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj)
1401	ELSE
1402	WRITE(narea+100,'(3(a,g11.3))')'Stable; transport contrib to ghamu/v: zsc_uw_1=',zsc_uw_1(ji,jj), ' zsc_vw_1=',zsc_vw_1(ji,jj), &
1403	&' zsc_uw_2=',zsc_uw_2(ji,jj)
1404	END IF
1405	WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm )
1406	WRITE(narea+100,*)
1407	FLUSH(narea+100)
1408	END IF
1409	#endif
1410
1411	IF(ln_dia_osm) THEN
1412	IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
1413	IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
1414	IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
1415	IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
1416	IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
1417	IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
1418	END IF
1419	!
1420	! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
1421
1422
1423	! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer.
1424
1425	DO jj = 2, jpjm1
1426	DO ji = 2, jpim1
1427	IF ( .not. lconv(ji,jj) ) THEN
1428	DO jk = 2, ibld(ji,jj)
1429	znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
1430	IF ( znd >= 0.0 ) THEN
1431	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
1432	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
1433	ELSE
1434	ghamu(ji,jj,jk) = 0._wp
1435	ghamv(ji,jj,jk) = 0._wp
1436	ENDIF
1437	END DO
1438	ENDIF
1439	END DO
1440	END DO
1441
1442	! pynocline contributions
1443	DO jj = 2, jpjm1
1444	DO ji = 2, jpim1
1445	IF ( .not. lconv(ji,jj) ) THEN
1446	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
1447	DO jk= 2, ibld(ji,jj)
1448	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
1449	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
1450	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
1451	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
1452	END DO
1453	END IF
1454	END IF
1455	END DO
1456	END DO
1457	IF(ln_dia_osm) THEN
1458	IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
1459	IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
1460	END IF
1461
1462	DO jj=2, jpjm1
1463	DO ji = 2, jpim1
1464	ghamt(ji,jj,ibld(ji,jj)) = 0._wp
1465	ghams(ji,jj,ibld(ji,jj)) = 0._wp
1466	ghamu(ji,jj,ibld(ji,jj)) = 0._wp
1467	ghamv(ji,jj,ibld(ji,jj)) = 0._wp
1468	END DO ! ji loop
1469	END DO ! jj loop
1470	#ifdef key_osm_debug
1471	IF(narea==nn_narea_db) THEN
1472	ji=iloc_db; jj=jloc_db
1473	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1474	WRITE(narea+100,'(a)')'Tweak gham[uv] to go to zero near surface, add pycnocline viscosity/diffusivity & set=0 at ibld'
1475	WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm )
1476	WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm )
1477	WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm )
1478	WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm )
1479	WRITE(narea+100,*)
1480	FLUSH(narea+100)
1481	END IF
1482	#endif
1483
1484	IF(ln_dia_osm) THEN
1485	IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
1486	IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
1487	IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
1488	IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
1489	IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
1490	END IF
1491	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
1492	! Need to put in code for contributions that are applied explicitly to
1493	! the prognostic variables
1494	! 1. Entrainment flux
1495	!
1496	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1497
1498
1499
1500	! rotate non-gradient velocity terms back to model reference frame
1501
1502	DO jj = 2, jpjm1
1503	DO ji = 2, jpim1
1504	DO jk = 2, ibld(ji,jj)
1505	ztemp = ghamu(ji,jj,jk)
1506	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
1507	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
1508	END DO
1509	END DO
1510	END DO
1511
1512	IF(ln_dia_osm) THEN
1513	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1514	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
1515	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
1516	END IF
1517
1518	! KPP-style Ri# mixing
1519	IF( ln_kpprimix) THEN
1520	DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form)
1521	DO jj = 1, jpjm1
1522	DO ji = 1, jpim1 ! vector opt.
1523	z3du(ji,jj,jk) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) &
1524	& * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) &
1525	& / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
1526	z3dv(ji,jj,jk) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) &
1527	& * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) &
1528	& / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
1529	END DO
1530	END DO
1531	END DO
1532	!
1533	DO jk = 2, jpkm1
1534	DO jj = 2, jpjm1
1535	DO ji = 2, jpim1 ! vector opt.
1536	! ! shear prod. at w-point weightened by mask
1537	zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
1538	& + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
1539	! ! local Richardson number
1540	zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
1541	zfri = MIN( zri / rn_riinfty , 1.0_wp )
1542	zfri = ( 1.0_wp - zfri * zfri )
1543	zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk)
1544	END DO
1545	END DO
1546	END DO
1547
1548	DO jj = 2, jpjm1
1549	DO ji = 2, jpim1
1550	DO jk = ibld(ji,jj) + 1, jpkm1
1551	zdiffut(ji,jj,jk) = MAX(zdiffut(ji,jj,jk), zrimix(ji,jj,jk)*rn_difri)
1552	zviscos(ji,jj,jk) = MAX(zviscos(ji,jj,jk), zrimix(ji,jj,jk)*rn_difri)
1553	END DO
1554	END DO
1555	END DO
1556
1557	END IF ! ln_kpprimix = .true.
1558
1559	! KPP-style set diffusivity large if unstable below BL
1560	IF( ln_convmix) THEN
1561	DO jj = 2, jpjm1
1562	DO ji = 2, jpim1
1563	DO jk = ibld(ji,jj) + 1, jpkm1
1564	IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = MAX(rn_difconv,zdiffut(ji,jj,jk))
1565	END DO
1566	END DO
1567	END DO
1568	END IF ! ln_convmix = .true.
1569	#ifdef key_osm_debug
1570	IF(narea==nn_narea_db) THEN
1571	ji=iloc_db; jj=jloc_db
1572	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1573	WRITE(narea+100,'(a)') ' After including KPP Ri# diffusivity & viscosity'
1574	WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm )
1575	WRITE(narea+100,'(a,*(g11.3))') ' zviscos[imld-1..ibld+2] =', ( zviscos(ji,jj,jk), jk=jl,jm )
1576	WRITE(narea+100,*)
1577	FLUSH(narea+100)
1578	END IF
1579	#endif
1580
1581
1582
1583	IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing
1584	DO jj = 2 , jpjm1
1585	DO ji = 2, jpim1
1586	IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer
1587	! Calculate MLE flux contribution from surface fluxes
1588	DO jk = 1, ibld(ji,jj)
1589	znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln)
1590	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - ( zwth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd )
1591	ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd )
1592	END DO
1593	DO jk = 1, mld_prof(ji,jj)
1594	znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1595	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + ( zwth0(ji,jj) - zrad0(ji,jj) + zradh(ji,jj) ) * ( 1.0 - znd )
1596	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd )
1597	END DO
1598	! Viscosity for MLEs
1599	DO jk = 1, mld_prof(ji,jj)
1600	znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1601	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 )
1602	END DO
1603	ELSE
1604	! Surface transports limited to OSBL.
1605	! Viscosity for MLEs
1606	DO jk = 1, mld_prof(ji,jj)
1607	znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
1608	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 )
1609	END DO
1610	ENDIF
1611	END DO
1612	END DO
1613	#ifdef key_osm_debug
1614	IF(narea==nn_narea_db) THEN
1615	ji=iloc_db; jj=jloc_db
1616	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1617	WRITE(narea+100,'(a)') ' After including FK diffusivity & non-local terms'
1618	WRITE(narea+100,'(a,*(g11.3))') ' zdiffut[imld-1..ibld+2] =', ( zdiffut(ji,jj,jk), jk=jl,jm )
1619	WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm )
1620	WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm )
1621	WRITE(narea+100,*)
1622	FLUSH(narea+100)
1623	END IF
1624	#endif
1625	ENDIF
1626
1627	IF(ln_dia_osm) THEN
1628	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1629	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
1630	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
1631	END IF
1632
1633
1634	! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
1635	!CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
1636
1637	! GN 25/8: need to change tmask --> wmask
1638
1639	DO jk = 2, jpkm1
1640	DO jj = 2, jpjm1
1641	DO ji = 2, jpim1
1642	p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1643	p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
1644	END DO
1645	END DO
1646	END DO
1647	! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids
1648	CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1., &
1649	& ghamu, 'W', 1. , ghamv, 'W', 1. )
1650	DO jk = 2, jpkm1
1651	DO jj = 2, jpjm1
1652	DO ji = 2, jpim1
1653	ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
1654	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1655
1656	ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
1657	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1658
1659	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk)
1660	ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk)
1661	END DO
1662	END DO
1663	END DO
1664	! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)
1665	CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
1666	! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged)
1667	! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged)
1668	CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., &
1669	& ghamu, 'U', -1. , ghamv, 'V', -1. )
1670	#ifdef key_osm_debug
1671	IF(narea==nn_narea_db) THEN
1672	ji=iloc_db; jj=jloc_db
1673	jl = imld(ji,jj) - 1; jm = MIN(ibld(ji,jj) + 2, mbkt(ji,jj) )
1674	WRITE(narea+100,'(a)') ' Final diffusivity & viscosity, & non-local terms'
1675	WRITE(narea+100,'(a,*(g11.3))') ' p_avt[imld-1..ibld+2] =', ( p_avt(ji,jj,jk), jk=jl,jm )
1676	WRITE(narea+100,'(a,*(g11.3))') ' p_avm[imld-1..ibld+2] =', ( p_avm(ji,jj,jk), jk=jl,jm )
1677	WRITE(narea+100,'(a,*(g11.3))') ' ghamt[imld-1..ibld+2] =', ( ghamt(ji,jj,jk), jk=jl,jm )
1678	WRITE(narea+100,'(a,*(g11.3))') ' ghams[imld-1..ibld+2] =', ( ghams(ji,jj,jk), jk=jl,jm )
1679	WRITE(narea+100,'(a,*(g11.3))') ' ghamu[imld-1..ibld+2] =', ( ghamu(ji,jj,jk), jk=jl,jm )
1680	WRITE(narea+100,'(a,*(g11.3))') ' ghamv[imld-1..ibld+2] =', ( ghamv(ji,jj,jk), jk=jl,jm )
1681	WRITE(narea+100,*)
1682	FLUSH(narea+100)
1683	END IF
1684	#endif
1685
1686	IF(ln_dia_osm) THEN
1687	SELECT CASE (nn_osm_wave)
1688	! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1).
1689	CASE(0:1)
1690	IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)zustkezcos_wind ) ! x surface Stokes drift
1691	IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)zustkezsin_wind ) ! y surface Stokes drift
1692	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2zustke )
1693	! Stokes drift read in from sbcwave (=2).
1694	CASE(2:3)
1695	IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift
1696	IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift
1697	IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period
1698	IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height
1699	IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.rpi1.026/(0.877grav) )wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum
1700	IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)wndm2tmask(:,:,1) ) ! significant wave height from NP spectrum
1701	IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10
1702	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2 &
1703	& SQRT(ut0sd2 + vt0sd2 ) )
1704	END SELECT
1705	IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL>
1706	IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL>
1707	IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL>
1708	IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL>
1709	IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! <Tw_0>
1710	IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0>
1711	IF ( iom_use("zwb0") ) CALL iom_put( "zwb0", tmask(:,:,1)*zwb0 ) ! <Sw_0>
1712	IF ( iom_use("zwbav") ) CALL iom_put( "zwbav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb buoyancy flux
1713	IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth
1714	IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k
1715	IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base
1716	IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base
1717	IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base
1718	IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base
1719	IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base
1720	IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth
1721	IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth
1722	IF ( iom_use("zdt_ml") ) CALL iom_put( "zdt_ml", tmask(:,:,1)*zdt_ml ) ! dt at ml base
1723	IF ( iom_use("zds_ml") ) CALL iom_put( "zds_ml", tmask(:,:,1)*zds_ml ) ! ds at ml base
1724	IF ( iom_use("zdb_ml") ) CALL iom_put( "zdb_ml", tmask(:,:,1)*zdb_ml ) ! db at ml base
1725	IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth
1726	IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points
1727	IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc ) ! convective velocity scale
1728	IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale
1729	IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale
1730	IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale
1731	IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir #
1732	IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.rau0tmask(:,:,1)zustar*3 ) ! BL depth internal to zdf_osm routine
1733	IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.rau0tmask(:,:,1)zustar2zustke )
1734	IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine
1735	IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine
1736	IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine
1737	IF ( iom_use("jp_ext") ) CALL iom_put( "jp_ext", tmask(:,:,1)*jp_ext ) ! =1 if pycnocline resolved internal to zdf_osm routine
1738	IF ( iom_use("j_ddh") ) CALL iom_put( "j_ddh", tmask(:,:,1)*j_ddh ) ! index forpyc thicknessh internal to zdf_osm routine
1739	IF ( iom_use("zshear") ) CALL iom_put( "zshear", tmask(:,:,1)*zshear ) ! shear production of TKE internal to zdf_osm routine
1740	IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine
1741	IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine
1742	IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux
1743	IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux
1744	IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux
1745	IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML
1746
1747	IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth
1748	IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth
1749	IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux
1750	IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML
1751	IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k
1752	IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt
1753	IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt
1754	IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt
1755	IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt
1756	IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt
1757	IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt
1758	IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1759	IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1760
1761	END IF
1762
1763	CONTAINS
1764	! subroutine code changed, needs syntax checking.
1765	SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
1766
1767	!!---------------------------------------------------------------------
1768	!! * ROUTINE zdf_osm_diffusivity_viscosity *
1769	!!
1770	!! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline.
1771	!!
1772	!! ** Method :
1773	!!
1774	!! !!----------------------------------------------------------------------
1775	REAL(wp), DIMENSION(:,:,:) :: zdiffut
1776	REAL(wp), DIMENSION(:,:,:) :: zviscos
1777	! local
1778
1779	! Scales used to calculate eddy diffusivity and viscosity profiles
1780	REAL(wp), DIMENSION(jpi,jpj) :: pdifml_sc, pvisml_sc
1781	REAL(wp), DIMENSION(jpi,jpj) :: pdifpyc_n_sc, pdifpyc_s_sc, pdifpyc_shr
1782	REAL(wp), DIMENSION(jpi,jpj) :: pvispyc_n_sc, pvispyc_s_sc,pvispyc_shr
1783	REAL(wp), DIMENSION(jpi,jpj) :: pbeta_d_sc, pbeta_v_sc
1784	REAL(wp), DIMENSION(jpi,jpj) :: pb_coup, pc_coup_vis, pc_coup_dif
1785	!
1786	REAL(wp) :: pvel_sc_pyc, pvel_sc_ml, pstab_fac, pz_b
1787	REAL(wp) :: pa_cubic, pb_cubic, pc_cubic, pd_cubic
1788	REAL(wp) :: pznd_ml, pznd_pyc
1789	REAL(wp) :: zmsku, zmskv
1790
1791	REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375
1792	REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142
1793	REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15
1794
1795	pb_coup(:,:) = 0._wp
1796
1797	DO jj = 2, jpjm1
1798	DO ji = 2, jpim1
1799	IF ( lconv(ji,jj) ) THEN
1800
1801	pvel_sc_pyc = ( 0.15 * zvstr(ji,jj)3 + zwstrc(ji,jj)3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird
1802	pvel_sc_ml = ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
1803	pstab_fac = ( zhml(ji,jj) / pvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )1.25 ) )2
1804
1805	pdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * pvel_sc_ml
1806	pvisml_sc(ji,jj) = rn_vis_ml * pdifml_sc(ji,jj)
1807	#ifdef key_osm_debug
1808	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1809	WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, lconv=T: zdifml_sc=',pdifml_sc(ji,jj),' zvisml_sc=',pvisml_sc(ji,jj)
1810	WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',pvel_sc_pyc,' zvel_sc_ml=',pvel_sc_ml,' pstab_fac=',pstab_fac
1811	FLUSH(narea+100)
1812	END IF
1813	#endif
1814	IF ( lpyc(ji,jj) ) THEN
1815	pdifpyc_n_sc(ji,jj) = rn_dif_pyc * pvel_sc_ml * zdh(ji,jj)
1816	pvispyc_n_sc(ji,jj) = 0.09 * pvel_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)
1817	pvispyc_n_sc(ji,jj) = rn_vis_pyc * pvel_sc_ml * zdh(ji,jj) + pvispyc_n_sc(ji,jj) * pstab_fac
1818	#ifdef key_osm_debug
1819	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1820	WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w/o shear contributions: zdifpyc_n_sc',pdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',pvispyc_n_sc(ji,jj)
1821	FLUSH(narea+100)
1822	END IF
1823	#endif
1824	IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN
1825	pdifpyc_n_sc(ji,jj) = pdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )*pthird zhbl(ji,jj)
1826	pvispyc_n_sc(ji,jj) = pvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )*pthird zhbl(ji,jj)
1827	ENDIF
1828	#ifdef key_osm_debug
1829	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1830	WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w shear contributions: zdifpyc_n_sc',pdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',pvispyc_n_sc(ji,jj)
1831	FLUSH(narea+100)
1832	END IF
1833	#endif
1834	pdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * pvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) )
1835	pvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * pvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) )
1836	#ifdef key_osm_debug
1837	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1838	WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',pdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',pvispyc_s_sc(ji,jj)
1839	FLUSH(narea+100)
1840	END IF
1841	#endif
1842	pdifpyc_s_sc(ji,jj) = 0.09 * pdifpyc_s_sc(ji,jj) * pstab_fac
1843	pvispyc_s_sc(ji,jj) = pvispyc_s_sc(ji,jj) * pstab_fac
1844	#ifdef key_osm_debug
1845	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1846	WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',pdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',pvispyc_s_sc(ji,jj)
1847	FLUSH(narea+100)
1848	END IF
1849	#endif
1850
1851	pdifpyc_s_sc(ji,jj) = MAX( pdifpyc_s_sc(ji,jj), -0.5 * pdifpyc_n_sc(ji,jj) )
1852	pvispyc_s_sc(ji,jj) = MAX( pvispyc_s_sc(ji,jj), -0.5 * pvispyc_n_sc(ji,jj) )
1853
1854	pbeta_d_sc(ji,jj) = 1.0 - ( ( pdifpyc_n_sc(ji,jj) + 1.4 * pdifpyc_s_sc(ji,jj) ) / ( pdifml_sc(ji,jj) + epsln ) )**p2third
1855	pbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( pvispyc_n_sc(ji,jj) + pvispyc_s_sc(ji,jj) ) / ( pvisml_sc(ji,jj) + epsln )
1856	ELSE
1857	pdifpyc_n_sc(ji,jj) = rn_dif_pyc * pvel_sc_ml * zdh(ji,jj) ! ag 19/03
1858	pdifpyc_s_sc(ji,jj) = 0._wp ! ag 19/03
1859	pvispyc_n_sc(ji,jj) = rn_vis_pyc * pvel_sc_ml * zdh(ji,jj) ! ag 19/03
1860	pvispyc_s_sc(ji,jj) = 0._wp ! ag 19/03
1861	IF(lcoup(ji,jj) ) THEN ! ag 19/03
1862	! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd\|ub\|
1863	! already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90
1864	! Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub)
1865	! wet-cell averaging ..
1866	zmsku = 0.5 * ( 2. - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
1867	zmskv = 0.5 * ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
1868	pb_coup(ji,jj) = 0.4 * SQRT(-rCdU_bot(ji,jj) * SQRT( ( zmsku( ub(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )*2 &
1869	& + ( zmskv( vb(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )*2 ) )
1870
1871	pz_b = -gdepw_n(ji,jj,mbkt(ji,jj) + 1 ) ! ag 19/03
1872	pc_coup_vis(ji,jj) = -0.5 * ( 0.5 * pvisml_sc(ji,jj) / zhml(ji,jj) - pb_coup(ji,jj) ) / ( zhml(ji,jj) + pz_b ) ! ag 19/03
1873	#ifdef key_osm_debug
1874	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1875	WRITE(narea+100,'(4(a,g11.3))')' lcoup = T; 1st pz_b= ', pz_b, ' pb_coup ', pb_coup(ji,jj), ' pc_coup_vis ', pc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj)
1876	WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv
1877	FLUSH(narea+100)
1878	END IF
1879	#endif
1880	! #ifdef key_osm_debug
1881	! WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), ' jjg= ', mjg(jj)
1882	! WRITE(narea+400,'(3(a,g11.3))') '1st pz_b= ', pz_b, 'pb_coup', pb_coup(ji,jj), ' pc_coup_vis', pc_coup_vis(ji,jj)
1883	! FLUSH(narea+400)
1884	! #endif
1885	pz_b = -zhml(ji,jj) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ! ag 19/03
1886	pbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( pb_coup(ji,jj) * pz_b + pc_coup_vis(ji,jj) * pz_b**2 ) / pvisml_sc(ji,jj) ! ag 19/03
1887	pbeta_d_sc(ji,jj) = 1.0 - ( ( pb_coup(ji,jj) * pz_b + pc_coup_vis(ji,jj) * pz_b2 ) / pdifml_sc(ji,jj) )p2third
1888	pc_coup_dif(ji,jj) = 0.5 * ( -pdifml_sc(ji,jj) / zhml(ji,jj) * ( 1.0 - pbeta_d_sc(ji,jj) )*1.5 + 1.5 (pdifml_sc(ji,jj) / zhml(ji,jj) )* pbeta_d_sc(ji,jj) * SQRT( 1.0 - pbeta_d_sc(ji,jj) ) -pb_coup(ji,jj) ) / pz_b ! ag 19/03
1889	#ifdef key_osm_debug
1890	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1891	WRITE(narea+100,'(2(a,g11.3))')' 2nd pz_b= ', pz_b, ' pc_coup_dif', pc_coup_dif(ji,jj)
1892	FLUSH(narea+100)
1893	END IF
1894	#endif
1895	! #ifdef key_osm_debug
1896	! WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_dif', pc_coup_dif(ji,jj)
1897	! FLUSH(narea+400)
1898	! #endif
1899	ELSE ! ag 19/03
1900	pbeta_d_sc(ji,jj) = 1.0 - ( ( pdifpyc_n_sc(ji,jj) + 1.4 * pdifpyc_s_sc(ji,jj) ) / ( pdifml_sc(ji,jj) + epsln ) )**p2third ! ag 19/03
1901	pbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( pvispyc_n_sc(ji,jj) + pvispyc_s_sc(ji,jj) ) / ( pvisml_sc(ji,jj) + epsln ) ! ag 19/03
1902	ENDIF ! ag 19/03
1903	ENDIF ! ag 19/03
1904	#ifdef key_osm_debug
1905	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1906	WRITE(narea+100,'(2(a,g11.3))')'lconv=T: zbeta_d_sc',pbeta_d_sc(ji,jj) ,' zbeta_v_sc=',pbeta_v_sc(ji,jj)
1907	WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',pdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',pvispyc_n_sc(ji,jj)
1908	WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',pdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',pvispyc_s_sc(ji,jj)
1909	FLUSH(narea+100)
1910	END IF
1911	#endif
1912	ELSE ! conv, stable
1913	pdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
1914	pvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
1915	#ifdef key_osm_debug
1916	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
1917	WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, lconv=F: zdifml_sc=',pdifml_sc(ji,jj),' zvisml_sc=',pvisml_sc(ji,jj)
1918	FLUSH(narea+100)
1919	END IF
1920	#endif
1921	END IF
1922
1923	END DO
1924	END DO
1925	!
1926	DO jj = 2, jpjm1
1927	DO ji = 2, jpim1
1928	IF ( lconv(ji,jj) ) THEN
1929	DO jk = 2, imld(ji,jj) ! mixed layer diffusivity
1930	pznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
1931	!
1932	zdiffut(ji,jj,jk) = pdifml_sc(ji,jj) * pznd_ml * ( 1.0 - pbeta_d_sc(ji,jj) * pznd_ml )**1.5
1933	!
1934	zviscos(ji,jj,jk) = pvisml_sc(ji,jj) * pznd_ml * ( 1.0 - pbeta_v_sc(ji,jj) * pznd_ml ) &
1935	& * ( 1.0 - 0.5 * pznd_ml**2 )
1936	END DO
1937
1938	! Coupling to bottom
1939
1940	IF ( lcoup(ji,jj) ) THEN ! ag 19/03
1941	DO jk = mbkt(ji,jj), imld(ji,jj), -1 ! ag 19/03
1942	pz_b = - ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mbkt(ji,jj) + 1 ) ) ! ag 19/03
1943	zviscos(ji,jj,jk) = pb_coup(ji,jj) * pz_b + pc_coup_vis(ji,jj) * pz_b**2 ! ag 19/03
1944	zdiffut(ji,jj,jk) = pb_coup(ji,jj) * pz_b + pc_coup_dif(ji,jj) * pz_b**2 ! ag 19/03
1945	END DO ! ag 19/03
1946	ENDIF ! ag 19/03
1947	! pycnocline
1948	IF ( lpyc(ji,jj) ) THEN
1949	! Diffusivity profile in the pycnocline given by cubic polynomial. Note, if lpyc TRUE can't be coupled to seabed.
1950	pa_cubic = 0.5
1951	pb_cubic = -1.75 * pdifpyc_s_sc(ji,jj) / pdifpyc_n_sc(ji,jj)
1952	pd_cubic = ( zdh(ji,jj) * pdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - pbeta_d_sc(ji,jj) ) * ( 2.5 * pbeta_d_sc(ji,jj) - 1.0 ) &
1953	& - 0.85 * pdifpyc_s_sc(ji,jj) ) / MAX(pdifpyc_n_sc(ji,jj), 1.e-8)
1954	pd_cubic = pd_cubic - pb_cubic - 2.0 * ( 1.0 - pa_cubic - pb_cubic )
1955	pc_cubic = 1.0 - pa_cubic - pb_cubic - pd_cubic
1956	DO jk = imld(ji,jj) , ibld(ji,jj)
1957	pznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
1958	!
1959	zdiffut(ji,jj,jk) = pdifpyc_n_sc(ji,jj) * ( pa_cubic + pb_cubic * pznd_pyc + pc_cubic * pznd_pyc*2 + pd_cubic pznd_pyc**3 )
1960
1961	zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + pdifpyc_s_sc(ji,jj) * ( 1.75 * pznd_pyc - 0.15 * pznd_pyc*2 - 0.2 pznd_pyc**3 )
1962	END DO
1963	! viscosity profiles.
1964	pa_cubic = 0.5
1965	pb_cubic = -1.75 * pvispyc_s_sc(ji,jj) / pvispyc_n_sc(ji,jj)
1966	pd_cubic = ( 0.5 * pvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * pvispyc_s_sc(ji,jj) ) / MAX(pvispyc_n_sc(ji,jj), 1.e-8)
1967	pd_cubic = pd_cubic - pb_cubic - 2.0 * ( 1.0 - pa_cubic - pb_cubic )
1968	pc_cubic = 1.0 - pa_cubic - pb_cubic - pd_cubic
1969	DO jk = imld(ji,jj) , ibld(ji,jj)
1970	pznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
1971	zviscos(ji,jj,jk) = pvispyc_n_sc(ji,jj) * ( pa_cubic + pb_cubic * pznd_pyc + pc_cubic * pznd_pyc*2 + pd_cubic pznd_pyc**3 )
1972	zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + pvispyc_s_sc(ji,jj) * ( 1.75 * pznd_pyc - 0.15 * pznd_pyc*2 -0.2 pznd_pyc**3 )
1973	END DO
1974	! IF ( zdhdt(ji,jj) > 0._wp ) THEN
1975	! zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 )
1976	! zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 )
1977	! ELSE
1978	! zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
1979	! zviscos(ji,jj,ibld(ji,jj)) = 0._wp
1980	! ENDIF
1981	ELSE
1982	! lpyc set false but not coupled to the bottom.
1983	IF ( .not. lcoup(ji,jj) ) THEN
1984	zdiffut(ji,jj,ibld(ji,jj)) = 0.5 * pdifpyc_n_sc(ji,jj)
1985	zviscos(ji,jj,ibld(ji,jj)) = 0.5 * pvispyc_n_sc(ji,jj)
1986	ENDIF
1987	ENDIF
1988	ELSE
1989	! stable conditions
1990	DO jk = 2, ibld(ji,jj)
1991	pznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
1992	zdiffut(ji,jj,jk) = 0.75 * pdifml_sc(ji,jj) * pznd_ml * ( 1.0 - pznd_ml )**1.5
1993	zviscos(ji,jj,jk) = 0.375 * pvisml_sc(ji,jj) * pznd_ml * (1.0 - pznd_ml) * ( 1.0 - pznd_ml**2 )
1994	END DO
1995
1996	IF ( zdhdt(ji,jj) > 0._wp ) THEN
1997	zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj))
1998	zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj))
1999	ENDIF
2000	ENDIF ! end if ( lconv )
2001	!
2002	END DO ! end of ji loop
2003	END DO ! end of jj loop
2004	IF ( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1)*pb_coup(:,:) ) ! BBL-coupling velocity scale
2005
2006	END SUBROUTINE zdf_osm_diffusivity_viscosity
2007
2008	SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear )
2009
2010	!!---------------------------------------------------------------------
2011	!! * ROUTINE zdf_osm_osbl_state *
2012	!!
2013	!! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number
2014	!!
2015	!! ** Method :
2016	!!
2017	!! !!----------------------------------------------------------------------
2018
2019	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.
2020
2021	LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear
2022
2023	REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer.
2024	REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline
2025
2026	! Local Variables
2027
2028	INTEGER :: jj, ji
2029
2030	REAL(wp), DIMENSION(jpi,jpj) :: zekman
2031	REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b ! Richardson numbers
2032	REAL(wp) :: zshear_u, zshear_v, zwb_shr
2033	REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
2034
2035	REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.8
2036	REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03
2037	REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15
2038	REAL, PARAMETER :: rn_ri_p_thresh = 27.0
2039	REAL, PARAMETER :: zri_c = 0.25
2040	REAL, PARAMETER :: zek = 4.0
2041	REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress.
2042
2043	! Determins stability and set flag lconv
2044	DO jj = 2, jpjm1
2045	DO ji = 2, jpim1
2046	IF ( zhol(ji,jj) < 0._wp ) THEN
2047	lconv(ji,jj) = .TRUE.
2048	ELSE
2049	lconv(ji,jj) = .FALSE.
2050	ENDIF
2051	END DO
2052	END DO
2053
2054	zekman(2:jpim1,2:jpjm1) = EXP( - zek * ABS( ff_t(2:jpim1,2:jpjm1) ) * zhbl(2:jpim1,2:jpjm1) / MAX(zustar(2:jpim1,2:jpjm1), 1.e-8 ) )
2055
2056	zshear(:,:) = zlarge
2057	zshear(2:jpim1,2:jpjm1) = 0._wp
2058	#ifdef key_osm_debug
2059	IF(narea==nn_narea_db) THEN
2060	ji=iloc_db; jj=jloc_db
2061	WRITE(narea+100,'(a,g11.3)') &
2062	& 'zdf_osm_osbl_state start: zekman=', zekman(ji,jj)
2063	FLUSH(narea+100)
2064	END IF
2065	#endif
2066	j_ddh(2:jpim1,2:jpjm1) = 1
2067
2068	DO jj = 2, jpjm1
2069	DO ji = 2, jpim1
2070	IF ( lconv(ji,jj) ) THEN
2071	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
2072	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 &
2073	& / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp )
2074
2075	IF ( ff_t(ji,jj) >= 0._wp ) THEN
2076	! Northern Hemisphere
2077	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 )
2078	ELSE
2079	! Southern Hemisphere
2080	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 )
2081	ENDIF
2082	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 ) )
2083	#ifdef key_osm_debug
2084	! IF(narea==nn_narea_db)THEN
2085	! WRITE(narea+100,'(2(a,i10.4))')'ji',ji,'jj',jj
2086	! WRITE(narea+100,'(2(a,i10.4))')'iloc_db',iloc_db,'jloc_db',jloc_db
2087	! WRITE(narea+100,'(2(a,i10.4))')'iloc_db+',mi0(nn_idb),'jloc_db+',mj0(nn_jdb)
2088	! FLUSH(narea+100)
2089	! END IF
2090	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2091	WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear=',zshear(ji,jj)
2092	WRITE(narea+100,'(2(a,g11.3))')'zdf_osm_osbl_state 1st zshear: zri_b=',zri_b(ji,jj),' zri_p=',zri_p(ji,jj)
2093	FLUSH(narea+100)
2094	END IF
2095	#endif
2096	! Stability Dependence
2097	zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75 * MAX(0._wp,( zri_b(ji,jj) - zri_c ) / zri_c ) )
2098	#ifdef key_osm_debug
2099	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2100	WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear inc ri part=',zshear(ji,jj)
2101	FLUSH(narea+100)
2102	END IF
2103	#endif
2104
2105	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2106	! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when !
2107	! full code available !
2108	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2109	IF ( zshear(ji,jj) > 1.e-10 ) THEN
2110	IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN(hu_n(ji,jj), hu_n(ji-1,jj), hv_n(ji,jj), hv_n(ji,jj-1))>100._wp ) THEN
2111	! Growing shear layer
2112	j_ddh(ji,jj) = 0
2113	lshear(ji,jj) = .TRUE.
2114	ELSE
2115	j_ddh(ji,jj) = 1
2116	! IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN
2117	! shear production large enough to determine layer charcteristics, but can't maintain a shear layer.
2118	lshear(ji,jj) = .TRUE.
2119	! ELSE
2120	ENDIF
2121	ELSE
2122	j_ddh(ji,jj) = 2
2123	lshear(ji,jj) = .FALSE.
2124	ENDIF
2125	! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline.
2126	! zshear(ji,jj) = 0.5 * zshear(ji,jj)
2127	! lshear(ji,jj) = .FALSE.
2128	! ENDIF
2129	ELSE ! zdb_bl test, note zshear set to zero
2130	j_ddh(ji,jj) = 2
2131	lshear(ji,jj) = .FALSE.
2132	ENDIF
2133	ENDIF
2134	END DO
2135	END DO
2136
2137	! Calculate entrainment buoyancy flux due to surface fluxes.
2138
2139	DO jj = 2, jpjm1
2140	DO ji = 2, jpim1
2141	IF ( lconv(ji,jj) ) THEN
2142	zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
2143	zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
2144	zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
2145	zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
2146	IF (nn_osm_SD_reduce > 0 ) THEN
2147	! Effective Stokes drift already reduced from surface value
2148	zr_stokes = 1.0_wp
2149	ELSE
2150	! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
2151	! requires further reduction where BL is deep
2152	zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
2153	& * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
2154	END IF
2155	zwb_ent(ji,jj) = - 2.0 * zalpha_c * zrf_conv * zwbav(ji,jj) &
2156	& - zalpha_s * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
2157	& + zr_stokes * ( zalpha_s * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
2158	& - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
2159	!
2160	#ifdef key_osm_debug
2161	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2162	WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state conv+shear0/lang: zwb_ent=',zwb_ent(ji,jj)
2163	FLUSH(narea+100)
2164	END IF
2165	#endif
2166	ENDIF
2167	END DO ! ji loop
2168	END DO ! jj loop
2169
2170	zwb_min(:,:) = zlarge
2171
2172	DO jj = 2, jpjm1
2173	DO ji = 2, jpim1
2174	IF ( lshear(ji,jj) ) THEN
2175	IF ( lconv(ji,jj) ) THEN
2176	! Unstable OSBL
2177	zwb_shr = -za_wb_s * zri_b(ji,jj) * zshear(ji,jj)
2178	#ifdef key_osm_debug
2179	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2180	WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zwb_shr: zwb_shr=',zwb_shr
2181	FLUSH(narea+100)
2182	END IF
2183	#endif
2184	IF ( j_ddh(ji,jj) == 0 ) THEN
2185
2186	! ! Developing shear layer, additional shear production possible.
2187
2188	! zshear_u = MAX( zustar(ji,jj)*2 MAX( zdu_ml(ji,jj), 0._wp ) / zhbl(ji,jj), 0._wp )
2189	! zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 )
2190	! zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )
2191
2192	! 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 )
2193	! zwb_shr = MAX( zwb_shr, -0.25 * zshear_u )
2194	#ifdef key_osm_debug
2195	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2196	WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state j_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, &
2197	& ' zshear=',zshear(ji,jj),' zshear_u=', zshear_u
2198	FLUSH(narea+100)
2199	END IF
2200	#endif
2201
2202	ENDIF
2203	zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
2204	! zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
2205	ELSE ! IF ( lconv ) THEN - ENDIF
2206	! Stable OSBL - shear production not coded for first attempt.
2207	ENDIF ! lconv
2208	ENDIF ! lshear
2209	IF ( lconv(ji,jj) ) THEN
2210	! Unstable OSBL
2211	zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2._wp * zwbav(ji,jj)
2212	ENDIF ! lconv
2213	#ifdef key_osm_debug
2214	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2215	WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',zwb_ent(ji,jj), &
2216	& ' zwb_min=',zwb_min(ji,jj), ' zwb0tot=', zwb0tot(ji,jj), ' zwbav= ', zwbav(ji,jj)
2217	FLUSH(narea+100)
2218	END IF
2219	#endif
2220	END DO ! ji
2221	END DO ! jj
2222	END SUBROUTINE zdf_osm_osbl_state
2223
2224
2225	SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv )
2226	!!---------------------------------------------------------------------
2227	!! * ROUTINE zdf_vertical_average *
2228	!!
2229	!! ** Purpose : Determines vertical averages from surface to jnlev.
2230	!!
2231	!! ** Method : Averages are calculated from the surface to jnlev.
2232	!! The external level used to calculate differences is ibld+ibld_ext
2233	!!
2234	!!----------------------------------------------------------------------
2235
2236	INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over.
2237	INTEGER, DIMENSION(jpi,jpj) :: jp_ext
2238
2239	! Alan: do we need zb?
2240	REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb ! Average temperature and salinity
2241	REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components
2242	REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
2243	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components.
2244
2245	INTEGER :: jk, ji, jj, ibld_ext
2246	REAL(wp) :: zthick, zthermal, zbeta
2247
2248
2249	zt(2:jpim1,2:jpjm1) = 0._wp
2250	zs(2:jpim1,2:jpjm1) = 0._wp
2251	zu(2:jpim1,2:jpjm1) = 0._wp
2252	zv(2:jpim1,2:jpjm1) = 0._wp
2253	DO jj = 2, jpjm1 ! Vertical slab
2254	DO ji = 2, jpim1
2255	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
2256	zbeta = rab_n(ji,jj,1,jp_sal)
2257	! average over depth of boundary layer
2258	zthick = epsln
2259	DO jk = 2, jnlev_av(ji,jj)
2260	zthick = zthick + e3t_n(ji,jj,jk)
2261	zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
2262	zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
2263	zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) &
2264	& * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
2265	& / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
2266	zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) &
2267	& * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
2268	& / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
2269	END DO
2270	zt(ji,jj) = zt(ji,jj) / zthick
2271	zs(ji,jj) = zs(ji,jj) / zthick
2272	zu(ji,jj) = zu(ji,jj) / zthick
2273	zv(ji,jj) = zv(ji,jj) / zthick
2274	zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj)
2275	ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj)
2276	IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN ! ag 09/03
2277	! Two external levels are available.
2278	zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem)
2279	zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal)
2280	zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) &
2281	& / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) )
2282	zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) &
2283	& / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) )
2284	zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
2285	ELSE
2286	zdt(ji,jj) = 0._wp
2287	zds(ji,jj) = 0._wp
2288	zdu(ji,jj) = 0._wp
2289	zdv(ji,jj) = 0._wp
2290	zdb(ji,jj) = 0._wp
2291	ENDIF
2292	END DO
2293	END DO
2294	END SUBROUTINE zdf_osm_vertical_average
2295
2296	SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
2297	!!---------------------------------------------------------------------
2298	!! * ROUTINE zdf_velocity_rotation *
2299	!!
2300	!! ** Purpose : Rotates frame of reference of averaged velocity components.
2301	!!
2302	!! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w
2303	!!
2304	!!----------------------------------------------------------------------
2305
2306	REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle
2307	REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current
2308	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline
2309
2310	INTEGER :: ji, jj
2311	REAL(wp) :: ztemp
2312
2313	DO jj = 2, jpjm1
2314	DO ji = 2, jpim1
2315	ztemp = zu(ji,jj)
2316	zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
2317	zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
2318	ztemp = zdu(ji,jj)
2319	zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
2320	zdv(ji,jj) = zdv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
2321	END DO
2322	END DO
2323	END SUBROUTINE zdf_osm_velocity_rotation
2324
2325	SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
2326	!!---------------------------------------------------------------------
2327	!! * ROUTINE zdf_osm_osbl_state_fk *
2328	!!
2329	!! ** 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.
2330	!! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined
2331	!! lflux :: determines whether effects of surface flux extend below the base of the OSBL
2332	!! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl.
2333	!!
2334	!! ** Method :
2335	!!
2336	!!
2337	!!----------------------------------------------------------------------
2338
2339	! Outputs
2340	LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle
2341	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk
2342	!
2343	REAL(wp), DIMENSION(jpi,jpj) :: znd_param
2344	REAL(wp) :: zbuoy, ztmp, zpe_mle_layer
2345	REAL(wp) :: zpe_mle_ref, zdbdz_mle_int
2346
2347	znd_param(2:jpim1,2:jpjm1) = 0._wp
2348
2349	DO jj = 2, jpjm1
2350	DO ji = 2, jpim1
2351	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
2352	zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj)
2353	END DO
2354	END DO
2355	DO jj = 2, jpjm1
2356	DO ji = 2, jpim1
2357	!
2358	IF ( lconv(ji,jj) ) THEN
2359	IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
2360	zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2361	zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2362	zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2363	zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
2364	! Calculate potential energies of actual profile and reference profile.
2365	zpe_mle_layer = 0._wp
2366	zpe_mle_ref = 0._wp
2367	zthermal = rab_n(ji,jj,1,jp_tem)
2368	zbeta = rab_n(ji,jj,1,jp_sal)
2369
2370	DO jk = ibld(ji,jj), mld_prof(ji,jj)
2371	zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) )
2372	zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
2373	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)
2374	END DO
2375	! Non-dimensional parameter to diagnose the presence of thermocline
2376
2377	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) )
2378	ENDIF
2379	ENDIF
2380	#ifdef key_osm_debug
2381	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2382	WRITE(narea+100,'(4(a,g11.3))')'start of zdf_osm_osbl_state_fk: zwb_fk=',zwb_fk(ji,jj), &
2383	& ' znd_param=',znd_param(ji,jj), ' zpe_mle_ref=', zpe_mle_ref, ' zpe_mle_layer=', zpe_mle_layer
2384	FLUSH(narea+100)
2385	END IF
2386	#endif
2387	END DO
2388	END DO
2389
2390	! Diagnosis
2391	DO jj = 2, jpjm1
2392	DO ji = 2, jpim1
2393	IF ( lconv(ji,jj) ) THEN
2394	IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent(ji,jj) > 0.5 ) THEN
2395	IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
2396	! MLE layer growing
2397	IF ( znd_param (ji,jj) > 100. ) THEN
2398	! Thermocline present
2399	lflux(ji,jj) = .FALSE.
2400	lmle(ji,jj) =.FALSE.
2401	ELSE
2402	! Thermocline not present
2403	lflux(ji,jj) = .TRUE.
2404	lmle(ji,jj) = .TRUE.
2405	ENDIF ! znd_param > 100
2406	!
2407	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
2408	lpyc(ji,jj) = .FALSE.
2409	ELSE
2410	lpyc(ji,jj) = .TRUE.
2411	ENDIF
2412	ELSE
2413	! MLE layer restricted to OSBL or just below.
2414	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
2415	! Weak stratification MLE layer can grow.
2416	lpyc(ji,jj) = .FALSE.
2417	lflux(ji,jj) = .TRUE.
2418	lmle(ji,jj) = .TRUE.
2419	ELSE
2420	! Strong stratification
2421	lpyc(ji,jj) = .TRUE.
2422	lflux(ji,jj) = .FALSE.
2423	lmle(ji,jj) = .FALSE.
2424	ENDIF ! zdb_bl < rn_mle_thresh_bl and
2425	ENDIF ! zhmle > 1.2 zhbl
2426	ELSE
2427	lpyc(ji,jj) = .TRUE.
2428	lflux(ji,jj) = .FALSE.
2429	lmle(ji,jj) = .FALSE.
2430	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
2431	ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5
2432	ELSE
2433	! Stable Boundary Layer
2434	lpyc(ji,jj) = .FALSE.
2435	lflux(ji,jj) = .FALSE.
2436	lmle(ji,jj) = .FALSE.
2437	ENDIF ! lconv
2438	#ifdef key_osm_debug
2439	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2440	WRITE(narea+100,'(3(a,g11.3),/,4(a,l2))')'end of zdf_osm_osbl_state_fk:zwb_ent=',zwb_ent(ji,jj), &
2441	& ' zhmle=',zhmle(ji,jj), ' zhbl=', zhbl(ji,jj), &
2442	& ' lpyc= ', lpyc(ji,jj), ' lflux= ', lflux(ji,jj), ' lmle= ', lmle(ji,jj), ' lconv= ', lconv(ji,jj)
2443	FLUSH(narea+100)
2444	END IF
2445	#endif
2446	END DO
2447	END DO
2448	END SUBROUTINE zdf_osm_osbl_state_fk
2449
2450	SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz )
2451	!!---------------------------------------------------------------------
2452	!! * ROUTINE zdf_osm_external_gradients *
2453	!!
2454	!! ** Purpose : Calculates the gradients below the OSBL
2455	!!
2456	!! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient.
2457	!!
2458	!!----------------------------------------------------------------------
2459
2460	INTEGER, DIMENSION(jpi,jpj) :: jbase
2461	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy.
2462
2463	INTEGER :: jj, ji, jkb, jkb1
2464	REAL(wp) :: zthermal, zbeta
2465
2466
2467	DO jj = 2, jpjm1
2468	DO ji = 2, jpim1
2469	IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
2470	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
2471	zbeta = rab_n(ji,jj,1,jp_sal)
2472	jkb = jbase(ji,jj)
2473	jkb1 = MIN(jkb + 1, mbkt(ji,jj))
2474	zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
2475	& / e3w_n(ji,jj,jkb1)
2476	zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
2477	& / e3w_n(ji,jj,jkb1)
2478	zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
2479	ELSE
2480	zdtdz(ji,jj) = 0._wp
2481	zdsdz(ji,jj) = 0._wp
2482	zdbdz(ji,jj) = 0._wp
2483	END IF
2484	END DO
2485	END DO
2486	END SUBROUTINE zdf_osm_external_gradients
2487
2488	SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha )
2489
2490	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline
2491	REAL(wp), DIMENSION(jpi,jpj) :: zalpha
2492
2493	INTEGER :: jk, jj, ji
2494	REAL(wp) :: ztgrad, zsgrad, zbgrad
2495	REAL(wp) :: zgamma_b_nd, znd
2496	REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc
2497	REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
2498
2499	DO jj = 2, jpjm1
2500	DO ji = 2, jpim1
2501	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
2502	IF ( lconv(ji,jj) ) THEN ! convective conditions
2503	IF ( lpyc(ji,jj) ) THEN
2504	zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
2505	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 ) )
2506	zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp )
2507
2508	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
2509	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2510	! Commented lines in this section are not needed in new code, once tested !
2511	! can be removed !
2512	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2513	! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
2514	! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
2515	zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)
2516	zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
2517	DO jk = 2, ibld(ji,jj)
2518	znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp
2519	IF ( znd <= zzeta_m ) THEN
2520	! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &
2521	! & EXP( -6.0 * ( znd -zzeta_m )**2 )
2522	! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &
2523	! & EXP( -6.0 * ( znd -zzeta_m )**2 )
2524	zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &
2525	& EXP( -6.0 * ( znd -zzeta_m )**2 )
2526	ELSE
2527	! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
2528	! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
2529	zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
2530	ENDIF
2531	END DO
2532	#ifdef key_osm_debug
2533	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2534	WRITE(narea+100,'(a,/,3(a,g11.3),/,2(a,g11.3),/)')'end of zdf_osm_pycnocline_scalar_profiles:lconv=lpyc=T',&
2535	& 'zzeta_m=', zzeta_m, ' zalpha=', zalpha(ji,jj), ' ztmp=', ztmp,&
2536	& ' zbgrad=', zbgrad, ' zgamma_b_nd=', zgamma_b_nd
2537	FLUSH(narea+100)
2538	END IF
2539	#endif
2540	ENDIF ! if no pycnocline pycnocline gradients set to zero
2541	ELSE
2542	! stable conditions
2543	! if pycnocline profile only defined when depth steady of increasing.
2544	IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady.
2545	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
2546	IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline
2547	ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
2548	ztgrad = zdt_bl(ji,jj) * ztmp
2549	zsgrad = zds_bl(ji,jj) * ztmp
2550	zbgrad = zdb_bl(ji,jj) * ztmp
2551	DO jk = 2, ibld(ji,jj)
2552	znd = gdepw_n(ji,jj,jk) * ztmp
2553	zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
2554	zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
2555	zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
2556	END DO
2557	ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
2558	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
2559	ztgrad = zdt_bl(ji,jj) * ztmp
2560	zsgrad = zds_bl(ji,jj) * ztmp
2561	zbgrad = zdb_bl(ji,jj) * ztmp
2562	DO jk = 2, ibld(ji,jj)
2563	znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp
2564	zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
2565	zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
2566	zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
2567	END DO
2568	ENDIF ! IF (zhol >=0.5)
2569	#ifdef key_osm_debug
2570	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2571	WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_pycnocline_scalar_profiles:lconv=F ztgrad=',&
2572	& ztgrad, ' zsgrad=', zsgrad, ' zbgrad=', zbgrad
2573	FLUSH(narea+100)
2574	END IF
2575	#endif
2576	ENDIF ! IF (zdb_bl> 0.)
2577	ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
2578	ENDIF ! IF (lconv)
2579	ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) )
2580	END DO
2581	END DO
2582
2583	END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
2584
2585	SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
2586	!!---------------------------------------------------------------------
2587	!! * ROUTINE zdf_osm_pycnocline_shear_profiles *
2588	!!
2589	!! ** Purpose : Calculates velocity shear in the pycnocline
2590	!!
2591	!! ** Method :
2592	!!
2593	!!----------------------------------------------------------------------
2594
2595	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
2596
2597	INTEGER :: jk, jj, ji
2598	REAL(wp) :: zugrad, zvgrad, znd
2599	REAL(wp) :: zzeta_v = 0.45
2600	!
2601	DO jj = 2, jpjm1
2602	DO ji = 2, jpim1
2603	!
2604	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
2605	IF ( lconv (ji,jj) ) THEN
2606	! Unstable conditions. Shouldn;t be needed with no pycnocline code.
2607	! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
2608	! & ( ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zhml(ji,jj) ) * &
2609	! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
2610	!Alan is this right?
2611	! zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
2612	! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
2613	! & ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird + epsln ) &
2614	! & )/ (zdh(ji,jj) + epsln )
2615	! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
2616	! znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
2617	! IF ( znd <= 0.0 ) THEN
2618	! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
2619	! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
2620	! ELSE
2621	! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
2622	! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
2623	! ENDIF
2624	! END DO
2625	ELSE
2626	! stable conditions
2627	zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
2628	zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
2629	DO jk = 2, ibld(ji,jj)
2630	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
2631	IF ( znd < 1.0 ) THEN
2632	zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
2633	ELSE
2634	zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
2635	ENDIF
2636	zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
2637	END DO
2638	ENDIF
2639	!
2640	END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
2641	END DO
2642	END DO
2643	END SUBROUTINE zdf_osm_pycnocline_shear_profiles
2644
2645	SUBROUTINE zdf_osm_calculate_dhdt( zdhdt )
2646	!!---------------------------------------------------------------------
2647	!! * ROUTINE zdf_osm_calculate_dhdt *
2648	!!
2649	!! ** Purpose : Calculates the rate at which hbl changes.
2650	!!
2651	!! ** Method :
2652	!!
2653	!!----------------------------------------------------------------------
2654
2655	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! Rate of change of hbl
2656
2657	INTEGER :: jj, ji
2658	REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi
2659	REAL(wp) :: zvel_max,zddhdt
2660	REAL(wp), PARAMETER :: zzeta_m = 0.3
2661	REAL(wp), PARAMETER :: zgamma_c = 2.0
2662	REAL(wp), PARAMETER :: zdhoh = 0.1
2663	REAL(wp), PARAMETER :: zalpha_b = 0.3
2664	REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
2665
2666	DO jj = 2, jpjm1
2667	DO ji = 2, jpim1
2668
2669	IF ( lshear(ji,jj) ) THEN
2670	IF ( lconv(ji,jj) ) THEN ! Convective
2671
2672	IF ( ln_osm_mle ) THEN
2673
2674	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
2675	! Fox-Kemper buoyancy flux average over OSBL
2676	zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * &
2677	(1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
2678	ELSE
2679	zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
2680	ENDIF
2681	zvel_max = ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
2682	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
2683	! OSBL is deepening, entrainment > restratification
2684	IF ( zdb_bl(ji,jj) > 1.0e-15 ) THEN
2685	zgamma_b_nd = MAX( zdbdz_bl_ext(ji,jj), 0._wp ) * zdh(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
2686	zpsi = ( 1.0 - 0.5 * zdh(ji,jj) / zhbl(ji,jj) ) * ( zwb0(ji,jj) - MIN( ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ), 0._wp ) ) * zdh(ji,jj) / zhbl(ji,jj)
2687	zpsi = zpsi + 1.75 * ( 1.0 - 0.5 * zdh(ji,jj) / zhbl(ji,jj) )( zdh(ji,jj) / zhbl(ji,jj) + zgamma_b_nd ) MIN( ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ), 0._wp )
2688	zpsi = zalpha_b * MAX ( zpsi, 0._wp )
2689	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 / ( zvel_max + MAX( zdb_bl(ji,jj), 1.e-15 ) )
2690	#ifdef key_osm_debug
2691	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2692	WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt, OSBL is deepening, entrainment > restratification: zdhdt=',zdhdt(ji,jj)
2693	WRITE(narea+100,'(3(a,g11.3))') ' zpsi=',zpsi, ' zgamma_b_nd=', zgamma_b_nd, ' zdh=', zdh(ji,jj)
2694	FLUSH(narea+100)
2695	END IF
2696	#endif
2697	IF ( j_ddh(ji,jj) == 1 ) THEN
2698	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
2699	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 )
2700	ELSE
2701	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 )
2702	ENDIF
2703	! Relaxation to dh_ref = zari * hbl
2704	zddhdt = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
2705	#ifdef key_osm_debug
2706	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2707	WRITE(narea+100,'(a,g11.3)')'Inside 1st major loop of zdf_osm_calculate_dhdt,j_ddh(ji,jj) == 1: zari=',zari
2708	FLUSH(narea+100)
2709	END IF
2710	#endif
2711
2712	ELSE IF ( j_ddh(ji,jj) == 0 ) THEN
2713	! Growing shear layer
2714	zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
2715	zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
2716	ELSE
2717	zddhdt = 0._wp
2718	ENDIF ! j_ddh
2719	zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( 1.0 -0.5 * zdh(ji,jj) / zhbl(ji,jj) ) * &
2720	& zdb_ml(ji,jj) * MAX(zddhdt,0._wp) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
2721	ELSE ! zdb_bl >0
2722	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)
2723	ENDIF
2724	ELSE ! zwb_min + 2*zwb_fk_b < 0
2725	! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
2726	zdhdt(ji,jj) = - MIN(zvel_mle(ji,jj), hbl(ji,jj)/10800.)
2727
2728
2729	ENDIF
2730
2731	ELSE
2732	! Fox-Kemper not used.
2733
2734	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) / &
2735	& MAX((zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln)
2736	zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
2737	! added ajgn 23 July as temporay fix
2738
2739	ENDIF ! ln_osm_mle
2740
2741	ELSE ! lconv - Stable
2742	zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
2743	IF ( zdhdt(ji,jj) < 0._wp ) THEN
2744	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
2745	zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
2746	ELSE
2747	zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
2748	ENDIF
2749	zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
2750	zdhdt(ji,jj) = MAX(zdhdt(ji,jj), -hbl(ji,jj)/5400.)
2751	ENDIF ! lconv
2752	ELSE ! lshear
2753	IF ( lconv(ji,jj) ) THEN ! Convective
2754
2755	IF ( ln_osm_mle ) THEN
2756
2757	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
2758	! Fox-Kemper buoyancy flux average over OSBL
2759	zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * &
2760	(1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
2761	ELSE
2762	zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
2763	ENDIF
2764	zvel_max = ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
2765	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
2766	! OSBL is deepening, entrainment > restratification
2767	IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
2768	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
2769	ELSE
2770	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)
2771	ENDIF
2772	ELSE
2773	! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
2774	zdhdt(ji,jj) = - MIN(zvel_mle(ji,jj), hbl(ji,jj)/10800.)
2775
2776
2777	ENDIF
2778
2779	ELSE
2780	! Fox-Kemper not used.
2781
2782	zvel_max = -zwb_ent(ji,jj) / &
2783	& MAX((zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln)
2784	zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
2785	! added ajgn 23 July as temporay fix
2786
2787	ENDIF ! ln_osm_mle
2788
2789	ELSE ! Stable
2790	zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
2791	IF ( zdhdt(ji,jj) < 0._wp ) THEN
2792	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
2793	zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj)
2794	ELSE
2795	zpert = MAX( zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
2796	ENDIF
2797	zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
2798	zdhdt(ji,jj) = MAX(zdhdt(ji,jj), -hbl(ji,jj)/5400.)
2799	ENDIF ! lconv
2800	ENDIF ! lshear
2801	#ifdef key_osm_debug
2802	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2803	WRITE(narea+100,'(4(a,g11.3))')'end of 1st major loop of zdf_osm_calculate_dhdt: zdhdt=',zdhdt(ji,jj), &
2804	& ' zpert=', zpert, ' zddhdt=', zddhdt, ' zvel_max=', zvel_max
2805
2806	IF ( ln_osm_mle ) THEN
2807	WRITE(narea+100,'(3(a,g11.3),/)') 'zvel_mle=',zvel_mle(ji,jj), ' zwb_fk_b=', zwb_fk_b(ji,jj), &
2808	& ' zwb_ent + 2zwb_fk_b =', zwb_ent(ji,jj) + 2.0 zwb_fk_b(ji,jj)
2809	FLUSH(narea+100)
2810	END IF
2811	END IF
2812	#endif
2813	END DO
2814	END DO
2815	END SUBROUTINE zdf_osm_calculate_dhdt
2816
2817	SUBROUTINE zdf_osm_timestep_hbl( zdhdt )
2818	!!---------------------------------------------------------------------
2819	!! * ROUTINE zdf_osm_timestep_hbl *
2820	!!
2821	!! ** Purpose : Increments hbl.
2822	!!
2823	!! ** Method : If thechange in hbl exceeds one model level the change is
2824	!! is calculated by moving down the grid, changing the buoyancy
2825	!! jump. This is to ensure that the change in hbl does not
2826	!! overshoot a stable layer.
2827	!!
2828	!!----------------------------------------------------------------------
2829
2830
2831	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! rates of change of hbl.
2832
2833	INTEGER :: jk, jj, ji, jm
2834	REAL(wp) :: zhbl_s, zvel_max, zdb
2835	REAL(wp) :: zthermal, zbeta
2836
2837	DO jj = 2, jpjm1
2838	DO ji = 2, jpim1
2839	#ifdef key_osm_debug
2840	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2841	WRITE(narea+100,'(2(a,i7))')'start of zdf_osm_timestep_hbl: old ibld=',imld(ji,jj),' trial ibld=', ibld(ji,jj)
2842	FLUSH(narea+100)
2843	END IF
2844	#endif
2845	IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
2846	!
2847	! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
2848	!
2849	zhbl_s = hbl(ji,jj)
2850	jm = imld(ji,jj)
2851	zthermal = rab_n(ji,jj,1,jp_tem)
2852	zbeta = rab_n(ji,jj,1,jp_sal)
2853
2854
2855	IF ( lconv(ji,jj) ) THEN
2856	!unstable
2857
2858	IF( ln_osm_mle ) THEN
2859	zvel_max = ( zwstrl(ji,jj)3 + zwstrc(ji,jj)3 )**p2third / hbl(ji,jj)
2860	ELSE
2861
2862	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) / &
2863	& ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
2864
2865	ENDIF
2866	#ifdef key_osm_debug
2867	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2868	WRITE(narea+100,'(a,g11.3)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=T: zvel_max=',zvel_max
2869	FLUSH(narea+100)
2870	END IF
2871	#endif
2872
2873	DO jk = imld(ji,jj), ibld(ji,jj)
2874	zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
2875	& - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
2876	& 0.0 ) + zvel_max
2877
2878
2879	IF ( ln_osm_mle ) THEN
2880	zhbl_s = zhbl_s + MIN( &
2881	& rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
2882	& e3w_n(ji,jj,jm) )
2883	ELSE
2884	zhbl_s = zhbl_s + MIN( &
2885	& rn_rdt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
2886	& e3w_n(ji,jj,jm) )
2887	ENDIF
2888
2889	! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2890	IF ( zhbl_s >= gdepw_n(ji,jj,mbkt(ji,jj) + 1) ) THEN
2891	zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2892	lpyc(ji,jj) = .FALSE.
2893	ENDIF
2894	IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
2895	#ifdef key_osm_debug
2896	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2897	WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm
2898	WRITE(narea+100,'(2(a,g11.3),a,l7)')'zdb=',zdb,' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj)
2899	FLUSH(narea+100)
2900	END IF
2901	#endif
2902	END DO
2903	hbl(ji,jj) = zhbl_s
2904	ibld(ji,jj) = jm
2905	ELSE
2906	! stable
2907	#ifdef key_osm_debug
2908	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2909	WRITE(narea+100,'(a)')'In zdf_osm_timestep_hbl, ibld - imld > 1, lconv=F'
2910	FLUSH(narea+100)
2911	END IF
2912	#endif
2913	DO jk = imld(ji,jj), ibld(ji,jj)
2914	zdb = MAX( &
2915	& grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&
2916	& - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
2917	& 0.0 ) + &
2918	& 2.0 * zvstr(ji,jj)**2 / zhbl_s
2919
2920	! Alan is thuis right? I have simply changed hbli to hbl
2921	zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
2922	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) ) ) * &
2923	& zustar(ji,jj)*3 / zhbl_s ) ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
2924	zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
2925	zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
2926
2927	! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2928	IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN
2929	zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
2930	lpyc(ji,jj) = .FALSE.
2931	ENDIF
2932	IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
2933	#ifdef key_osm_debug
2934	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2935	WRITE(narea+100,'(2(a,i7))')' jk=',jk,' jm=', jm
2936	WRITE(narea+100,'(4(a,g11.3),a,l7)')'zdb=',zdb,' zhol',zhol(ji,jj),' zdhdt',zdhdt(ji,jj),' zhbl_s=', zhbl_s,' lpyc=',lpyc(ji,jj)
2937	FLUSH(narea+100)
2938	END IF
2939	#endif
2940	END DO
2941	ENDIF ! IF ( lconv )
2942	hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
2943	ibld(ji,jj) = MAX(jm, 4 )
2944	ELSE
2945	! change zero or one model level.
2946	hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
2947	ENDIF
2948	zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
2949	#ifdef key_osm_debug
2950	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
2951	WRITE(narea+100,'(2(a,g11.3),a,i7,/)')'end of zdf_osm_timestep_hbl: hbl=', hbl(ji,jj),' zhbl=', zhbl(ji,jj),' ibld=', ibld(ji,jj)
2952	FLUSH(narea+100)
2953	END IF
2954	#endif
2955	END DO
2956	END DO
2957
2958	END SUBROUTINE zdf_osm_timestep_hbl
2959
2960	SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
2961	!!---------------------------------------------------------------------
2962	!! * ROUTINE zdf_osm_pycnocline_thickness *
2963	!!
2964	!! ** Purpose : Calculates thickness of the pycnocline
2965	!!
2966	!! ** Method : The thickness is calculated from a prognostic equation
2967	!! that relaxes the pycnocine thickness to a diagnostic
2968	!! value. The time change is calculated assuming the
2969	!! thickness relaxes exponentially. This is done to deal
2970	!! with large timesteps.
2971	!!
2972	!!----------------------------------------------------------------------
2973
2974	REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness.
2975	!
2976	INTEGER :: jj, ji
2977	INTEGER :: inhml
2978	REAL(wp) :: zari, ztau, zdh_ref, zddhdt, zvel_max
2979	REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
2980
2981	DO jj = 2, jpjm1
2982	DO ji = 2, jpim1
2983
2984	IF ( lshear(ji,jj) ) THEN
2985	IF ( lconv(ji,jj) ) THEN
2986	IF ( zdb_bl(ji,jj) > 1.0e-15) THEN
2987	IF ( j_ddh(ji,jj) == 0 ) THEN
2988	zvel_max = ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
2989	! ddhdt for pycnocline determined in osm_calculate_dhdt
2990	zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
2991	zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
2992	! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer
2993	dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) )
2994	ELSE
2995	! Need to recalculate because hbl has been updated.
2996	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
2997	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 )
2998	ELSE
2999	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 )
3000	ENDIF
3001	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 )
3002	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
3003	IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)
3004	ENDIF
3005	ELSE
3006	ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln), 2.0 * rn_rdt )
3007	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.2 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
3008	IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2 * hbl(ji,jj)
3009	ENDIF
3010	ELSE ! lconv
3011	! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL
3012
3013	ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
3014	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
3015	! boundary layer deepening
3016	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
3017	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
3018	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
3019	& / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )
3020	zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
3021	ELSE
3022	zdh_ref = 0.2 * hbl(ji,jj)
3023	ENDIF
3024	ELSE ! IF(dhdt < 0)
3025	zdh_ref = 0.2 * hbl(ji,jj)
3026	ENDIF ! IF (dhdt >= 0)
3027	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
3028	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
3029	ENDIF
3030
3031	ELSE ! lshear
3032	! for lshear = .FALSE. calculate ddhdt here
3033
3034	IF ( lconv(ji,jj) ) THEN
3035
3036	IF( ln_osm_mle ) THEN
3037	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN
3038	! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
3039	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
3040	IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability
3041	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 )
3042	ELSE ! unstable
3043	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 )
3044	ENDIF
3045	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
3046	zdh_ref = zari * hbl(ji,jj)
3047	ELSE
3048	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
3049	zdh_ref = 0.2 * hbl(ji,jj)
3050	ENDIF
3051	ELSE
3052	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
3053	zdh_ref = 0.2 * hbl(ji,jj)
3054	ENDIF
3055	ELSE ! ln_osm_mle
3056	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
3057	IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability
3058	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 )
3059	ELSE ! unstable
3060	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 )
3061	ENDIF
3062	ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
3063	zdh_ref = zari * hbl(ji,jj)
3064	ELSE
3065	ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
3066	zdh_ref = 0.2 * hbl(ji,jj)
3067	ENDIF
3068
3069	END IF ! ln_osm_mle
3070
3071	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
3072	! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
3073	IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
3074	! Alan: this hml is never defined or used
3075	ELSE ! IF (lconv)
3076	ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
3077	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
3078	! boundary layer deepening
3079	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
3080	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
3081	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
3082	& / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )
3083	zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
3084	ELSE
3085	zdh_ref = 0.2 * hbl(ji,jj)
3086	ENDIF
3087	ELSE ! IF(dhdt < 0)
3088	zdh_ref = 0.2 * hbl(ji,jj)
3089	ENDIF ! IF (dhdt >= 0)
3090	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
3091	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
3092	ENDIF ! IF (lconv)
3093	ENDIF ! lshear
3094
3095	hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
3096	inhml = MAX( INT( dh(ji,jj) / MAX(e3t_n(ji,jj,ibld(ji,jj)-1), 1.e-3) ) , 1 )
3097	imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
3098	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
3099	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
3100	#ifdef key_osm_debug
3101	IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
3102	WRITE(narea+100,'(4(a,g11.3),2(a,i7),/,5(a,g11.3),/)') 'end of zdf_osm_pycnocline_thickness:hml=',hml(ji,jj), &
3103	& ' zhml=',zhml(ji,jj),' zdh=', zdh(ji,jj), ' dh=', dh(ji,jj), ' imld=', imld(ji,jj), ' inhml=', inhml, &
3104	& 'zvel_max=', zvel_max, ' ztau=', ztau,' zdh_ref=', zdh_ref, ' zar=', zari, ' zddhdt=', zddhdt
3105	FLUSH(narea+100)
3106	END IF
3107	#endif
3108
3109	END DO
3110	END DO
3111
3112	END SUBROUTINE zdf_osm_pycnocline_thickness
3113
3114
3115	SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
3116	!!----------------------------------------------------------------------
3117	!! * ROUTINE zdf_osm_horizontal_gradients *
3118	!!
3119	!! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.
3120	!!
3121	!! ** Method :
3122	!!
3123	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
3124	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
3125
3126
3127	REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points
3128	REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient.
3129	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == !
3130	REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy
3131
3132	INTEGER :: ji, jj, jk ! dummy loop indices
3133	INTEGER :: ii, ij, ik, ikmax ! local integers
3134	REAL(wp) :: zc
3135	REAL(wp) :: zN2_c ! local buoyancy difference from 10m value
3136	REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH
3137	REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv
3138	REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv
3139	!!----------------------------------------------------------------------
3140	!
3141	! !== MLD used for MLE ==!
3142
3143	mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point
3144	zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2
3145	zN2_c = grav * rn_osm_mle_rho_c * r1_rau0 ! convert density criteria into N^2 criteria
3146	DO jk = nlb10, jpkm1
3147	DO jj = 1, jpj ! Mixed layer level: w-level
3148	DO ji = 1, jpi
3149	ikt = mbkt(ji,jj)
3150	zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
3151	IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
3152	END DO
3153	END DO
3154	END DO
3155	DO jj = 1, jpj
3156	DO ji = 1, jpi
3157	mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
3158	zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
3159	END DO
3160	END DO
3161	! ensure mld_prof .ge. ibld
3162	!
3163	ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation
3164	!
3165	ztm(:,:) = 0._wp
3166	zsm(:,:) = 0._wp
3167	DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer
3168	DO jj = 1, jpj
3169	DO ji = 1, jpi
3170	zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points
3171	ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem)
3172	zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal)
3173	END DO
3174	END DO
3175	END DO
3176	! average temperature and salinity.
3177	ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
3178	zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
3179	! calculate horizontal gradients at u & v points
3180
3181	DO jj = 2, jpjm1
3182	DO ji = 1, jpim1
3183	zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
3184	zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
3185	zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))
3186	ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )
3187	ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )
3188	END DO
3189	END DO
3190
3191	DO jj = 1, jpjm1
3192	DO ji = 2, jpim1
3193	zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
3194	zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
3195	zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))
3196	ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )
3197	ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )
3198	END DO
3199	END DO
3200
3201	! ensure salinity > 0 in unset values so EOS doesn't give FP error with fpe0 on
3202	ztsm_midu(:,jpj,jp_sal) = 10.
3203	ztsm_midv(jpi,:,jp_sal) = 10.
3204
3205	CALL eos_rab(ztsm_midu, zmld_midu, zabu)
3206	CALL eos_rab(ztsm_midv, zmld_midv, zabv)
3207
3208	DO jj = 2, jpjm1
3209	DO ji = 1, jpim1
3210	dbdx_mle(ji,jj) = grav(zdtdx(ji,jj)zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
3211	END DO
3212	END DO
3213	DO jj = 1, jpjm1
3214	DO ji = 2, jpim1
3215	dbdy_mle(ji,jj) = grav(zdtdy(ji,jj)zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
3216	END DO
3217	END DO
3218
3219	DO jj = 2, jpjm1
3220	DO ji = 2, jpim1
3221	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
3222	zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &
3223	& + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
3224	END DO
3225	END DO
3226
3227	END SUBROUTINE zdf_osm_zmld_horizontal_gradients
3228	SUBROUTINE zdf_osm_mle_parameters( zmld, mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
3229	!!----------------------------------------------------------------------
3230	!! * ROUTINE zdf_osm_mle_parameters *
3231	!!
3232	!! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.
3233	!!
3234	!! ** Method :
3235	!!
3236	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
3237	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
3238
3239	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == !
3240	INTEGER, DIMENSION(jpi,jpj) :: mld_prof
3241	REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle
3242	INTEGER :: ji, jj, jk ! dummy loop indices
3243	INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers
3244	INTEGER , DIMENSION(jpi,jpj) :: inml_mle
3245	REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle
3246
3247	! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
3248
3249	DO jj = 2, jpjm1
3250	DO ji = 2, jpim1
3251	IF ( lconv(ji,jj) ) THEN
3252	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
3253	! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
3254	zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1)
3255	zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2
3256	ENDIF
3257	END DO
3258	END DO
3259	! Timestep mixed layer eddy depth.
3260	DO jj = 2, jpjm1
3261	DO ji = 2, jpim1
3262	IF ( lmle(ji,jj) ) THEN ! MLE layer growing.
3263	! Buoyancy gradient at base of MLE layer.
3264	zthermal = rab_n(ji,jj,1,jp_tem)
3265	zbeta = rab_n(ji,jj,1,jp_sal)
3266	jkb = mld_prof(ji,jj)
3267	jkb1 = MIN(jkb + 1, mbkt(ji,jj))
3268	!
3269	zbuoy = grav * ( zthermal * tsn(ji,jj,mld_prof(ji,jj)+2,jp_tem) - zbeta * tsn(ji,jj,mld_prof(ji,jj)+2,jp_sal) )
3270	zdb_mle = zb_bl(ji,jj) - zbuoy
3271	! Timestep hmle.
3272	hmle(ji,jj) = hmle(ji,jj) + zwb0tot(ji,jj) * rn_rdt / zdb_mle
3273	ELSE
3274	IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN
3275	hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau
3276	ELSE
3277	hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt /rn_osm_mle_tau
3278	ENDIF
3279	ENDIF
3280	hmle(ji,jj) = MAX(MIN(hmle(ji,jj), ht_n(ji,jj)), gdepw_n(ji,jj,4))
3281	IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), rn_osm_hmle_limit*hbl(ji,jj) )
3282	! For now try just set hmle to zmld
3283	hmle(ji,jj) = zmld(ji,jj)
3284	END DO
3285	END DO
3286
3287	mld_prof = 4
3288	DO jk = 5, jpkm1
3289	DO jj = 2, jpjm1
3290	DO ji = 2, jpim1
3291	IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
3292	END DO
3293	END DO
3294	END DO
3295	DO jj = 2, jpjm1
3296	DO ji = 2, jpim1
3297	zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj))
3298	END DO
3299	END DO
3300	END SUBROUTINE zdf_osm_mle_parameters
3301
3302	END SUBROUTINE zdf_osm
3303
3304
3305	SUBROUTINE zdf_osm_init
3306	!!----------------------------------------------------------------------
3307	!! * ROUTINE zdf_osm_init *
3308	!!
3309	!! ** Purpose : Initialization of the vertical eddy diffivity and
3310	!! viscosity when using a osm turbulent closure scheme
3311	!!
3312	!! ** Method : Read the namosm namelist and check the parameters
3313	!! called at the first timestep (nit000)
3314	!!
3315	!! ** input : Namlist namosm
3316	!!----------------------------------------------------------------------
3317	INTEGER :: ios ! local integer
3318	INTEGER :: ji, jj, jk ! dummy loop indices
3319	REAL z1_t2
3320	REAL(wp) :: zlarge = -1.e10_wp, zero = 0._wp
3321	!!
3322	#ifdef key_osm_debug
3323	NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
3324	& ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd &
3325	& ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
3326	& ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter &
3327	& ,nn_idb, nn_jdb, nn_kdb, nn_narea_db
3328	#else
3329	NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
3330	& ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd &
3331	& ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
3332	& ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter
3333	#endif
3334	! Namelist for Fox-Kemper parametrization.
3335	NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
3336	& rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit
3337
3338	!!----------------------------------------------------------------------
3339	!
3340	REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model
3341	READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
3342	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
3343
3344	REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
3345	READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
3346	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
3347	IF(lwm) WRITE ( numond, namzdf_osm )
3348
3349	IF(lwp) THEN ! Control print
3350	WRITE(numout,*)
3351	WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
3352	WRITE(numout,*) '~~~~~~~~~~~~'
3353	WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters'
3354	WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la
3355	WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle
3356	WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la
3357	WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd
3358	WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0
3359	WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes
3360	WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave
3361	WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave
3362	SELECT CASE (nn_osm_wave)
3363	CASE(0)
3364	WRITE(numout,*) ' calculated assuming constant La#=0.3'
3365	CASE(1)
3366	WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves'
3367	CASE(2)
3368	WRITE(numout,*) ' calculated from ECMWF wave fields'
3369	END SELECT
3370	WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce', nn_osm_SD_reduce
3371	WRITE(numout,*) ' fraction of hbl to average SD over/fit'
3372	WRITE(numout,*) ' exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac
3373	SELECT CASE (nn_osm_SD_reduce)
3374	CASE(0)
3375	WRITE(numout,*) ' No reduction'
3376	CASE(1)
3377	WRITE(numout,*) ' Average SD over upper rn_osm_hblfrac of BL'
3378	CASE(2)
3379	WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL'
3380	END SELECT
3381	WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter
3382	WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm
3383	WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s'
3384	WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix
3385	WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
3386	WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri
3387	WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix
3388	WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv
3389	#ifdef key_osm_debug
3390	WRITE(numout,*) 'nn_idb', nn_idb, 'nn_jdb', nn_jdb, 'nn_kdb', nn_kdb, 'nn_narea_db', nn_narea_db
3391
3392	iloc_db = mi0(nn_idb)
3393	jloc_db = mj0(nn_jdb)
3394	WRITE(numout,*) 'iloc_db ', iloc_db , 'jloc_db', jloc_db
3395	#endif
3396	ENDIF
3397
3398
3399	! ! Check wave coupling settings !
3400	! ! Further work needed - see ticket #2447 !
3401	IF( nn_osm_wave == 2 ) THEN
3402	IF (.NOT. ( ln_wave .AND. ln_sdw )) &
3403	& CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' )
3404	END IF
3405
3406	! ! allocate zdfosm arrays
3407	IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
3408
3409
3410	IF( ln_osm_mle ) THEN
3411	! Initialise Fox-Kemper parametrization
3412	REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme
3413	READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
3414	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist')
3415
3416	REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme
3417	READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
3418	904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist')
3419	IF(lwm) WRITE ( numond, namosm_mle )
3420
3421	IF(lwp) THEN ! Namelist print
3422	WRITE(numout,*)
3423	WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
3424	WRITE(numout,*) '~~~~~~~~~~~~~'
3425	WRITE(numout,*) ' Namelist namosm_mle : '
3426	WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle
3427	WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce
3428	WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm'
3429	WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's'
3430	WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg'
3431	WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c
3432	WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s'
3433	WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's'
3434	WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit
3435	WRITE(numout,*) ' fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T. rn_osm_hmle_limit = ', rn_osm_hmle_limit
3436	ENDIF !
3437	ENDIF
3438	!
3439	IF(lwp) THEN
3440	WRITE(numout,*)
3441	IF( ln_osm_mle ) THEN
3442	WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
3443	IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation'
3444	IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation'
3445	ELSE
3446	WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
3447	ENDIF
3448	ENDIF
3449	!
3450	IF( ln_osm_mle ) THEN ! MLE initialisation
3451	!
3452	rb_c = grav * rn_osm_mle_rho_c /rau0 ! Mixed Layer buoyancy criteria
3453	IF(lwp) WRITE(numout,*)
3454	IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 '
3455	IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
3456	!
3457	IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation !
3458	!
3459	ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation
3460	rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) )
3461	!
3462	ENDIF
3463	! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
3464	z1_t2 = 2.e-5
3465	do jj=1,jpj
3466	do ji = 1,jpi
3467	r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
3468	end do
3469	end do
3470	! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
3471	! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 )
3472	!
3473	ENDIF
3474
3475	call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle
3476
3477
3478	IF( ln_zdfddm) THEN
3479	IF(lwp) THEN
3480	WRITE(numout,*)
3481	WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
3482	WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm '
3483	ENDIF
3484	ENDIF
3485
3486
3487	!set constants not in namelist
3488	!-----------------------------
3489
3490	IF(lwp) THEN
3491	WRITE(numout,*)
3492	ENDIF
3493
3494	dstokes(:,:) = zlarge
3495	IF (nn_osm_wave == 0) THEN
3496	dstokes(:,:) = rn_osm_dstokes
3497	END IF
3498
3499	! Horizontal average : initialization of weighting arrays
3500	! -------------------
3501
3502	SELECT CASE ( nn_ave )
3503
3504	CASE ( 0 ) ! no horizontal average
3505	IF(lwp) WRITE(numout,*) ' no horizontal average on avt'
3506	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
3507	! weighting mean arrays etmean
3508	! ( 1 1 )
3509	! avt = 1/4 ( 1 1 )
3510	!
3511	etmean(:,:,:) = 0.e0
3512
3513	DO jk = 1, jpkm1
3514	DO jj = 2, jpjm1
3515	DO ji = 2, jpim1 ! vector opt.
3516	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
3517	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
3518	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
3519	END DO
3520	END DO
3521	END DO
3522
3523	CASE ( 1 ) ! horizontal average
3524	IF(lwp) WRITE(numout,*) ' horizontal average on avt'
3525	! weighting mean arrays etmean
3526	! ( 1/2 1 1/2 )
3527	! avt = 1/8 ( 1 2 1 )
3528	! ( 1/2 1 1/2 )
3529	etmean(:,:,:) = 0.e0
3530
3531	DO jk = 1, jpkm1
3532	DO jj = 2, jpjm1
3533	DO ji = 2, jpim1 ! vector opt.
3534	etmean(ji,jj,jk) = tmask(ji, jj,jk) &
3535	& / MAX( 1., 2.* tmask(ji,jj,jk) &
3536	& +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) &
3537	& +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
3538	& +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) &
3539	& +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) )
3540	END DO
3541	END DO
3542	END DO
3543
3544	CASE DEFAULT
3545	WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave
3546	CALL ctl_stop( ctmp1 )
3547
3548	END SELECT
3549
3550	! Initialization of vertical eddy coef. to the background value
3551	! -------------------------------------------------------------
3552	DO jk = 1, jpk
3553	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
3554	END DO
3555
3556	! zero the surface flux for non local term and osm mixed layer depth
3557	! ------------------------------------------------------------------
3558	ghamt(:,:,:) = 0.
3559	ghams(:,:,:) = 0.
3560	ghamu(:,:,:) = 0.
3561	ghamv(:,:,:) = 0.
3562	!
3563	IF( lwxios ) THEN
3564	CALL iom_set_rstw_var_active('wn')
3565	CALL iom_set_rstw_var_active('hbl')
3566	CALL iom_set_rstw_var_active('dh')
3567	IF( ln_osm_mle ) THEN
3568	CALL iom_set_rstw_var_active('hmle')
3569	END IF
3570	ENDIF
3571	END SUBROUTINE zdf_osm_init
3572
3573
3574	SUBROUTINE osm_rst( kt, cdrw )
3575	!!---------------------------------------------------------------------
3576	!! * ROUTINE osm_rst *
3577	!!
3578	!! ** Purpose : Read or write BL fields in restart file
3579	!!
3580	!! ** Method : use of IOM library. If the restart does not contain
3581	!! required fields, they are recomputed from stratification
3582	!!----------------------------------------------------------------------
3583
3584	INTEGER, INTENT(in) :: kt
3585	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
3586
3587	INTEGER :: id1, id2, id3 ! iom enquiry index
3588	INTEGER :: ji, jj, jk ! dummy loop indices
3589	INTEGER :: iiki, ikt ! local integer
3590	REAL(wp) :: zhbf ! tempory scalars
3591	REAL(wp) :: zN2_c ! local scalar
3592	REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth
3593	INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
3594	!!----------------------------------------------------------------------
3595	!
3596	!!-----------------------------------------------------------------------------
3597	! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
3598	!!-----------------------------------------------------------------------------
3599	IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
3600	id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. )
3601	IF( id1 > 0 ) THEN ! 'wn' exists; read
3602	CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios )
3603	WRITE(numout,*) ' ===>>>> : wn read from restart file'
3604	ELSE
3605	wn(:,:,:) = 0._wp
3606	WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
3607	END IF
3608
3609	id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. )
3610	id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. )
3611	IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return
3612	CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
3613	CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios )
3614	hml (:,:) = hbl(:,:) - dh(:,:) ! Initialise ML depth.
3615	WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file'
3616	IF( ln_osm_mle ) THEN
3617	id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. )
3618	IF( id3 > 0) THEN
3619	CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios )
3620	WRITE(numout,*) ' ===>>>> : hmle read from restart file'
3621	ELSE
3622	WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl'
3623	hmle(:,:) = hbl(:,:) ! Initialise MLE depth.
3624	END IF
3625	END IF
3626	RETURN
3627	ELSE ! 'hbl' & 'dh' not in restart file, recalculate
3628	WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
3629	END IF
3630	END IF
3631
3632	!!-----------------------------------------------------------------------------
3633	! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
3634	!!-----------------------------------------------------------------------------
3635	IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return
3636	IF(lwp) WRITE(numout,*) '---- osm-rst ----'
3637	CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn, ldxios = lwxios )
3638	CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios )
3639	CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios )
3640	IF( ln_osm_mle ) THEN
3641	CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios )
3642	END IF
3643	RETURN
3644	END IF
3645
3646	!!-----------------------------------------------------------------------------
3647	! Getting hbl, no restart file with hbl, so calculate from surface stratification
3648	!!-----------------------------------------------------------------------------
3649	IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
3650	! w-level of the mixing and mixed layers
3651	CALL eos_rab( tsn, rab_n )
3652	CALL bn2(tsn, rab_n, rn2)
3653	imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point
3654	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
3655	zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
3656	!
3657	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
3658	DO jk = 1, jpkm1
3659	DO jj = 1, jpj ! Mixed layer level: w-level
3660	DO ji = 1, jpi
3661	ikt = mbkt(ji,jj)
3662	hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
3663	IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
3664	END DO
3665	END DO
3666	END DO
3667	!
3668	DO jj = 1, jpj
3669	DO ji = 1, jpi
3670	iiki = MAX(4,imld_rst(ji,jj))
3671	hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth
3672	dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth
3673	hml (ji,jj) = hbl(ji,jj) - dh(ji,jj)
3674	END DO
3675	END DO
3676
3677	WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
3678
3679	IF( ln_osm_mle ) THEN
3680	hmle(:,:) = hbl(:,:) ! Initialise MLE depth.
3681	WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
3682	END IF
3683
3684	wn(:,:,:) = 0._wp
3685	WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
3686	END SUBROUTINE osm_rst
3687
3688
3689	SUBROUTINE tra_osm( kt )
3690	!!----------------------------------------------------------------------
3691	!! * ROUTINE tra_osm *
3692	!!
3693	!! ** Purpose : compute and add to the tracer trend the non-local tracer flux
3694	!!
3695	!! ** Method : ???
3696	!!----------------------------------------------------------------------
3697	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace
3698	!!----------------------------------------------------------------------
3699	INTEGER, INTENT(in) :: kt
3700	INTEGER :: ji, jj, jk
3701	!
3702	IF( kt == nit000 ) THEN
3703	IF(lwp) WRITE(numout,*)
3704	IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
3705	IF(lwp) WRITE(numout,*) '~~~~~~~ '
3706	ENDIF
3707
3708	IF( l_trdtra ) THEN !* Save ta and sa trends
3709	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
3710	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)
3711	ENDIF
3712
3713	DO jk = 1, jpkm1
3714	DO jj = 2, jpjm1
3715	DO ji = 2, jpim1
3716	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &
3717	& - ( ghamt(ji,jj,jk ) &
3718	& - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
3719	tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) &
3720	& - ( ghams(ji,jj,jk ) &
3721	& - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
3722	END DO
3723	END DO
3724	END DO
3725
3726	! save the non-local tracer flux trends for diagnostics
3727	IF( l_trdtra ) THEN
3728	ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
3729	ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
3730
3731	CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt )
3732	CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds )
3733	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
3734	ENDIF
3735
3736	IF(ln_ctl) THEN
3737	CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm - Ta: ', mask1=tmask, &
3738	& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
3739	ENDIF
3740	!
3741	END SUBROUTINE tra_osm
3742
3743
3744	SUBROUTINE trc_osm( kt ) ! Dummy routine
3745	!!----------------------------------------------------------------------
3746	!! * ROUTINE trc_osm *
3747	!!
3748	!! ** Purpose : compute and add to the passive tracer trend the non-local
3749	!! passive tracer flux
3750	!!
3751	!!
3752	!! ** Method : ???
3753	!!----------------------------------------------------------------------
3754	!
3755	!!----------------------------------------------------------------------
3756	INTEGER, INTENT(in) :: kt
3757	WRITE(,) 'trc_osm: Not written yet', kt
3758	END SUBROUTINE trc_osm
3759
3760
3761	SUBROUTINE dyn_osm( kt )
3762	!!----------------------------------------------------------------------
3763	!! * ROUTINE dyn_osm *
3764	!!
3765	!! ** Purpose : compute and add to the velocity trend the non-local flux
3766	!! copied/modified from tra_osm
3767	!!
3768	!! ** Method : ???
3769	!!----------------------------------------------------------------------
3770	INTEGER, INTENT(in) :: kt !
3771	!
3772	INTEGER :: ji, jj, jk ! dummy loop indices
3773	!!----------------------------------------------------------------------
3774	!
3775	IF( kt == nit000 ) THEN
3776	IF(lwp) WRITE(numout,*)
3777	IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
3778	IF(lwp) WRITE(numout,*) '~~~~~~~ '
3779	ENDIF
3780	!code saving tracer trends removed, replace with trdmxl_oce
3781
3782	DO jk = 1, jpkm1 ! add non-local u and v fluxes
3783	DO jj = 2, jpjm1
3784	DO ji = 2, jpim1
3785	ua(ji,jj,jk) = ua(ji,jj,jk) &
3786	& - ( ghamu(ji,jj,jk ) &
3787	& - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
3788	va(ji,jj,jk) = va(ji,jj,jk) &
3789	& - ( ghamv(ji,jj,jk ) &
3790	& - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
3791	END DO
3792	END DO
3793	END DO
3794	!
3795	! code for saving tracer trends removed
3796	!
3797	END SUBROUTINE dyn_osm
3798
3799	!!======================================================================
3800
3801	END MODULE zdfosm

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_FKOSM_m11715/src/OCE/ZDF/zdfosm.F90 @ 14774

Download in other formats: