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

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

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

Context Navigation

source: NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90 @ 12928

Last change on this file since 12928 was 12928, checked in by smueller, 4 years ago
Synchronizing with ^{/NEMO/trunk@12925 (ticket #2170)}
Property svn:keywords set to `Id`
File size: 122.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	!!----------------------------------------------------------------------
37
38	!!----------------------------------------------------------------------
39	!! 'ln_zdfosm' OSMOSIS scheme
40	!!----------------------------------------------------------------------
41	!! zdf_osm : update momentum and tracer Kz from osm scheme
42	!! zdf_osm_init : initialization, namelist read, and parameters control
43	!! osm_rst : read (or initialize) and write osmosis restart fields
44	!! tra_osm : compute and add to the T & S trend the non-local flux
45	!! trc_osm : compute and add to the passive tracer trend the non-local flux (TBD)
46	!! dyn_osm : compute and add to u & v trensd the non-local flux
47	!!
48	!! Subroutines in revised code.
49	!!----------------------------------------------------------------------
50	USE oce ! ocean dynamics and active tracers
51	! uses ww from previous time step (which is now wb) to calculate hbl
52	USE dom_oce ! ocean space and time domain
53	USE zdf_oce ! ocean vertical physics
54	USE sbc_oce ! surface boundary condition: ocean
55	USE sbcwave ! surface wave parameters
56	USE phycst ! physical constants
57	USE eosbn2 ! equation of state
58	USE traqsr ! details of solar radiation absorption
59	USE zdfddm ! double diffusion mixing (avs array)
60	USE iom ! I/O library
61	USE lib_mpp ! MPP library
62	USE trd_oce ! ocean trends definition
63	USE trdtra ! tracers trends
64	!
65	USE in_out_manager ! I/O manager
66	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
67	USE prtctl ! Print control
68	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
69
70	IMPLICIT NONE
71	PRIVATE
72
73	PUBLIC zdf_osm ! routine called by step.F90
74	PUBLIC zdf_osm_init ! routine called by nemogcm.F90
75	PUBLIC osm_rst ! routine called by step.F90
76	PUBLIC tra_osm ! routine called by step.F90
77	PUBLIC trc_osm ! routine called by trcstp.F90
78	PUBLIC dyn_osm ! routine called by step.F90
79
80	PUBLIC ln_osm_mle ! logical needed by tra_mle_init in tramle.F90
81
82	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamu !: non-local u-momentum flux
83	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamv !: non-local v-momentum flux
84	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghamt !: non-local temperature flux (gamma/<ws>o)
85	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ghams !: non-local salinity flux (gamma/<ws>o)
86	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: etmean !: averaging operator for avt
87	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbl !: boundary layer depth
88	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dh ! depth of pycnocline
89	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hml ! ML depth
90	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dstokes !: penetration depth of the Stokes drift.
91
92	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_ft ! inverse of the modified Coriolis parameter at t-pts
93	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmle ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization
94	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdx_mle ! zonal buoyancy gradient in ML
95	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: dbdy_mle ! meridional buoyancy gradient in ML
96	INTEGER, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mld_prof ! level of base of MLE layer.
97
98	! !! Namelist namzdf_osm
99	LOGICAL :: ln_use_osm_la ! Use namelist rn_osm_la
100
101	LOGICAL :: ln_osm_mle !: flag to activate the Mixed Layer Eddy (MLE) parameterisation
102
103	REAL(wp) :: rn_osm_la ! Turbulent Langmuir number
104	REAL(wp) :: rn_osm_dstokes ! Depth scale of Stokes drift
105	REAL(wp) :: rn_osm_hbl0 = 10._wp ! Initial value of hbl for 1D runs
106	INTEGER :: nn_ave ! = 0/1 flag for horizontal average on avt
107	INTEGER :: nn_osm_wave = 0 ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
108	LOGICAL :: ln_dia_osm ! Use namelist rn_osm_la
109
110
111	LOGICAL :: ln_kpprimix = .true. ! Shear instability mixing
112	REAL(wp) :: rn_riinfty = 0.7 ! local Richardson Number limit for shear instability
113	REAL(wp) :: rn_difri = 0.005 ! maximum shear mixing at Rig = 0 (m2/s)
114	LOGICAL :: ln_convmix = .true. ! Convective instability mixing
115	REAL(wp) :: rn_difconv = 1._wp ! diffusivity when unstable below BL (m2/s)
116
117	! OSMOSIS mixed layer eddy parametrization constants
118	INTEGER :: nn_osm_mle ! = 0/1 flag for horizontal average on avt
119	REAL(wp) :: rn_osm_mle_ce ! MLE coefficient
120	! ! parameters used in nn_osm_mle = 0 case
121	REAL(wp) :: rn_osm_mle_lf ! typical scale of mixed layer front
122	REAL(wp) :: rn_osm_mle_time ! time scale for mixing momentum across the mixed layer
123	! ! parameters used in nn_osm_mle = 1 case
124	REAL(wp) :: rn_osm_mle_lat ! reference latitude for a 5 km scale of ML front
125	REAL(wp) :: rn_osm_mle_rho_c ! Density criterion for definition of MLD used by FK
126	REAL(wp) :: r5_21 = 5.e0 / 21.e0 ! factor used in mle streamfunction computation
127	REAL(wp) :: rb_c ! ML buoyancy criteria = g rho_c /rau0 where rho_c is defined in zdfmld
128	REAL(wp) :: rc_f ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case
129	REAL(wp) :: rn_osm_mle_thresh ! Threshold buoyancy for deepening of MLE layer below OSBL base.
130	REAL(wp) :: rn_osm_mle_tau ! Adjustment timescale for MLE.
131
132
133	! !!! General constants
134	REAL(wp) :: epsln = 1.0e-20_wp ! a small positive number to ensure no div by zero
135	REAL(wp) :: depth_tol = 1.0e-6_wp ! a small-ish positive number to give a hbl slightly shallower than gdepw
136	REAL(wp) :: pthird = 1._wp/3._wp ! 1/3
137	REAL(wp) :: p2third = 2._wp/3._wp ! 2/3
138
139	INTEGER :: idebug = 236
140	INTEGER :: jdebug = 228
141	!! * Substitutions
142	# include "do_loop_substitute.h90"
143	!!----------------------------------------------------------------------
144	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
145	!! $Id$
146	!! Software governed by the CeCILL license (see ./LICENSE)
147	!!----------------------------------------------------------------------
148	CONTAINS
149
150	INTEGER FUNCTION zdf_osm_alloc()
151	!!----------------------------------------------------------------------
152	!! * FUNCTION zdf_osm_alloc *
153	!!----------------------------------------------------------------------
154	ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &
155	& hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
156	& etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
157
158	ALLOCATE( hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), &
159	& mld_prof(jpi,jpj), STAT= zdf_osm_alloc )
160
161	! ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), & ! would ths be better ?
162	! & hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
163	! & etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
164	! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
165	! IF ( ln_osm_mle ) THEN
166	! Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc )
167	! IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays')
168	! ENDIF
169
170	IF( zdf_osm_alloc /= 0 ) CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
171	CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
172	END FUNCTION zdf_osm_alloc
173
174
175	SUBROUTINE zdf_osm( kt, Kbb, Kmm, Krhs, p_avm, p_avt )
176	!!----------------------------------------------------------------------
177	!! * ROUTINE zdf_osm *
178	!!
179	!! ** Purpose : Compute the vertical eddy viscosity and diffusivity
180	!! coefficients and non local mixing using the OSMOSIS scheme
181	!!
182	!! ** Method : The boundary layer depth hosm is diagnosed at tracer points
183	!! from profiles of buoyancy, and shear, and the surface forcing.
184	!! Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
185	!!
186	!! Kx = hosm Wx(sigma) G(sigma)
187	!!
188	!! and the non local term ghamt = Cs / Ws(sigma) / hosm
189	!! Below hosm the coefficients are the sum of mixing due to internal waves
190	!! shear instability and double diffusion.
191	!!
192	!! -1- Compute the now interior vertical mixing coefficients at all depths.
193	!! -2- Diagnose the boundary layer depth.
194	!! -3- Compute the now boundary layer vertical mixing coefficients.
195	!! -4- Compute the now vertical eddy vicosity and diffusivity.
196	!! -5- Smoothing
197	!!
198	!! N.B. The computation is done from jk=2 to jpkm1
199	!! Surface value of avt are set once a time to zero
200	!! in routine zdf_osm_init.
201	!!
202	!! ** Action : update the non-local terms ghamts
203	!! update avt (before vertical eddy coef.)
204	!!
205	!! References : Large W.G., Mc Williams J.C. and Doney S.C.
206	!! Reviews of Geophysics, 32, 4, November 1994
207	!! Comments in the code refer to this paper, particularly
208	!! the equation number. (LMD94, here after)
209	!!----------------------------------------------------------------------
210	INTEGER , INTENT(in ) :: kt ! ocean time step
211	INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
212	REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points)
213	!!
214	INTEGER :: ji, jj, jk ! dummy loop indices
215
216	INTEGER :: jl ! dummy loop indices
217
218	INTEGER :: ikbot, jkmax, jkm1, jkp2 !
219
220	REAL(wp) :: ztx, zty, zflageos, zstabl, zbuofdep,zucube !
221	REAL(wp) :: zbeta, zthermal !
222	REAL(wp) :: zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
223	REAL(wp) :: zwsun, zwmun, zcons, zconm, zwcons, zwconm !
224
225	REAL(wp) :: zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed ! In situ density
226	INTEGER :: jm ! dummy loop indices
227	REAL(wp) :: zr1, zr2, zr3, zr4, zrhop ! Compression terms
228	REAL(wp) :: zflag, zrn2, zdep21, zdep32, zdep43
229	REAL(wp) :: zesh2, zri, zfri ! Interior richardson mixing
230	REAL(wp) :: zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
231	REAL(wp) :: zt,zs,zu,zv,zrh ! variables used in constructing averages
232	! Scales
233	REAL(wp), DIMENSION(jpi,jpj) :: zrad0 ! Surface solar temperature flux (deg m/s)
234	REAL(wp), DIMENSION(jpi,jpj) :: zradh ! Radiative flux at bl base (Buoyancy units)
235	REAL(wp), DIMENSION(jpi,jpj) :: zradav ! Radiative flux, bl average (Buoyancy Units)
236	REAL(wp), DIMENSION(jpi,jpj) :: zustar ! friction velocity
237	REAL(wp), DIMENSION(jpi,jpj) :: zwstrl ! Langmuir velocity scale
238	REAL(wp), DIMENSION(jpi,jpj) :: zvstr ! Velocity scale that ends to zustar for large Langmuir number.
239	REAL(wp), DIMENSION(jpi,jpj) :: zwstrc ! Convective velocity scale
240	REAL(wp), DIMENSION(jpi,jpj) :: zuw0 ! Surface u-momentum flux
241	REAL(wp), DIMENSION(jpi,jpj) :: zvw0 ! Surface v-momentum flux
242	REAL(wp), DIMENSION(jpi,jpj) :: zwth0 ! Surface heat flux (Kinematic)
243	REAL(wp), DIMENSION(jpi,jpj) :: zws0 ! Surface freshwater flux
244	REAL(wp), DIMENSION(jpi,jpj) :: zwb0 ! Surface buoyancy flux
245	REAL(wp), DIMENSION(jpi,jpj) :: zwthav ! Heat flux - bl average
246	REAL(wp), DIMENSION(jpi,jpj) :: zwsav ! freshwater flux - bl average
247	REAL(wp), DIMENSION(jpi,jpj) :: zwbav ! Buoyancy flux - bl average
248	REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent ! Buoyancy entrainment flux
249
250
251	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b ! MLE buoyancy flux averaged over OSBL
252	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk ! max MLE buoyancy flux
253	REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff
254	REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle ! velocity scale for dhdt with stable ML and FK
255
256	REAL(wp), DIMENSION(jpi,jpj) :: zustke ! Surface Stokes drift
257	REAL(wp), DIMENSION(jpi,jpj) :: zla ! Trubulent Langmuir number
258	REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
259	REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
260	REAL(wp), DIMENSION(jpi,jpj) :: zhol ! Stability parameter for boundary layer
261	LOGICAL, DIMENSION(jpi,jpj) :: lconv ! unstable/stable bl
262
263	! mixed-layer variables
264
265	INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
266	INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
267
268	REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
269	REAL(wp) :: zugrad,zvgrad ! temporary variables for calculating pycnocline shear
270
271	REAL(wp), DIMENSION(jpi,jpj) :: zhbl ! bl depth - grid
272	REAL(wp), DIMENSION(jpi,jpj) :: zhml ! ml depth - grid
273
274	REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid
275	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! ML depth on grid
276
277	REAL(wp), DIMENSION(jpi,jpj) :: zdh ! pycnocline depth - grid
278	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
279	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt_2 ! correction to dhdt due to internal structure.
280	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_ext ! external temperature/salinity and buoyancy gradients
281
282	REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy ! horizontal gradients for Fox-Kemper parametrization.
283
284	REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl ! averages over the depth of the blayer
285	REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml ! averages over the depth of the mixed layer
286	REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer
287	REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
288	REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
289	REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse ! momentum fluxes at the top of the pycnocline
290	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc ! parametrized gradient of temperature in pycnocline
291	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc ! parametrised gradient of salinity in pycnocline
292	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc ! parametrised gradient of buoyancy in the pycnocline
293	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc ! u-shear across the pycnocline
294	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc ! v-shear across the pycnocline
295
296	! Flux-gradient relationship variables
297
298	REAL(wp) :: zl_c,zl_l,zl_eps ! Used to calculate turbulence length scale.
299
300	REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity
301	REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
302	REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms.
303	REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
304
305	! For calculating Ri#-dependent mixing
306	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du ! u-shear^2
307	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv ! v-shear^2
308	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
309
310	! Temporary variables
311	INTEGER :: inhml
312	REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
313	REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb ! temporary variables
314	REAL(wp) :: zthick, zz0, zz1 ! temporary variables
315	REAL(wp) :: zvel_max, zhbl_s ! temporary variables
316	REAL(wp) :: zfac ! temporary variable
317	REAL(wp) :: zus_x, zus_y ! temporary Stokes drift
318	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
319	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
320
321	INTEGER :: ibld_ext=0 ! does not have to be zero for modified scheme
322	REAL(wp) :: zwb_min, zgamma_b_nd, zgamma_b, zdhoh, ztau, zddhdt
323	REAL(wp) :: zzeta_s = 0._wp
324	REAL(wp) :: zzeta_v = 0.46
325	REAL(wp) :: zabsstke
326
327	! For debugging
328	INTEGER :: ikt
329	!!--------------------------------------------------------------------
330	!
331	ibld(:,:) = 0 ; imld(:,:) = 0
332	zrad0(:,:) = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:) = 0._wp ; zustar(:,:) = 0._wp
333	zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:) = 0._wp ; zuw0(:,:) = 0._wp
334	zvw0(:,:) = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:) = 0._wp ; zwb0(:,:) = 0._wp
335	zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:) = 0._wp ; zwb_ent(:,:) = 0._wp
336	zustke(:,:) = 0._wp ; zla(:,:) = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
337	zhol(:,:) = 0._wp
338	lconv(:,:) = .FALSE.
339	! mixed layer
340	! no initialization of zhbl or zhml (or zdh?)
341	zhbl(:,:) = 1._wp ; zhml(:,:) = 1._wp ; zdh(:,:) = 1._wp ; zdhdt(:,:) = 0._wp
342	zt_bl(:,:) = 0._wp ; zs_bl(:,:) = 0._wp ; zu_bl(:,:) = 0._wp ; zv_bl(:,:) = 0._wp
343	zrh_bl(:,:) = 0._wp ; zt_ml(:,:) = 0._wp ; zs_ml(:,:) = 0._wp ; zu_ml(:,:) = 0._wp
344
345	zv_ml(:,:) = 0._wp ; zrh_ml(:,:) = 0._wp ; zdt_bl(:,:) = 0._wp ; zds_bl(:,:) = 0._wp
346	zdu_bl(:,:) = 0._wp ; zdv_bl(:,:) = 0._wp ; zdrh_bl(:,:) = 0._wp ; zdb_bl(:,:) = 0._wp
347	zdt_ml(:,:) = 0._wp ; zds_ml(:,:) = 0._wp ; zdu_ml(:,:) = 0._wp ; zdv_ml(:,:) = 0._wp
348	zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:) = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp
349	zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp
350	!
351	zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
352	zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
353	!
354	zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_ext(:,:) = 0._wp
355
356	IF ( ln_osm_mle ) THEN ! only initialise arrays if needed
357	zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp
358	zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp
359	zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp
360	zhmle(:,:) = 0._wp ; zmld(:,:) = 0._wp
361	ENDIF
362	zwb_fk_b(:,:) = 0._wp ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt
363
364	! Flux-Gradient arrays.
365	zdifml_sc(:,:) = 0._wp ; zvisml_sc(:,:) = 0._wp ; zdifpyc_sc(:,:) = 0._wp
366	zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp
367	zsc_wth_1(:,:) = 0._wp ; zsc_ws_1(:,:) = 0._wp ; zsc_uw_1(:,:) = 0._wp
368	zsc_uw_2(:,:) = 0._wp ; zsc_vw_1(:,:) = 0._wp ; zsc_vw_2(:,:) = 0._wp
369	zhbl_t(:,:) = 0._wp ; zdhdt(:,:) = 0._wp
370
371	zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
372	ghams(:,:,:) = 0._wp ; ghamu(:,:,:) = 0._wp ; ghamv(:,:,:) = 0._wp
373
374	zdhdt_2(:,:) = 0._wp
375	! hbl = MAX(hbl,epsln)
376	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
377	! Calculate boundary layer scales
378	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
379
380	! Assume two-band radiation model for depth of OSBL
381	zz0 = rn_abs ! surface equi-partition in 2-bands
382	zz1 = 1. - rn_abs
383	DO_2D_00_00
384	! Surface downward irradiance (so always +ve)
385	zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_rcp
386	! Downwards irradiance at base of boundary layer
387	zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
388	! Downwards irradiance averaged over depth of the OSBL
389	zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
390	& + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
391	END_2D
392	! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
393	DO_2D_00_00
394	zthermal = rab_n(ji,jj,1,jp_tem)
395	zbeta = rab_n(ji,jj,1,jp_sal)
396	! Upwards surface Temperature flux for non-local term
397	zwth0(ji,jj) = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1)
398	! Upwards surface salinity flux for non-local term
399	zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)
400	! Non radiative upwards surface buoyancy flux
401	zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj)
402	! turbulent heat flux averaged over depth of OSBL
403	zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
404	! turbulent salinity flux averaged over depth of the OBSL
405	zwsav(ji,jj) = 0.5 * zws0(ji,jj)
406	! turbulent buoyancy flux averaged over the depth of the OBSBL
407	zwbav(ji,jj) = grav * zthermal * zwthav(ji,jj) - grav * zbeta * zwsav(ji,jj)
408	! Surface upward velocity fluxes
409	zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
410	zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
411	! Friction velocity (zustar), at T-point : LMD94 eq. 2
412	zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
413	zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
414	zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
415	END_2D
416	! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
417	SELECT CASE (nn_osm_wave)
418	! Assume constant La#=0.3
419	CASE(0)
420	DO_2D_00_00
421	zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
422	zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
423	zustke(ji,jj) = MAX ( SQRT( zus_xzus_x + zus_yzus_y), 1.0e-8 )
424	! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
425	END_2D
426	! Assume Pierson-Moskovitz wind-wave spectrum
427	CASE(1)
428	DO_2D_00_00
429	! Use wind speed wndm included in sbc_oce module
430	zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
431	dstokes(ji,jj) = MAX( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
432	END_2D
433	! Use ECMWF wave fields as output from SBCWAVE
434	CASE(2)
435	zfac = 2.0_wp * rpi / 16.0_wp
436	DO_2D_00_00
437	! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
438	! The coefficient 0.8 gives La=0.3 in this situation.
439	! It could represent the effects of the spread of wave directions
440	! around the mean wind. The effect of this adjustment needs to be tested.
441	zabsstke = SQRT(ut0sd(ji,jj)2 + vt0sd(ji,jj)2)
442	zustke(ji,jj) = MAX (0.8 * ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj) * vt0sd(ji,jj) ), 1.0e-8)
443	dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)hsw(ji,jj) / ( MAX(zabsstkewmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
444	END_2D
445	END SELECT
446
447	! Langmuir velocity scale (zwstrl), La # (zla)
448	! mixed scale (zvstr), convective velocity scale (zwstrc)
449	DO_2D_00_00
450	! Langmuir velocity scale (zwstrl), at T-point
451	zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
452	zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)
453	IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))
454	! Velocity scale that tends to zustar for large Langmuir numbers
455	zvstr(ji,jj) = ( zwstrl(ji,jj)**3 + &
456	& ( 1.0 - EXP( -0.5 * zla(ji,jj)*2 ) ) zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
457
458	! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
459	! Note zustke and zwstrl are not amended.
460	!
461	! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
462	IF ( zwbav(ji,jj) > 0.0) THEN
463	zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
464	zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
465	lconv(ji,jj) = .TRUE.
466	ELSE
467	zhol(ji,jj) = -hbl(ji,jj) * 2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3 + epsln )
468	lconv(ji,jj) = .FALSE.
469	ENDIF
470	END_2D
471
472	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
473	! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
474	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
475	! BL must be always 4 levels deep.
476	! For calculation of lateral buoyancy gradients for FK in
477	! zdf_osm_zmld_horizontal_gradients need halo values for ibld, so must
478	! previously exist for hbl also.
479	hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) )
480	ibld(:,:) = 4
481	DO_3D_11_11( 5, jpkm1 )
482	IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
483	ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
484	ENDIF
485	END_3D
486
487	DO_2D_00_00
488	zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
489	imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 ))
490	zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
491	END_2D
492	! Averages over well-mixed and boundary layer
493	! Alan: do we need zb_nl?, zb_ml?
494	CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, Kmm, Kbb )
495	CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, Kmm, Kbb )
496	! External gradient
497	CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext, Kmm )
498
499
500	IF ( ln_osm_mle ) THEN
501	CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, Kmm )
502	CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle, Kmm )
503	ENDIF
504
505	! Rate of change of hbl
506	CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
507	! Calculate averages over depth of boundary layer
508	DO_2D_00_00
509	zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wn here, so subtract it
510	! adjustment to represent limiting by ocean bottom
511	zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)! ht_n(:,:))
512	END_2D
513
514	imld(:,:) = ibld(:,:) ! use imld to hold previous blayer index
515	ibld(:,:) = 4
516
517
518	DO_3D_00_00( 4, jpkm1 )
519	IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
520	ibld(ji,jj) = jk
521	ENDIF
522	END_3D
523
524	!
525	! Step through model levels taking account of buoyancy change to determine the effect on dhdt
526	!
527	CALL zdf_osm_timestep_hbl( zdhdt, zdhdt_2, Kmm )
528	! Alan: do we need zb_ml?
529	CALL zdf_osm_vertical_average( ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, Kmm, Kbb )
530	!
531	!
532	CALL zdf_osm_pycnocline_thickness( dh, zdh, Kmm )
533	dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. ) ! Limit delta for shallow boundary layers for calculating flux-gradient terms.
534	!
535	! Average over the depth of the mixed layer in the convective boundary layer
536	! Alan: do we need zb_ml?
537	CALL zdf_osm_vertical_average( imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml, Kmm, Kbb )
538	! rotate mean currents and changes onto wind align co-ordinates
539	!
540	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
541	CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
542	zuw_bse = 0._wp
543	zvw_bse = 0._wp
544	zwth_ent = 0._wp
545	zws_ent = 0._wp
546	DO_2D_00_00
547	IF ( ibld(ji,jj) < mbkt(ji,jj) ) THEN
548	IF ( lconv(ji,jj) ) THEN
549	zuw_bse(ji,jj) = -0.0075((zvstr(ji,jj)3+0.5zwstrc(ji,jj)3)pthird*zdu_ml(ji,jj) + &
550	& 1.5zustar(ji,jj)2(zhbl(ji,jj)-zhml(ji,jj)) )/ &
551	& ( zhml(ji,jj)MIN(zla(ji,jj)*(8./3.),1.) + epsln)
552	zvw_bse(ji,jj) = 0.01(-(zvstr(ji,jj)3+0.5zwstrc(ji,jj)3)pthird*zdv_ml(ji,jj)+ &
553	& 2.0ff_t(ji,jj)zustke(ji,jj)dstokes(ji,jj)zla(ji,jj))
554	IF ( zdb_ml(ji,jj) > 0._wp ) THEN
555	zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
556	zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
557	ENDIF
558	ELSE
559	zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )
560	zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )
561	ENDIF
562	ENDIF
563	END_2D
564
565	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
566	! Pycnocline gradients for scalars and velocity
567	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
568
569	CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext, Kmm )
570	CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, Kmm )
571	CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc, Kmm )
572	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
573	! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
574	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
575
576	DO_2D_00_00
577	IF ( lconv(ji,jj) ) THEN
578	zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
579	zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
580	zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zdh(ji,jj)
581	zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zdh(ji,jj)
582	zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third
583	zbeta_v_sc(ji,jj) = 1.0 - 2.0 * (0.142 /0.375) * zdh(ji,jj) / ( zhml(ji,jj) + epsln )
584	ELSE
585	zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
586	zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
587	END IF
588	END_2D
589	!
590	DO_2D_00_00
591	IF ( lconv(ji,jj) ) THEN
592	DO jk = 2, imld(ji,jj) ! mixed layer diffusivity
593	zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
594	!
595	zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
596	!
597	zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) &
598	& * ( 1.0 - 0.5 * zznd_ml**2 )
599	END DO
600	! pycnocline - if present linear profile
601	IF ( zdh(ji,jj) > 0._wp ) THEN
602	zgamma_b = 6.0
603	DO jk = imld(ji,jj)+1 , ibld(ji,jj)
604	zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
605	!
606	zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
607	!
608	zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * EXP( zgamma_b * zznd_pyc )
609	END DO
610	IF ( ibld_ext == 0 ) THEN
611	zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
612	zviscos(ji,jj,ibld(ji,jj)) = 0._wp
613	ELSE
614	zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw(ji, jj, ibld(ji,jj)-1, Kmm) )
615	zviscos(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj) * ( hbl(ji,jj) - gdepw(ji, jj, ibld(ji,jj)-1, Kmm) )
616	ENDIF
617	ENDIF
618	! Temporary fix to ensure zdiffut is +ve; won't be necessary with ww taken out
619	zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t(ji,jj,ibld(ji,jj), Kmm), 1.e-6)
620	zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t(ji,jj,ibld(ji,jj), Kmm), 1.e-5)
621	! could be taken out, take account of entrainment represents as a diffusivity
622	! should remove w from here, represents entrainment
623	ELSE
624	! stable conditions
625	DO jk = 2, ibld(ji,jj)
626	zznd_ml = gdepw(ji,jj,jk, Kmm) / zhbl(ji,jj)
627	zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
628	zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
629	END DO
630
631	IF ( ibld_ext == 0 ) THEN
632	zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
633	zviscos(ji,jj,ibld(ji,jj)) = 0._wp
634	ELSE
635	zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w(ji, jj, ibld(ji,jj), Kmm)
636	zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 0._wp) * e3w(ji, jj, ibld(ji,jj), Kmm)
637	ENDIF
638	ENDIF ! end if ( lconv )
639	!
640	END_2D
641
642	!
643	! calculate non-gradient components of the flux-gradient relationships
644	!
645	! Stokes term in scalar flux, flux-gradient relationship
646	WHERE ( lconv )
647	zsc_wth_1 = zwstrl*3 zwth0 / ( zvstr*3 + 0.5 zwstrc**3 + epsln)
648	!
649	zsc_ws_1 = zwstrl*3 zws0 / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
650	ELSEWHERE
651	zsc_wth_1 = 2.0 * zwthav
652	!
653	zsc_ws_1 = 2.0 * zwsav
654	ENDWHERE
655
656
657	DO_2D_00_00
658	IF ( lconv(ji,jj) ) THEN
659	DO jk = 2, imld(ji,jj)
660	zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
661	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
662	!
663	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
664	END DO ! end jk loop
665	ELSE ! else for if (lconv)
666	! Stable conditions
667	DO jk = 2, ibld(ji,jj)
668	zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
669	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
670	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
671	!
672	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
673	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
674	END DO
675	ENDIF ! endif for check on lconv
676
677	END_2D
678
679	! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero)
680	WHERE ( lconv )
681	zsc_uw_1 = ( zwstrl*3 + 0.5 zwstrc3 )pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 )
682	zsc_uw_2 = ( zwstrl*3 + 0.5 zwstrc3 )pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 )
683	zsc_vw_1 = ff_t * zhml * zustke*3 MIN( zla(8.0/3.0), 0.12 ) / ( ( zvstr3 + 0.5 * zwstrc3 )(2.0/3.0) + epsln )
684	ELSEWHERE
685	zsc_uw_1 = zustar**2
686	zsc_vw_1 = ff_t * zhbl * zustke*3 MIN( zla(8.0/3.0), 0.12 ) / (zvstr2 + epsln)
687	ENDWHERE
688	IF(ln_dia_osm) THEN
689	IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
690	IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
691	END IF
692	DO_2D_00_00
693	IF ( lconv(ji,jj) ) THEN
694	DO jk = 2, imld(ji,jj)
695	zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
696	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) &
697	& + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) &
698	& * ( 1.0 - EXP ( -2.0 * zznd_d ) )
699	!
700	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( - zznd_d ) &
701	& * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
702	END DO ! end jk loop
703	ELSE
704	! Stable conditions
705	DO jk = 2, ibld(ji,jj) ! corrected to ibld
706	zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
707	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) &
708	& * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
709	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
710	END DO ! end jk loop
711	ENDIF
712	END_2D
713
714	! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
715
716	WHERE ( lconv )
717	zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
718	zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
719	ELSEWHERE
720	zsc_wth_1 = 0._wp
721	zsc_ws_1 = 0._wp
722	ENDWHERE
723
724	DO_2D_00_00
725	IF (lconv(ji,jj) ) THEN
726	DO jk = 2, imld(ji,jj)
727	zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
728	! calculate turbulent length scale
729	zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) &
730	& * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) )
731	zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) &
732	& * ( 1.0 - EXP ( - 5.0 * ( 1.0 - zznd_ml ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
733	zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
734	! non-gradient buoyancy terms
735	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
736	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 * zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
737	END DO
738	ELSE
739	DO jk = 2, ibld(ji,jj)
740	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
741	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj)
742	END DO
743	ENDIF
744	END_2D
745
746	WHERE ( lconv )
747	zsc_uw_1 = -zwb0 * zustar*2 zhml / ( zvstr*3 + 0.5 zwstrc**3 + epsln )
748	zsc_uw_2 = zwb0 * zustke * zhml / ( zvstr*3 + 0.5 zwstrc3 + epsln )(2.0/3.0)
749	zsc_vw_1 = 0._wp
750	ELSEWHERE
751	zsc_uw_1 = 0._wp
752	zsc_vw_1 = 0._wp
753	ENDWHERE
754
755	DO_2D_00_00
756	IF ( lconv(ji,jj) ) THEN
757	DO jk = 2 , imld(ji,jj)
758	zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
759	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) &
760	& * ( 1.0 - EXP( -0.5 * zznd_d ) ) &
761	& * zsc_uw_2(ji,jj) )
762	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
763	END DO ! jk loop
764	ELSE
765	! stable conditions
766	DO jk = 2, ibld(ji,jj)
767	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
768	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
769	END DO
770	ENDIF
771	END_2D
772
773	IF(ln_dia_osm) THEN
774	IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
775	IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
776	END IF
777	! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
778
779	WHERE ( lconv )
780	zsc_wth_1 = zwth0
781	zsc_ws_1 = zws0
782	ELSEWHERE
783	zsc_wth_1 = 2.0 * zwthav
784	zsc_ws_1 = zws0
785	ENDWHERE
786
787	DO_2D_00_00
788	IF ( lconv(ji,jj) ) THEN
789	DO jk = 2, imld(ji,jj)
790	zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
791	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) &
792	& * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) &
793	& - EXP( - 6.0 * zznd_ml ) ) ) &
794	& * ( 1.0 - EXP( - 15.0 * ( 1.0 - zznd_ml ) ) )
795	!
796	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) &
797	& * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) &
798	& - EXP( - 6.0 * zznd_ml ) ) ) &
799	& * ( 1.0 - EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
800	END DO
801	ELSE
802	DO jk = 2, ibld(ji,jj)
803	zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
804	znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
805	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
806	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
807	ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
808	& 7.5 * EXP ( -10.0 * ( 0.95 - znd )*2 ) ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
809	END DO
810	ENDIF
811	END_2D
812
813	WHERE ( lconv )
814	zsc_uw_1 = zustar**2
815	zsc_vw_1 = ff_t * zustke * zhml
816	ELSEWHERE
817	zsc_uw_1 = zustar**2
818	zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
819	zsc_vw_1 = ff_t * zustke * zhbl
820	zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )*2 ) zsc_vw_1
821	ENDWHERE
822
823	DO_2D_00_00
824	IF ( lconv(ji,jj) ) THEN
825	DO jk = 2, imld(ji,jj)
826	zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
827	zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
828	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
829	& + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml*4 ) - EXP ( -8.0 zznd_ml ) ) * zsc_uw_1(ji,jj)
830	!
831	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
832	& + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
833	END DO
834	ELSE
835	DO jk = 2, ibld(ji,jj)
836	znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
837	zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
838	IF ( zznd_d <= 2.0 ) THEN
839	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
840	&* ( 2.25 - 3.0 * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
841	!
842	ELSE
843	ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
844	& + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
845	!
846	ENDIF
847
848	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
849	& + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d*2 ) zsc_vw_1(ji,jj)
850	ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
851	& + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
852	END DO
853	ENDIF
854	END_2D
855
856	IF(ln_dia_osm) THEN
857	IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
858	IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
859	IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
860	IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
861	IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
862	IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
863	END IF
864	!
865	! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
866
867	DO_2D_00_00
868	IF ( lconv(ji,jj) ) THEN
869	DO jk = 2, ibld(ji,jj)
870	znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
871	IF ( znd >= 0.0 ) THEN
872	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
873	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
874	ELSE
875	ghamu(ji,jj,jk) = 0._wp
876	ghamv(ji,jj,jk) = 0._wp
877	ENDIF
878	END DO
879	ELSE
880	DO jk = 2, ibld(ji,jj)
881	znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
882	IF ( znd >= 0.0 ) THEN
883	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
884	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
885	ELSE
886	ghamu(ji,jj,jk) = 0._wp
887	ghamv(ji,jj,jk) = 0._wp
888	ENDIF
889	END DO
890	ENDIF
891	END_2D
892
893	IF(ln_dia_osm) THEN
894	IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
895	IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
896	END IF
897	! pynocline contributions
898	! Temporary fix to avoid instabilities when zdb_bl becomes very very small
899	zsc_uw_1 = 0._wp ! 50.0 * zla*(8.0/3.0) zustar*2 zhbl / ( zdb_bl + epsln )
900	DO_2D_00_00
901	IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
902	DO jk= 2, ibld(ji,jj)
903	znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
904	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
905	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
906	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
907	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )*(7.0/4.0) zdbdz_pyc(ji,jj,jk)
908	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
909	END DO
910	END IF
911	END_2D
912
913	! Entrainment contribution.
914
915	DO_2D_00_00
916	IF ( lconv(ji,jj) .AND. ibld(ji,jj) + ibld_ext < mbkt(ji,jj)) THEN
917	DO jk = 1, imld(ji,jj) - 1
918	znd=gdepw(ji,jj,jk,kmm) / zhml(ji,jj)
919	! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
920	! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
921	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd
922	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd
923	END DO
924	DO jk = imld(ji,jj), ibld(ji,jj)
925	znd = -( gdepw(ji,jj,jk,kmm) - zhml(ji,jj) ) / zdh(ji,jj)
926	! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
927	! ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
928	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd )
929	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd )
930	END DO
931	ENDIF
932
933	ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
934	ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
935	ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
936	ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
937	END_2D
938
939	IF(ln_dia_osm) THEN
940	IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
941	IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
942	IF ( iom_use("zuw_bse") ) CALL iom_put( "zuw_bse", tmask(:,:,1)*zuw_bse )
943	IF ( iom_use("zvw_bse") ) CALL iom_put( "zvw_bse", tmask(:,:,1)*zvw_bse )
944	IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
945	IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
946	IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
947	END IF
948	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
949	! Need to put in code for contributions that are applied explicitly to
950	! the prognostic variables
951	! 1. Entrainment flux
952	!
953	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
954
955
956
957	! rotate non-gradient velocity terms back to model reference frame
958
959	DO_2D_00_00
960	DO jk = 2, ibld(ji,jj)
961	ztemp = ghamu(ji,jj,jk)
962	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
963	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
964	END DO
965	END_2D
966
967	IF(ln_dia_osm) THEN
968	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
969	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
970	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
971	END IF
972
973	! KPP-style Ri# mixing
974	IF( ln_kpprimix) THEN
975	DO_3D_10_10( 2, jpkm1 )
976	z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) &
977	& * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) &
978	& / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) )
979	z3dv(ji,jj,jk) = 0.5 * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj ,jk,Kmm) ) &
980	& * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj ,jk,Kbb) ) * wvmask(ji,jj,jk) &
981	& / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) )
982	END_3D
983	!
984	DO_3D_00_00( 2, jpkm1 )
985	! ! shear prod. at w-point weightened by mask
986	zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
987	& + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
988	! ! local Richardson number
989	zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
990	zfri = MIN( zri / rn_riinfty , 1.0_wp )
991	zfri = ( 1.0_wp - zfri * zfri )
992	zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk)
993	END_3D
994
995	DO_2D_00_00
996	DO jk = ibld(ji,jj) + 1, jpkm1
997	zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
998	zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
999	END DO
1000	END_2D
1001
1002	END IF ! ln_kpprimix = .true.
1003
1004	! KPP-style set diffusivity large if unstable below BL
1005	IF( ln_convmix) THEN
1006	DO_2D_00_00
1007	DO jk = ibld(ji,jj) + 1, jpkm1
1008	IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
1009	END DO
1010	END_2D
1011	END IF ! ln_convmix = .true.
1012
1013
1014
1015	IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing
1016	DO_2D_00_00
1017	IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL.
1018	! Calculate MLE flux profiles
1019	! DO jk = 1, mld_prof(ji,jj)
1020	! znd = - gdepw(ji,jj,jk, Kmm ) / MAX(zhmle(ji,jj),epsln)
1021	! ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + &
1022	! & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )*2 ) ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 )
1023	! ghams(ji,jj,jk) = ghams(ji,jj,jk) + &
1024	! & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )*2 ) ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 )
1025	! END DO
1026	! Calculate MLE flux contribution from surface fluxes
1027	DO jk = 1, ibld(ji,jj)
1028	znd = gdepw(ji,jj,jk, Kmm ) / MAX(zhbl(ji,jj),epsln)
1029	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd )
1030	ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd )
1031	END DO
1032	DO jk = 1, mld_prof(ji,jj)
1033	znd = gdepw(ji,jj,jk, Kmm ) / MAX(zhmle(ji,jj),epsln)
1034	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd )
1035	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd )
1036	END DO
1037	! Viscosity for MLEs
1038	DO jk = ibld(ji,jj), mld_prof(ji,jj)
1039	zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) )
1040	END DO
1041	! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj)
1042	jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj))
1043	jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t(ji,jj,jl, Kmm ) ), 2 ) + mld_prof(ji,jj), mbkt(ji,jj))
1044	DO jk = mld_prof(ji,jj), jl
1045	zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * &
1046	& ( gdepw(ji,jj,jk, Kmm ) - gdepw(ji,jj,jl, Kmm ) ) / &
1047	& ( gdepw(ji,jj,mld_prof(ji,jj), Kmm ) - gdepw(ji,jj,jl, Kmm ) - epsln))
1048	END DO
1049	ENDIF
1050	END_2D
1051	ENDIF
1052
1053	IF(ln_dia_osm) THEN
1054	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
1055	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
1056	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
1057	END IF
1058
1059
1060	! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
1061	!CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
1062
1063	! GN 25/8: need to change tmask --> wmask
1064
1065	DO_3D_00_00( 2, jpkm1 )
1066	p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
1067	p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
1068	END_3D
1069	! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids
1070	CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1., &
1071	& ghamu, 'W', 1. , ghamv, 'W', 1. )
1072	DO_3D_00_00( 2, jpkm1 )
1073	ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
1074	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
1075
1076	ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
1077	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
1078
1079	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk)
1080	ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk)
1081	END_3D
1082	! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)
1083	CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
1084	! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged)
1085	! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged)
1086	CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., &
1087	& ghamu, 'U', -1. , ghamv, 'V', -1. )
1088
1089	IF(ln_dia_osm) THEN
1090	SELECT CASE (nn_osm_wave)
1091	! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1).
1092	CASE(0:1)
1093	IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)zustkezcos_wind ) ! x surface Stokes drift
1094	IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)zustkezsin_wind ) ! y surface Stokes drift
1095	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rho0tmask(:,:,1)zustar2zustke )
1096	! Stokes drift read in from sbcwave (=2).
1097	CASE(2)
1098	IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift
1099	IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift
1100	IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period
1101	IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height
1102	IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.rpi1.026/(0.877grav) )wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum
1103	IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)wndm2tmask(:,:,1) ) ! significant wave height from NP spectrum
1104	IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10
1105	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rho0tmask(:,:,1)zustar2 &
1106	& SQRT(ut0sd2 + vt0sd2 ) )
1107	END SELECT
1108	IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL>
1109	IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL>
1110	IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL>
1111	IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL>
1112	IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! <Tw_0>
1113	IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0>
1114	IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth
1115	IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k
1116	IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base
1117	IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base
1118	IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base
1119	IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base
1120	IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base
1121	IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth
1122	IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth
1123	IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth
1124	IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points
1125	IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc ) ! convective velocity scale
1126	IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale
1127	IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale
1128	IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale
1129	IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir #
1130	IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.rho0tmask(:,:,1)zustar*3 ) ! BL depth internal to zdf_osm routine
1131	IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.rho0tmask(:,:,1)zustar2zustke )
1132	IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine
1133	IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine
1134	IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine
1135	IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine
1136	IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine
1137	IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux
1138	IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux
1139	IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux
1140	IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux
1141	IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML
1142
1143	IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth
1144	IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth
1145	IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux
1146	IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML
1147	IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k
1148	IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt
1149	IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt
1150	IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt
1151	IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt
1152	IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt
1153	IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt
1154	IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1155	IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
1156
1157	END IF
1158
1159	CONTAINS
1160
1161
1162	! Alan: do we need zb?
1163	SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, zu, zv, zdt, zds, zdb, zdu, zdv, Kmm, Kbb )
1164	!!---------------------------------------------------------------------
1165	!! * ROUTINE zdf_vertical_average *
1166	!!
1167	!! ** Purpose : Determines vertical averages from surface to jnlev.
1168	!!
1169	!! ** Method : Averages are calculated from the surface to jnlev.
1170	!! The external level used to calculate differences is ibld+ibld_ext
1171	!!
1172	!!----------------------------------------------------------------------
1173
1174	INTEGER, INTENT(in) :: Kmm, Kbb ! time-level indices
1175	INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over.
1176
1177	! Alan: do we need zb?
1178	REAL(wp), DIMENSION(jpi,jpj) :: zt, zs ! Average temperature and salinity
1179	REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components
1180	REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
1181	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components.
1182
1183	INTEGER :: jk, ji, jj
1184	REAL(wp) :: zthick, zthermal, zbeta
1185
1186
1187	zt = 0._wp
1188	zs = 0._wp
1189	zu = 0._wp
1190	zv = 0._wp
1191	DO_2D_00_00
1192	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1193	zbeta = rab_n(ji,jj,1,jp_sal)
1194	! average over depth of boundary layer
1195	zthick = epsln
1196	DO jk = 2, jnlev_av(ji,jj)
1197	zthick = zthick + e3t(ji,jj,jk, Kmm )
1198	zt(ji,jj) = zt(ji,jj) + e3t(ji,jj,jk, Kmm ) * ts(ji,jj,jk,jp_tem, Kmm )
1199	zs(ji,jj) = zs(ji,jj) + e3t(ji,jj,jk, Kmm ) * ts(ji,jj,jk,jp_sal, Kmm )
1200	zu(ji,jj) = zu(ji,jj) + e3t(ji,jj,jk, Kmm ) &
1201	& * ( uu(ji,jj,jk, Kbb ) + uu(ji - 1,jj,jk, Kbb ) ) &
1202	& / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
1203	zv(ji,jj) = zv(ji,jj) + e3t(ji,jj,jk, Kmm ) &
1204	& * ( vv(ji,jj,jk, Kbb ) + vv(ji,jj - 1,jk, Kbb ) ) &
1205	& / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
1206	END DO
1207	zt(ji,jj) = zt(ji,jj) / zthick
1208	zs(ji,jj) = zs(ji,jj) / zthick
1209	zu(ji,jj) = zu(ji,jj) / zthick
1210	zv(ji,jj) = zv(ji,jj) / zthick
1211	! Alan: do we need zb?
1212	zdt(ji,jj) = zt(ji,jj) - ts(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem, Kmm )
1213	zds(ji,jj) = zs(ji,jj) - ts(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal, Kmm )
1214	zdu(ji,jj) = zu(ji,jj) - ( uu(ji,jj,ibld(ji,jj)+ibld_ext, Kbb ) + uu(ji-1,jj,ibld(ji,jj)+ibld_ext, Kbb ) ) &
1215	& / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) )
1216	zdv(ji,jj) = zv(ji,jj) - ( vv(ji,jj,ibld(ji,jj)+ibld_ext, Kbb ) + vv(ji,jj-1,ibld(ji,jj)+ibld_ext, Kbb ) ) &
1217	& / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) )
1218	zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
1219	END_2D
1220	END SUBROUTINE zdf_osm_vertical_average
1221
1222	SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
1223	!!---------------------------------------------------------------------
1224	!! * ROUTINE zdf_velocity_rotation *
1225	!!
1226	!! ** Purpose : Rotates frame of reference of averaged velocity components.
1227	!!
1228	!! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w
1229	!!
1230	!!----------------------------------------------------------------------
1231
1232	REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle
1233	REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current
1234	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline
1235
1236	INTEGER :: ji, jj
1237	REAL(wp) :: ztemp
1238
1239	DO_2D_00_00
1240	ztemp = zu(ji,jj)
1241	zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
1242	zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
1243	ztemp = zdu(ji,jj)
1244	zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
1245	zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj)
1246	END_2D
1247	END SUBROUTINE zdf_osm_velocity_rotation
1248
1249	SUBROUTINE zdf_osm_external_gradients( zdtdz, zdsdz, zdbdz, Kmm )
1250	!!---------------------------------------------------------------------
1251	!! * ROUTINE zdf_osm_external_gradients *
1252	!!
1253	!! ** Purpose : Calculates the gradients below the OSBL
1254	!!
1255	!! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient.
1256	!!
1257	!!----------------------------------------------------------------------
1258
1259	INTEGER, INTENT(in) :: Kmm ! time-level index
1260	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy.
1261
1262	INTEGER :: jj, ji, jkb, jkb1
1263	REAL(wp) :: zthermal, zbeta
1264
1265
1266	DO_2D_00_00
1267	IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1268	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
1269	zbeta = rab_n(ji,jj,1,jp_sal)
1270	jkb = ibld(ji,jj) + ibld_ext
1271	jkb1 = MIN(jkb + 1, mbkt(ji,jj))
1272	zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem, Kmm ) - ts(ji,jj,jkb,jp_tem, Kmm ) ) &
1273	& / e3t(ji,jj,ibld(ji,jj), Kmm )
1274	zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal, Kmm ) - ts(ji,jj,jkb,jp_sal, Kmm ) ) &
1275	& / e3t(ji,jj,ibld(ji,jj), Kmm )
1276	zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
1277	END IF
1278	END_2D
1279	END SUBROUTINE zdf_osm_external_gradients
1280
1281	SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, Kmm )
1282
1283	INTEGER, INTENT(in) :: Kmm ! time-level index
1284	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline
1285
1286	INTEGER :: jk, jj, ji
1287	REAL(wp) :: ztgrad, zsgrad, zbgrad
1288	REAL(wp) :: zgamma_b_nd, zgamma_c, znd
1289	REAL(wp) :: zzeta_s=0.3, ztmp
1290
1291	DO_2D_00_00
1292	IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1293	IF ( lconv(ji,jj) ) THEN ! convective conditions
1294	IF ( zdbdz_ext(ji,jj) > 0._wp .AND. &
1295	& (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh &
1296	& .OR. zdb_bl(ji,jj) > 0._wp)) THEN ! zdhdt could be <0 due to FK, hence check zdhdt>0
1297	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
1298	ztgrad = 0.5 * zdt_ml(ji,jj) * ztmp + zdtdz_ext(ji,jj)
1299	zsgrad = 0.5 * zds_ml(ji,jj) * ztmp + zdsdz_ext(ji,jj)
1300	zbgrad = 0.5 * zdb_ml(ji,jj) * ztmp + zdbdz_ext(ji,jj)
1301	zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
1302	zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
1303	& ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
1304	DO jk = 2, ibld(ji,jj)+ibld_ext
1305	znd = -( gdepw(ji,jj,jk, Kmm ) - zhbl(ji,jj) ) * ztmp
1306	IF ( znd <= zzeta_s ) THEN
1307	zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) * ztmp * &
1308	& EXP( -6.0 * ( znd -zzeta_s )**2 )
1309	zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) * ztmp * &
1310	& EXP( -6.0 * ( znd -zzeta_s )**2 )
1311	zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) * ztmp * &
1312	& EXP( -6.0 * ( znd -zzeta_s )**2 )
1313	ELSE
1314	zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1315	zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1316	zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
1317	ENDIF
1318	END DO
1319	ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing
1320	ELSE
1321	! stable conditions
1322	! if pycnocline profile only defined when depth steady of increasing.
1323	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! Depth increasing, or steady.
1324	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1325	IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline
1326	ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
1327	ztgrad = zdt_bl(ji,jj) * ztmp
1328	zsgrad = zds_bl(ji,jj) * ztmp
1329	zbgrad = zdb_bl(ji,jj) * ztmp
1330	DO jk = 2, ibld(ji,jj)
1331	znd = gdepw(ji,jj,jk, Kmm ) * ztmp
1332	zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1333	zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1334	zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
1335	END DO
1336	ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
1337	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
1338	ztgrad = zdt_bl(ji,jj) * ztmp
1339	zsgrad = zds_bl(ji,jj) * ztmp
1340	zbgrad = zdb_bl(ji,jj) * ztmp
1341	DO jk = 2, ibld(ji,jj)
1342	znd = -( gdepw(ji,jj,jk, Kmm ) - zhml(ji,jj) ) * ztmp
1343	zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1344	zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1345	zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
1346	END DO
1347	ENDIF ! IF (zhol >=0.5)
1348	ENDIF ! IF (zdb_bl> 0.)
1349	ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
1350	ENDIF ! IF (lconv)
1351	END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1352	END_2D
1353
1354	END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
1355
1356	SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz, Kmm )
1357	!!---------------------------------------------------------------------
1358	!! * ROUTINE zdf_osm_pycnocline_shear_profiles *
1359	!!
1360	!! ** Purpose : Calculates velocity shear in the pycnocline
1361	!!
1362	!! ** Method :
1363	!!
1364	!!----------------------------------------------------------------------
1365
1366	INTEGER, INTENT(in) :: Kmm ! time-level index
1367	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
1368
1369	INTEGER :: jk, jj, ji
1370	REAL(wp) :: zugrad, zvgrad, znd
1371	REAL(wp) :: zzeta_v = 0.45
1372	!
1373	DO_2D_00_00
1374	!
1375	IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
1376	IF ( lconv (ji,jj) ) THEN
1377	! Unstable conditions
1378	zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
1379	& ( ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zhml(ji,jj) ) * &
1380	& MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
1381	!Alan is this right?
1382	zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
1383	& 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
1384	& ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird + epsln ) &
1385	& )/ (zdh(ji,jj) + epsln )
1386	DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
1387	znd = -( gdepw(ji,jj,jk, Kmm ) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
1388	IF ( znd <= 0.0 ) THEN
1389	zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
1390	zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
1391	ELSE
1392	zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
1393	zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
1394	ENDIF
1395	END DO
1396	ELSE
1397	! stable conditions
1398	zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
1399	zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
1400	DO jk = 2, ibld(ji,jj)
1401	znd = gdepw(ji,jj,jk, Kmm ) / zhbl(ji,jj)
1402	IF ( znd < 1.0 ) THEN
1403	zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
1404	ELSE
1405	zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
1406	ENDIF
1407	zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
1408	END DO
1409	ENDIF
1410	!
1411	END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
1412	END_2D
1413	END SUBROUTINE zdf_osm_pycnocline_shear_profiles
1414
1415	SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
1416	!!---------------------------------------------------------------------
1417	!! * ROUTINE zdf_osm_calculate_dhdt *
1418	!!
1419	!! ** Purpose : Calculates the rate at which hbl changes.
1420	!!
1421	!! ** Method :
1422	!!
1423	!!----------------------------------------------------------------------
1424
1425	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! Rate of change of hbl
1426
1427	INTEGER :: jj, ji
1428	REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert
1429	REAL(wp) :: zvel_max, zwb_min
1430	REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
1431	REAL(wp) :: zzeta_m = 0.3
1432	REAL(wp) :: zgamma_c = 2.0
1433	REAL(wp) :: zdhoh = 0.1
1434	REAL(wp) :: alpha_bc = 0.5
1435
1436	DO_2D_00_00
1437	IF ( lconv(ji,jj) ) THEN ! Convective
1438	! Alan is this right? Yes, it's a bit different from the previous relationship
1439	! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) &
1440	! & - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)*3 + 0.03 zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1441	zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
1442	zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
1443	zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
1444	zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
1445	zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
1446	& * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
1447
1448	zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
1449	& - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
1450	& + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
1451	& - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
1452	!
1453	zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj)
1454
1455	IF ( ln_osm_mle ) THEN
1456	! zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + &
1457	! & ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 ) ! Fox-Kemper buoyancy flux average over OSBL
1458	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
1459	zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * &
1460	(1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
1461	ELSE
1462	zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
1463	ENDIF
1464	zvel_max = ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
1465	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
1466	! OSBL is deepening, entrainment > restratification
1467	IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1468	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1469	ELSE
1470	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)
1471	ENDIF
1472	ELSE
1473	! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
1474	zdhdt(ji,jj) = - zvel_mle(ji,jj)
1475
1476
1477	ENDIF
1478
1479	ELSE
1480	! Fox-Kemper not used.
1481
1482	zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
1483	& ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
1484	zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
1485	! added ajgn 23 July as temporay fix
1486
1487	ENDIF
1488
1489	zdhdt_2(ji,jj) = 0._wp
1490
1491	! commented out ajgn 23 July as temporay fix
1492	! IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
1493	! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary.
1494	! zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1495	! zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj)
1496	! zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
1497	! zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj)
1498	! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min )
1499	! ! Alan no idea what this should be?
1500	! zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) &
1501	! & + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) &
1502	! & * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj)
1503	! zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
1504	! IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN
1505	! zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj)
1506	! ENDIF
1507	ELSE ! Stable
1508	zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
1509	zdhdt_2(ji,jj) = 0._wp
1510	IF ( zdhdt(ji,jj) < 0._wp ) THEN
1511	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
1512	zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
1513	ELSE
1514	zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
1515	ENDIF
1516	zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert
1517	ENDIF
1518	END_2D
1519	END SUBROUTINE zdf_osm_calculate_dhdt
1520
1521	SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2, Kmm )
1522	!!---------------------------------------------------------------------
1523	!! * ROUTINE zdf_osm_timestep_hbl *
1524	!!
1525	!! ** Purpose : Increments hbl.
1526	!!
1527	!! ** Method : If thechange in hbl exceeds one model level the change is
1528	!! is calculated by moving down the grid, changing the buoyancy
1529	!! jump. This is to ensure that the change in hbl does not
1530	!! overshoot a stable layer.
1531	!!
1532	!!----------------------------------------------------------------------
1533
1534
1535	INTEGER, INTENT(in) :: Kmm ! time-level index
1536	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2 ! rates of change of hbl.
1537
1538	INTEGER :: jk, jj, ji, jm
1539	REAL(wp) :: zhbl_s, zvel_max, zdb
1540	REAL(wp) :: zthermal, zbeta
1541
1542	DO_2D_00_00
1543	IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
1544	!
1545	! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
1546	!
1547	zhbl_s = hbl(ji,jj)
1548	jm = imld(ji,jj)
1549	zthermal = rab_n(ji,jj,1,jp_tem)
1550	zbeta = rab_n(ji,jj,1,jp_sal)
1551
1552
1553	IF ( lconv(ji,jj) ) THEN
1554	!unstable
1555
1556	IF( ln_osm_mle ) THEN
1557	zvel_max = ( zwstrl(ji,jj)3 + zwstrc(ji,jj)3 )**p2third / hbl(ji,jj)
1558	ELSE
1559
1560	zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
1561	& ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
1562
1563	ENDIF
1564
1565	DO jk = imld(ji,jj), ibld(ji,jj)
1566	zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) &
1567	& - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), &
1568	& 0.0 ) + zvel_max
1569
1570
1571	IF ( ln_osm_mle ) THEN
1572	zhbl_s = zhbl_s + MIN( &
1573	& rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
1574	& e3w(ji,jj,jm, Kmm ) )
1575	ELSE
1576	zhbl_s = zhbl_s + MIN( &
1577	& rn_Dt * ( -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
1578	& e3w(ji,jj,jm, Kmm ) )
1579	ENDIF
1580
1581	zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm ) - depth_tol)
1582
1583	IF ( zhbl_s >= gdepw(ji,jj,jm+1, Kmm) ) jm = jm + 1
1584	END DO
1585	hbl(ji,jj) = zhbl_s
1586	ibld(ji,jj) = jm
1587	ELSE
1588	! stable
1589	DO jk = imld(ji,jj), ibld(ji,jj)
1590	zdb = MAX( &
1591	& grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem, Kmm) )&
1592	& - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal, Kmm) ) ),&
1593	& 0.0 ) + &
1594	& 2.0 * zvstr(ji,jj)**2 / zhbl_s
1595
1596	! Alan is thuis right? I have simply changed hbli to hbl
1597	zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
1598	zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)*3 / zhbl_s - 0.15 / 2.0 ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * &
1599	& zustar(ji,jj)*3 / zhbl_s ) ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
1600	zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
1601	zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm, Kmm ) )
1602
1603	zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm ) - depth_tol)
1604	IF ( zhbl_s >= gdepw(ji,jj,jm, Kmm ) ) jm = jm + 1
1605	END DO
1606	ENDIF ! IF ( lconv )
1607	hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,4, Kmm ) )
1608	ibld(ji,jj) = MAX(jm, 4 )
1609	ELSE
1610	! change zero or one model level.
1611	hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4, Kmm ) )
1612	ENDIF
1613	zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj), Kmm )
1614	END_2D
1615
1616	END SUBROUTINE zdf_osm_timestep_hbl
1617
1618	SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh, Kmm )
1619	!!---------------------------------------------------------------------
1620	!! * ROUTINE zdf_osm_pycnocline_thickness *
1621	!!
1622	!! ** Purpose : Calculates thickness of the pycnocline
1623	!!
1624	!! ** Method : The thickness is calculated from a prognostic equation
1625	!! that relaxes the pycnocine thickness to a diagnostic
1626	!! value. The time change is calculated assuming the
1627	!! thickness relaxes exponentially. This is done to deal
1628	!! with large timesteps.
1629	!!
1630	!!----------------------------------------------------------------------
1631
1632	INTEGER, INTENT(in) :: Kmm ! time-level index
1633	REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness.
1634	!
1635	INTEGER :: jj, ji
1636	INTEGER :: inhml
1637	REAL(wp) :: zari, ztau, zddhdt
1638
1639
1640	DO_2D_00_00
1641	IF ( lconv(ji,jj) ) THEN
1642
1643	IF( ln_osm_mle ) THEN
1644	IF ( ( zwb_ent(ji,jj) + zwb_fk_b(ji,jj) ) < 0._wp ) THEN
1645	! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
1646	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1647	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability
1648	zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1649	(1.0 + zdb_bl(ji,jj)*2 / ( 4.5 zvstr(ji,jj)*2 zdbdz_ext(ji,jj) ) ), 0.2 )
1650	ELSE ! unstable
1651	zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1652	(1.0 + zdb_bl(ji,jj)*2 / ( 4.5 zwstrc(ji,jj)*2 zdbdz_ext(ji,jj) ) ), 0.2 )
1653	ENDIF
1654	ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird
1655	zddhdt = zari * hbl(ji,jj)
1656	ELSE
1657	ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
1658	zddhdt = 0.2 * hbl(ji,jj)
1659	ENDIF
1660	ELSE
1661	ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
1662	zddhdt = 0.2 * hbl(ji,jj)
1663	ENDIF
1664	ELSE ! ln_osm_mle
1665	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
1666	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN ! near neutral stability
1667	zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1668	(1.0 + zdb_bl(ji,jj)*2 / ( 4.5 zvstr(ji,jj)*2 zdbdz_ext(ji,jj) ) ), 0.2 )
1669	ELSE ! unstable
1670	zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
1671	(1.0 + zdb_bl(ji,jj)*2 / ( 4.5 zwstrc(ji,jj)*2 zdbdz_ext(ji,jj) ) ), 0.2 )
1672	ENDIF
1673	ztau = hbl(ji,jj) / (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird
1674	zddhdt = zari * hbl(ji,jj)
1675	ELSE
1676	ztau = hbl(ji,jj) / ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
1677	zddhdt = 0.2 * hbl(ji,jj)
1678	ENDIF
1679
1680	END IF ! ln_osm_mle
1681
1682	dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_Dt / ztau ) )
1683	! Alan: this hml is never defined or used
1684	ELSE ! IF (lconv)
1685	ztau = hbl(ji,jj) / zvstr(ji,jj)
1686	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
1687	! boundary layer deepening
1688	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
1689	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
1690	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
1691	& / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 )
1692	zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj)
1693	ELSE
1694	zddhdt = 0.2 * hbl(ji,jj)
1695	ENDIF
1696	ELSE ! IF(dhdt < 0)
1697	zddhdt = 0.2 * hbl(ji,jj)
1698	ENDIF ! IF (dhdt >= 0)
1699	dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_Dt / ztau ) )
1700	IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt ! can be a problem with dh>hbl for rapid collapse
1701	! Alan: this hml is never defined or used -- do we need it?
1702	ENDIF ! IF (lconv)
1703
1704	hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
1705	inhml = MAX( INT( dh(ji,jj) / e3t(ji,jj,ibld(ji,jj), Kmm) ) , 1 )
1706	imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
1707	zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj), Kmm)
1708	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
1709	END_2D
1710
1711	END SUBROUTINE zdf_osm_pycnocline_thickness
1712
1713
1714	SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, Kmm )
1715	!!----------------------------------------------------------------------
1716	!! * ROUTINE zdf_osm_horizontal_gradients *
1717	!!
1718	!! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.
1719	!!
1720	!! ** Method :
1721	!!
1722	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
1723	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
1724
1725	INTEGER, INTENT(in) :: Kmm ! time-level index
1726	REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points
1727	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == !
1728	REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy
1729
1730	INTEGER :: ji, jj, jk ! dummy loop indices
1731	INTEGER :: ii, ij, ik, ikmax ! local integers
1732	REAL(wp) :: zc
1733	REAL(wp) :: zN2_c ! local buoyancy difference from 10m value
1734	REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH
1735	REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv
1736	REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv
1737	!!----------------------------------------------------------------------
1738	!
1739	! !== MLD used for MLE ==!
1740
1741	mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point
1742	zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2
1743	zN2_c = grav * rn_osm_mle_rho_c * r1_rho0 ! convert density criteria into N^2 criteria
1744	DO_3D_11_11( nlb10, jpkm1 ) ! Mixed layer level: w-level
1745	ikt = mbkt(ji,jj)
1746	zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk, Kmm )
1747	IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
1748	END_3D
1749	DO_2D_11_11
1750	mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
1751	zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj), Kmm )
1752	END_2D
1753	! ensure mld_prof .ge. ibld
1754	!
1755	ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation
1756	!
1757	ztm(:,:) = 0._wp
1758	zsm(:,:) = 0._wp
1759	DO_3D_11_11( 1, ikmax ) ! MLD and mean buoyancy and N2 over the mixed layer
1760	zc = e3t(ji,jj,jk, Kmm ) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points
1761	ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem, Kmm )
1762	zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal, Kmm )
1763	END_3D
1764	! average temperature and salinity.
1765	ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1, Kmm), zmld(:,:) )
1766	zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1, Kmm), zmld(:,:) )
1767	! calculate horizontal gradients at u & v points
1768
1769	DO_2D_00_10
1770	zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
1771	zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
1772	zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))
1773	ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )
1774	ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )
1775	END_2D
1776
1777	DO_2D_10_00
1778	zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
1779	zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
1780	zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))
1781	ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )
1782	ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )
1783	END_2D
1784
1785	CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm )
1786	CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm )
1787
1788	DO_2D_00_10
1789	dbdx_mle(ji,jj) = grav(zdtdx(ji,jj)zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
1790	END_2D
1791	DO_2D_10_00
1792	dbdy_mle(ji,jj) = grav(zdtdy(ji,jj)zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
1793	END_2D
1794
1795	END SUBROUTINE zdf_osm_zmld_horizontal_gradients
1796	SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle, Kmm )
1797	!!----------------------------------------------------------------------
1798	!! * ROUTINE zdf_osm_mle_parameters *
1799	!!
1800	!! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.
1801	!!
1802	!! ** Method :
1803	!!
1804	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
1805	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
1806
1807	INTEGER, INTENT(in) :: Kmm ! time-level index
1808	REAL(wp), DIMENSION(jpi,jpj) :: hmle, zwb_fk, zvel_mle, zdiff_mle
1809	INTEGER :: ji, jj, jk ! dummy loop indices
1810	INTEGER :: ii, ij, ik ! local integers
1811	INTEGER , DIMENSION(jpi,jpj) :: inml_mle
1812	REAL(wp) :: zdb_mle, ztmp, zdbds_mle
1813
1814	mld_prof(:,:) = 4
1815	DO_3D_00_00( 5, jpkm1 )
1816	IF ( hmle(ji,jj) >= gdepw(ji,jj,jk, Kmm ) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
1817	END_3D
1818	! DO_2D_00_10
1819	! zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj), Kmm )
1820	! END_2D
1821	! Timestep mixed layer eddy depth.
1822	DO_2D_00_00
1823	zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rho0 ! check ALMG
1824	IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN
1825	hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_Dt / MAX( zdb_mle, rn_osm_mle_thresh ) ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy.
1826	ELSE
1827	! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10
1828	! IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN
1829	! hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_Dt / rn_osm_mle_tau
1830	! ELSE
1831	! hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_Dt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD
1832	! ENDIF
1833	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
1834	hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau
1835	ELSE
1836	hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD
1837	ENDIF
1838	ENDIF
1839	hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj))
1840	END_2D
1841
1842	mld_prof = 4
1843	DO_3D_00_00( 5, jpkm1 )
1844	IF ( hmle(ji,jj) >= gdepw(ji,jj,jk, Kmm ) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
1845	END_3D
1846	DO_2D_00_00
1847	zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj), Kmm )
1848	END_2D
1849	! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
1850
1851	DO_2D_00_00
1852	IF ( lconv(ji,jj) ) THEN
1853	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
1854	! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) &
1855	! & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) )
1856	! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) &
1857	! & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) )
1858	zdbds_mle = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &
1859	& + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
1860	zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) ztmp zdbds_mle * zdbds_mle
1861	! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
1862	zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1)
1863	zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3 / rn_osm_mle_lf
1864	ENDIF
1865	END_2D
1866	END SUBROUTINE zdf_osm_mle_parameters
1867
1868	END SUBROUTINE zdf_osm
1869
1870
1871	SUBROUTINE zdf_osm_init( Kmm )
1872	!!----------------------------------------------------------------------
1873	!! * ROUTINE zdf_osm_init *
1874	!!
1875	!! ** Purpose : Initialization of the vertical eddy diffivity and
1876	!! viscosity when using a osm turbulent closure scheme
1877	!!
1878	!! ** Method : Read the namosm namelist and check the parameters
1879	!! called at the first timestep (nit000)
1880	!!
1881	!! ** input : Namlist namosm
1882	!!----------------------------------------------------------------------
1883	INTEGER, INTENT(in) :: Kmm ! time level index (middle)
1884	!
1885	INTEGER :: ios ! local integer
1886	INTEGER :: ji, jj, jk ! dummy loop indices
1887	REAL z1_t2
1888	!!
1889	NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
1890	& ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
1891	& ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
1892	& ,ln_osm_mle
1893	! Namelist for Fox-Kemper parametrization.
1894	NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
1895	& rn_osm_mle_rho_c,rn_osm_mle_thresh,rn_osm_mle_tau
1896
1897	!!----------------------------------------------------------------------
1898	!
1899	READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
1900	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
1901
1902	READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
1903	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
1904	IF(lwm) WRITE ( numond, namzdf_osm )
1905
1906	IF(lwp) THEN ! Control print
1907	WRITE(numout,*)
1908	WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
1909	WRITE(numout,*) '~~~~~~~~~~~~'
1910	WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters'
1911	WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la
1912	WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle
1913	WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la
1914	WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0
1915	WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes
1916	WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave
1917	WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave
1918	SELECT CASE (nn_osm_wave)
1919	CASE(0)
1920	WRITE(numout,*) ' calculated assuming constant La#=0.3'
1921	CASE(1)
1922	WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves'
1923	CASE(2)
1924	WRITE(numout,*) ' calculated from ECMWF wave fields'
1925	END SELECT
1926	WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm
1927	WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix
1928	WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
1929	WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri
1930	WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix
1931	WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv
1932	ENDIF
1933
1934	! ! allocate zdfosm arrays
1935	IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
1936
1937
1938	IF( ln_osm_mle ) THEN
1939	! Initialise Fox-Kemper parametrization
1940	! Namelist namosm_mle in reference namelist : Tracer advection scheme
1941	READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
1942	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist')
1943
1944	! Namelist namosm_mle in configuration namelist : Tracer advection scheme
1945	READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
1946	904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist')
1947	IF(lwm) WRITE ( numond, namosm_mle )
1948
1949	IF(lwp) THEN ! Namelist print
1950	WRITE(numout,*)
1951	WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
1952	WRITE(numout,*) '~~~~~~~~~~~~~'
1953	WRITE(numout,*) ' Namelist namosm_mle : '
1954	WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle
1955	WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce
1956	WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm'
1957	WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's'
1958	WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg'
1959	WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c
1960	WRITE(numout,*) ' Threshold used to define ML for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s'
1961	WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's'
1962	ENDIF !
1963	ENDIF
1964	!
1965	IF(lwp) THEN
1966	WRITE(numout,*)
1967	IF( ln_osm_mle ) THEN
1968	WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
1969	IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation'
1970	IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation'
1971	ELSE
1972	WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
1973	ENDIF
1974	ENDIF
1975	!
1976	IF( ln_osm_mle ) THEN ! MLE initialisation
1977	!
1978	rb_c = grav * rn_osm_mle_rho_c /rho0 ! Mixed Layer buoyancy criteria
1979	IF(lwp) WRITE(numout,*)
1980	IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 '
1981	IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
1982	!
1983	IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation !
1984	!
1985	ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation
1986	rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) )
1987	!
1988	ENDIF
1989	! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
1990	z1_t2 = 2.e-5
1991	do jj=1,jpj
1992	do ji = 1,jpi
1993	r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
1994	end do
1995	end do
1996	! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
1997	! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 )
1998	!
1999	ENDIF
2000
2001	call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl, dh, hmle
2002
2003
2004	IF( ln_zdfddm) THEN
2005	IF(lwp) THEN
2006	WRITE(numout,*)
2007	WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
2008	WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm '
2009	ENDIF
2010	ENDIF
2011
2012
2013	!set constants not in namelist
2014	!-----------------------------
2015
2016	IF(lwp) THEN
2017	WRITE(numout,*)
2018	ENDIF
2019
2020	IF (nn_osm_wave == 0) THEN
2021	dstokes(:,:) = rn_osm_dstokes
2022	END IF
2023
2024	! Horizontal average : initialization of weighting arrays
2025	! -------------------
2026
2027	SELECT CASE ( nn_ave )
2028
2029	CASE ( 0 ) ! no horizontal average
2030	IF(lwp) WRITE(numout,*) ' no horizontal average on avt'
2031	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
2032	! weighting mean arrays etmean
2033	! ( 1 1 )
2034	! avt = 1/4 ( 1 1 )
2035	!
2036	etmean(:,:,:) = 0.e0
2037
2038	DO_3D_00_00( 1, jpkm1 )
2039	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
2040	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
2041	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
2042	END_3D
2043
2044	CASE ( 1 ) ! horizontal average
2045	IF(lwp) WRITE(numout,*) ' horizontal average on avt'
2046	! weighting mean arrays etmean
2047	! ( 1/2 1 1/2 )
2048	! avt = 1/8 ( 1 2 1 )
2049	! ( 1/2 1 1/2 )
2050	etmean(:,:,:) = 0.e0
2051
2052	DO_3D_00_00( 1, jpkm1 )
2053	etmean(ji,jj,jk) = tmask(ji, jj,jk) &
2054	& / MAX( 1., 2.* tmask(ji,jj,jk) &
2055	& +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) &
2056	& +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
2057	& +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) &
2058	& +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) )
2059	END_3D
2060
2061	CASE DEFAULT
2062	WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave
2063	CALL ctl_stop( ctmp1 )
2064
2065	END SELECT
2066
2067	! Initialization of vertical eddy coef. to the background value
2068	! -------------------------------------------------------------
2069	DO jk = 1, jpk
2070	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
2071	END DO
2072
2073	! zero the surface flux for non local term and osm mixed layer depth
2074	! ------------------------------------------------------------------
2075	ghamt(:,:,:) = 0.
2076	ghams(:,:,:) = 0.
2077	ghamu(:,:,:) = 0.
2078	ghamv(:,:,:) = 0.
2079	!
2080	IF( lwxios ) THEN
2081	CALL iom_set_rstw_var_active('wn')
2082	CALL iom_set_rstw_var_active('hbl')
2083	CALL iom_set_rstw_var_active('dh')
2084	IF( ln_osm_mle ) THEN
2085	CALL iom_set_rstw_var_active('hmle')
2086	END IF
2087	ENDIF
2088	END SUBROUTINE zdf_osm_init
2089
2090
2091	SUBROUTINE osm_rst( kt, Kmm, cdrw )
2092	!!---------------------------------------------------------------------
2093	!! * ROUTINE osm_rst *
2094	!!
2095	!! ** Purpose : Read or write BL fields in restart file
2096	!!
2097	!! ** Method : use of IOM library. If the restart does not contain
2098	!! required fields, they are recomputed from stratification
2099	!!----------------------------------------------------------------------
2100
2101	INTEGER , INTENT(in) :: kt ! ocean time step index
2102	INTEGER , INTENT(in) :: Kmm ! ocean time level index (middle)
2103	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
2104
2105	INTEGER :: id1, id2, id3 ! iom enquiry index
2106	INTEGER :: ji, jj, jk ! dummy loop indices
2107	INTEGER :: iiki, ikt ! local integer
2108	REAL(wp) :: zhbf ! tempory scalars
2109	REAL(wp) :: zN2_c ! local scalar
2110	REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth
2111	INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
2112	!!----------------------------------------------------------------------
2113	!
2114	!!-----------------------------------------------------------------------------
2115	! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
2116	!!-----------------------------------------------------------------------------
2117	IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
2118	id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. )
2119	IF( id1 > 0 ) THEN ! 'wn' exists; read
2120	CALL iom_get( numror, jpdom_autoglo, 'wn', ww, ldxios = lrxios )
2121	WRITE(numout,*) ' ===>>>> : ww read from restart file'
2122	ELSE
2123	ww(:,:,:) = 0._wp
2124	WRITE(numout,*) ' ===>>>> : ww not in restart file, set to zero initially'
2125	END IF
2126
2127	id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. )
2128	id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. )
2129	IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return
2130	CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
2131	CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios )
2132	WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file'
2133	IF( ln_osm_mle ) THEN
2134	id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. )
2135	IF( id3 > 0) THEN
2136	CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios )
2137	WRITE(numout,*) ' ===>>>> : hmle read from restart file'
2138	ELSE
2139	WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl'
2140	hmle(:,:) = hbl(:,:) ! Initialise MLE depth.
2141	END IF
2142	END IF
2143	RETURN
2144	ELSE ! 'hbl' & 'dh' not in restart file, recalculate
2145	WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
2146	END IF
2147	END IF
2148
2149	!!-----------------------------------------------------------------------------
2150	! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
2151	!!-----------------------------------------------------------------------------
2152	IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return
2153	IF(lwp) WRITE(numout,*) '---- osm-rst ----'
2154	CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww , ldxios = lwxios )
2155	CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios )
2156	CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios )
2157	IF( ln_osm_mle ) THEN
2158	CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios )
2159	END IF
2160	RETURN
2161	END IF
2162
2163	!!-----------------------------------------------------------------------------
2164	! Getting hbl, no restart file with hbl, so calculate from surface stratification
2165	!!-----------------------------------------------------------------------------
2166	IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
2167	! w-level of the mixing and mixed layers
2168	CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )
2169	CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2, Kmm)
2170	imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point
2171	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
2172	zN2_c = grav * rho_c * r1_rho0 ! convert density criteria into N^2 criteria
2173	!
2174	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
2175	DO_3D_11_11( 1, jpkm1 )
2176	ikt = mbkt(ji,jj)
2177	hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
2178	IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
2179	END_3D
2180	!
2181	DO_2D_11_11
2182	iiki = MAX(4,imld_rst(ji,jj))
2183	hbl (ji,jj) = gdepw(ji,jj,iiki, Kmm ) ! Turbocline depth
2184	dh (ji,jj) = e3t(ji,jj,iiki-1, Kmm ) ! Turbocline depth
2185	END_2D
2186
2187	WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
2188
2189	IF( ln_osm_mle ) THEN
2190	hmle(:,:) = hbl(:,:) ! Initialise MLE depth.
2191	WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
2192	END IF
2193
2194	ww(:,:,:) = 0._wp
2195	WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
2196	END SUBROUTINE osm_rst
2197
2198
2199	SUBROUTINE tra_osm( kt, Kmm, pts, Krhs )
2200	!!----------------------------------------------------------------------
2201	!! * ROUTINE tra_osm *
2202	!!
2203	!! ** Purpose : compute and add to the tracer trend the non-local tracer flux
2204	!!
2205	!! ** Method : ???
2206	!!----------------------------------------------------------------------
2207	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace
2208	!!----------------------------------------------------------------------
2209	INTEGER , INTENT(in) :: kt ! time step index
2210	INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices
2211	REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation
2212	!
2213	INTEGER :: ji, jj, jk
2214	!
2215	IF( kt == nit000 ) THEN
2216	IF(lwp) WRITE(numout,*)
2217	IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
2218	IF(lwp) WRITE(numout,*) '~~~~~~~ '
2219	ENDIF
2220
2221	IF( l_trdtra ) THEN !* Save ta and sa trends
2222	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
2223	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
2224	ENDIF
2225
2226	DO_3D_00_00( 1, jpkm1 )
2227	pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) &
2228	& - ( ghamt(ji,jj,jk ) &
2229	& - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm)
2230	pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) &
2231	& - ( ghams(ji,jj,jk ) &
2232	& - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm)
2233	END_3D
2234
2235	! save the non-local tracer flux trends for diagnostics
2236	IF( l_trdtra ) THEN
2237	ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
2238	ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
2239	CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt )
2240	CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds )
2241	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
2242	ENDIF
2243
2244	IF(sn_cfctl%l_prtctl) THEN
2245	CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, &
2246	& tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
2247	ENDIF
2248	!
2249	END SUBROUTINE tra_osm
2250
2251
2252	SUBROUTINE trc_osm( kt ) ! Dummy routine
2253	!!----------------------------------------------------------------------
2254	!! * ROUTINE trc_osm *
2255	!!
2256	!! ** Purpose : compute and add to the passive tracer trend the non-local
2257	!! passive tracer flux
2258	!!
2259	!!
2260	!! ** Method : ???
2261	!!----------------------------------------------------------------------
2262	!
2263	!!----------------------------------------------------------------------
2264	INTEGER, INTENT(in) :: kt
2265	WRITE(,) 'trc_osm: Not written yet', kt
2266	END SUBROUTINE trc_osm
2267
2268
2269	SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs )
2270	!!----------------------------------------------------------------------
2271	!! * ROUTINE dyn_osm *
2272	!!
2273	!! ** Purpose : compute and add to the velocity trend the non-local flux
2274	!! copied/modified from tra_osm
2275	!!
2276	!! ** Method : ???
2277	!!----------------------------------------------------------------------
2278	INTEGER , INTENT( in ) :: kt ! ocean time step index
2279	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
2280	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
2281	!
2282	INTEGER :: ji, jj, jk ! dummy loop indices
2283	!!----------------------------------------------------------------------
2284	!
2285	IF( kt == nit000 ) THEN
2286	IF(lwp) WRITE(numout,*)
2287	IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
2288	IF(lwp) WRITE(numout,*) '~~~~~~~ '
2289	ENDIF
2290	!code saving tracer trends removed, replace with trdmxl_oce
2291
2292	DO_3D_00_00( 1, jpkm1 )
2293	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) &
2294	& - ( ghamu(ji,jj,jk ) &
2295	& - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm)
2296	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) &
2297	& - ( ghamv(ji,jj,jk ) &
2298	& - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm)
2299	END_3D
2300	!
2301	! code for saving tracer trends removed
2302	!
2303	END SUBROUTINE dyn_osm
2304
2305	!!======================================================================
2306
2307	END MODULE zdfosm

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

Download in other formats: