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

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

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

Context Navigation

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

Last change on this file since 13859 was 13859, checked in by agn, 3 years ago
Added Alan's commit c46eae7f7fd788 + corr to zznd+pyc def as email 1613 23 Nov 20
Property svn:keywords set to `Id`
File size: 162.3 KB

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

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

Download in other formats: