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 @ 13858

Last change on this file since 13858 was 13858, checked in by agn, 3 years ago
first try of Alan's new code
Property svn:keywords set to `Id`
File size: 159.4 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
	1051	IF(ln_dia_osm) THEN
	1052	IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
	1053	IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
	1054	IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
	1055	IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
	1056	IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
	1057	IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
	1058	END IF
[8930]	1059	!
	1060	! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
	1061
[13858]	1062
	1063	! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer.
	1064
[8930]	1065	DO jj = 2, jpjm1
	1066	DO ji = 2, jpim1
[13858]	1067	IF ( .not. lconv(ji,jj) ) THEN
[8930]	1068	DO jk = 2, ibld(ji,jj)
[13858]	1069	znd = ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
[8930]	1070	IF ( znd >= 0.0 ) THEN
	1071	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
	1072	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
	1073	ELSE
	1074	ghamu(ji,jj,jk) = 0._wp
	1075	ghamv(ji,jj,jk) = 0._wp
	1076	ENDIF
	1077	END DO
	1078	ENDIF
	1079	END DO
	1080	END DO
	1081
	1082	! pynocline contributions
	1083	DO jj = 2, jpjm1
	1084	DO ji = 2, jpim1
[13858]	1085	IF ( .not. lconv(ji,jj) ) THEN
	1086	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
[12216]	1087	DO jk= 2, ibld(ji,jj)
	1088	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	1089	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
	1090	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
	1091	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
	1092	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
	1093	END DO
	1094	END IF
[13858]	1095	END IF
[12216]	1096	END DO
[8930]	1097	END DO
[13858]	1098	IF(ln_dia_osm) THEN
	1099	IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
	1100	IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
	1101	END IF
[8930]	1102
	1103	DO jj=2, jpjm1
	1104	DO ji = 2, jpim1
[12216]	1105	ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
	1106	ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
	1107	ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
	1108	ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
[8930]	1109	END DO ! ji loop
	1110	END DO ! jj loop
	1111
[12216]	1112	IF(ln_dia_osm) THEN
	1113	IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
	1114	IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
	1115	IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
	1116	IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
	1117	IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
	1118	END IF
[8930]	1119	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	1120	! Need to put in code for contributions that are applied explicitly to
	1121	! the prognostic variables
	1122	! 1. Entrainment flux
	1123	!
	1124	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	1125
	1126
	1127
	1128	! rotate non-gradient velocity terms back to model reference frame
	1129
	1130	DO jj = 2, jpjm1
	1131	DO ji = 2, jpim1
	1132	DO jk = 2, ibld(ji,jj)
	1133	ztemp = ghamu(ji,jj,jk)
	1134	ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
	1135	ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
	1136	END DO
	1137	END DO
	1138	END DO
	1139
[12312]	1140	IF(ln_dia_osm) THEN
	1141	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
	1142	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
	1143	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
	1144	END IF
	1145
[8930]	1146	! KPP-style Ri# mixing
	1147	IF( ln_kpprimix) THEN
	1148	DO jk = 2, jpkm1 !* Shear production at uw- and vw-points (energy conserving form)
	1149	DO jj = 1, jpjm1
	1150	DO ji = 1, jpim1 ! vector opt.
	1151	z3du(ji,jj,jk) = 0.5 * ( un(ji,jj,jk-1) - un(ji ,jj,jk) ) &
	1152	& * ( ub(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) &
	1153	& / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
	1154	z3dv(ji,jj,jk) = 0.5 * ( vn(ji,jj,jk-1) - vn(ji,jj ,jk) ) &
	1155	& * ( vb(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) &
	1156	& / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
	1157	END DO
	1158	END DO
	1159	END DO
	1160	!
	1161	DO jk = 2, jpkm1
	1162	DO jj = 2, jpjm1
	1163	DO ji = 2, jpim1 ! vector opt.
	1164	! ! shear prod. at w-point weightened by mask
	1165	zesh2 = ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) &
	1166	& + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
	1167	! ! local Richardson number
	1168	zri = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
	1169	zfri = MIN( zri / rn_riinfty , 1.0_wp )
	1170	zfri = ( 1.0_wp - zfri * zfri )
	1171	zrimix(ji,jj,jk) = zfri * zfri * zfri * wmask(ji, jj, jk)
	1172	END DO
	1173	END DO
	1174	END DO
	1175
	1176	DO jj = 2, jpjm1
	1177	DO ji = 2, jpim1
	1178	DO jk = ibld(ji,jj) + 1, jpkm1
	1179	zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
	1180	zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
	1181	END DO
	1182	END DO
	1183	END DO
	1184
	1185	END IF ! ln_kpprimix = .true.
	1186
	1187	! KPP-style set diffusivity large if unstable below BL
	1188	IF( ln_convmix) THEN
	1189	DO jj = 2, jpjm1
	1190	DO ji = 2, jpim1
	1191	DO jk = ibld(ji,jj) + 1, jpkm1
	1192	IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
	1193	END DO
	1194	END DO
	1195	END DO
	1196	END IF ! ln_convmix = .true.
	1197
[12218]	1198
[13398]	1199
[12218]	1200	IF ( ln_osm_mle ) THEN ! set up diffusivity and non-gradient mixing
	1201	DO jj = 2 , jpjm1
	1202	DO ji = 2, jpim1
[13858]	1203	IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer
	1204	! Calculate MLE flux contribution from surface fluxes
[12218]	1205	DO jk = 1, ibld(ji,jj)
	1206	znd = gdepw_n(ji,jj,jk) / MAX(zhbl(ji,jj),epsln)
	1207	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd )
	1208	ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd )
	1209	END DO
	1210	DO jk = 1, mld_prof(ji,jj)
	1211	znd = gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
	1212	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd )
	1213	ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd )
	1214	END DO
	1215	! Viscosity for MLEs
[13858]	1216	DO jk = 1, mld_prof(ji,jj)
	1217	znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
	1218	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]	1219	END DO
[13858]	1220	ELSE
	1221	! Surface transports limited to OSBL.
	1222	! Viscosity for MLEs
	1223	DO jk = 1, mld_prof(ji,jj)
	1224	znd = -gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
	1225	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]	1226	END DO
[13858]	1227	ENDIF
[12218]	1228	END DO
	1229	END DO
	1230	ENDIF
	1231
	1232	IF(ln_dia_osm) THEN
	1233	IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
	1234	IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
	1235	IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
	1236	END IF
	1237
	1238
[8930]	1239	! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
[12216]	1240	!CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
[8930]	1241
	1242	! GN 25/8: need to change tmask --> wmask
	1243
	1244	DO jk = 2, jpkm1
	1245	DO jj = 2, jpjm1
	1246	DO ji = 2, jpim1
	1247	p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
	1248	p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
	1249	END DO
	1250	END DO
	1251	END DO
	1252	! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and v grids
[10425]	1253	CALL lbc_lnk_multi( 'zdfosm', p_avt, 'W', 1. , p_avm, 'W', 1., &
[9104]	1254	& ghamu, 'W', 1. , ghamv, 'W', 1. )
[8930]	1255	DO jk = 2, jpkm1
	1256	DO jj = 2, jpjm1
	1257	DO ji = 2, jpim1
	1258	ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
	1259	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
	1260
	1261	ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
	1262	& / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
	1263
	1264	ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk)
	1265	ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk)
	1266	END DO
	1267	END DO
	1268	END DO
[12319]	1269	! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)
	1270	CALL lbc_lnk_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
[8930]	1271	! Lateral boundary conditions on final outputs for gham[ts], on W-grid (sign unchanged)
[9104]	1272	! Lateral boundary conditions on final outputs for gham[uv], on [UV]-grid (sign unchanged)
[10425]	1273	CALL lbc_lnk_multi( 'zdfosm', ghamt, 'W', 1. , ghams, 'W', 1., &
[12216]	1274	& ghamu, 'U', -1. , ghamv, 'V', -1. )
[8930]	1275
[12216]	1276	IF(ln_dia_osm) THEN
[8930]	1277	SELECT CASE (nn_osm_wave)
	1278	! Stokes drift set by assumimg onstant La#=0.3(=0) or Pierson-Moskovitz spectrum (=1).
	1279	CASE(0:1)
	1280	IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)zustkezcos_wind ) ! x surface Stokes drift
	1281	IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)zustkezsin_wind ) ! y surface Stokes drift
	1282	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2zustke )
	1283	! Stokes drift read in from sbcwave (=2).
[13400]	1284	CASE(2:3)
[12216]	1285	IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) ) ! x surface Stokes drift
	1286	IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) ) ! y surface Stokes drift
	1287	IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) ) ! wave mean period
	1288	IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) ) ! significant wave height
	1289	IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", (2.rpi1.026/(0.877grav) )wndm*tmask(:,:,1) ) ! wave mean period from NP spectrum
	1290	IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)wndm2tmask(:,:,1) ) ! significant wave height from NP spectrum
	1291	IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) ) ! U_10
[8930]	1292	IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.rau0tmask(:,:,1)zustar2 &
	1293	& SQRT(ut0sd2 + vt0sd2 ) )
	1294	END SELECT
	1295	IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt ) ! <Tw_NL>
	1296	IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams ) ! <Sw_NL>
	1297	IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu ) ! <uw_NL>
	1298	IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv ) ! <vw_NL>
	1299	IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 ) ! <Tw_0>
	1300	IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 ) ! <Sw_0>
	1301	IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl ) ! boundary-layer depth
[12216]	1302	IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld ) ! boundary-layer max k
	1303	IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl ) ! dt at ml base
	1304	IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl ) ! ds at ml base
	1305	IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl ) ! db at ml base
	1306	IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl ) ! du at ml base
	1307	IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl ) ! dv at ml base
	1308	IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh ) ! Initial boundary-layer depth
	1309	IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml ) ! Initial boundary-layer depth
[8930]	1310	IF ( iom_use("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes ) ! Stokes drift penetration depth
	1311	IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke ) ! Stokes drift magnitude at T-points
	1312	IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc ) ! convective velocity scale
	1313	IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl ) ! Langmuir velocity scale
	1314	IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar ) ! friction velocity scale
[12216]	1315	IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr ) ! mixed velocity scale
	1316	IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla ) ! langmuir #
[8930]	1317	IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.rau0tmask(:,:,1)zustar*3 ) ! BL depth internal to zdf_osm routine
	1318	IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.rau0tmask(:,:,1)zustar2zustke )
	1319	IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl ) ! BL depth internal to zdf_osm routine
	1320	IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml ) ! ML depth internal to zdf_osm routine
[12216]	1321	IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld ) ! index for ML depth internal to zdf_osm routine
	1322	IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh ) ! pyc thicknessh internal to zdf_osm routine
[8930]	1323	IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol ) ! ML depth internal to zdf_osm routine
[12218]	1324	IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav ) ! upward BL-avged turb temp flux
	1325	IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent ) ! upward turb temp entrainment flux
	1326	IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent ) ! upward turb buoyancy entrainment flux
	1327	IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent ) ! upward turb salinity entrainment flux
	1328	IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml ) ! average T in ML
	1329
	1330	IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle ) ! FK layer depth
	1331	IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld ) ! FK target layer depth
	1332	IF ( iom_use("zwb_fk") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk ) ! FK b flux
	1333	IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b ) ! FK b flux averaged over ML
	1334	IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )! FK layer max k
	1335	IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx ) ! FK dtdx at u-pt
	1336	IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy ) ! FK dtdy at v-pt
	1337	IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx ) ! FK dtdx at u-pt
	1338	IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy ) ! FK dsdy at v-pt
	1339	IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle ) ! FK dbdx at u-pt
	1340	IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle ) ! FK dbdy at v-pt
	1341	IF ( iom_use("zdiff_mle") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
	1342	IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
	1343
[8946]	1344	END IF
[12218]	1345
	1346	CONTAINS
[13858]	1347	! subroutine code changed, needs syntax checking.
	1348	SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
[12218]	1349
[13858]	1350	!!---------------------------------------------------------------------
	1351	!! * ROUTINE zdf_osm_diffusivity_viscosity *
	1352	!!
	1353	!! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline.
	1354	!!
	1355	!! ** Method :
	1356	!!
	1357	!! !!----------------------------------------------------------------------
	1358	REAL(wp), DIMENSION(:,:,:) :: zdiffut
	1359	REAL(wp), DIMENSION(:,:,:) :: zviscos
	1360	! local
[12218]	1361
[13858]	1362	! Scales used to calculate eddy diffusivity and viscosity profiles
	1363	REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc
	1364	REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr
	1365	REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr
	1366	REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc
	1367	!
	1368	REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac
	1369
	1370	REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375
	1371	REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142
	1372	REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15
	1373
	1374	DO jj = 2, jpjm1
	1375	DO ji = 2, jpim1
	1376	IF ( lconv(ji,jj) ) THEN
	1377
	1378	zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)3 + zwstrc(ji,jj)3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird
	1379	zvel_sc_ml = ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
	1380	zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )1.25 ) )2
	1381
	1382	zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml
	1383	zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj)
	1384
	1385	IF ( lpyc(ji,jj) ) THEN
	1386	zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj)
	1387
	1388	IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN
	1389	zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )*pthird zhbl(ji,jj)
	1390	ENDIF
	1391
	1392	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) )
	1393	zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac
	1394	zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) )
	1395
	1396	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)
	1397	zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac
	1398	IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN
	1399	zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )*pthird zhbl(ji,jj)
	1400	ENDIF
	1401
	1402	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) ) )
	1403	zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac
	1404	zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) )
	1405
	1406	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
	1407	zbeta_v_sc(ji,jj) = 1.0 - 2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln )
	1408	ELSE
	1409	zbeta_d_sc = 1.0
	1410	zbeta_v_sc = 1.0
	1411	ENDIF
	1412	ELSE
	1413	zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
	1414	zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
	1415	END IF
	1416	END DO
	1417	END DO
	1418	!
	1419	DO jj = 2, jpjm1
	1420	DO ji = 2, jpim1
	1421	IF ( lconv(ji,jj) ) THEN
	1422	DO jk = 2, imld(ji,jj) ! mixed layer diffusivity
	1423	zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
	1424	!
	1425	zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
	1426	!
	1427	zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) &
	1428	& * ( 1.0 - 0.5 * zznd_ml**2 )
	1429	END DO
	1430	! pycnocline
	1431	IF ( lpyc(ji,jj) ) THEN
	1432	! Diffusivity profile in the pycnocline given by cubic polynomial.
	1433	za_cubic = 0.5
	1434	zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj)
	1435	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 ) &
	1436	& - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8)
	1437	zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic )
	1438	zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
	1439	DO jk = imld(ji,jj) , ibld(ji,jj)
	1440	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
	1441	!
	1442	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 )
	1443
	1444	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 )
	1445	END DO
	1446	! viscosity profiles.
	1447	za_cubic = 0.5
	1448	zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)
	1449	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)
	1450	zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic )
	1451	zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
	1452	DO jk = imld(ji,jj) , ibld(ji,jj)
	1453	zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
	1454	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 )
	1455	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 )
	1456	END DO
	1457	IF ( zdhdt(ji,jj) > 0._wp ) THEN
	1458	zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 )
	1459	zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w_n(ji,jj,ibld(ji,jj)+1), 1.0e-6 )
	1460	ELSE
	1461	zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
	1462	zviscos(ji,jj,ibld(ji,jj)) = 0._wp
	1463	ENDIF
	1464	ENDIF
	1465	ELSE
	1466	! stable conditions
	1467	DO jk = 2, ibld(ji,jj)
	1468	zznd_ml = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	1469	zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
	1470	zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
	1471	END DO
	1472
	1473	IF ( zdhdt(ji,jj) > 0._wp ) THEN
	1474	zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj))
	1475	zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w_n(ji, jj, ibld(ji,jj))
	1476	ENDIF
	1477	ENDIF ! end if ( lconv )
	1478	!
	1479	END DO ! end of ji loop
	1480	END DO ! end of jj loop
	1481
	1482	END SUBROUTINE zdf_osm_diffusivity_viscosity
	1483
	1484	SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i )
	1485
	1486	!!---------------------------------------------------------------------
	1487	!! * ROUTINE zdf_osm_osbl_state *
	1488	!!
	1489	!! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number
	1490	!!
	1491	!! ** Method :
	1492	!!
	1493	!! !!----------------------------------------------------------------------
	1494
	1495	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.
	1496
	1497	LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear
	1498
	1499	REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer.
	1500	REAL(wp), DIMENSION(jpi,jpj) :: zshear ! production of TKE due to shear across the pycnocline
	1501	REAL(wp), DIMENSION(jpi,jpj) :: zri_i ! Interfacial Richardson Number
	1502
	1503	! Local Variables
	1504
	1505	INTEGER :: jj, ji
	1506
	1507	REAL(wp), DIMENSION(jpi,jpj) :: zekman
	1508	REAL(wp) :: zri_p, zri_b ! Richardson numbers
	1509	REAL(wp) :: zshear_u, zshear_v, zwb_shr
	1510	REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
	1511
	1512	REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1
	1513	REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59
	1514	REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04
	1515	REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15
	1516	REAL, PARAMETER :: rn_ri_p_thresh = 27.0
	1517	REAL, PARAMETER :: zrot=0._wp ! dummy rotation rate of surface stress.
	1518
	1519	! Determins stability and set flag lconv
	1520	DO jj = 2, jpjm1
	1521	DO ji = 2, jpim1
	1522	IF ( zhol(ji,jj) < 0._wp ) THEN
	1523	lconv(ji,jj) = .TRUE.
	1524	ELSE
	1525	lconv(ji,jj) = .FALSE.
	1526	ENDIF
	1527	END DO
	1528	END DO
	1529
	1530	zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) )
	1531
	1532	WHERE ( lconv )
	1533	zri_i = zdb_ml * zhml2 / MAX( ( zvstr3 + 0.5 * zwstrc3 )p2third * zdh, 1.e-12 )
	1534	END WHERE
	1535
	1536	zshear(:,:) = 0._wp
	1537	j_ddh(:,:) = 1
	1538
	1539	DO jj = 2, jpjm1
	1540	DO ji = 2, jpim1
	1541	IF ( lconv(ji,jj) ) THEN
	1542	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
	1543	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 &
	1544	& / MAX( zekman(ji,jj), 1.e-6 ) , 5._wp )
	1545
	1546	zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)2 + zdv_ml(ji,jj)2, 1.e-8 )
	1547
	1548	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 ) )
	1549	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1550	! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when !
	1551	! full code available !
	1552	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1553	IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN
	1554	! Growing shear layer
	1555	j_ddh(ji,jj) = 0
	1556	lshear(ji,jj) = .TRUE.
	1557	ELSE
	1558	j_ddh(ji,jj) = 1
	1559	IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN
	1560	! shear production large enough to determine layer charcteristics, but can't maintain a shear layer.
	1561	lshear(ji,jj) = .TRUE.
	1562	ELSE
	1563	! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline.
	1564	zshear(ji,jj) = 0.5 * zshear(ji,jj)
	1565	lshear(ji,jj) = .FALSE.
	1566	ENDIF
	1567	ENDIF
	1568	ELSE ! zdb_bl test, note zshear set to zero
	1569	j_ddh(ji,jj) = 2
	1570	lshear(ji,jj) = .FALSE.
	1571	ENDIF
	1572	ENDIF
	1573	END DO
	1574	END DO
	1575
	1576	! Calculate entrainment buoyancy flux due to surface fluxes.
	1577
	1578	DO jj = 2, jpjm1
	1579	DO ji = 2, jpim1
	1580	IF ( lconv(ji,jj) ) THEN
	1581	zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
	1582	zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
	1583	zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
	1584	zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
	1585	IF (nn_osm_SD_reduce > 0 ) THEN
	1586	! Effective Stokes drift already reduced from surface value
	1587	zr_stokes = 1.0_wp
	1588	ELSE
	1589	! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
	1590	! requires further reduction where BL is deep
	1591	zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
	1592	& * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
	1593	END IF
	1594	zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
	1595	& - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
	1596	& + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
	1597	& - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
	1598	!
	1599	ENDIF
	1600	END DO ! ji loop
	1601	END DO ! jj loop
	1602
	1603	zwb_min(:,:) = 0._wp
	1604
	1605	DO jj = 2, jpjm1
	1606	DO ji = 2, jpjm1
	1607	IF ( lshear(ji,jj) ) THEN
	1608	IF ( lconv(ji,jj) ) THEN
	1609	! Unstable OSBL
	1610	zwb_shr = -za_wb_s * zshear(ji,jj)
	1611	IF ( j_ddh(ji,jj) == 0 ) THEN
	1612
	1613	! Developing shear layer, additional shear production possible.
	1614
	1615	zshear_u = MAX( zustar(ji,jj)*2 zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp )
	1616	zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) )
	1617	zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )
	1618
	1619	zwb_shr = -za_wb_s * zshear(ji,jj)
	1620
	1621	ENDIF
	1622	zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
	1623	zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
	1624	ELSE ! IF ( lconv ) THEN - ENDIF
	1625	! Stable OSBL - shear production not coded for first attempt.
	1626	ENDIF ! lconv
	1627	ELSE ! lshear
	1628	IF ( lconv(ji,jj) ) THEN
	1629	! Unstable OSBL
	1630	zwb_shr = -za_wb_s * zshear(ji,jj)
	1631	zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
	1632	zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
	1633	ENDIF ! lconv
	1634	ENDIF ! lshear
	1635	END DO ! ji
	1636	END DO ! jj
	1637	END SUBROUTINE zdf_osm_osbl_state
	1638
	1639
	1640	SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv )
[12218]	1641	!!---------------------------------------------------------------------
	1642	!! * ROUTINE zdf_vertical_average *
	1643	!!
	1644	!! ** Purpose : Determines vertical averages from surface to jnlev.
	1645	!!
	1646	!! ** Method : Averages are calculated from the surface to jnlev.
	1647	!! The external level used to calculate differences is ibld+ibld_ext
	1648	!!
	1649	!!----------------------------------------------------------------------
	1650
	1651	INTEGER, DIMENSION(jpi,jpj) :: jnlev_av ! Number of levels to average over.
[13858]	1652	INTEGER, DIMENSION(jpi,jpj) :: jp_ext
[12218]	1653
	1654	! Alan: do we need zb?
[13858]	1655	REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb ! Average temperature and salinity
[12218]	1656	REAL(wp), DIMENSION(jpi,jpj) :: zu,zv ! Average current components
	1657	REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
	1658	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Difference for velocity components.
	1659
[13858]	1660	INTEGER :: jk, ji, jj, ibld_ext
[12218]	1661	REAL(wp) :: zthick, zthermal, zbeta
	1662
	1663
	1664	zt = 0._wp
	1665	zs = 0._wp
	1666	zu = 0._wp
	1667	zv = 0._wp
	1668	DO jj = 2, jpjm1 ! Vertical slab
	1669	DO ji = 2, jpim1
	1670	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
	1671	zbeta = rab_n(ji,jj,1,jp_sal)
	1672	! average over depth of boundary layer
	1673	zthick = epsln
	1674	DO jk = 2, jnlev_av(ji,jj)
	1675	zthick = zthick + e3t_n(ji,jj,jk)
	1676	zt(ji,jj) = zt(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
	1677	zs(ji,jj) = zs(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
	1678	zu(ji,jj) = zu(ji,jj) + e3t_n(ji,jj,jk) &
	1679	& * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
	1680	& / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
	1681	zv(ji,jj) = zv(ji,jj) + e3t_n(ji,jj,jk) &
	1682	& * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
	1683	& / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
	1684	END DO
	1685	zt(ji,jj) = zt(ji,jj) / zthick
	1686	zs(ji,jj) = zs(ji,jj) / zthick
	1687	zu(ji,jj) = zu(ji,jj) / zthick
	1688	zv(ji,jj) = zv(ji,jj) / zthick
[13858]	1689	zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj)
	1690	ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj)
	1691	IF ( ibld_ext < mbkt(ji,jj) ) THEN
	1692	zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld_ext,jp_tem)
	1693	zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld_ext,jp_sal)
	1694	zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld_ext) + ub(ji-1,jj,ibld_ext ) ) &
	1695	& / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) )
	1696	zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld_ext) + vb(ji,jj-1,ibld_ext ) ) &
	1697	& / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) )
	1698	zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
	1699	ELSE
	1700	zdt(ji,jj) = 0._wp
	1701	zds(ji,jj) = 0._wp
	1702	zdu(ji,jj) = 0._wp
	1703	zdv(ji,jj) = 0._wp
	1704	zdb(ji,jj) = 0._wp
	1705	ENDIF
[12218]	1706	END DO
	1707	END DO
	1708	END SUBROUTINE zdf_osm_vertical_average
	1709
	1710	SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
	1711	!!---------------------------------------------------------------------
	1712	!! * ROUTINE zdf_velocity_rotation *
	1713	!!
	1714	!! ** Purpose : Rotates frame of reference of averaged velocity components.
	1715	!!
	1716	!! ** Method : The velocity components are rotated into frame specified by zcos_w and zsin_w
	1717	!!
	1718	!!----------------------------------------------------------------------
	1719
	1720	REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w ! Cos and Sin of rotation angle
	1721	REAL(wp), DIMENSION(jpi,jpj) :: zu, zv ! Components of current
	1722	REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv ! Change in velocity components across pycnocline
	1723
	1724	INTEGER :: ji, jj
	1725	REAL(wp) :: ztemp
	1726
	1727	DO jj = 2, jpjm1
	1728	DO ji = 2, jpim1
	1729	ztemp = zu(ji,jj)
	1730	zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
	1731	zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
	1732	ztemp = zdu(ji,jj)
	1733	zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
	1734	zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj)
	1735	END DO
	1736	END DO
	1737	END SUBROUTINE zdf_osm_velocity_rotation
	1738
[13858]	1739	SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
[12218]	1740	!!---------------------------------------------------------------------
[13858]	1741	!! * ROUTINE zdf_osm_osbl_state_fk *
	1742	!!
	1743	!! ** 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.
	1744	!! lpyc :: determines whether pycnocline flux-grad relationship needs to be determined
	1745	!! lflux :: determines whether effects of surface flux extend below the base of the OSBL
	1746	!! lmle :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl.
	1747	!!
	1748	!! ** Method :
	1749	!!
	1750	!!
	1751	!!----------------------------------------------------------------------
	1752
	1753	! Outputs
	1754	LOGICAL, DIMENSION(jpi,jpj) :: lpyc, lflux, lmle
	1755	REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk
	1756	!
	1757	REAL(wp), DIMENSION(jpi,jpj) :: znd_param
	1758	REAL(wp) :: zbuoy, ztmp, zpe_mle_layer
	1759	REAL(wp) :: zpe_mle_ref, zwb_ent, zdbdz_mle_int
	1760
	1761	znd_param(:,:) = 0._wp
	1762
	1763	DO jj = 2, jpjm1
	1764	DO ji = 2, jpim1
	1765	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
	1766	zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj)
	1767	END DO
	1768	END DO
	1769	DO jj = 2, jpjm1
	1770	DO ji = 2, jpim1
	1771	!
	1772	IF ( lconv(ji,jj) ) THEN
	1773	IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
	1774	zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
	1775	zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
	1776	zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
	1777	zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
	1778	! Calculate potential energies of actual profile and reference profile.
	1779	zpe_mle_layer = 0._wp
	1780	zpe_mle_ref = 0._wp
	1781	DO jk = ibld(ji,jj), mld_prof(ji,jj)
	1782	zbuoy = grav * ( zthermal * tsn(ji,jj,jk,jp_tem) - zbeta * tsn(ji,jj,jk,jp_sal) )
	1783	zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw_n(ji,jj,jk) * e3w_n(ji,jj,jk)
	1784	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)
	1785	END DO
	1786	! Non-dimensional parameter to diagnose the presence of thermocline
	1787
	1788	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) )
	1789	ENDIF
	1790	ENDIF
	1791	END DO
	1792	END DO
	1793
	1794	! Diagnosis
	1795	DO jj = 2, jpjm1
	1796	DO ji = 2, jpim1
	1797	IF ( lconv(ji,jj) ) THEN
	1798	zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) &
	1799	& - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) &
	1800	& + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 &
	1801	& - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
	1802	IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN
	1803	IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
	1804	! MLE layer growing
	1805	IF ( znd_param (ji,jj) > 100. ) THEN
	1806	! Thermocline present
	1807	lflux(ji,jj) = .FALSE.
	1808	lmle(ji,jj) =.FALSE.
	1809	ELSE
	1810	! Thermocline not present
	1811	lflux(ji,jj) = .TRUE.
	1812	lmle(ji,jj) = .TRUE.
	1813	ENDIF ! znd_param > 100
	1814	!
	1815	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
	1816	lpyc(ji,jj) = .FALSE.
	1817	ELSE
	1818	lpyc = .TRUE.
	1819	ENDIF
	1820	ELSE
	1821	! MLE layer restricted to OSBL or just below.
	1822	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
	1823	! Weak stratification MLE layer can grow.
	1824	lpyc(ji,jj) = .FALSE.
	1825	lflux(ji,jj) = .TRUE.
	1826	lmle(ji,jj) = .TRUE.
	1827	ELSE
	1828	! Strong stratification
	1829	lpyc(ji,jj) = .TRUE.
	1830	lflux(ji,jj) = .FALSE.
	1831	lmle(ji,jj) = .FALSE.
	1832	ENDIF ! zdb_bl < rn_mle_thresh_bl and
	1833	ENDIF ! zhmle > 1.2 zhbl
	1834	ELSE
	1835	lpyc(ji,jj) = .TRUE.
	1836	lflux(ji,jj) = .FALSE.
	1837	lmle(ji,jj) = .FALSE.
	1838	IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
	1839	ENDIF ! -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5
	1840	ELSE
	1841	! Stable Boundary Layer
	1842	lpyc(ji,jj) = .FALSE.
	1843	lflux(ji,jj) = .FALSE.
	1844	lmle(ji,jj) = .FALSE.
	1845	ENDIF ! lconv
	1846	END DO
	1847	END DO
	1848	END SUBROUTINE zdf_osm_osbl_state_fk
	1849
	1850	SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz )
	1851	!!---------------------------------------------------------------------
[12218]	1852	!! * ROUTINE zdf_osm_external_gradients *
	1853	!!
	1854	!! ** Purpose : Calculates the gradients below the OSBL
	1855	!!
	1856	!! ** Method : Uses ibld and ibld_ext to determine levels to calculate the gradient.
	1857	!!
	1858	!!----------------------------------------------------------------------
	1859
[13858]	1860	INTEGER, DIMENSION(jpi,jpj) :: jbase
[12218]	1861	REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz ! External gradients of temperature, salinity and buoyancy.
	1862
	1863	INTEGER :: jj, ji, jkb, jkb1
	1864	REAL(wp) :: zthermal, zbeta
	1865
	1866
	1867	DO jj = 2, jpjm1
	1868	DO ji = 2, jpim1
[13858]	1869	IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
[12218]	1870	zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
	1871	zbeta = rab_n(ji,jj,1,jp_sal)
[13858]	1872	jkb = jbase(ji,jj)
[12218]	1873	jkb1 = MIN(jkb + 1, mbkt(ji,jj))
	1874	zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
	1875	& / e3t_n(ji,jj,ibld(ji,jj))
	1876	zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
	1877	& / e3t_n(ji,jj,ibld(ji,jj))
	1878	zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
[13858]	1879	ELSE
	1880	zdtdz(ji,jj) = 0._wp
	1881	zdsdz(ji,jj) = 0._wp
	1882	zdbdz(ji,jj) = 0._wp
[12218]	1883	END IF
	1884	END DO
	1885	END DO
	1886	END SUBROUTINE zdf_osm_external_gradients
	1887
[13858]	1888	SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha )
[12218]	1889
	1890	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz ! gradients in the pycnocline
[13858]	1891	REAL(wp), DIMENSION(jpi,jpj) :: zalpha
[12218]	1892
	1893	INTEGER :: jk, jj, ji
	1894	REAL(wp) :: ztgrad, zsgrad, zbgrad
[13858]	1895	REAL(wp) :: zgamma_b_nd, znd
	1896	REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc
	1897	REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
[12218]	1898
	1899	DO jj = 2, jpjm1
	1900	DO ji = 2, jpim1
[13858]	1901	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
[12218]	1902	IF ( lconv(ji,jj) ) THEN ! convective conditions
[13858]	1903	IF ( lpyc(ji,jj) ) THEN
	1904	zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
	1905	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 ) )
	1906	zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp )
	1907
	1908	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
	1909	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1910	! Commented lines in this section are not needed in new code, once tested !
	1911	! can be removed !
	1912	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	1913	! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
	1914	! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
	1915	zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)
	1916	zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
	1917	DO jk = 2, ibld(ji,jj)+ibld_ext
	1918	znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) * ztmp
	1919	IF ( znd <= zzeta_m ) THEN
	1920	! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &
	1921	! & EXP( -6.0 * ( znd -zzeta_m )**2 )
	1922	! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &
	1923	! & EXP( -6.0 * ( znd -zzeta_m )**2 )
	1924	zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &
	1925	& EXP( -6.0 * ( znd -zzeta_m )**2 )
	1926	ELSE
	1927	! zdtdz(ji,jj,jk) = ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
	1928	! zdsdz(ji,jj,jk) = zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
	1929	zdbdz(ji,jj,jk) = zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
	1930	ENDIF
	1931	END DO
	1932	ENDIF ! if no pycnocline pycnocline gradients set to zero
[12218]	1933	ELSE
	1934	! stable conditions
	1935	! if pycnocline profile only defined when depth steady of increasing.
[13858]	1936	IF ( zdhdt(ji,jj) > 0.0 ) THEN ! Depth increasing, or steady.
[12218]	1937	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
	1938	IF ( zhol(ji,jj) >= 0.5 ) THEN ! Very stable - 'thick' pycnocline
[12384]	1939	ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
	1940	ztgrad = zdt_bl(ji,jj) * ztmp
	1941	zsgrad = zds_bl(ji,jj) * ztmp
	1942	zbgrad = zdb_bl(ji,jj) * ztmp
[12218]	1943	DO jk = 2, ibld(ji,jj)
[12384]	1944	znd = gdepw_n(ji,jj,jk) * ztmp
[12218]	1945	zdtdz(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
	1946	zdbdz(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
	1947	zdsdz(ji,jj,jk) = zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
	1948	END DO
	1949	ELSE ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
[12384]	1950	ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
	1951	ztgrad = zdt_bl(ji,jj) * ztmp
	1952	zsgrad = zds_bl(ji,jj) * ztmp
	1953	zbgrad = zdb_bl(ji,jj) * ztmp
[12218]	1954	DO jk = 2, ibld(ji,jj)
[12384]	1955	znd = -( gdepw_n(ji,jj,jk) - zhml(ji,jj) ) * ztmp
[12218]	1956	zdtdz(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	1957	zdbdz(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	1958	zdsdz(ji,jj,jk) = zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
	1959	END DO
	1960	ENDIF ! IF (zhol >=0.5)
	1961	ENDIF ! IF (zdb_bl> 0.)
	1962	ENDIF ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
	1963	ENDIF ! IF (lconv)
[13858]	1964	ENDIF ! IF ( ibld(ji,jj) < mbkt(ji,jj) )
[12218]	1965	END DO
	1966	END DO
	1967
	1968	END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
	1969
	1970	SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
	1971	!!---------------------------------------------------------------------
	1972	!! * ROUTINE zdf_osm_pycnocline_shear_profiles *
	1973	!!
	1974	!! ** Purpose : Calculates velocity shear in the pycnocline
	1975	!!
	1976	!! ** Method :
	1977	!!
	1978	!!----------------------------------------------------------------------
	1979
	1980	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
	1981
	1982	INTEGER :: jk, jj, ji
	1983	REAL(wp) :: zugrad, zvgrad, znd
	1984	REAL(wp) :: zzeta_v = 0.45
[8946]	1985	!
[12218]	1986	DO jj = 2, jpjm1
	1987	DO ji = 2, jpim1
	1988	!
[13858]	1989	IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
[12218]	1990	IF ( lconv (ji,jj) ) THEN
[13858]	1991	! Unstable conditions. Shouldn;t be needed with no pycnocline code.
	1992	! zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
	1993	! & ( ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird * zhml(ji,jj) ) * &
	1994	! & MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
[12218]	1995	!Alan is this right?
[13858]	1996	! zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
	1997	! & 2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
	1998	! & ( ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird + epsln ) &
	1999	! & )/ (zdh(ji,jj) + epsln )
	2000	! DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
	2001	! znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
	2002	! IF ( znd <= 0.0 ) THEN
	2003	! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
	2004	! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
	2005	! ELSE
	2006	! zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
	2007	! zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
	2008	! ENDIF
	2009	! END DO
[12218]	2010	ELSE
	2011	! stable conditions
	2012	zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
	2013	zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
	2014	DO jk = 2, ibld(ji,jj)
	2015	znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
	2016	IF ( znd < 1.0 ) THEN
	2017	zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
	2018	ELSE
	2019	zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
	2020	ENDIF
	2021	zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
	2022	END DO
	2023	ENDIF
	2024	!
	2025	END IF ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
	2026	END DO
	2027	END DO
	2028	END SUBROUTINE zdf_osm_pycnocline_shear_profiles
[8930]	2029
[13858]	2030	SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt )
[12218]	2031	!!---------------------------------------------------------------------
	2032	!! * ROUTINE zdf_osm_calculate_dhdt *
	2033	!!
	2034	!! ** Purpose : Calculates the rate at which hbl changes.
	2035	!!
	2036	!! ** Method :
	2037	!!
	2038	!!----------------------------------------------------------------------
[8946]	2039
[13858]	2040	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt ! Rate of change of hbl
[12218]	2041
	2042	INTEGER :: jj, ji
[13858]	2043	REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi
	2044	REAL(wp) :: zvel_max!, zwb_min
[12218]	2045	REAL(wp) :: zzeta_m = 0.3
	2046	REAL(wp) :: zgamma_c = 2.0
	2047	REAL(wp) :: zdhoh = 0.1
	2048	REAL(wp) :: alpha_bc = 0.5
[13858]	2049	REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
	2050
	2051	DO jj = 2, jpjm1
	2052	DO ji = 2, jpim1
	2053
	2054	IF ( lshear(ji,jj) ) THEN
	2055	IF ( lconv(ji,jj) ) THEN ! Convective
[12218]	2056
[13858]	2057	IF ( ln_osm_mle ) THEN
	2058
	2059	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
	2060	! Fox-Kemper buoyancy flux average over OSBL
	2061	zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * &
	2062	(1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
	2063	ELSE
	2064	zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
	2065	ENDIF
	2066	zvel_max = ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
	2067	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
	2068	! OSBL is deepening, entrainment > restratification
	2069	IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
	2070	! *** Used for shear Needs to be changed to work stabily
	2071	! zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml
	2072	! zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd )
	2073	! zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) )
	2074	! za_1 = 1.0 / zgamma_b**2 - 0.017
	2075	! za_2 = 1.0 / zgamma_b**3 - 0.0025
	2076	! zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl )
	2077	zpsi = 0._wp
	2078	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
	2079	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)
	2080	IF ( j_ddh(ji,jj) == 1 ) THEN
	2081	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
	2082	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 )
	2083	ELSE
	2084	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 )
	2085	ENDIF
	2086	! Relaxation to dh_ref = zari * hbl
	2087	zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
	2088
	2089	ELSE ! j_ddh == 0
	2090	! Growing shear layer
	2091	zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
	2092	ENDIF ! j_ddh
	2093	zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj)
	2094	ELSE ! zdb_bl >0
	2095	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)
	2096	ENDIF
	2097	ELSE ! zwb_min + 2*zwb_fk_b < 0
	2098	! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
	2099	zdhdt(ji,jj) = - zvel_mle(ji,jj)
	2100
	2101
	2102	ENDIF
	2103
[13400]	2104	ELSE
[13858]	2105	! Fox-Kemper not used.
[12218]	2106
[13858]	2107	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) / &
	2108	& MAX((zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln)
	2109	zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
	2110	! added ajgn 23 July as temporay fix
[12218]	2111
[13858]	2112	ENDIF ! ln_osm_mle
	2113
	2114	ELSE ! lconv - Stable
	2115	zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
	2116	IF ( zdhdt(ji,jj) < 0._wp ) THEN
	2117	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
	2118	zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
	2119	ELSE
	2120	zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
	2121	ENDIF
	2122	zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
	2123	ENDIF ! lconv
	2124	ELSE ! lshear
	2125	IF ( lconv(ji,jj) ) THEN ! Convective
	2126
[12218]	2127	IF ( ln_osm_mle ) THEN
[13858]	2128
[12218]	2129	IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
[13858]	2130	! Fox-Kemper buoyancy flux average over OSBL
[12218]	2131	zwb_fk_b(ji,jj) = zwb_fk(ji,jj) * &
	2132	(1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
	2133	ELSE
	2134	zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
	2135	ENDIF
	2136	zvel_max = ( zwstrl(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )p2third / hbl(ji,jj)
	2137	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
	2138	! OSBL is deepening, entrainment > restratification
[13858]	2139	IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
[12218]	2140	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
	2141	ELSE
	2142	zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / MAX( zvel_max, 1.0e-15)
	2143	ENDIF
	2144	ELSE
	2145	! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
	2146	zdhdt(ji,jj) = - zvel_mle(ji,jj)
	2147
	2148
	2149	ENDIF
	2150
	2151	ELSE
	2152	! Fox-Kemper not used.
	2153
[13858]	2154	zvel_max = -zwb_ent(ji,jj) / &
[13396]	2155	& MAX((zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird, epsln)
[12218]	2156	zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
	2157	! added ajgn 23 July as temporay fix
	2158
[13858]	2159	ENDIF ! ln_osm_mle
[12218]	2160
	2161	ELSE ! Stable
	2162	zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
	2163	IF ( zdhdt(ji,jj) < 0._wp ) THEN
	2164	! For long timsteps factor in brackets slows the rapid collapse of the OSBL
[13858]	2165	zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj)
[12218]	2166	ELSE
[13858]	2167	zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
[12218]	2168	ENDIF
[13396]	2169	zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
[13858]	2170	ENDIF ! lconv
	2171	ENDIF ! lshear
	2172	END DO
	2173	END DO
[12218]	2174	END SUBROUTINE zdf_osm_calculate_dhdt
	2175
[13858]	2176	SUBROUTINE zdf_osm_timestep_hbl( zdhdt )
[12218]	2177	!!---------------------------------------------------------------------
	2178	!! * ROUTINE zdf_osm_timestep_hbl *
	2179	!!
	2180	!! ** Purpose : Increments hbl.
	2181	!!
	2182	!! ** Method : If thechange in hbl exceeds one model level the change is
	2183	!! is calculated by moving down the grid, changing the buoyancy
	2184	!! jump. This is to ensure that the change in hbl does not
	2185	!! overshoot a stable layer.
	2186	!!
	2187	!!----------------------------------------------------------------------
	2188
	2189
[13858]	2190	REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! rates of change of hbl.
[12218]	2191
	2192	INTEGER :: jk, jj, ji, jm
	2193	REAL(wp) :: zhbl_s, zvel_max, zdb
	2194	REAL(wp) :: zthermal, zbeta
	2195
	2196	DO jj = 2, jpjm1
	2197	DO ji = 2, jpim1
	2198	IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
	2199	!
	2200	! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
	2201	!
	2202	zhbl_s = hbl(ji,jj)
	2203	jm = imld(ji,jj)
	2204	zthermal = rab_n(ji,jj,1,jp_tem)
	2205	zbeta = rab_n(ji,jj,1,jp_sal)
	2206
	2207
	2208	IF ( lconv(ji,jj) ) THEN
	2209	!unstable
	2210
	2211	IF( ln_osm_mle ) THEN
	2212	zvel_max = ( zwstrl(ji,jj)3 + zwstrc(ji,jj)3 )**p2third / hbl(ji,jj)
	2213	ELSE
	2214
	2215	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) / &
	2216	& ( zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3 )pthird
	2217
	2218	ENDIF
	2219
	2220	DO jk = imld(ji,jj), ibld(ji,jj)
	2221	zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) ) &
	2222	& - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
	2223	& 0.0 ) + zvel_max
	2224
	2225
	2226	IF ( ln_osm_mle ) THEN
	2227	zhbl_s = zhbl_s + MIN( &
[13858]	2228	& rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
[12218]	2229	& e3w_n(ji,jj,jm) )
	2230	ELSE
	2231	zhbl_s = zhbl_s + MIN( &
[13858]	2232	& rn_rdt * ( -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
[12218]	2233	& e3w_n(ji,jj,jm) )
	2234	ENDIF
	2235
[13858]	2236	! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
	2237	IF ( zhbl_s >= gdepw_n(ji,jj,mbkt(ji,jj) + 1) ) THEN
	2238	zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
	2239	lpyc(ji,jj) = .FALSE.
	2240	ENDIF
[12218]	2241	IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) jm = jm + 1
	2242	END DO
	2243	hbl(ji,jj) = zhbl_s
	2244	ibld(ji,jj) = jm
	2245	ELSE
	2246	! stable
	2247	DO jk = imld(ji,jj), ibld(ji,jj)
	2248	zdb = MAX( &
	2249	& grav * ( zthermal * ( zt_bl(ji,jj) - tsn(ji,jj,jm,jp_tem) )&
	2250	& - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
	2251	& 0.0 ) + &
	2252	& 2.0 * zvstr(ji,jj)**2 / zhbl_s
	2253
	2254	! Alan is thuis right? I have simply changed hbli to hbl
	2255	zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
	2256	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) ) ) * &
	2257	& zustar(ji,jj)*3 / zhbl_s ) ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
	2258	zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
	2259	zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
	2260
[13858]	2261	! zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
	2262	IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN
	2263	zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
	2264	lpyc(ji,jj) = .FALSE.
	2265	ENDIF
[12218]	2266	IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
	2267	END DO
	2268	ENDIF ! IF ( lconv )
	2269	hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
	2270	ibld(ji,jj) = MAX(jm, 4 )
	2271	ELSE
	2272	! change zero or one model level.
	2273	hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
	2274	ENDIF
	2275	zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
	2276	END DO
	2277	END DO
	2278
	2279	END SUBROUTINE zdf_osm_timestep_hbl
	2280
	2281	SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
	2282	!!---------------------------------------------------------------------
	2283	!! * ROUTINE zdf_osm_pycnocline_thickness *
	2284	!!
	2285	!! ** Purpose : Calculates thickness of the pycnocline
	2286	!!
	2287	!! ** Method : The thickness is calculated from a prognostic equation
	2288	!! that relaxes the pycnocine thickness to a diagnostic
	2289	!! value. The time change is calculated assuming the
	2290	!! thickness relaxes exponentially. This is done to deal
	2291	!! with large timesteps.
	2292	!!
	2293	!!----------------------------------------------------------------------
	2294
	2295	REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh ! pycnocline thickness.
[13858]	2296	!
[12218]	2297	INTEGER :: jj, ji
	2298	INTEGER :: inhml
[13858]	2299	REAL(wp) :: zari, ztau, zdh_ref
	2300	REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth
[12218]	2301
[13858]	2302	DO jj = 2, jpjm1
	2303	DO ji = 2, jpim1
[12218]	2304
[13858]	2305	IF ( lshear(ji,jj) ) THEN
[12218]	2306	IF ( lconv(ji,jj) ) THEN
[13858]	2307	IF ( j_ddh(ji,jj) == 0 ) THEN
	2308	! ddhdt for pycnocline determined in osm_calculate_dhdt
	2309	dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_rdt
	2310	ELSE
	2311	! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt
	2312	IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
	2313	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 )
	2314	ELSE
	2315	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 )
	2316	ENDIF
	2317	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 )
	2318	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
	2319	IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)
	2320	ENDIF
	2321
	2322	ELSE ! lconv
	2323	! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL
[12218]	2324
[13858]	2325	ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
	2326	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
	2327	! boundary layer deepening
	2328	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
	2329	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
	2330	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
	2331	& / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )
	2332	zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
	2333	ELSE
	2334	zdh_ref = 0.2 * hbl(ji,jj)
	2335	ENDIF
	2336	ELSE ! IF(dhdt < 0)
	2337	zdh_ref = 0.2 * hbl(ji,jj)
	2338	ENDIF ! IF (dhdt >= 0)
	2339	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
	2340	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
	2341	! Alan: this hml is never defined or used -- do we need it?
	2342	ENDIF
	2343
	2344	ELSE ! lshear
	2345	! for lshear = .FALSE. calculate ddhdt here
	2346
	2347	IF ( lconv(ji,jj) ) THEN
	2348
[12218]	2349	IF( ln_osm_mle ) THEN
[13858]	2350	IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN
[12218]	2351	! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
[13858]	2352	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
[13396]	2353	IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability
[13858]	2354	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]	2355	ELSE ! unstable
[13858]	2356	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]	2357	ENDIF
[13396]	2358	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
[13858]	2359	zdh_ref = zari * hbl(ji,jj)
[12218]	2360	ELSE
[13396]	2361	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
[13858]	2362	zdh_ref = 0.2 * hbl(ji,jj)
[12218]	2363	ENDIF
	2364	ELSE
[13396]	2365	ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
[13858]	2366	zdh_ref = 0.2 * hbl(ji,jj)
[12218]	2367	ENDIF
	2368	ELSE ! ln_osm_mle
[13858]	2369	IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
[13396]	2370	IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN ! near neutral stability
[13858]	2371	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]	2372	ELSE ! unstable
[13858]	2373	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]	2374	ENDIF
[13396]	2375	ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
[13858]	2376	zdh_ref = zari * hbl(ji,jj)
[12218]	2377	ELSE
[13396]	2378	ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)*3 + 0.5 zwstrc(ji,jj)3)pthird)
[13858]	2379	zdh_ref = 0.2 * hbl(ji,jj)
[12218]	2380	ENDIF
	2381
	2382	END IF ! ln_osm_mle
	2383
[13858]	2384	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
	2385	! IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
	2386	IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
[12218]	2387	! Alan: this hml is never defined or used
	2388	ELSE ! IF (lconv)
[13396]	2389	ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
[12218]	2390	IF ( zdhdt(ji,jj) >= 0.0 ) THEN ! probably shouldn't include wm here
	2391	! boundary layer deepening
	2392	IF ( zdb_bl(ji,jj) > 0._wp ) THEN
	2393	! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
	2394	zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
[13396]	2395	& / MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01 , 0.2 )
[13858]	2396	zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
[12218]	2397	ELSE
[13858]	2398	zdh_ref = 0.2 * hbl(ji,jj)
[12218]	2399	ENDIF
	2400	ELSE ! IF(dhdt < 0)
[13858]	2401	zdh_ref = 0.2 * hbl(ji,jj)
[12218]	2402	ENDIF ! IF (dhdt >= 0)
[13858]	2403	dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_rdt / ztau ) )
	2404	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]	2405	ENDIF ! IF (lconv)
[13858]	2406	ENDIF ! lshear
	2407
	2408	hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
	2409	inhml = MAX( INT( dh(ji,jj) / MAX(e3t_n(ji,jj,ibld(ji,jj)), 1.e-3) ) , 1 )
	2410	imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
	2411	zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
	2412	zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
	2413	END DO
	2414	END DO
[12218]	2415
	2416	END SUBROUTINE zdf_osm_pycnocline_thickness
	2417
	2418
[13858]	2419	SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
[12218]	2420	!!----------------------------------------------------------------------
	2421	!! * ROUTINE zdf_osm_horizontal_gradients *
	2422	!!
	2423	!! ** Purpose : Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.
	2424	!!
	2425	!! ** Method :
	2426	!!
	2427	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
	2428	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
	2429
	2430
	2431	REAL(wp), DIMENSION(jpi,jpj) :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points
[13858]	2432	REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle ! Magnitude of horizontal buoyancy gradient.
[12218]	2433	REAL(wp), DIMENSION(jpi,jpj) :: zmld ! == estimated FK BLD used for MLE horiz gradients == !
	2434	REAL(wp), DIMENSION(jpi,jpj) :: zdtdx, zdtdy, zdsdx, zdsdy
	2435
	2436	INTEGER :: ji, jj, jk ! dummy loop indices
	2437	INTEGER :: ii, ij, ik, ikmax ! local integers
	2438	REAL(wp) :: zc
	2439	REAL(wp) :: zN2_c ! local buoyancy difference from 10m value
	2440	REAL(wp), DIMENSION(jpi,jpj) :: ztm, zsm, zLf_NH, zLf_MH
	2441	REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv
	2442	REAL(wp), DIMENSION(jpi,jpj) :: zmld_midu, zmld_midv
	2443	!!----------------------------------------------------------------------
	2444	!
	2445	! !== MLD used for MLE ==!
	2446
	2447	mld_prof(:,:) = nlb10 ! Initialization to the number of w ocean point
	2448	zmld(:,:) = 0._wp ! here hmlp used as a dummy variable, integrating vertically N^2
	2449	zN2_c = grav * rn_osm_mle_rho_c * r1_rau0 ! convert density criteria into N^2 criteria
	2450	DO jk = nlb10, jpkm1
[13398]	2451	DO jj = 1, jpj ! Mixed layer level: w-level
[12218]	2452	DO ji = 1, jpi
	2453	ikt = mbkt(ji,jj)
	2454	zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
	2455	IF( zmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
	2456	END DO
	2457	END DO
	2458	END DO
	2459	DO jj = 1, jpj
	2460	DO ji = 1, jpi
	2461	mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
	2462	zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
	2463	END DO
	2464	END DO
	2465	! ensure mld_prof .ge. ibld
	2466	!
	2467	ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 ) ! max level of the computation
	2468	!
	2469	ztm(:,:) = 0._wp
	2470	zsm(:,:) = 0._wp
	2471	DO jk = 1, ikmax ! MLD and mean buoyancy and N2 over the mixed layer
	2472	DO jj = 1, jpj
	2473	DO ji = 1, jpi
	2474	zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1 ) ) ! zc being 0 outside the ML t-points
	2475	ztm(ji,jj) = ztm(ji,jj) + zc * tsn(ji,jj,jk,jp_tem)
	2476	zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal)
	2477	END DO
	2478	END DO
	2479	END DO
	2480	! average temperature and salinity.
	2481	ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
	2482	zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
	2483	! calculate horizontal gradients at u & v points
	2484
	2485	DO jj = 2, jpjm1
	2486	DO ji = 1, jpim1
	2487	zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
	2488	zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) ) * umask(ji,jj,1) / e1u(ji,jj)
	2489	zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))
	2490	ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )
	2491	ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )
	2492	END DO
	2493	END DO
	2494
	2495	DO jj = 1, jpjm1
	2496	DO ji = 2, jpim1
	2497	zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
	2498	zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
	2499	zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))
	2500	ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )
	2501	ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )
	2502	END DO
	2503	END DO
	2504
	2505	CALL eos_rab(ztsm_midu, zmld_midu, zabu)
	2506	CALL eos_rab(ztsm_midv, zmld_midv, zabv)
	2507
	2508	DO jj = 2, jpjm1
	2509	DO ji = 1, jpim1
	2510	dbdx_mle(ji,jj) = grav(zdtdx(ji,jj)zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
	2511	END DO
	2512	END DO
	2513	DO jj = 1, jpjm1
	2514	DO ji = 2, jpim1
	2515	dbdy_mle(ji,jj) = grav(zdtdy(ji,jj)zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
	2516	END DO
	2517	END DO
	2518
[13858]	2519	DO jj = 2, jpjm1
	2520	DO ji = 2, jpim1
	2521	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
	2522	zdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,jj) * dbdx_mle(ji,jj) + dbdy_mle(ji,jj) * dbdy_mle(ji,jj) &
	2523	& + dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
	2524	END DO
	2525	END DO
	2526
[12218]	2527	END SUBROUTINE zdf_osm_zmld_horizontal_gradients
[13858]	2528	SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
[12218]	2529	!!----------------------------------------------------------------------
	2530	!! * ROUTINE zdf_osm_mle_parameters *
	2531	!!
	2532	!! ** Purpose : Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.
	2533	!!
	2534	!! ** Method :
	2535	!!
	2536	!! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
	2537	!! Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
	2538
[13858]	2539	INTEGER, DIMENSION(jpi,jpj) :: mld_prof
	2540	REAL(wp), DIMENSION(jpi,jpj) :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle
[12218]	2541	INTEGER :: ji, jj, jk ! dummy loop indices
[13858]	2542	INTEGER :: ii, ij, ik, jkb, jkb1 ! local integers
[12218]	2543	INTEGER , DIMENSION(jpi,jpj) :: inml_mle
[13858]	2544	REAL(wp) :: ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle
[12218]	2545
[13858]	2546	! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
	2547
	2548	DO jj = 2, jpjm1
	2549	DO ji = 2, jpim1
	2550	IF ( lconv(ji,jj) ) THEN
	2551	ztmp = r1_ft(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
	2552	! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
	2553	zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1)
	2554	zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2
	2555	ENDIF
[12218]	2556	END DO
	2557	END DO
	2558	! Timestep mixed layer eddy depth.
	2559	DO jj = 2, jpjm1
	2560	DO ji = 2, jpim1
[13858]	2561	IF ( lmle(ji,jj) ) THEN ! MLE layer growing.
	2562	! Buoyancy gradient at base of MLE layer.
	2563	zthermal = rab_n(ji,jj,1,jp_tem)
	2564	zbeta = rab_n(ji,jj,1,jp_sal)
	2565	jkb = mld_prof(ji,jj)
	2566	jkb1 = MIN(jkb + 1, mbkt(ji,jj))
	2567	!
	2568	zbuoy = grav * ( zthermal * tsn(ji,jj,mld_prof(ji,jj)+2,jp_tem) - zbeta * tsn(ji,jj,mld_prof(ji,jj)+2,jp_sal) )
	2569	zdb_mle = zb_bl(ji,jj) - zbuoy
	2570	! Timestep hmle.
	2571	hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / zdb_mle
	2572	ELSE
	2573	IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN
	2574	hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt / rn_osm_mle_tau
	2575	ELSE
	2576	hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_rdt /rn_osm_mle_tau
	2577	ENDIF
	2578	ENDIF
	2579	hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj))
	2580	IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) )
[12218]	2581	END DO
	2582	END DO
	2583
	2584	mld_prof = 4
	2585	DO jk = 5, jpkm1
	2586	DO jj = 2, jpjm1
	2587	DO ji = 2, jpim1
	2588	IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
	2589	END DO
	2590	END DO
	2591	END DO
	2592	DO jj = 2, jpjm1
	2593	DO ji = 2, jpim1
	2594	zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj))
	2595	END DO
	2596	END DO
	2597	END SUBROUTINE zdf_osm_mle_parameters
	2598
	2599	END SUBROUTINE zdf_osm
	2600
	2601
[8930]	2602	SUBROUTINE zdf_osm_init
	2603	!!----------------------------------------------------------------------
	2604	!! * ROUTINE zdf_osm_init *
	2605	!!
	2606	!! ** Purpose : Initialization of the vertical eddy diffivity and
	2607	!! viscosity when using a osm turbulent closure scheme
	2608	!!
	2609	!! ** Method : Read the namosm namelist and check the parameters
	2610	!! called at the first timestep (nit000)
	2611	!!
	2612	!! ** input : Namlist namosm
	2613	!!----------------------------------------------------------------------
	2614	INTEGER :: ios ! local integer
	2615	INTEGER :: ji, jj, jk ! dummy loop indices
[12216]	2616	REAL z1_t2
[8930]	2617	!!
	2618	NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
[13400]	2619	& ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd &
[12218]	2620	& ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
[13858]	2621	& ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter
[12218]	2622	! Namelist for Fox-Kemper parametrization.
	2623	NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
[13398]	2624	& rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit
[12218]	2625
[8930]	2626	!!----------------------------------------------------------------------
	2627	!
	2628	REWIND( numnam_ref ) ! Namelist namzdf_osm in reference namelist : Osmosis ML model
	2629	READ ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
[12178]	2630	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist' )
[8930]	2631
	2632	REWIND( numnam_cfg ) ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
	2633	READ ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
[12178]	2634	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist' )
[8930]	2635	IF(lwm) WRITE ( numond, namzdf_osm )
	2636
	2637	IF(lwp) THEN ! Control print
	2638	WRITE(numout,*)
	2639	WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
	2640	WRITE(numout,*) '~~~~~~~~~~~~'
[12218]	2641	WRITE(numout,*) ' Namelist namzdf_osm : set osm mixing parameters'
	2642	WRITE(numout,*) ' Use rn_osm_la ln_use_osm_la = ', ln_use_osm_la
	2643	WRITE(numout,*) ' Use MLE in OBL, i.e. Fox-Kemper param ln_osm_mle = ', ln_osm_mle
[8930]	2644	WRITE(numout,*) ' Turbulent Langmuir number rn_osm_la = ', rn_osm_la
[13400]	2645	WRITE(numout,*) ' Stokes drift reduction factor rn_zdfosm_adjust_sd = ', rn_zdfosm_adjust_sd
[8930]	2646	WRITE(numout,*) ' Initial hbl for 1D runs rn_osm_hbl0 = ', rn_osm_hbl0
[12216]	2647	WRITE(numout,*) ' Depth scale of Stokes drift rn_osm_dstokes = ', rn_osm_dstokes
[8930]	2648	WRITE(numout,*) ' horizontal average flag nn_ave = ', nn_ave
	2649	WRITE(numout,*) ' Stokes drift nn_osm_wave = ', nn_osm_wave
	2650	SELECT CASE (nn_osm_wave)
	2651	CASE(0)
	2652	WRITE(numout,*) ' calculated assuming constant La#=0.3'
	2653	CASE(1)
	2654	WRITE(numout,*) ' calculated from Pierson Moskowitz wind-waves'
	2655	CASE(2)
	2656	WRITE(numout,*) ' calculated from ECMWF wave fields'
[13400]	2657	END SELECT
	2658	WRITE(numout,*) ' Stokes drift reduction nn_osm_SD_reduce', nn_osm_SD_reduce
	2659	WRITE(numout,*) ' fraction of hbl to average SD over/fit'
	2660	WRITE(numout,*) ' exponential with nn_osm_SD_reduce = 1 or 2 rn_osm_hblfrac = ', rn_osm_hblfrac
	2661	SELECT CASE (nn_osm_SD_reduce)
	2662	CASE(0)
	2663	WRITE(numout,*) ' No reduction'
	2664	CASE(1)
	2665	WRITE(numout,*) ' Average SD over upper rn_osm_hblfrac of BL'
	2666	CASE(2)
	2667	WRITE(numout,*) ' Fit exponential to slope rn_osm_hblfrac of BL'
[8930]	2668	END SELECT
[13400]	2669	WRITE(numout,*) ' reduce surface SD and depth scale under ice ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter
[8930]	2670	WRITE(numout,*) ' Output osm diagnostics ln_dia_osm = ', ln_dia_osm
[13858]	2671	WRITE(numout,*) ' Threshold used to define BL rn_osm_bl_thresh = ', rn_osm_bl_thresh, 'm^2/s'
[8930]	2672	WRITE(numout,*) ' Use KPP-style shear instability mixing ln_kpprimix = ', ln_kpprimix
	2673	WRITE(numout,*) ' local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
	2674	WRITE(numout,*) ' maximum shear diffusivity at Rig = 0 (m2/s) rn_difri = ', rn_difri
	2675	WRITE(numout,*) ' Use large mixing below BL when unstable ln_convmix = ', ln_convmix
	2676	WRITE(numout,*) ' diffusivity when unstable below BL (m2/s) rn_difconv = ', rn_difconv
	2677	ENDIF
	2678
[13858]	2679
	2680	! ! Check wave coupling settings !
	2681	! ! Further work needed - see ticket #2447 !
	2682	IF( nn_osm_wave == 2 ) THEN
	2683	IF (.NOT. ( ln_wave .AND. ln_sdw )) &
	2684	& CALL ctl_stop( 'zdf_osm_init : ln_zdfosm and nn_osm_wave=2, ln_wave and ln_sdw must be true' )
	2685	END IF
	2686
[8930]	2687	! ! allocate zdfosm arrays
	2688	IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
	2689
	2690
[12218]	2691	IF( ln_osm_mle ) THEN
	2692	! Initialise Fox-Kemper parametrization
	2693	REWIND( numnam_ref ) ! Namelist namosm_mle in reference namelist : Tracer advection scheme
	2694	READ ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
	2695	903 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namosm_mle in reference namelist')
	2696
	2697	REWIND( numnam_cfg ) ! Namelist namosm_mle in configuration namelist : Tracer advection scheme
	2698	READ ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
	2699	904 IF( ios > 0 ) CALL ctl_nam ( ios , 'namosm_mle in configuration namelist')
	2700	IF(lwm) WRITE ( numond, namosm_mle )
	2701
	2702	IF(lwp) THEN ! Namelist print
	2703	WRITE(numout,*)
	2704	WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
	2705	WRITE(numout,*) '~~~~~~~~~~~~~'
	2706	WRITE(numout,*) ' Namelist namosm_mle : '
	2707	WRITE(numout,*) ' MLE type: =0 standard Fox-Kemper ; =1 new formulation nn_osm_mle = ', nn_osm_mle
	2708	WRITE(numout,*) ' magnitude of the MLE (typical value: 0.06 to 0.08) rn_osm_mle_ce = ', rn_osm_mle_ce
	2709	WRITE(numout,*) ' scale of ML front (ML radius of deformation) (nn_osm_mle=0) rn_osm_mle_lf = ', rn_osm_mle_lf, 'm'
	2710	WRITE(numout,*) ' maximum time scale of MLE (nn_osm_mle=0) rn_osm_mle_time = ', rn_osm_mle_time, 's'
	2711	WRITE(numout,*) ' reference latitude (degrees) of MLE coef. (nn_osm_mle=1) rn_osm_mle_lat = ', rn_osm_mle_lat, 'deg'
	2712	WRITE(numout,*) ' Density difference used to define ML for FK rn_osm_mle_rho_c = ', rn_osm_mle_rho_c
[13858]	2713	WRITE(numout,*) ' Threshold used to define MLE for FK rn_osm_mle_thresh = ', rn_osm_mle_thresh, 'm^2/s'
[12218]	2714	WRITE(numout,*) ' Timescale for OSM-FK rn_osm_mle_tau = ', rn_osm_mle_tau, 's'
[13398]	2715	WRITE(numout,*) ' switch to limit hmle ln_osm_hmle_limit = ', ln_osm_hmle_limit
	2716	WRITE(numout,*) ' fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T. rn_osm_hmle_limit = ', rn_osm_hmle_limit
[12218]	2717	ENDIF !
	2718	ENDIF
	2719	!
	2720	IF(lwp) THEN
	2721	WRITE(numout,*)
	2722	IF( ln_osm_mle ) THEN
	2723	WRITE(numout,*) ' ==>>> Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
	2724	IF( nn_osm_mle == 0 ) WRITE(numout,*) ' Fox-Kemper et al 2010 formulation'
	2725	IF( nn_osm_mle == 1 ) WRITE(numout,*) ' New formulation'
	2726	ELSE
	2727	WRITE(numout,*) ' ==>>> Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
	2728	ENDIF
	2729	ENDIF
	2730	!
	2731	IF( ln_osm_mle ) THEN ! MLE initialisation
	2732	!
	2733	rb_c = grav * rn_osm_mle_rho_c /rau0 ! Mixed Layer buoyancy criteria
	2734	IF(lwp) WRITE(numout,*)
	2735	IF(lwp) WRITE(numout,*) ' ML buoyancy criteria = ', rb_c, ' m/s2 '
	2736	IF(lwp) WRITE(numout,*) ' associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
	2737	!
	2738	IF( nn_osm_mle == 0 ) THEN ! MLE array allocation & initialisation !
	2739	!
	2740	ELSEIF( nn_osm_mle == 1 ) THEN ! MLE array allocation & initialisation
	2741	rc_f = rn_osm_mle_ce/ ( 5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat ) )
	2742	!
	2743	ENDIF
	2744	! ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
	2745	z1_t2 = 2.e-5
	2746	do jj=1,jpj
	2747	do ji = 1,jpi
	2748	r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
	2749	end do
	2750	end do
	2751	! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
	2752	! r1_ft(:,:) = 1._wp / SQRT( ff_t(:,:) * ff_t(:,:) + z1_t2 )
	2753	!
	2754	ENDIF
	2755
	2756	call osm_rst( nit000, 'READ' ) !* read or initialize hbl, dh, hmle
	2757
	2758
[8930]	2759	IF( ln_zdfddm) THEN
	2760	IF(lwp) THEN
	2761	WRITE(numout,*)
	2762	WRITE(numout,*) ' Double diffusion mixing on temperature and salinity '
	2763	WRITE(numout,*) ' CAUTION : done in routine zdfosm, not in routine zdfddm '
	2764	ENDIF
	2765	ENDIF
	2766
	2767
	2768	!set constants not in namelist
	2769	!-----------------------------
	2770
	2771	IF(lwp) THEN
	2772	WRITE(numout,*)
	2773	ENDIF
	2774
	2775	IF (nn_osm_wave == 0) THEN
	2776	dstokes(:,:) = rn_osm_dstokes
	2777	END IF
	2778
	2779	! Horizontal average : initialization of weighting arrays
	2780	! -------------------
	2781
	2782	SELECT CASE ( nn_ave )
	2783
	2784	CASE ( 0 ) ! no horizontal average
	2785	IF(lwp) WRITE(numout,*) ' no horizontal average on avt'
	2786	IF(lwp) WRITE(numout,*) ' only in very high horizontal resolution !'
	2787	! weighting mean arrays etmean
	2788	! ( 1 1 )
	2789	! avt = 1/4 ( 1 1 )
	2790	!
	2791	etmean(:,:,:) = 0.e0
	2792
	2793	DO jk = 1, jpkm1
	2794	DO jj = 2, jpjm1
	2795	DO ji = 2, jpim1 ! vector opt.
	2796	etmean(ji,jj,jk) = tmask(ji,jj,jk) &
	2797	& / MAX( 1., umask(ji-1,jj ,jk) + umask(ji,jj,jk) &
	2798	& + vmask(ji ,jj-1,jk) + vmask(ji,jj,jk) )
	2799	END DO
	2800	END DO
	2801	END DO
	2802
	2803	CASE ( 1 ) ! horizontal average
	2804	IF(lwp) WRITE(numout,*) ' horizontal average on avt'
	2805	! weighting mean arrays etmean
	2806	! ( 1/2 1 1/2 )
	2807	! avt = 1/8 ( 1 2 1 )
	2808	! ( 1/2 1 1/2 )
	2809	etmean(:,:,:) = 0.e0
	2810
	2811	DO jk = 1, jpkm1
	2812	DO jj = 2, jpjm1
	2813	DO ji = 2, jpim1 ! vector opt.
	2814	etmean(ji,jj,jk) = tmask(ji, jj,jk) &
	2815	& / MAX( 1., 2.* tmask(ji,jj,jk) &
	2816	& +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk) &
	2817	& +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
	2818	& +1. * ( tmask(ji-1,jj ,jk) + tmask(ji ,jj+1,jk) &
	2819	& +tmask(ji ,jj-1,jk) + tmask(ji+1,jj ,jk) ) )
	2820	END DO
	2821	END DO
	2822	END DO
	2823
	2824	CASE DEFAULT
	2825	WRITE(ctmp1,*) ' bad flag value for nn_ave = ', nn_ave
	2826	CALL ctl_stop( ctmp1 )
	2827
	2828	END SELECT
	2829
	2830	! Initialization of vertical eddy coef. to the background value
	2831	! -------------------------------------------------------------
	2832	DO jk = 1, jpk
	2833	avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
	2834	END DO
	2835
	2836	! zero the surface flux for non local term and osm mixed layer depth
	2837	! ------------------------------------------------------------------
	2838	ghamt(:,:,:) = 0.
	2839	ghams(:,:,:) = 0.
	2840	ghamu(:,:,:) = 0.
	2841	ghamv(:,:,:) = 0.
	2842	!
[9367]	2843	IF( lwxios ) THEN
	2844	CALL iom_set_rstw_var_active('wn')
	2845	CALL iom_set_rstw_var_active('hbl')
[12319]	2846	CALL iom_set_rstw_var_active('dh')
	2847	IF( ln_osm_mle ) THEN
	2848	CALL iom_set_rstw_var_active('hmle')
	2849	END IF
[9367]	2850	ENDIF
[8930]	2851	END SUBROUTINE zdf_osm_init
	2852
[8946]	2853
[8930]	2854	SUBROUTINE osm_rst( kt, cdrw )
	2855	!!---------------------------------------------------------------------
	2856	!! * ROUTINE osm_rst *
	2857	!!
	2858	!! ** Purpose : Read or write BL fields in restart file
	2859	!!
	2860	!! ** Method : use of IOM library. If the restart does not contain
	2861	!! required fields, they are recomputed from stratification
	2862	!!----------------------------------------------------------------------
	2863
	2864	INTEGER, INTENT(in) :: kt
	2865	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
	2866
[12317]	2867	INTEGER :: id1, id2, id3 ! iom enquiry index
[8930]	2868	INTEGER :: ji, jj, jk ! dummy loop indices
	2869	INTEGER :: iiki, ikt ! local integer
	2870	REAL(wp) :: zhbf ! tempory scalars
	2871	REAL(wp) :: zN2_c ! local scalar
	2872	REAL(wp) :: rho_c = 0.01_wp !: density criterion for mixed layer depth
[12216]	2873	INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
[8930]	2874	!!----------------------------------------------------------------------
	2875	!
	2876	!!-----------------------------------------------------------------------------
	2877	! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
	2878	!!-----------------------------------------------------------------------------
	2879	IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
	2880	id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. )
	2881	IF( id1 > 0 ) THEN ! 'wn' exists; read
[9367]	2882	CALL iom_get( numror, jpdom_autoglo, 'wn', wn, ldxios = lrxios )
[8930]	2883	WRITE(numout,*) ' ===>>>> : wn read from restart file'
	2884	ELSE
	2885	wn(:,:,:) = 0._wp
	2886	WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
	2887	END IF
[12216]	2888
[8930]	2889	id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. )
[12216]	2890	id2 = iom_varid( numror, 'dh' , ldstop = .FALSE. )
[8930]	2891	IF( id1 > 0 .AND. id2 > 0) THEN ! 'hbl' exists; read and return
[9367]	2892	CALL iom_get( numror, jpdom_autoglo, 'hbl' , hbl , ldxios = lrxios )
[12216]	2893	CALL iom_get( numror, jpdom_autoglo, 'dh', dh, ldxios = lrxios )
	2894	WRITE(numout,*) ' ===>>>> : hbl & dh read from restart file'
[12317]	2895	IF( ln_osm_mle ) THEN
	2896	id3 = iom_varid( numror, 'hmle' , ldstop = .FALSE. )
	2897	IF( id3 > 0) THEN
	2898	CALL iom_get( numror, jpdom_autoglo, 'hmle' , hmle , ldxios = lrxios )
	2899	WRITE(numout,*) ' ===>>>> : hmle read from restart file'
	2900	ELSE
	2901	WRITE(numout,*) ' ===>>>> : hmle not found, set to hbl'
	2902	hmle(:,:) = hbl(:,:) ! Initialise MLE depth.
	2903	END IF
	2904	END IF
[8930]	2905	RETURN
[12216]	2906	ELSE ! 'hbl' & 'dh' not in restart file, recalculate
[8930]	2907	WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
	2908	END IF
	2909	END IF
	2910
	2911	!!-----------------------------------------------------------------------------
	2912	! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
	2913	!!-----------------------------------------------------------------------------
	2914	IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return
	2915	IF(lwp) WRITE(numout,*) '---- osm-rst ----'
[12317]	2916	CALL iom_rstput( kt, nitrst, numrow, 'wn' , wn, ldxios = lwxios )
	2917	CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl, ldxios = lwxios )
	2918	CALL iom_rstput( kt, nitrst, numrow, 'dh' , dh, ldxios = lwxios )
	2919	IF( ln_osm_mle ) THEN
	2920	CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios )
	2921	END IF
[8930]	2922	RETURN
	2923	END IF
	2924
	2925	!!-----------------------------------------------------------------------------
	2926	! Getting hbl, no restart file with hbl, so calculate from surface stratification
	2927	!!-----------------------------------------------------------------------------
	2928	IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
	2929	! w-level of the mixing and mixed layers
	2930	CALL eos_rab( tsn, rab_n )
	2931	CALL bn2(tsn, rab_n, rn2)
	2932	imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point
	2933	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
	2934	zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
	2935	!
	2936	hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2
	2937	DO jk = 1, jpkm1
	2938	DO jj = 1, jpj ! Mixed layer level: w-level
	2939	DO ji = 1, jpi
	2940	ikt = mbkt(ji,jj)
	2941	hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
	2942	IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level
	2943	END DO
	2944	END DO
	2945	END DO
	2946	!
	2947	DO jj = 1, jpj
	2948	DO ji = 1, jpi
[12216]	2949	iiki = MAX(4,imld_rst(ji,jj))
	2950	hbl (ji,jj) = gdepw_n(ji,jj,iiki ) ! Turbocline depth
	2951	dh (ji,jj) = e3t_n(ji,jj,iiki-1 ) ! Turbocline depth
[8930]	2952	END DO
	2953	END DO
[12218]	2954
[12317]	2955	WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
[12218]	2956
[12317]	2957	IF( ln_osm_mle ) THEN
	2958	hmle(:,:) = hbl(:,:) ! Initialise MLE depth.
	2959	WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
	2960	END IF
	2961
[12216]	2962	wn(:,:,:) = 0._wp
	2963	WRITE(numout,*) ' ===>>>> : wn not in restart file, set to zero initially'
[8930]	2964	END SUBROUTINE osm_rst
	2965
[8946]	2966
[8930]	2967	SUBROUTINE tra_osm( kt )
	2968	!!----------------------------------------------------------------------
	2969	!! * ROUTINE tra_osm *
	2970	!!
	2971	!! ** Purpose : compute and add to the tracer trend the non-local tracer flux
	2972	!!
	2973	!! ** Method : ???
	2974	!!----------------------------------------------------------------------
	2975	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace
	2976	!!----------------------------------------------------------------------
	2977	INTEGER, INTENT(in) :: kt
	2978	INTEGER :: ji, jj, jk
	2979	!
	2980	IF( kt == nit000 ) THEN
	2981	IF(lwp) WRITE(numout,*)
	2982	IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
	2983	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	2984	ENDIF
	2985
	2986	IF( l_trdtra ) THEN !* Save ta and sa trends
	2987	ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
	2988	ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)
	2989	ENDIF
	2990
	2991	DO jk = 1, jpkm1
	2992	DO jj = 2, jpjm1
	2993	DO ji = 2, jpim1
	2994	tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &
	2995	& - ( ghamt(ji,jj,jk ) &
	2996	& - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
	2997	tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) &
	2998	& - ( ghams(ji,jj,jk ) &
	2999	& - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
	3000	END DO
	3001	END DO
	3002	END DO
	3003
[12216]	3004	! save the non-local tracer flux trends for diagnostics
[8930]	3005	IF( l_trdtra ) THEN
	3006	ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
	3007	ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
[12216]	3008
	3009	CALL trd_tra( kt, 'TRA', jp_tem, jptra_osm, ztrdt )
	3010	CALL trd_tra( kt, 'TRA', jp_sal, jptra_osm, ztrds )
[8930]	3011	DEALLOCATE( ztrdt ) ; DEALLOCATE( ztrds )
	3012	ENDIF
	3013
	3014	IF(ln_ctl) THEN
	3015	CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm - Ta: ', mask1=tmask, &
	3016	& tab3d_2=tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
	3017	ENDIF
	3018	!
	3019	END SUBROUTINE tra_osm
	3020
[8946]	3021
[8930]	3022	SUBROUTINE trc_osm( kt ) ! Dummy routine
	3023	!!----------------------------------------------------------------------
	3024	!! * ROUTINE trc_osm *
	3025	!!
	3026	!! ** Purpose : compute and add to the passive tracer trend the non-local
	3027	!! passive tracer flux
	3028	!!
	3029	!!
	3030	!! ** Method : ???
	3031	!!----------------------------------------------------------------------
[8946]	3032	!
[8930]	3033	!!----------------------------------------------------------------------
	3034	INTEGER, INTENT(in) :: kt
	3035	WRITE(,) 'trc_osm: Not written yet', kt
	3036	END SUBROUTINE trc_osm
	3037
[8946]	3038
[8930]	3039	SUBROUTINE dyn_osm( kt )
	3040	!!----------------------------------------------------------------------
	3041	!! * ROUTINE dyn_osm *
	3042	!!
	3043	!! ** Purpose : compute and add to the velocity trend the non-local flux
	3044	!! copied/modified from tra_osm
	3045	!!
	3046	!! ** Method : ???
	3047	!!----------------------------------------------------------------------
[8946]	3048	INTEGER, INTENT(in) :: kt !
	3049	!
	3050	INTEGER :: ji, jj, jk ! dummy loop indices
[8930]	3051	!!----------------------------------------------------------------------
	3052	!
	3053	IF( kt == nit000 ) THEN
	3054	IF(lwp) WRITE(numout,*)
	3055	IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
	3056	IF(lwp) WRITE(numout,*) '~~~~~~~ '
	3057	ENDIF
	3058	!code saving tracer trends removed, replace with trdmxl_oce
	3059
[8946]	3060	DO jk = 1, jpkm1 ! add non-local u and v fluxes
[8930]	3061	DO jj = 2, jpjm1
	3062	DO ji = 2, jpim1
	3063	ua(ji,jj,jk) = ua(ji,jj,jk) &
	3064	& - ( ghamu(ji,jj,jk ) &
	3065	& - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
	3066	va(ji,jj,jk) = va(ji,jj,jk) &
	3067	& - ( ghamv(ji,jj,jk ) &
	3068	& - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
	3069	END DO
	3070	END DO
	3071	END DO
[9089]	3072	!
[8930]	3073	! code for saving tracer trends removed
	3074	!
	3075	END SUBROUTINE dyn_osm
	3076
[8946]	3077	!!======================================================================
[12216]	3078
[8930]	3079	END MODULE zdfosm

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

Download in other formats: