Context Navigation

zdfosm.F90 @ 15540

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

Download in other formats: